Belle II Software  release-05-01-25
PXDDigitizerModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2011 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Martin Ritter, Benjamin Schwenker *
7  * *
8  **************************************************************************/
9 
10 #include <pxd/modules/pxdSimulation/PXDDigitizerModule.h>
11 #include <vxd/geometry/GeoCache.h>
12 #include <vxd/geometry/GeoTools.h>
13 #include <pxd/reconstruction/PXDGainCalibrator.h>
14 #include <pxd/reconstruction/PXDPixelMasker.h>
15 
16 #include <framework/logging/Logger.h>
17 #include <framework/gearbox/Unit.h>
18 #include <framework/datastore/DataStore.h>
19 #include <framework/datastore/StoreArray.h>
20 #include <framework/datastore/RelationArray.h>
21 #include <framework/datastore/RelationIndex.h>
22 
23 #include <mdst/dataobjects/MCParticle.h>
24 #include <pxd/dataobjects/PXDTrueHit.h>
25 #include <pxd/dataobjects/PXDDigit.h>
26 #include <cmath>
27 
28 #include <TRandom.h>
29 
30 using namespace std;
31 using namespace Belle2;
32 using namespace Belle2::PXD;
33 
34 //-----------------------------------------------------------------
35 // Register the Module
36 //-----------------------------------------------------------------
37 REG_MODULE(PXDDigitizer)
38 
39 //-----------------------------------------------------------------
40 // Implementation
41 //-----------------------------------------------------------------
42 
44  , m_noiseFraction(0), m_eToADU(0), m_chargeThreshold(0), m_chargeThresholdElectrons(0)
45  , m_currentHit(nullptr), m_currentParticle(0), m_currentTrueHit(0), m_currentSensor(nullptr)
46  , m_currentSensorInfo(nullptr), m_nGates(0), m_pxdIntegrationTime(0), m_timePerGate(0)
47  , m_triggerGate(0), m_gated(false)
48 {
49  //Set module properties
50  setDescription("Digitize PXDSimHits");
51  setPropertyFlags(c_ParallelProcessingCertified);
52  addParam("ElectronicEffects", m_applyNoise, "Apply electronic effects?",
53  true);
54  addParam("ElectronicNoise", m_elNoise,
55  "Noise added by the electronics, set in ENC", 150.0);
56 
57  addParam("MCParticles", m_storeMCParticlesName,
58  "MCParticle collection name", string(""));
59  addParam("Digits", m_storeDigitsName, "Digits collection name", string(""));
60  addParam("SimHits", m_storeSimHitsName, "SimHit collection name",
61  string(""));
62  addParam("TrueHits", m_storeTrueHitsName, "TrueHit collection name",
63  string(""));
64 
65  addParam("PoissonSmearing", m_applyPoisson,
66  "Apply Poisson smearing of electrons collected on pixels?", true);
67  addParam("IntegrationWindow", m_applyWindow, "Use integration window?",
68  true);
69 
70  addParam("SegmentLength", m_segmentLength, "Maximum segment length (in mm)",
71  0.005);
72  addParam("ElectronGroupSize", m_elGroupSize,
73  "Split Signalpoints in smaller groups of N electrons (in e)", 100);
74  addParam("ElectronStepTime", m_elStepTime,
75  "Time step for tracking electron groups in readout plane (in ns)",
76  1.0);
77  addParam("ElectronMaxSteps", m_elMaxSteps,
78  "Maximum number of steps when propagating electrons", 200);
79  addParam("Gq", m_gq, "Gq of a pixel in nA/electron", 0.6);
80  addParam("ADCUnit", m_ADCUnit, "Slope of the linear ADC transfer curve in nA/ADU", 130.0);
81  addParam("PedestalMean", m_pedestalMean, "Mean of pedestals in ADU", 100.0);
82  addParam("PedestalRMS", m_pedestalRMS, "RMS of pedestals in ADU", 30.0);
83  addParam("GatingTime", m_gatingTime,
84  "Time window during which the PXD is not collecting charge in nano seconds", 400.0);
85  addParam("GatingWithoutReadout", m_gatingWithoutReadout,
86  "Digits from gated rows not sent to DHH ", true);
87  addParam("GatingWithoutReadoutTime", m_gatingWithoutReadoutTime,
88  "Time window during which digits from gated rows are not sent to DHH in nano seconds", 1400.0);
89  addParam("HardwareDelay", m_hwdelay,
90  "Constant time delay between bunch crossing and switching on triggergate in nano seconds", 0.0);
91 
92 }
93 
94 void PXDDigitizerModule::initialize()
95 {
96  //Register output collections
97  StoreArray<PXDDigit> storeDigits(m_storeDigitsName);
98  storeDigits.registerInDataStore();
99  StoreArray<MCParticle> storeParticles(m_storeMCParticlesName);
100  storeDigits.registerRelationTo(storeParticles);
101  StoreArray<PXDTrueHit> storeTrueHits(m_storeTrueHitsName);
102  storeDigits.registerRelationTo(storeTrueHits);
103 
104  m_storePXDTiming.isOptional();
105 
106  //Set default names for the relations
107  m_relMCParticleSimHitName = DataStore::relationName(
108  DataStore::arrayName<MCParticle>(m_storeMCParticlesName),
109  DataStore::arrayName<PXDSimHit>(m_storeSimHitsName));
110  m_relTrueHitSimHitName = DataStore::relationName(
111  DataStore::arrayName<PXDTrueHit>(m_storeTrueHitsName),
112  DataStore::arrayName<PXDSimHit>(m_storeSimHitsName));
113  m_relDigitMCParticleName = DataStore::relationName(
114  DataStore::arrayName<PXDDigit>(m_storeDigitsName),
115  DataStore::arrayName<MCParticle>(m_storeMCParticlesName));
116  m_relDigitTrueHitName = DataStore::relationName(
117  DataStore::arrayName<PXDDigit>(m_storeDigitsName),
118  DataStore::arrayName<PXDTrueHit>(m_storeTrueHitsName));
119 
120  //Convert parameters to correct units
121  m_elNoise *= Unit::e;
122  m_eToADU = m_ADCUnit / m_gq;
123  m_elStepTime *= Unit::ns;
124  m_segmentLength *= Unit::mm;
125 
126  // Get number of PXD gates
127  auto gTools = VXD::GeoCache::getInstance().getGeoTools();
128  m_nGates = gTools->getNumberOfPXDReadoutGates();
129 
130  // Active integration time (per gate) in ns
131  m_pxdIntegrationTime = 20000.0;
132  // Readout time per gate (sample - clear - cycle of rolling shutter) in ns
133  m_timePerGate = m_pxdIntegrationTime / m_nGates;
134 
135  B2DEBUG(20, "PXDDigitizer Parameters (in system units, *=calculated +=set in xml):");
136  B2DEBUG(20, " --> ElectronicEffects: " << (m_applyNoise ? "true" : "false"));
137  B2DEBUG(20, " --> ElectronicNoise: " << m_elNoise << " e-");
138  B2DEBUG(20, " --> +ChargeThreshold: " << "set in xml by sensor, nominal 4 ADU");
139  B2DEBUG(20, " --> *NoiseFraction: " << "set in xml by sensor, nominal 1.0e-5");
140  B2DEBUG(20,
141  " --> MCParticles: " << DataStore::arrayName<MCParticle>(m_storeMCParticlesName));
142  B2DEBUG(20,
143  " --> SimHits: " << DataStore::arrayName<PXDSimHit>(m_storeSimHitsName));
144  B2DEBUG(20,
145  " --> Digits: " << DataStore::arrayName<PXDDigit>(m_storeDigitsName));
146  B2DEBUG(20,
147  " --> TrueHits: " << DataStore::arrayName<PXDTrueHit>(m_storeTrueHitsName));
148  B2DEBUG(20, " --> MCSimHitRel: " << m_relMCParticleSimHitName);
149  B2DEBUG(20, " --> DigitMCRel: " << m_relDigitMCParticleName);
150  B2DEBUG(20, " --> TrueSimRel: " << m_relTrueHitSimHitName);
151  B2DEBUG(20, " --> DigitTrueRel: " << m_relDigitTrueHitName);
152  B2DEBUG(20, " --> PoissonSmearing: " << (m_applyPoisson ? "true" : "false"));
153  B2DEBUG(20, " --> +IntegrationWindow: " << (m_applyWindow ? "true" : "false") << ", size defined in xml");
154  B2DEBUG(20, " --> SegmentLength: " << m_segmentLength << " cm");
155  B2DEBUG(20, " --> ElectronGroupSize: " << m_elGroupSize << " e-");
156  B2DEBUG(20, " --> ElectronStepTime: " << m_elStepTime << " ns");
157  B2DEBUG(20, " --> ElectronMaxSteps: " << m_elMaxSteps);
158  B2DEBUG(20, " --> ADU unit: " << m_eToADU << " e-/ADU");
159 }
160 
161 void PXDDigitizerModule::beginRun()
162 {
163  //Fill map with all possible sensors This is too slow to be done every event so
164  //we fill it once and only clear the content of the sensors per event, not
165  //the whole map.
166  // NOTE: Some VXDIDs will be added to the map on the way if we have multiple
167  // frames, but they will stay and will be cleared appropriately, so this is not
168  // a problem, and after a few events the performance will stabilize.
169  m_sensors.clear();
170  VXD::GeoCache& geo = VXD::GeoCache::getInstance();
171  for (VxdID layer : geo.getLayers(SensorInfo::PXD)) {
172  for (VxdID ladder : geo.getLadders(layer)) {
173  for (VxdID sensor : geo.getSensors(ladder)) {
174  m_sensors[sensor] = Sensor();
175  }
176  }
177  }
178 }
179 
180 void PXDDigitizerModule::event()
181 {
182 
183  //Check for injection vetos
184  if (m_storePXDTiming.isValid()) {
185  m_gated = (*m_storePXDTiming).isGated();
186  m_triggerGate = (*m_storePXDTiming).getTriggerGate();
187  m_gatingStartTimes = (*m_storePXDTiming).getGatingStartTimes();
188  m_gatedChannelIntervals.clear();
189 
190  B2DEBUG(20, "Reading from PXDInjectionBGTiming: Event has TriggerGate " << m_triggerGate << " and is gated " << m_gated);
191 
192  if (m_gatingWithoutReadout) {
193  for (auto gatingStartTime : m_gatingStartTimes) {
194  B2DEBUG(20, "Gating start time " << gatingStartTime << "ns ");
195 
196  // Gated readout can only come from gating in the time period from t=0us to t=20us.
197  if (gatingStartTime >= 0) {
198  // Compute the gate where gating started. Note the wrap around of the rolling shutter
199  int firstGated = (m_triggerGate + int(gatingStartTime / m_timePerGate)) % m_nGates ;
200  // Compute the gate where gating stopped. Note the wrap around of the rolling shutter
201  int lastGated = (m_triggerGate + int((gatingStartTime + m_gatingWithoutReadoutTime) / m_timePerGate)) % m_nGates ;
202 
203  if (lastGated >= firstGated) {
204  // Gated channels in the same frame
205  B2DEBUG(20, "Gated channel interval from " << firstGated << " to " << lastGated);
206  m_gatedChannelIntervals.push_back(pair<int, int>(firstGated, lastGated));
207  } else {
208  // Gated channels in two consecutive frames
209  B2DEBUG(20, "Gated channel interval from " << firstGated << " to " << m_nGates);
210  B2DEBUG(20, "Gated channel interval from " << 0 << " to " << lastGated);
211  m_gatedChannelIntervals.push_back(pair<int, int>(firstGated, m_nGates));
212  m_gatedChannelIntervals.push_back(pair<int, int>(0, lastGated));
213  }
214  }
215  }
216  }
217  } else {
218  m_gated = false;
219  m_triggerGate = gRandom->Integer(m_nGates);
220  m_gatingStartTimes.clear();
221  m_gatedChannelIntervals.clear();
222  }
223 
224  //Clear sensors and process SimHits
225  for (Sensors::value_type& sensor : m_sensors) {
226  sensor.second.clear();
227  }
228  m_currentSensor = 0;
229  m_currentSensorInfo = 0;
230  //Clear return value
231  setReturnValue(0);
232 
233  StoreArray<MCParticle> storeMCParticles(m_storeMCParticlesName);
234  StoreArray<PXDSimHit> storeSimHits(m_storeSimHitsName);
235  StoreArray<PXDTrueHit> storeTrueHits(m_storeTrueHitsName);
236 
237  RelationArray mcParticlesToSimHits(storeMCParticles, storeSimHits,
238  m_relMCParticleSimHitName);
239  RelationArray mcParticlesToTrueHits(storeMCParticles, storeTrueHits); // not used here
240  RelationArray trueHitsToSimHits(storeTrueHits, storeSimHits,
241  m_relTrueHitSimHitName);
242 
243  StoreArray<PXDDigit> storeDigits(m_storeDigitsName);
244  storeDigits.clear();
245 
246  RelationArray relDigitMCParticle(storeDigits, storeMCParticles,
247  m_relDigitMCParticleName);
248  relDigitMCParticle.clear();
249 
250  RelationArray relDigitTrueHit(storeDigits, storeTrueHits,
251  m_relDigitTrueHitName);
252  relDigitTrueHit.clear();
253 
254  unsigned int nSimHits = storeSimHits.getEntries();
255  if (nSimHits == 0)
256  return;
257 
258  RelationIndex<MCParticle, PXDSimHit> relMCParticleSimHit(storeMCParticles,
259  storeSimHits, m_relMCParticleSimHitName);
260  RelationIndex<PXDTrueHit, PXDSimHit> relTrueHitSimHit(storeTrueHits,
261  storeSimHits, m_relTrueHitSimHitName);
262 
263  //Check sensor info and set pointers to current sensor
264  for (unsigned int i = 0; i < nSimHits; ++i) {
265  m_currentHit = storeSimHits[i];
267  relMCParticleSimHit.getFirstElementTo(m_currentHit);
268  if (mcRel) {
269  m_currentParticle = mcRel->indexFrom;
270  if (mcRel->weight < 0) {
271  //This simhit is from a particle which was not saved by the simulation
272  //so we do not take it into account for relations. Otherwise we might
273  //end up adding positive and negative weights together
274  m_currentParticle = -1;
275  }
276  } else {
277  // Don't warn if this is a background SimHit
278  if (m_currentHit->getBackgroundTag() == BackgroundMetaData::bg_none)
279  B2WARNING(
280  "Could not find MCParticle which produced PXDSimhit " << i);
281  m_currentParticle = -1;
282  }
284  relTrueHitSimHit.getFirstElementTo(m_currentHit);
285  //We only care about true hits from particles which have not been ignored
286  if (trueRel && trueRel->weight > 0) {
287  m_currentTrueHit = trueRel->indexFrom;
288  } else {
289  m_currentTrueHit = -1;
290  }
291 
292  VxdID sensorID = m_currentHit->getSensorID();
293  if (!m_currentSensorInfo || sensorID != m_currentSensorInfo->getID()) {
294  m_currentSensorInfo =
295  dynamic_cast<const SensorInfo*>(&VXD::GeoCache::get(
296  sensorID));
297  if (!m_currentSensorInfo)
298  B2FATAL(
299  "SensorInformation for Sensor " << sensorID << " not found, make sure that the geometry is set up correctly");
300  m_currentSensor = &m_sensors[sensorID];
301  B2DEBUG(20, "Sensor Parameters for Sensor " << sensorID << ": " << endl
302  << " --> Width: " << m_currentSensorInfo->getWidth() << endl
303  << " --> Length: " << m_currentSensorInfo->getLength() << endl
304  << " --> uPitch: " << m_currentSensorInfo->getUPitch() << endl
305  << " --> vPitch " << m_currentSensorInfo->getVPitch(-m_currentSensorInfo->getLength() / 2.0)
306  << ", " << m_currentSensorInfo->getVPitch(m_currentSensorInfo->getLength() / 2.0) << endl
307  << " --> Thickness: " << m_currentSensorInfo->getThickness() << endl
308  << " --> BulkDoping: " << m_currentSensorInfo->getBulkDoping() << endl
309  << " --> BackVoltage: " << m_currentSensorInfo->getBackVoltage() << endl
310  << " --> TopVoltage: " << m_currentSensorInfo->getTopVoltage() << endl
311  << " --> SourceBorder: " << m_currentSensorInfo->getSourceBorder(-0.4 * m_currentSensorInfo->getLength())
312  << ", " << m_currentSensorInfo->getSourceBorder(+0.4 * m_currentSensorInfo->getLength()) << endl
313  << " --> ClearBorder: " << m_currentSensorInfo->getClearBorder(-0.4 * m_currentSensorInfo->getLength())
314  << ", " << m_currentSensorInfo->getClearBorder(+0.4 * m_currentSensorInfo->getLength()) << endl
315  << " --> DrainBorder: " << m_currentSensorInfo->getDrainBorder(-0.4 * m_currentSensorInfo->getLength())
316  << ", " << m_currentSensorInfo->getDrainBorder(+0.4 * m_currentSensorInfo->getLength()) << endl
317  << " --> GateDepth: " << m_currentSensorInfo->getGateDepth() << endl
318  << " --> DoublePixel: " << m_currentSensorInfo->getDoublePixel() << endl);
319 
320  }
321  B2DEBUG(10,
322  "Processing hit " << i << " in Sensor " << sensorID << ", related to MCParticle " << m_currentParticle);
323  processHit();
324  }
325 
326  addNoiseDigits();
327  saveDigits();
328  //Return number of created digits
329  setReturnValue(storeDigits.getEntries());
330 }
331 
332 void PXDDigitizerModule::processHit()
333 {
334  if (m_applyWindow) {
335  //Ignore a hit which is outside of the active frame time of the hit gate
336  float currentHitTime = m_currentHit->getGlobalTime();
337 
338  // First we have to now the current hit gate
339  const TVector3 startPoint = m_currentHit->getPosIn();
340  const TVector3 stopPoint = m_currentHit->getPosOut();
341  auto row = m_currentSensorInfo->getVCellID(0.5 * (stopPoint.y() + startPoint.y()));
342 
343  if (m_currentSensorInfo->getID().getLayerNumber() == 1)
344  row = 767 - row;
345  int gate = row / 4;
346 
347  // Compute the distance in unit of gates to the trigger gate. Do it in
348  // direction of PXD rolling shutter.
349  int distanceToTriggerGate = 0;
350  if (gate >= m_triggerGate) distanceToTriggerGate = gate - m_triggerGate;
351  else distanceToTriggerGate = 192 - m_triggerGate + gate;
352 
353  // The triggergate is switched ON at t0=0ns. Compute the integration time window for the current hit gate.
354  // Time intervals of subsequent readout gates are shifted by timePerGate.
355  float gateMinTime = -m_pxdIntegrationTime + m_timePerGate + m_hwdelay + distanceToTriggerGate * m_timePerGate ;
356  float gateMaxTime = gateMinTime + m_pxdIntegrationTime;
357 
358  B2DEBUG(30,
359  "Checking if hit is in timeframe " << gateMinTime << " <= " << currentHitTime <<
360  " <= " << gateMaxTime);
361 
362  if (currentHitTime
363  < gateMinTime
364  || currentHitTime
365  > gateMaxTime)
366  return;
367 
368  if (m_gated) {
369  //Ignore simHits when the PXD is gated, i.e. PXD does not collect any charge.
370  for (auto gatingStartTime : m_gatingStartTimes) {
371  B2DEBUG(30,
372  "Checking if hit is in gated timeframe " << gatingStartTime << " <= " << currentHitTime <<
373  " <= " << gatingStartTime + m_gatingTime);
374 
375  if (currentHitTime
376  > gatingStartTime
377  && currentHitTime
378  < gatingStartTime + m_gatingTime)
379  return;
380  }
381  }
382  }
383 
384  //Get step length and direction
385  const TVector3 startPoint = m_currentHit->getPosIn();
386  TVector3 stopPoint = m_currentHit->getPosOut();
387  double dx = stopPoint.x() - startPoint.x();
388  double dy = stopPoint.y() - startPoint.y();
389  double dz = stopPoint.z() - startPoint.z();
390  double trackLength2 = dx * dx + dy * dy + dz * dz;
391 
392  // Set magnetic field to save calls to getBField()
393  m_currentBField = m_currentSensorInfo->getBField(0.5 * (startPoint + stopPoint));
394 
395  if (m_currentHit->getPDGcode() == 22 || trackLength2 <= 0.01 * Unit::um * Unit::um) {
396  //Photons deposit the energy at the end of their step
397  driftCharge(stopPoint, m_currentHit->getElectrons());
398  } else {
399  //Otherwise, split into segments of (default) max. 5µm and
400  //drift the charges from the center of each segment
401  auto segments = m_currentHit->getElectronsConstantDistance(m_segmentLength);
402  double lastFraction {0};
403  double lastElectrons {0};
404 
405  for (auto& segment : segments) {
406  //Simhit returns step fraction and cumulative electrons. We want the
407  //center of these steps and electrons in this step
408  const double f = (segment.first + lastFraction) * 0.5;
409  const double e = segment.second - lastElectrons;
410  //Update last values
411  std::tie(lastFraction, lastElectrons) = segment;
412 
413  //And drift charge from that position
414  stopPoint.SetXYZ(startPoint.x() + f * dx, startPoint.y() + f * dy, startPoint.z() + f * dz);
415  driftCharge(stopPoint, e);
416  }
417  }
418 }
419 
420 inline double pol3(double x, const double* c)
421 {
422  return c[0] + x * (c[1] + x * (c[2] + x * c[3]));
423 };
424 
425 void PXDDigitizerModule::driftCharge(const TVector3& r, double electrons)
426 {
427  //Get references to current sensor/info for ease of use
428  const SensorInfo& info = *m_currentSensorInfo;
429  Sensor& sensor = *m_currentSensor;
430 
431  static const double cx[] = {0, 2.611995e-01, -1.948316e+01, 9.167064e+02};
432  static const double cu[] = {4.223817e-04, -1.041434e-03, 5.139596e-02, -1.180229e+00};
433  static const double cv[] = {4.085681e-04, -8.615459e-04, 5.706439e-02, -1.774319e+00};
434  static const double ct[] = {2.823743e-04, -1.138627e-06, 4.310888e-09, -1.559542e-11}; // temp range [250, 350] degrees with m_elStepTime = 1 ns
435 
436  double dz = 0.5 * info.getThickness() - info.getGateDepth() - r.Z(), adz = fabs(dz);
437 
438  double dx = fabs(m_currentBField.y()) > Unit::T ? copysign(pol3(adz, cx),
439  dz) : 0; // two configurations: with default B field (B~1.5T) and B=0T
440  double sigmaDrift_u = pol3(adz, cu);
441  double sigmaDrift_v = pol3(adz, cv);
442  double sigmaDiffus = pol3(info.getTemperature() - 300.0,
443  ct); // sqrt(2 * Const::uTherm * info.getElectronMobility(0) * m_elStepTime);
444 
445  int nGroups = (int)(electrons / m_elGroupSize) + 1;
446  double groupCharge = electrons / nGroups;
447  while (nGroups--) {
448  double du = gRandom->Gaus(0.0, sigmaDrift_u), dv = gRandom->Gaus(0.0, sigmaDrift_v);
449  double uPos = r.x() + du + dx;
450  double vPos = r.y() + dv;
451  int step = 0;
452  for (; step < m_elMaxSteps; ++step) {
453  int id = info.getTrappedID(uPos, vPos);
454  //Check if cloud inside of IG region
455  if (id >= 0) {
456  // Trapped in the internal gate region
457  sensor[Digit(id % 250, id / 250)].add(groupCharge, m_currentParticle, m_currentTrueHit);
458  break;
459  }
460  //Random walk with drift
461  dv = gRandom->Gaus(0.0, sigmaDiffus), du = gRandom->Gaus(0.0, sigmaDiffus);
462  uPos += du;
463  vPos += dv;
464  }
465  if (step == m_elMaxSteps) {
466  // cloud is not trapped but save it anyway into the closest pixel
467  int iu = info.getUCellID(uPos, vPos, 1), iv = info.getVCellID(vPos, 1);
468  sensor[Digit(iu, iv)].add(groupCharge, m_currentParticle, m_currentTrueHit);
469  }
470  }
471 }
472 
473 
474 double PXDDigitizerModule::addNoise(double charge)
475 {
476  if (charge <= 0) {
477  //Noise Pixel, add chargeThreshold electrons;
478  charge = 1.1 * m_chargeThresholdElectrons;
479  } else {
480  if (m_applyPoisson) {
481  // For big charge assume Gaussian distr.
482  if (charge > (1000. * Unit::e))
483  charge = gRandom->Gaus(charge, sqrt(charge));
484  else
485  // Otherwise Poisson distr.
486  charge = gRandom->Poisson(charge);
487  }
488  if (m_applyNoise) {
489  charge += gRandom->Gaus(0., m_elNoise);
490  }
491  }
492  return charge;
493 }
494 
495 void PXDDigitizerModule::addNoiseDigits()
496 {
497  if (!m_applyNoise)
498  return;
499 
500  for (Sensors::value_type& sensor : m_sensors) {
501  Sensor& s = sensor.second;
502 
503  //Calculate the number of pixels on an empty sensor which will exceed the noise cut
504  const SensorInfo& info =
505  dynamic_cast<const SensorInfo&>(VXD::GeoCache::get(sensor.first));
506  int nU = info.getUCells();
507  int nV = info.getVCells();
508  int nPixels = gRandom->Poisson(info.getNoiseFraction() * nU * nV);
509 
510  // Existing digits will get noise in PXDDIgitizer::AddNoise.
511  // Here we add zero charge to nPixels digits.
512  // In part, this will create some new zero-charge digits.
513  // If we happen to select an existing (fired) digit, it remains unchanged.
514  for (int i = 0; i < nPixels; ++i) {
515  Digit d(gRandom->Integer(nU), gRandom->Integer(nV));
516  s[d].add(0.0);
517  }
518  }
519 }
520 
521 bool PXDDigitizerModule::checkIfGated(int gate)
522 {
523  return std::any_of(m_gatedChannelIntervals.begin(), m_gatedChannelIntervals.end(),
524  [gate](auto interval) {return gate >= interval.first && gate < interval.second;});
525 }
526 
527 void PXDDigitizerModule::saveDigits()
528 {
529  StoreArray<MCParticle> storeMCParticles(m_storeMCParticlesName);
530  StoreArray<PXDTrueHit> storeTrueHits(m_storeTrueHitsName);
531  StoreArray<PXDDigit> storeDigits(m_storeDigitsName);
532  RelationArray relDigitMCParticle(storeDigits, storeMCParticles,
533  m_relDigitMCParticleName);
534  RelationArray relDigitTrueHit(storeDigits, storeTrueHits,
535  m_relDigitTrueHitName);
536 
537 
538  for (Sensors::value_type& sensor : m_sensors) {
539  auto sensorID = sensor.first;
540  const SensorInfo& info =
541  dynamic_cast<const SensorInfo&>(VXD::GeoCache::get(sensorID));
542  m_chargeThreshold = info.getChargeThreshold();
543  m_chargeThresholdElectrons = m_chargeThreshold * m_eToADU;
544  for (Sensor::value_type& digitAndValue : sensor.second) {
545  const Digit& d = digitAndValue.first;
546  const DigitValue& v = digitAndValue.second;
547 
548  //Add Noise where applicable
549  double charge = addNoise(v.charge());
550 
551  // Draw a pedestal value
552  double pedestal = std::max(gRandom->Gaus(m_pedestalMean, m_pedestalRMS), 0.0);
553 
554  // Get gain correction
555  double gain = PXDGainCalibrator::getInstance().getGainCorrection(sensorID, d.u(), d.v());
556 
557  // Add pedestal to charge
558  charge = round(charge * gain / m_eToADU + pedestal);
559 
560  // Clipping of ADC codes at 255
561  charge = std::min(charge, 255.0);
562 
563  // Subtraction of integer pedestal in DHP
564  charge = charge - round(pedestal);
565 
566  // Zero Suppression in DHP
567  if (charge <= m_chargeThreshold)
568  continue;
569 
570  // Check if the readout digits is coming from a gated row
571  if (m_gatingWithoutReadout) {
572  int row = d.v();
573  if (sensorID.getLayerNumber() == 1)
574  row = 767 - row;
575  int gate = row / 4;
576  if (checkIfGated(gate)) continue;
577  }
578 
579  // Check if the readout digit is coming from a masked or dead area
580  if (PXD::PXDPixelMasker::getInstance().pixelDead(sensorID, d.u(), d.v())
581  || !PXD::PXDPixelMasker::getInstance().pixelOK(sensorID, d.u(), d.v())) {
582  continue;
583  }
584 
585  //Add the digit to datastore
586  int digIndex = storeDigits.getEntries();
587  storeDigits.appendNew(
588  PXDDigit(sensorID, d.u(), d.v(), charge));
589 
590  //If the digit has any relations to MCParticles, add the Relation
591  if (v.particles().size() > 0) {
592  relDigitMCParticle.add(digIndex, v.particles().begin(),
593  v.particles().end());
594  }
595  //If the digit has any truehits to TrueHit, add the Relation
596  if (v.truehits().size() > 0) {
597  relDigitTrueHit.add(digIndex, v.truehits().begin(),
598  v.truehits().end());
599  }
600  }
601  }
602 }
603 
604 
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:256
Belle2::RelationArray::clear
void clear() override
Clear all elements from the relation.
Definition: RelationArray.h:279
Belle2::RelationIndexContainer::Element::weight
RelationElement::weight_type weight
weight of the relation.
Definition: RelationIndexContainer.h:82
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::PXD::DigitValue
Class representing the charge and particle contributions for one pixel.
Definition: PXDDigitizerModule.h:57
Belle2::StoreArray::registerRelationTo
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:150
Belle2::RelationIndex::getFirstElementTo
const Element * getFirstElementTo(const TO &to) const
Return a pointer to the first relation Element of the given object.
Definition: RelationIndex.h:222
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::PXD::Digit
Class to represent the coordinates of one pixel.
Definition: PXDDigitizerModule.h:37
Belle2::RelationIndexContainer::Element::indexFrom
RelationElement::index_type indexFrom
index of the element from which the relation points.
Definition: RelationIndexContainer.h:70
Belle2::RelationIndex
Provides access to fast ( O(log n) ) bi-directional lookups on a specified relation.
Definition: RelationIndex.h:84
Belle2::VxdID::getID
baseType getID() const
Get the unique id.
Definition: VxdID.h:104
Belle2::VXD::GeoCache::getLayers
const std::set< Belle2::VxdID > getLayers(SensorInfoBase::SensorType sensortype=SensorInfoBase::VXD)
Return a set of all known Layers.
Definition: GeoCache.cc:177
Belle2::PXDDigit
The PXD digit class.
Definition: PXDDigit.h:38
Belle2::VXD::GeoCache::getSensors
const std::set< Belle2::VxdID > & getSensors(Belle2::VxdID ladder) const
Return a set of all sensor IDs belonging to a given ladder.
Definition: GeoCache.cc:205
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::RelationIndexContainer::Element
Element type for the index.
Definition: RelationIndexContainer.h:62
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::StoreArray::clear
void clear() override
Delete all entries in this array.
Definition: StoreArray.h:217
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::RelationArray::add
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
Definition: RelationArray.h:291
Belle2::PXD::Sensor
std::map< Digit, DigitValue > Sensor
Map of all hits in one Sensor.
Definition: PXDDigitizerModule.h:89
Belle2::VXD::GeoCache
Class to faciliate easy access to sensor information of the VXD like coordinate transformations or pi...
Definition: GeoCache.h:41
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::VXD::GeoCache::getLadders
const std::set< Belle2::VxdID > & getLadders(Belle2::VxdID layer) const
Return a set of all ladder IDs belonging to a given layer.
Definition: GeoCache.cc:194
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
Belle2::PXD::PXDDigitizerModule
The PXD Digitizer module.
Definition: PXDDigitizerModule.h:101