Belle II Software  release-05-02-19
PXDBackgroundModule.cc
1 #include <pxd/modules/pxdBackground/PXDBackgroundModule.h>
2 
3 #include <framework/datastore/DataStore.h>
4 #include <framework/datastore/StoreObjPtr.h>
5 #include <framework/datastore/StoreArray.h>
6 #include <framework/datastore/RelationArray.h>
7 
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>
18 #include <cmath>
19 #include <fstream>
20 #include <boost/format.hpp>
21 
22 #include "TVector3.h"
23 
24 using namespace std;
25 using boost::format;
26 using namespace Belle2;
27 using namespace Belle2::PXD;
28 
29 
30 
31 //-----------------------------------------------------------------
32 // Register the Module
33 //-----------------------------------------------------------------
34 REG_MODULE(PXDBackground)
35 
36 
37 //-----------------------------------------------------------------
38 // Implementation
39 //-----------------------------------------------------------------
40 
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))
51 {
52  //Set module properties
53  setDescription("PXD background module");
54  setPropertyFlags(c_ParallelProcessingCertified); // specify this flag if you need parallel processing
55  // FIXME: This information can in principle be extracted from bg files, though not trivially.
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);
65 }
66 
67 const TVector3& PXDBackgroundModule::pointToGlobal(VxdID sensorID, const TVector3& local)
68 {
69  static TVector3 result(0, 0, 0);
70 
71  const PXD::SensorInfo& info = getInfo(sensorID);
72  result = info.pointToGlobal(local);
73  return result;
74 }
75 
76 const TVector3& PXDBackgroundModule::vectorToGlobal(VxdID sensorID, const TVector3& local)
77 {
78  static TVector3 result(0, 0, 0);
79 
80  const PXD::SensorInfo& info = getInfo(sensorID);
81  result = info.vectorToGlobal(local);
82  return result;
83 }
84 
85 PXDBackgroundModule::~PXDBackgroundModule()
86 {
87 }
88 
89 void PXDBackgroundModule::initialize()
90 {
91  //Register collections
92  StoreObjPtr<FileMetaData> storeFileMetaData(m_storeFileMetaDataName, DataStore::c_Persistent);
93  StoreObjPtr<BackgroundMetaData> storeBgMetaData(m_storeBgMetaDataName, DataStore::c_Persistent);
94  StoreArray<MCParticle> storeMCParticles(m_storeMCParticlesName);
95  StoreArray<PXDSimHit> storeSimHits(m_storeSimHitsName);
96  StoreArray<PXDTrueHit> storeTrueHits(m_storeTrueHitsName);
97  StoreArray<PXDDigit> storeDigits(m_storeDigitsName);
98  StoreArray<PXDCluster> storeClusters(m_storeClustersName);
99 
100  RelationArray relDigitsMCParticles(storeDigits, storeMCParticles);
101  RelationArray relDigitsTrueHits(storeDigits, storeTrueHits);
102  RelationArray relMCParticlesTrueHits(storeMCParticles, storeTrueHits);
103  RelationArray relTrueHitsSimHits(storeTrueHits, storeSimHits);
104 
105  // Add two new StoreArrays
106  StoreArray<PXDEnergyDepositionEvent> storeEnergyDeposits(m_storeEnergyDepositsName);
107  storeEnergyDeposits.registerInDataStore();
108  StoreArray<PXDNeutronFluxEvent> storeNeutronFluxes(m_storeNeutronFluxesName);
109  storeNeutronFluxes.registerInDataStore();
110  StoreArray<PXDOccupancyEvent> storeOccupancyEvents(m_storeOccupancyEventsName);
111  storeOccupancyEvents.registerInDataStore();
112 
113  //Store names to speed up creation later
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();
126 
127  m_componentTime *= Unit::us;
128  m_integrationTime *= Unit::us;
129 }
130 
131 void PXDBackgroundModule::beginRun()
132 {
133 }
134 
135 void PXDBackgroundModule::event()
136 {
137  //Register collections
138  const StoreObjPtr<FileMetaData> storeFileMetaData(m_storeFileMetaDataName, DataStore::c_Persistent);
139  const StoreObjPtr<BackgroundMetaData> storeBgMetaData(m_storeBgMetaDataName, DataStore::c_Persistent);
140  const StoreArray<MCParticle> storeMCParticles(m_storeMCParticlesName);
141  const StoreArray<PXDSimHit> storeSimHits(m_storeSimHitsName);
142  const StoreArray<PXDTrueHit> storeTrueHits(m_storeTrueHitsName);
143  const StoreArray<PXDDigit> storeDigits(m_storeDigitsName);
144  const StoreArray<PXDCluster> storeClsuters(m_storeClustersName);
145 
146  // Add two new StoreArrays
147  StoreArray<PXDEnergyDepositionEvent> storeEnergyDeposits(m_storeEnergyDepositsName);
148  StoreArray<PXDNeutronFluxEvent> storeNeutronFluxes(m_storeNeutronFluxesName);
149  StoreArray<PXDOccupancyEvent> storeOccupancyEvents(m_storeOccupancyEventsName);
150 
151  // Relations
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);
156 
157  // unsigned long numberOfEvents = storeFileMetaData->getNEvents();
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);
163 
164  VxdID currentSensorID(0);
165  double currentSensorThickness(0);
166  double currentSensorMass(0);
167  double currentSensorArea(0);
168 
169  // Exposition and dose
170  if (m_doseReportingLevel > c_reportNone) {
171  B2DEBUG(100, "Expo and dose");
172  currentSensorID.setID(0);
173  for (const PXDSimHit& hit : storeSimHits) {
174  // Update if we have a new sensor
175  VxdID sensorID = hit.getSensorID();
176  if (sensorID != currentSensorID) {
177  currentSensorID = sensorID;
178  currentSensorThickness = getSensorThickness(currentSensorID);
179  currentSensorMass = getSensorMass(currentSensorID);
180  currentSensorArea = getSensorArea(currentSensorID);
181  }
182  double hitEnergy = hit.getElectrons() * Const::ehEnergy;
183  // Dose in Gy/smy, normalize by sensor mass
184  m_sensorData[currentSensorID].m_dose +=
185  (hitEnergy / Unit::J) / (currentSensorMass / 1000) * (c_smy / currentComponentTime);
186  // Exposition in GeV/cm2/s
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);
193  storeEnergyDeposits.appendNew(
194  sensorID.getLayerNumber(), sensorID.getLadderNumber(), sensorID.getSensorNumber(),
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)
199  );
200  }
201  }
202  }
203 
204  // Neutron flux
205  if (m_nfluxReportingLevel > c_reportNone) {
206  B2DEBUG(100, "Neutron flux");
207  currentSensorID.setID(0);
208  for (const PXDTrueHit& hit : storeTrueHits) {
209  VxdID sensorID = hit.getSensorID();
210  // Update if we are on a new sensor
211  if (sensorID != currentSensorID) {
212  currentSensorID = sensorID;
213  currentSensorThickness = getSensorThickness(currentSensorID);
214  currentSensorMass = getSensorMass(currentSensorID);
215  currentSensorArea = getSensorArea(currentSensorID);
216  }
217  // J(TrueHit) = abs(step)/thickness * correctionFactor;
218  TVector3 entryPos(hit.getEntryU(), hit.getEntryV(), hit.getEntryW());
219  TVector3 exitPos(hit.getExitU(), hit.getExitV(), hit.getExitW());
220  double stepLength = (exitPos - entryPos).Mag();
221  // Identify what particle we've got. We need type and kinetic energy.
222  // TODO: TrueHit must carry pdg or SimHit must carry energy.
223  // NOTE: MCParticles may get remapped, then SimHits still carry correct pdg.
224  const PXDSimHit* simhit = hit.getRelatedTo<PXDSimHit>();
225  if (!simhit) { //either something is very wrong, or we just don't have the relation. Try to find an appropriate SimHit manually.
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;
231  simhit = &related;
232  }
233  }
234  }
235  // FIXME: Is there a difference between positrons and electrons wrt. NIEL?
236  // We fill neutronFluxBars with summary NIEL deposit for all kinds of particles by layer and component.
237  // Fluency plots are by component and are deposition histograms for a particular type of particle and compoonent.
238  // Special treatment of corrupt p's in TrueHits:
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);
243  int pdg = abs(simhit->getPDGcode());
244  double kineticEnergy(0.0);
245  double nielWeight(0.0);
246  if (pdg == 2112) {
247  double m0 = 0.940;
248  kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
249  nielWeight = m_nielNeutrons->getNielFactor(kineticEnergy / Unit::MeV);
250  }
251  if (pdg == 2212) {
252  double m0 = 0.938;
253  kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
254  nielWeight = m_nielProtons->getNielFactor(kineticEnergy / Unit::MeV);
255  }
256  if (pdg == 111 || pdg == 211) {
257  double m0 = 0.135;
258  kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
259  nielWeight = m_nielPions->getNielFactor(kineticEnergy / Unit::MeV);
260  }
261  if (pdg == 11) {
262  double m0 = 0.511e-3;
263  kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
264  nielWeight = m_nielElectrons->getNielFactor(kineticEnergy / Unit::MeV);
265  }
266  if (pdg == 22) {
267  double m0 = 0.0;
268  kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
269  }
270 
271  // Only set weight for supported particles
272  nielWeight = std::isfinite(nielWeight) ? nielWeight : 0.0;
273  m_sensorData[currentSensorID].m_neutronFlux += nielWeight * stepLength / currentSensorThickness / currentSensorArea /
274  currentComponentTime * c_smy;
275 
276  // Store data in a PXDNeutronFluxEvent object
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);
286  storeNeutronFluxes.appendNew(
287  sensorID.getLayerNumber(), sensorID.getLadderNumber(), sensorID.getSensorNumber(),
288  simhit->getPDGcode(), simhit->getGlobalTime(),
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)
293  );
294  }
295  }
296  }
297 
298  // Occupancy
299  if (m_occupancyReportingLevel > c_reportNone) {
300  B2DEBUG(100, "Fired pixels");
301  currentSensorID.setID(0);
302  double currentSensorCut = 0;
303  // Store fired pixels: count number of digits over threshold
304  std::map<VxdID, std::vector<float> > firedPixels;
305  for (const PXDDigit& storeDigit : storeDigits) {
306  // Filter out digits with signals below zero-suppression threshold
307  // ARE THRE SUCH DIGITS?
308  VxdID sensorID = storeDigit.getSensorID();
309  if (sensorID != currentSensorID) {
310  currentSensorID = sensorID;
311  auto info = getInfo(sensorID);
312  currentSensorCut = info.getChargeThreshold();
313  }
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());
318  }
319  // Process the map
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;
326  }
327 
328  B2DEBUG(100, "Occupancy");
329  currentSensorID.setID(0);
330  int nPixels = 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();
337  }
338 
339  double w_acceptance = m_integrationTime / currentComponentTime;
340  double occupancy = 1.0 / nPixels * cluster.getSize();
341  m_sensorData[sensorID].m_occupancy += w_acceptance * occupancy;
342 
343  if (m_occupancyReportingLevel == c_reportNTuple) {
344  storeOccupancyEvents.appendNew(
345  sensorID.getLayerNumber(), sensorID.getLadderNumber(),
346  sensorID.getSensorNumber(),
347  cluster.getU(), cluster.getV(), cluster.getSize(),
348  cluster.getCharge(), occupancy
349  );
350  }
351  }
352  }
353 }
354 
355 void PXDBackgroundModule::endRun()
356 {
357 }
358 
359 
360 void PXDBackgroundModule::terminate()
361 {
362  // Write out m_data
363  ofstream outfile;
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"
368  << "layer\t"
369  << "ladder\t"
370  << "sensor\t"
371  << "dose\t"
372  << "expo\t"
373  << "neutronFlux\t"
374  << "fired\t"
375  << "occupancy"
376  << endl;
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
389  << endl;
390  }
391  outfile << endl;
392 }
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:256
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:72
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:43
Belle2::PXDTrueHit
Class PXDTrueHit - Records of tracks that either enter or leave the sensitive volume.
Definition: PXDTrueHit.h:35
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::VxdID::getLadderNumber
baseType getLadderNumber() const
Get the ladder id.
Definition: VxdID.h:108
Belle2::StoreAccessorBase::getName
const std::string & getName() const
Return name under which the object is saved in the DataStore.
Definition: StoreAccessorBase.h:130
Belle2::PXDDigit
The PXD digit class.
Definition: PXDDigit.h:38
TNiel
TNiel - the class providing values for NIEL factors.
Definition: niel_fun.h:10
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::PXD::SensorInfo
Specific implementation of SensorInfo for PXD Sensors which provides additional pixel specific inform...
Definition: SensorInfo.h:34
Belle2::PXD
Namespace to encapsulate code needed for simulation and reconstrucion of the PXD.
Definition: PXDCalibrationUtilities.h:28
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Belle2::VXDSimHit::getGlobalTime
float getGlobalTime() const override
Return the time of the electron deposition.
Definition: VXDSimHit.h:80
Belle2::PXD::PXDBackgroundModule
PXD Background module.
Definition: PXDBackgroundModule.h:56
Belle2::VxdID::getSensorNumber
baseType getSensorNumber() const
Get the sensor id.
Definition: VxdID.h:110
Belle2::StoreArray< MCParticle >
Belle2::VxdID::getLayerNumber
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:106
Belle2::VXDSimHit::getPDGcode
int getPDGcode() const
Return the PDG code of the particle causing the electron deposition.
Definition: VXDSimHit.h:70
Belle2::VxdID::setID
void setID(baseType id)
Set the unique id.
Definition: VxdID.h:115
Belle2::PXDSimHit
Class PXDSimHit - Geant4 simulated hit for the PXD.
Definition: PXDSimHit.h:28