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>
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>
23 #include <mdst/dataobjects/MCParticle.h>
24 #include <pxd/dataobjects/PXDTrueHit.h>
25 #include <pxd/dataobjects/PXDDigit.h>
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)
50 setDescription(
"Digitize PXDSimHits");
51 setPropertyFlags(c_ParallelProcessingCertified);
52 addParam(
"ElectronicEffects", m_applyNoise,
"Apply electronic effects?",
54 addParam(
"ElectronicNoise", m_elNoise,
55 "Noise added by the electronics, set in ENC", 150.0);
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",
62 addParam(
"TrueHits", m_storeTrueHitsName,
"TrueHit collection name",
65 addParam(
"PoissonSmearing", m_applyPoisson,
66 "Apply Poisson smearing of electrons collected on pixels?",
true);
67 addParam(
"IntegrationWindow", m_applyWindow,
"Use integration window?",
70 addParam(
"SegmentLength", m_segmentLength,
"Maximum segment length (in mm)",
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)",
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);
94 void PXDDigitizerModule::initialize()
98 storeDigits.registerInDataStore();
104 m_storePXDTiming.isOptional();
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));
121 m_elNoise *= Unit::e;
122 m_eToADU = m_ADCUnit / m_gq;
123 m_elStepTime *= Unit::ns;
124 m_segmentLength *= Unit::mm;
127 auto gTools = VXD::GeoCache::getInstance().getGeoTools();
128 m_nGates = gTools->getNumberOfPXDReadoutGates();
131 m_pxdIntegrationTime = 20000.0;
133 m_timePerGate = m_pxdIntegrationTime / m_nGates;
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");
141 " --> MCParticles: " << DataStore::arrayName<MCParticle>(m_storeMCParticlesName));
143 " --> SimHits: " << DataStore::arrayName<PXDSimHit>(m_storeSimHitsName));
145 " --> Digits: " << DataStore::arrayName<PXDDigit>(m_storeDigitsName));
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");
161 void PXDDigitizerModule::beginRun()
174 m_sensors[sensor] =
Sensor();
180 void PXDDigitizerModule::event()
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();
190 B2DEBUG(20,
"Reading from PXDInjectionBGTiming: Event has TriggerGate " << m_triggerGate <<
" and is gated " << m_gated);
192 if (m_gatingWithoutReadout) {
193 for (
auto gatingStartTime : m_gatingStartTimes) {
194 B2DEBUG(20,
"Gating start time " << gatingStartTime <<
"ns ");
197 if (gatingStartTime >= 0) {
199 int firstGated = (m_triggerGate + int(gatingStartTime / m_timePerGate)) % m_nGates ;
201 int lastGated = (m_triggerGate + int((gatingStartTime + m_gatingWithoutReadoutTime) / m_timePerGate)) % m_nGates ;
203 if (lastGated >= firstGated) {
205 B2DEBUG(20,
"Gated channel interval from " << firstGated <<
" to " << lastGated);
206 m_gatedChannelIntervals.push_back(pair<int, int>(firstGated, lastGated));
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));
219 m_triggerGate = gRandom->Integer(m_nGates);
220 m_gatingStartTimes.clear();
221 m_gatedChannelIntervals.clear();
225 for (Sensors::value_type& sensor : m_sensors) {
226 sensor.second.clear();
229 m_currentSensorInfo = 0;
237 RelationArray mcParticlesToSimHits(storeMCParticles, storeSimHits,
238 m_relMCParticleSimHitName);
239 RelationArray mcParticlesToTrueHits(storeMCParticles, storeTrueHits);
241 m_relTrueHitSimHitName);
246 RelationArray relDigitMCParticle(storeDigits, storeMCParticles,
247 m_relDigitMCParticleName);
248 relDigitMCParticle.
clear();
251 m_relDigitTrueHitName);
252 relDigitTrueHit.
clear();
254 unsigned int nSimHits = storeSimHits.
getEntries();
259 storeSimHits, m_relMCParticleSimHitName);
261 storeSimHits, m_relTrueHitSimHitName);
264 for (
unsigned int i = 0; i < nSimHits; ++i) {
265 m_currentHit = storeSimHits[i];
274 m_currentParticle = -1;
278 if (m_currentHit->getBackgroundTag() == BackgroundMetaData::bg_none)
280 "Could not find MCParticle which produced PXDSimhit " << i);
281 m_currentParticle = -1;
286 if (trueRel && trueRel->
weight > 0) {
289 m_currentTrueHit = -1;
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(
297 if (!m_currentSensorInfo)
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);
322 "Processing hit " << i <<
" in Sensor " << sensorID <<
", related to MCParticle " << m_currentParticle);
332 void PXDDigitizerModule::processHit()
336 float currentHitTime = m_currentHit->getGlobalTime();
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()));
343 if (m_currentSensorInfo->getID().getLayerNumber() == 1)
349 int distanceToTriggerGate = 0;
350 if (gate >= m_triggerGate) distanceToTriggerGate = gate - m_triggerGate;
351 else distanceToTriggerGate = 192 - m_triggerGate + gate;
355 float gateMinTime = -m_pxdIntegrationTime + m_timePerGate + m_hwdelay + distanceToTriggerGate * m_timePerGate ;
356 float gateMaxTime = gateMinTime + m_pxdIntegrationTime;
359 "Checking if hit is in timeframe " << gateMinTime <<
" <= " << currentHitTime <<
360 " <= " << gateMaxTime);
370 for (
auto gatingStartTime : m_gatingStartTimes) {
372 "Checking if hit is in gated timeframe " << gatingStartTime <<
" <= " << currentHitTime <<
373 " <= " << gatingStartTime + m_gatingTime);
378 < gatingStartTime + m_gatingTime)
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;
393 m_currentBField = m_currentSensorInfo->getBField(0.5 * (startPoint + stopPoint));
395 if (m_currentHit->getPDGcode() == 22 || trackLength2 <= 0.01 * Unit::um * Unit::um) {
397 driftCharge(stopPoint, m_currentHit->getElectrons());
401 auto segments = m_currentHit->getElectronsConstantDistance(m_segmentLength);
402 double lastFraction {0};
403 double lastElectrons {0};
405 for (
auto& segment : segments) {
408 const double f = (segment.first + lastFraction) * 0.5;
409 const double e = segment.second - lastElectrons;
411 std::tie(lastFraction, lastElectrons) = segment;
414 stopPoint.SetXYZ(startPoint.x() + f * dx, startPoint.y() + f * dy, startPoint.z() + f * dz);
415 driftCharge(stopPoint, e);
420 inline double pol3(
double x,
const double* c)
422 return c[0] + x * (c[1] + x * (c[2] + x * c[3]));
425 void PXDDigitizerModule::driftCharge(
const TVector3& r,
double electrons)
428 const SensorInfo& info = *m_currentSensorInfo;
429 Sensor& sensor = *m_currentSensor;
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};
436 double dz = 0.5 * info.getThickness() - info.getGateDepth() - r.Z(), adz = fabs(dz);
438 double dx = fabs(m_currentBField.y()) > Unit::T ? copysign(pol3(adz, cx),
440 double sigmaDrift_u = pol3(adz, cu);
441 double sigmaDrift_v = pol3(adz, cv);
442 double sigmaDiffus = pol3(info.getTemperature() - 300.0,
445 int nGroups = (int)(electrons / m_elGroupSize) + 1;
446 double groupCharge = electrons / 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;
452 for (; step < m_elMaxSteps; ++step) {
453 int id = info.getTrappedID(uPos, vPos);
457 sensor[
Digit(
id % 250,
id / 250)].add(groupCharge, m_currentParticle, m_currentTrueHit);
461 dv = gRandom->Gaus(0.0, sigmaDiffus), du = gRandom->Gaus(0.0, sigmaDiffus);
465 if (step == m_elMaxSteps) {
467 int iu = info.getUCellID(uPos, vPos, 1), iv = info.getVCellID(vPos, 1);
468 sensor[
Digit(iu, iv)].add(groupCharge, m_currentParticle, m_currentTrueHit);
474 double PXDDigitizerModule::addNoise(
double charge)
478 charge = 1.1 * m_chargeThresholdElectrons;
480 if (m_applyPoisson) {
482 if (charge > (1000. * Unit::e))
483 charge = gRandom->Gaus(charge, sqrt(charge));
486 charge = gRandom->Poisson(charge);
489 charge += gRandom->Gaus(0., m_elNoise);
495 void PXDDigitizerModule::addNoiseDigits()
500 for (Sensors::value_type& sensor : m_sensors) {
501 Sensor& s = sensor.second;
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);
514 for (
int i = 0; i < nPixels; ++i) {
515 Digit d(gRandom->Integer(nU), gRandom->Integer(nV));
521 bool PXDDigitizerModule::checkIfGated(
int gate)
523 return std::any_of(m_gatedChannelIntervals.begin(), m_gatedChannelIntervals.end(),
524 [gate](
auto interval) {return gate >= interval.first && gate < interval.second;});
527 void PXDDigitizerModule::saveDigits()
532 RelationArray relDigitMCParticle(storeDigits, storeMCParticles,
533 m_relDigitMCParticleName);
535 m_relDigitTrueHitName);
538 for (Sensors::value_type& sensor : m_sensors) {
539 auto sensorID = sensor.first;
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;
549 double charge = addNoise(v.charge());
552 double pedestal = std::max(gRandom->Gaus(m_pedestalMean, m_pedestalRMS), 0.0);
555 double gain = PXDGainCalibrator::getInstance().getGainCorrection(sensorID, d.u(), d.v());
558 charge = round(charge * gain / m_eToADU + pedestal);
561 charge = std::min(charge, 255.0);
564 charge = charge - round(pedestal);
567 if (charge <= m_chargeThreshold)
571 if (m_gatingWithoutReadout) {
573 if (sensorID.getLayerNumber() == 1)
576 if (checkIfGated(gate))
continue;
580 if (PXD::PXDPixelMasker::getInstance().pixelDead(sensorID, d.u(), d.v())
581 || !PXD::PXDPixelMasker::getInstance().pixelOK(sensorID, d.u(), d.v())) {
588 PXDDigit(sensorID, d.u(), d.v(), charge));
591 if (v.particles().size() > 0) {
592 relDigitMCParticle.
add(digIndex, v.particles().begin(),
593 v.particles().end());
596 if (v.truehits().size() > 0) {
597 relDigitTrueHit.
add(digIndex, v.truehits().begin(),