1 #include <pxd/modules/pxdBackground/PXDBackgroundModule.h>
3 #include <framework/datastore/DataStore.h>
4 #include <framework/datastore/StoreObjPtr.h>
5 #include <framework/datastore/StoreArray.h>
6 #include <framework/datastore/RelationArray.h>
8 #include <framework/dataobjects/FileMetaData.h>
9 #include <framework/dataobjects/BackgroundMetaData.h>
10 #include <mdst/dataobjects/MCParticle.h>
11 #include <pxd/dataobjects/PXDSimHit.h>
12 #include <pxd/dataobjects/PXDTrueHit.h>
13 #include <pxd/dataobjects/PXDDigit.h>
14 #include <pxd/dataobjects/PXDCluster.h>
15 #include <pxd/dataobjects/PXDEnergyDepositionEvent.h>
16 #include <pxd/dataobjects/PXDNeutronFluxEvent.h>
17 #include <pxd/dataobjects/PXDOccupancyEvent.h>
20 #include <boost/format.hpp>
42 Module(), m_outputDirectoryName(""),
43 m_doseReportingLevel(c_reportNTuple),
44 m_nfluxReportingLevel(c_reportNTuple),
45 m_occupancyReportingLevel(c_reportNTuple),
46 m_componentName(""), m_componentTime(0), m_integrationTime(20),
47 m_nielNeutrons(new
TNiel(c_niel_neutronFile)),
48 m_nielProtons(new
TNiel(c_niel_protonFile)),
49 m_nielPions(new
TNiel(c_niel_pionFile)),
50 m_nielElectrons(new
TNiel(c_niel_electronFile))
53 setDescription(
"PXD background module");
54 setPropertyFlags(c_ParallelProcessingCertified);
56 addParam(
"componentName", m_componentName,
"Background component name to process", m_componentName);
57 addParam(
"componentTime", m_componentTime,
"Background component time", m_componentTime);
58 addParam(
"integrationTime", m_integrationTime,
"PXD integration time", m_integrationTime);
59 addParam(
"doseReportingLevel", m_doseReportingLevel,
"0 - no data, 1 - summary only, 2 - summary + ntuple", m_doseReportingLevel);
60 addParam(
"nfluxReportingLevel", m_nfluxReportingLevel,
"0 - no data, 1 - summary only, 2 - summary + ntuple",
61 m_nfluxReportingLevel);
62 addParam(
"occupancyReportingLevel", m_occupancyReportingLevel,
"0 - no data, 1 - summary only, 2 - summary + ntuple",
63 m_occupancyReportingLevel);
64 addParam(
"outputDirectory", m_outputDirectoryName,
"Name of output directory", m_outputDirectoryName);
67 const TVector3& PXDBackgroundModule::pointToGlobal(
VxdID sensorID,
const TVector3& local)
69 static TVector3 result(0, 0, 0);
72 result = info.pointToGlobal(local);
76 const TVector3& PXDBackgroundModule::vectorToGlobal(
VxdID sensorID,
const TVector3& local)
78 static TVector3 result(0, 0, 0);
81 result = info.vectorToGlobal(local);
85 PXDBackgroundModule::~PXDBackgroundModule()
89 void PXDBackgroundModule::initialize()
100 RelationArray relDigitsMCParticles(storeDigits, storeMCParticles);
102 RelationArray relMCParticlesTrueHits(storeMCParticles, storeTrueHits);
103 RelationArray relTrueHitsSimHits(storeTrueHits, storeSimHits);
107 storeEnergyDeposits.registerInDataStore();
109 storeNeutronFluxes.registerInDataStore();
111 storeOccupancyEvents.registerInDataStore();
114 m_storeFileMetaDataName = storeFileMetaData.getName();
115 m_storeBgMetaDataName = storeBgMetaData.getName();
116 m_storeMCParticlesName = storeMCParticles.getName();
117 m_storeSimHitsName = storeSimHits.getName();
118 m_storeTrueHitsName = storeTrueHits.getName();
119 m_storeDigitsName = storeDigits.getName();
120 m_relDigitsMCParticlesName = relDigitsMCParticles.
getName();
121 m_relDigitsTrueHitsName = relDigitsTrueHits.
getName();
122 m_relParticlesTrueHitsName = relMCParticlesTrueHits.
getName();
123 m_relTrueHitsSimHitsName = relTrueHitsSimHits.
getName();
124 m_storeEnergyDepositsName = storeEnergyDeposits.getName();
125 m_storeNeutronFluxesName = storeNeutronFluxes.getName();
127 m_componentTime *= Unit::us;
128 m_integrationTime *= Unit::us;
131 void PXDBackgroundModule::beginRun()
135 void PXDBackgroundModule::event()
152 RelationArray relDigitsMCParticles(storeDigits, storeMCParticles, m_relDigitsMCParticlesName);
153 RelationArray relDigitsTrueHits(storeDigits, storeTrueHits, m_relDigitsTrueHitsName);
154 RelationArray relTrueHitsSimHits(storeTrueHits, storeSimHits, m_relTrueHitsSimHitsName);
155 RelationArray relTrueHitsMCParticles(storeMCParticles, storeTrueHits, m_relParticlesTrueHitsName);
158 double currentComponentTime = storeBgMetaData->getRealTime();
159 if (currentComponentTime != m_componentTime)
160 B2FATAL(
"Mismatch in component times:\n"
161 <<
"Steering file: " << m_componentTime <<
"\n"
162 <<
"Background file: " << currentComponentTime);
164 VxdID currentSensorID(0);
165 double currentSensorThickness(0);
166 double currentSensorMass(0);
167 double currentSensorArea(0);
170 if (m_doseReportingLevel > c_reportNone) {
171 B2DEBUG(100,
"Expo and dose");
172 currentSensorID.
setID(0);
173 for (
const PXDSimHit& hit : storeSimHits) {
175 VxdID sensorID = hit.getSensorID();
176 if (sensorID != currentSensorID) {
177 currentSensorID = sensorID;
178 currentSensorThickness = getSensorThickness(currentSensorID);
179 currentSensorMass = getSensorMass(currentSensorID);
180 currentSensorArea = getSensorArea(currentSensorID);
182 double hitEnergy = hit.getElectrons() * Const::ehEnergy;
184 m_sensorData[currentSensorID].m_dose +=
185 (hitEnergy / Unit::J) / (currentSensorMass / 1000) * (c_smy / currentComponentTime);
187 m_sensorData[currentSensorID].m_expo += hitEnergy / currentSensorArea / (currentComponentTime / Unit::s);
188 if (m_doseReportingLevel == c_reportNTuple) {
189 const TVector3 localPos = hit.getPosIn();
190 const TVector3 globalPos = pointToGlobal(currentSensorID, localPos);
191 float globalPosXYZ[3];
192 globalPos.GetXYZ(globalPosXYZ);
195 hit.getPDGcode(), hit.getGlobalTime(),
196 localPos.X(), localPos.Y(), globalPosXYZ, hitEnergy,
197 (hitEnergy / Unit::J) / (currentSensorMass / 1000) / (currentComponentTime / Unit::s),
198 (hitEnergy / Unit::J) / currentSensorArea / (currentComponentTime / Unit::s)
205 if (m_nfluxReportingLevel > c_reportNone) {
206 B2DEBUG(100,
"Neutron flux");
207 currentSensorID.
setID(0);
209 VxdID sensorID = hit.getSensorID();
211 if (sensorID != currentSensorID) {
212 currentSensorID = sensorID;
213 currentSensorThickness = getSensorThickness(currentSensorID);
214 currentSensorMass = getSensorMass(currentSensorID);
215 currentSensorArea = getSensorArea(currentSensorID);
218 TVector3 entryPos(hit.getEntryU(), hit.getEntryV(), hit.getEntryW());
219 TVector3 exitPos(hit.getExitU(), hit.getExitV(), hit.getExitW());
220 double stepLength = (exitPos - entryPos).Mag();
226 double minDistance = 1.0e10;
227 for (
const PXDSimHit& related : storeSimHits) {
228 double distance = (entryPos - related.getPosIn()).Mag();
229 if (distance < minDistance) {
230 minDistance = distance;
239 TVector3 hitMomentum(hit.getMomentum());
240 hitMomentum.SetX(std::isfinite(hitMomentum.X()) ? hitMomentum.X() : 0.0);
241 hitMomentum.SetY(std::isfinite(hitMomentum.Y()) ? hitMomentum.Y() : 0.0);
242 hitMomentum.SetZ(std::isfinite(hitMomentum.Z()) ? hitMomentum.Z() : 0.0);
244 double kineticEnergy(0.0);
245 double nielWeight(0.0);
248 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
249 nielWeight = m_nielNeutrons->getNielFactor(kineticEnergy / Unit::MeV);
253 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
254 nielWeight = m_nielProtons->getNielFactor(kineticEnergy / Unit::MeV);
256 if (pdg == 111 || pdg == 211) {
258 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
259 nielWeight = m_nielPions->getNielFactor(kineticEnergy / Unit::MeV);
262 double m0 = 0.511e-3;
263 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
264 nielWeight = m_nielElectrons->getNielFactor(kineticEnergy / Unit::MeV);
268 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
272 nielWeight = std::isfinite(nielWeight) ? nielWeight : 0.0;
273 m_sensorData[currentSensorID].m_neutronFlux += nielWeight * stepLength / currentSensorThickness / currentSensorArea /
274 currentComponentTime * c_smy;
277 if (m_nfluxReportingLevel == c_reportNTuple) {
278 TVector3 localPos(hit.getU(), hit.getV(), hit.getW());
279 const TVector3 globalPos = pointToGlobal(currentSensorID, localPos);
280 float globalPosXYZ[3];
281 globalPos.GetXYZ(globalPosXYZ);
282 TVector3 localMom = hit.getMomentum();
283 const TVector3 globalMom = vectorToGlobal(currentSensorID, localMom);
284 float globalMomXYZ[3];
285 globalMom.GetXYZ(globalMomXYZ);
289 hit.getU(), hit.getV(), globalPosXYZ, globalMomXYZ, kineticEnergy,
290 stepLength, nielWeight,
291 stepLength / currentSensorThickness / currentSensorArea / (currentComponentTime / Unit::s),
292 nielWeight * stepLength / currentSensorThickness / currentSensorArea / (currentComponentTime / Unit::s)
299 if (m_occupancyReportingLevel > c_reportNone) {
300 B2DEBUG(100,
"Fired pixels");
301 currentSensorID.
setID(0);
302 double currentSensorCut = 0;
304 std::map<VxdID, std::vector<float> > firedPixels;
305 for (
const PXDDigit& storeDigit : storeDigits) {
308 VxdID sensorID = storeDigit.getSensorID();
309 if (sensorID != currentSensorID) {
310 currentSensorID = sensorID;
311 auto info = getInfo(sensorID);
312 currentSensorCut = info.getChargeThreshold();
314 B2DEBUG(30,
"Digit charge: " << storeDigit.getCharge() <<
" threshold: " << currentSensorCut);
315 if (storeDigit.getCharge() < currentSensorCut)
continue;
316 B2DEBUG(30,
"Passed.");
317 firedPixels[sensorID].push_back(storeDigit.getCharge());
320 for (
auto idAndSet : firedPixels) {
321 VxdID sensorID = idAndSet.first;
322 double sensorArea = getSensorArea(sensorID);
323 int nFired = idAndSet.second.size();
324 double fired = nFired / (currentComponentTime / Unit::s) / sensorArea;
325 m_sensorData[sensorID].m_fired += fired;
328 B2DEBUG(100,
"Occupancy");
329 currentSensorID.
setID(0);
331 for (
auto cluster : storeClsuters) {
332 VxdID sensorID = cluster.getSensorID();
333 if (currentSensorID != sensorID) {
334 currentSensorID = sensorID;
335 auto info = getInfo(sensorID);
336 nPixels = info.getUCells() * info.getVCells();
339 double w_acceptance = m_integrationTime / currentComponentTime;
340 double occupancy = 1.0 / nPixels * cluster.getSize();
341 m_sensorData[sensorID].m_occupancy += w_acceptance * occupancy;
343 if (m_occupancyReportingLevel == c_reportNTuple) {
347 cluster.getU(), cluster.getV(), cluster.getSize(),
348 cluster.getCharge(), occupancy
355 void PXDBackgroundModule::endRun()
360 void PXDBackgroundModule::terminate()
364 string outfileName(m_outputDirectoryName + m_componentName +
"_summary.txt");
365 outfile.open(outfileName.c_str(), ios::out | ios::trunc);
366 outfile <<
"component_name\t"
367 <<
"component_time\t"
377 double componentTime = m_componentTime / Unit::us;
378 for (
auto vxdSensor : m_sensorData) {
379 outfile << m_componentName.c_str() <<
"\t"
380 << componentTime <<
"\t"
381 << vxdSensor.first.getLayerNumber() <<
"\t"
382 << vxdSensor.first.getLadderNumber() <<
"\t"
383 << vxdSensor.first.getSensorNumber() <<
"\t"
384 << vxdSensor.second.m_dose <<
"\t"
385 << vxdSensor.second.m_expo <<
"\t"
386 << vxdSensor.second.m_neutronFlux <<
"\t"
387 << vxdSensor.second.m_fired <<
"\t"
388 << vxdSensor.second.m_occupancy