11 #include <cdc/modules/cdcDigitizer/CDCDigitizerModule.h>
12 #include <cdc/utilities/ClosestApproach.h>
14 #include <framework/datastore/RelationArray.h>
16 #include <framework/gearbox/Unit.h>
17 #include <framework/logging/Logger.h>
31 m_cdcgp(), m_gcp(), m_aCDCSimHit(), m_posFlag(0),
32 m_driftLength(0.0), m_flightTime(0.0), m_globalTime(0.0),
33 m_tdcBinWidth(1.0), m_tdcBinWidthInv(1.0),
34 m_tdcResol(0.9825), m_driftV(4.0e-3),
35 m_driftVInv(250.0), m_propSpeedInv(27.25), m_align(true)
38 setDescription(
"Creates CDCHits from CDCSimHits.");
39 setPropertyFlags(c_ParallelProcessingCertified);
43 addParam(
"InputCDCSimHitsName", m_inputCDCSimHitsName,
"Name of input array. Should consist of CDCSimHits.",
string(
""));
44 addParam(
"OutputCDCHitsName", m_outputCDCHitsName,
"Name of output array. Will consist of CDCHits.",
string(
""));
45 addParam(
"OutputCDCHitsName4Trg", m_outputCDCHitsName4Trg,
46 "Name of output array for trigger. Can contain several hits per wire, "
47 "if they correspond to different time windows of 32ns.",
48 string(
"CDCHits4Trg"));
51 addParam(
"MCParticlesToCDCSimHitsName", m_MCParticlesToSimHitsName,
52 "Name of relation between MCParticles and CDCSimHits used",
string(
""));
53 addParam(
"CDCSimHistToCDCHitsName", m_SimHitsTOCDCHitsName,
54 "Name of relation between the CDCSimHits and the CDCHits used",
string(
""));
58 addParam(
"UseSimpleDigitization", m_useSimpleDigitization,
59 "If true, a simple x-t with a constant velocity is used for the drift-length to -time conversion",
false);
62 addParam(
"Fraction", m_fraction,
"Fraction of first Gaussian used to smear drift length in cm", 1.0);
63 addParam(
"Mean1", m_mean1,
"Mean value of first Gaussian used to smear drift length in cm", 0.0000);
64 addParam(
"Resolution1", m_resolution1,
"Resolution of first Gaussian used to smear drift length in cm", 0.0130);
65 addParam(
"Mean2", m_mean2,
"Mean value of second Gaussian used to smear drift length in cm", 0.0000);
66 addParam(
"Resolution2", m_resolution2,
"Resolution of second Gaussian used to smear drift length in cm", 0.0000);
69 addParam(
"DoSmearing", m_doSmearing,
70 "If false, drift length will not be smeared.",
true);
72 addParam(
"TrigTimeJitter", m_trigTimeJitter,
73 "Magnitude (w) of trigger timing jitter (ns). The trigger timing is randuminzed uniformly in a time window of [-w/2, +w/2].",
76 addParam(
"AddTimeWalk", m_addTimeWalk,
"A switch for time-walk (pulse-heght dep. delay); true: on; false: off",
true);
77 addParam(
"AddInWirePropagationDelay", m_addInWirePropagationDelay,
78 "A switch used to control adding propagation delay in the wire into the final drift time or not; this is for signal hits.",
true);
79 addParam(
"AddInWirePropagationDelay4Bg", m_addInWirePropagationDelay4Bg,
80 "The same switch but for beam bg. hits.",
true);
81 addParam(
"AddTimeOfFlight", m_addTimeOfFlight,
82 "A switch used to control adding time of flight into the final drift time or not; this is for signal hits.",
true);
83 addParam(
"AddTimeOfFlight4Bg", m_addTimeOfFlight4Bg,
84 "The same switch but for beam bg. hits.",
true);
85 addParam(
"OutputNegativeDriftTime", m_outputNegativeDriftTime,
"Output hits with negative drift time",
true);
86 addParam(
"Output2ndHit", m_output2ndHit,
87 "Output the 2nd hit if exists in the time window. Note that it is not well-simulated at all, partly because no cross-talk betw. channels is simulated.",
90 addParam(
"CorrectForWireSag", m_correctForWireSag,
91 "A switch for sense wire sag effect; true: drift-time is calculated with the sag taken into account; false: not. Here, sag means the perturbative part which corresponds to alignment in case of wire-position. The main part (corresponding to design+displacement in wire-position) is taken into account in FullSim; you can control it via CDCJobCntlParModifier.",
95 addParam(
"TDCThreshold4Outer", m_tdcThreshold4Outer,
96 "TDC threshold (dE in eV) for Layers#8-56. The value corresponds to He-C2H6 gas; for the gas+wire (MaterialDefinitionMode=0) case, (this value)*f will be used, where f is specified by GasToGasWire",
98 addParam(
"TDCThreshold4Inner", m_tdcThreshold4Inner,
99 "Same as TDCThreshold4Outer but for Layers#0-7,", 25.0);
100 addParam(
"GasToGasWire", m_gasToGasWire,
101 "(Approximate) ratio of dE in He/C2H6-gas to dE in gas+wire, where dE is energy deposit.", 1. / 1.478);
104 addParam(
"ADCThreshold", m_adcThreshold,
105 "Threshold for ADC-count (in unit of count). ADC-count < threshold is treated as count=0.", 2);
106 addParam(
"tMin", m_tMin,
"Lower edge of time window in ns", -100.);
107 addParam(
"tMaxOuter", m_tMaxOuter,
"Upper edge of time window in ns for the normal-cell layers", 500.);
108 addParam(
"tMaxInner", m_tMaxInner,
"Upper edge of time window in ns for the small-cell layers", 300.);
115 addParam(
"UseDB4FEE", m_useDB4FEE,
"Fetch and use FEE params. from database or not",
true);
116 addParam(
"UseDB4EDepToADC", m_useDB4EDepToADC,
"Uuse edep-to-ADC conversion params. from database or not",
true);
117 addParam(
"UseDB4RunGain", m_useDB4RunGain,
"Fetch and use run gain from database or not",
true);
118 addParam(
"OverallGainFactor", m_overallGainFactor,
"Overall gain factor for adjustment", 1.0);
121 addParam(
"TDCThresholdOffset", m_tdcThresholdOffset,
"Offset for TDC (digital) threshold (mV)", 3820.);
122 addParam(
"AnalogGain", m_analogGain,
"Analog gain (V/pC)", 1.1);
123 addParam(
"DigitalGain", m_digitalGain,
"Digital gain (V/pC)", 15.);
124 addParam(
"ADCBinWidth", m_adcBinWidth,
"ADC bin width (mV)", 2.);
126 addParam(
"AddFudgeFactorForSigma", m_addFudgeFactorForSigma,
127 "Additional fudge factor for space resol. (common to all cells)", 1.);
128 addParam(
"SpaceChargeEffect", m_spaceChargeEffect,
"Switch for space charge effect",
true);
130 addParam(
"AddXTalk", m_addXTalk,
"A switch for crosstalk; true: on; false: off",
true);
131 addParam(
"Issue2ndHitWarning", m_issue2ndHitWarning,
"=true: issue a warning when a 2nd TDC hit is found.",
true);
132 addParam(
"IncludeEarlyXTalks", m_includeEarlyXTalks,
"=true: include earlier x-talks as well than the signal hit in question.",
134 addParam(
"DebugLevel4XTalk", m_debugLevel4XTalk,
"Debug level fr crosstalk; 20-29 are usable.", 20);
137 #if defined(CDC_DEBUG)
139 cout <<
"CDCDigitizer constructor" << endl;
143 void CDCDigitizerModule::initialize()
145 m_simHits.isRequired(m_inputCDCSimHitsName);
148 m_cdcHits.registerInDataStore(m_outputCDCHitsName);
149 m_simHits.registerRelationTo(m_cdcHits);
150 m_mcParticles.registerRelationTo(m_cdcHits);
152 m_cdcHits4Trg.registerInDataStore(m_outputCDCHitsName4Trg);
153 m_simHits.registerRelationTo(m_cdcHits4Trg);
154 m_mcParticles.registerRelationTo(m_cdcHits4Trg);
156 m_cdcgp = &(CDCGeometryPar::Instance());
159 m_tdcBinWidthInv = 1. / m_tdcBinWidth;
160 m_tdcResol = m_tdcBinWidth / sqrt(12.);
162 m_driftVInv = 1. / m_driftV;
165 if (m_cdcgp->getMaterialDefinitionMode() == 0) {
166 m_scaleFac = m_gasToGasWire;
168 m_scaleFac *= Unit::GeV;
169 m_gcp = &(CDCGeoControlPar::getInstance());
170 m_totalFudgeFactor = m_cdcgp->getFudgeFactorForSigma(2);
171 m_totalFudgeFactor *= m_addFudgeFactorForSigma;
172 B2DEBUG(29,
"totalFugeF in Digi= " << m_totalFudgeFactor);
183 if ((*m_fEElectronicsFromDB).isValid()) {
184 (*m_fEElectronicsFromDB).
addCallback(
this, &CDCDigitizerModule::setFEElectronics);
187 B2FATAL(
"CDCDigitizer:: CDCFEElectronics not valid !");
203 if (m_useDB4RunGain) {
205 if ((*m_runGainFromDB).isValid()) {
206 (*m_runGainFromDB).
addCallback(
this, &CDCDigitizerModule::setRunGain);
209 B2FATAL(
"CDCDedxRunGain invalid!");
212 B2DEBUG(29,
"run-gain = " << m_runGain);
216 if ((*m_xTalkFromDB).isValid()) {
218 B2FATAL(
"CDCCrossTalkLibrary invalid!");
222 #if defined(CDC_DEBUG)
224 cout <<
"CDCDigitizer initialize" << endl;
226 cout <<
"m_tdcBinWidth= " << m_tdcBinWidth << endl;
227 cout <<
"m_tdcResol= " << m_tdcResol << endl;
228 cout <<
"m_driftV= " << m_driftV << endl;
229 cout <<
"m_driftVInv= " << m_driftVInv << endl;
230 cout <<
"m_propSpeedInv= " << m_propSpeedInv << endl;
240 if (m_useDB4EDepToADC) {
241 if (m_cdcgp->getEDepToADCMainFactor(0, 0) == 0.) {
242 B2FATAL(
"CDCEDepToADCConversion payloads are unavailable!");
248 void CDCDigitizerModule::event()
251 RelationArray mcParticlesToCDCSimHits(m_mcParticles, m_simHits);
258 map<WireID, SignalInfo> signalMap;
259 map<WireID, SignalInfo>::iterator iterSignalMap;
261 map<WireID, unsigned short> adcMap;
262 map<WireID, unsigned short>::iterator iterADCMap;
267 map<pair<WireID, unsigned>,
SignalInfo> signalMapTrg;
268 map<pair<WireID, unsigned>,
SignalInfo>::iterator iterSignalMapTrg;
271 double trigTiming = m_trigTimeJitter == 0. ? 0. : m_trigTimeJitter * (gRandom->Uniform() - 0.5);
274 int nHits = m_simHits.getEntries();
275 B2DEBUG(29,
"Number of CDCSimHits in the current event: " << nHits);
276 for (
int iHits = 0; iHits < nHits; ++iHits) {
278 m_aCDCSimHit = m_simHits[iHits];
281 m_wireID = m_aCDCSimHit->getWireID();
284 m_posFlag = m_aCDCSimHit->getLeftRightPassageRaw();
285 m_boardID = m_cdcgp->getBoardID(m_wireID);
287 m_posWire = m_aCDCSimHit->getPosWire();
288 m_posTrack = m_aCDCSimHit->getPosTrack();
289 m_momentum = m_aCDCSimHit->getMomentum();
290 m_flightTime = m_aCDCSimHit->getFlightTime();
291 m_globalTime = m_aCDCSimHit->getGlobalTime();
292 m_driftLength = m_aCDCSimHit->getDriftLength() * Unit::cm;
298 TVector3 bwpAlign = m_cdcgp->wireBackwardPosition(m_wireID, CDCGeometryPar::c_Aligned);
299 TVector3 fwpAlign = m_cdcgp->wireForwardPosition(m_wireID, CDCGeometryPar::c_Aligned);
301 TVector3 bwp = m_cdcgp->wireBackwardPosition(m_wireID);
302 TVector3 fwp = m_cdcgp->wireForwardPosition(m_wireID);
305 if ((bwpAlign - bwp).Mag() == 0. && (fwpAlign - fwp).Mag() == 0.) m_align =
false;
308 if (m_align || m_correctForWireSag) {
313 if (m_correctForWireSag) {
314 double zpos = m_posWire.z();
315 double bckYSag = bwp.y();
316 double forYSag = fwp.y();
321 const int layerID = m_wireID.getICLayer();
322 const int wireID = m_wireID.getIWire();
323 m_cdcgp->getWireSagEffect(set, layerID, wireID, zpos, bckYSag, forYSag);
328 const TVector3 L = 5. * m_momentum.Unit();
329 TVector3 posIn = m_posTrack - L;
330 TVector3 posOut = m_posTrack + L;
331 TVector3 posTrack = m_posTrack;
332 TVector3 posWire = m_posWire;
335 m_driftLength =
ClosestApproach(bwp, fwp, posIn, posOut, posTrack, posWire);
337 m_posTrack = posTrack;
340 double deltaTime = 0.;
342 m_flightTime += deltaTime;
343 m_globalTime += deltaTime;
344 m_posFlag = m_cdcgp->getNewLeftRightRaw(m_posWire, m_posTrack, m_momentum);
349 double hitDriftLength = m_driftLength;
350 double dDdt = getdDdt(hitDriftLength);
352 hitDriftLength = smearDriftLength(hitDriftLength, dDdt);
356 bool addTof = m_addTimeOfFlight4Bg;
357 bool addDelay = m_addInWirePropagationDelay4Bg;
358 if (m_aCDCSimHit->getBackgroundTag() == 0) {
359 addTof = m_addTimeOfFlight;
360 addDelay = m_addInWirePropagationDelay;
362 double hitDriftTime = getDriftTime(hitDriftLength, addTof, addDelay);
365 if (m_aCDCSimHit->getBackgroundTag() != 0) {
366 hitDriftTime += m_globalTime - m_flightTime;
370 hitDriftTime += trigTiming;
373 double tMin = m_tMin;
374 double tMax = m_tMaxOuter;
375 if (m_wireID.getISuperLayer() == 0) tMax = m_tMaxInner;
377 tMin = m_lowEdgeOfTimeWindow[m_boardID];
378 tMax = m_uprEdgeOfTimeWindow[m_boardID];
380 if (hitDriftTime < tMin || hitDriftTime > tMax)
continue;
383 const double stepLength = m_aCDCSimHit->getStepLength() * Unit::cm;
384 const double costh = m_momentum.z() / m_momentum.Mag();
385 const double hitdE = m_scaleFac * m_aCDCSimHit->getEnergyDep();
387 unsigned short layerID = m_wireID.getICLayer();
388 unsigned short cellID = m_wireID.getIWire();
390 unsigned short adcCount = getADCCount(m_wireID, hitdE, stepLength, costh);
391 const unsigned short adcTh = m_useDB4FEE ? m_adcThresh[m_boardID] : m_adcThreshold;
393 if (adcCount < adcTh) adcCount = 0;
394 iterADCMap = adcMap.find(m_wireID);
395 if (iterADCMap == adcMap.end()) {
396 adcMap.insert(make_pair(m_wireID, adcCount));
399 iterADCMap->second += adcCount;
406 double dEThreshold = 0.;
407 if (m_useDB4FEE && m_useDB4EDepToADC) {
408 dEThreshold = (m_tdcThresh[m_boardID] / m_cdcgp->getEDepToADCMainFactor(layerID, cellID)) * Unit::keV;
410 dEThreshold = (m_wireID.getISuperLayer() == 0) ? m_tdcThreshold4Inner : m_tdcThreshold4Outer;
411 dEThreshold *= Unit::eV;
413 B2DEBUG(29,
"hitdE,dEThreshold,driftLength " << hitdE <<
" " << dEThreshold <<
" " << hitDriftLength);
415 if (hitdE < dEThreshold) {
416 B2DEBUG(29,
"Below Ethreshold: " << hitdE <<
" " << dEThreshold);
421 unsigned short trigWindow = floor((hitDriftTime - m_tMin) * m_tdcBinWidthInv / 32);
422 iterSignalMapTrg = signalMapTrg.find(make_pair(m_wireID, trigWindow));
423 if (iterSignalMapTrg == signalMapTrg.end()) {
426 signalMapTrg.insert(make_pair(make_pair(m_wireID, trigWindow),
429 if (hitDriftTime < iterSignalMapTrg->second.m_driftTime) {
430 iterSignalMapTrg->second.m_driftTime = hitDriftTime;
431 iterSignalMapTrg->second.m_simHitIndex = iHits;
434 iterSignalMapTrg->second.m_charge += adcCount;
438 if (m_cdcgp->isBadWire(m_wireID)) {
444 if (m_cdcgp->isDeadWire(m_wireID, eff)) {
446 if (eff < gRandom->Uniform())
continue;
450 const double a = bwpAlign.X();
451 const double b = bwpAlign.Y();
452 const double c = bwpAlign.Z();
453 const TVector3 fmbAlign = fwpAlign - bwpAlign;
454 const double lmn = 1. / fmbAlign.Mag();
455 const double l = fmbAlign.X() * lmn;
456 const double m = fmbAlign.Y() * lmn;
457 const double n = fmbAlign.Z() * lmn;
459 double dx = m_aCDCSimHit->getPosIn().X() - a;
460 double dy = m_aCDCSimHit->getPosIn().Y() - b;
461 double dz = m_aCDCSimHit->getPosIn().Z() - c;
462 double sub = l * dx + m * dy + n * dz;
463 const double driftLFromIn = sqrt(dx * dx + dy * dy + dz * dz - sub * sub);
465 dx = m_aCDCSimHit->getPosOut().X() - a;
466 dy = m_aCDCSimHit->getPosOut().Y() - b;
467 dz = m_aCDCSimHit->getPosOut().Z() - c;
468 sub = l * dx + m * dy + n * dz;
469 const double driftLFromOut = sqrt(dx * dx + dy * dy + dz * dz - sub * sub);
471 const double maxDriftL = std::max(driftLFromIn, driftLFromOut);
472 const double minDriftL = m_driftLength;
473 B2DEBUG(29,
"driftLFromIn= " << driftLFromIn <<
" driftLFromOut= " << driftLFromOut <<
" minDriftL= " << minDriftL <<
" maxDriftL= "
475 maxDriftL <<
"m_driftLength= " << m_driftLength);
477 iterSignalMap = signalMap.find(m_wireID);
479 if (iterSignalMap == signalMap.end()) {
482 signalMap.insert(make_pair(m_wireID,
SignalInfo(iHits, hitDriftTime, adcCount, maxDriftL, minDriftL)));
483 B2DEBUG(29,
"Creating new Signal with encoded wire number: " << m_wireID);
486 if (hitDriftTime < iterSignalMap->second.m_driftTime) {
487 iterSignalMap->second.m_driftTime3 = iterSignalMap->second.m_driftTime2;
488 iterSignalMap->second.m_simHitIndex3 = iterSignalMap->second.m_simHitIndex2;
489 iterSignalMap->second.m_driftTime2 = iterSignalMap->second.m_driftTime;
490 iterSignalMap->second.m_simHitIndex2 = iterSignalMap->second.m_simHitIndex;
491 iterSignalMap->second.m_driftTime = hitDriftTime;
492 iterSignalMap->second.m_simHitIndex = iHits;
493 B2DEBUG(29,
"hitDriftTime of current Signal: " << hitDriftTime <<
", hitDriftLength: " << hitDriftLength);
494 }
else if (hitDriftTime < iterSignalMap->second.m_driftTime2) {
495 iterSignalMap->second.m_driftTime3 = iterSignalMap->second.m_driftTime2;
496 iterSignalMap->second.m_simHitIndex3 = iterSignalMap->second.m_simHitIndex2;
497 iterSignalMap->second.m_driftTime2 = hitDriftTime;
498 iterSignalMap->second.m_simHitIndex2 = iHits;
499 }
else if (hitDriftTime < iterSignalMap->second.m_driftTime3) {
500 iterSignalMap->second.m_driftTime3 = hitDriftTime;
501 iterSignalMap->second.m_simHitIndex3 = iHits;
505 iterSignalMap->second.m_charge += adcCount;
508 if (iterSignalMap->second.m_maxDriftL < maxDriftL) iterSignalMap->second.m_maxDriftL = maxDriftL;
509 if (iterSignalMap->second.m_minDriftL > minDriftL) iterSignalMap->second.m_minDriftL = minDriftL;
510 B2DEBUG(29,
"maxDriftL in struct= " << iterSignalMap->second.m_maxDriftL <<
"minDriftL in struct= " <<
511 iterSignalMap->second.m_minDriftL);
519 unsigned int iCDCHits = 0;
521 RelationArray mcParticlesToCDCHits(m_mcParticles, m_cdcHits);
523 for (iterSignalMap = signalMap.begin(); iterSignalMap != signalMap.end(); ++iterSignalMap) {
528 iterADCMap = adcMap.find(iterSignalMap->first);
529 unsigned short adcCount = iterADCMap != adcMap.end() ? iterADCMap->second : 0;
542 B2DEBUG(29,
"timewalk= " << m_cdcgp->getTimeWalk(iterSignalMap->first, adcCount));
543 iterSignalMap->second.m_driftTime += m_cdcgp->getTimeWalk(iterSignalMap->first, adcCount);
547 if (!m_outputNegativeDriftTime &&
548 iterSignalMap->second.m_driftTime < 0.) {
553 unsigned short tdcCount =
static_cast<unsigned short>((m_cdcgp->getT0(iterSignalMap->first) - iterSignalMap->second.m_driftTime) *
554 m_tdcBinWidthInv + 0.5);
557 double deltaDL = iterSignalMap->second.m_maxDriftL - iterSignalMap->second.m_minDriftL;
559 B2DEBUG(29,
"negative deltaDL= " << deltaDL);
562 const unsigned short boardID = m_cdcgp->getBoardID(iterSignalMap->first);
563 unsigned short tot = std::min(std::round(5.92749 * deltaDL + 2.59706),
static_cast<double>(m_widthOfTimeWindow[boardID]));
564 if (m_adcThresh[boardID] > 0) {
565 tot = std::min(
static_cast<int>(tot),
static_cast<int>(adcCount / m_adcThresh[boardID]));
568 CDCHit* firstHit = m_cdcHits.appendNew(tdcCount, adcCount, iterSignalMap->first, 0, tot);
571 cdcSimHitsToCDCHits.
add(iterSignalMap->second.m_simHitIndex, iCDCHits);
575 if (rels.
size() != 0) {
578 double weight = rels.
weight(0);
583 if (m_output2ndHit && iterSignalMap->second.m_simHitIndex2 >= 0) {
584 unsigned short tdcCount2 =
static_cast<unsigned short>((m_cdcgp->getT0(iterSignalMap->first) - iterSignalMap->second.m_driftTime2) *
585 m_tdcBinWidthInv + 0.5);
586 if (tdcCount2 != tdcCount) {
587 CDCHit* secondHit = m_cdcHits.appendNew(tdcCount2, adcCount, iterSignalMap->first, 0, tot);
600 cdcSimHitsToCDCHits.
add(iterSignalMap->second.m_simHitIndex2, iCDCHits);
604 rels = m_simHits[iterSignalMap->second.m_simHitIndex2]->getRelationsFrom<
MCParticle>();
605 if (rels.
size() != 0) {
608 double weight = rels.
weight(0);
613 if (iterSignalMap->second.m_simHitIndex3 >= 0) {
614 unsigned short tdcCount3 =
static_cast<unsigned short>((m_cdcgp->getT0(iterSignalMap->first) - iterSignalMap->second.m_driftTime3) *
615 m_tdcBinWidthInv + 0.5);
617 if (tdcCount3 != tdcCount) {
618 CDCHit* secondHit = m_cdcHits.appendNew(tdcCount3, adcCount, iterSignalMap->first, 0, tot);
627 cdcSimHitsToCDCHits.
add(iterSignalMap->second.m_simHitIndex3, iCDCHits);
631 rels = m_simHits[iterSignalMap->second.m_simHitIndex3]->getRelationsFrom<
MCParticle>();
632 if (rels.
size() != 0) {
635 double weight = rels.
weight(0);
652 if (m_addXTalk) addXTalk();
656 for (iterSignalMapTrg = signalMapTrg.begin(); iterSignalMapTrg != signalMapTrg.end(); ++iterSignalMapTrg) {
665 unsigned short adcCount = iterSignalMapTrg->second.m_charge;
666 unsigned short tdcCount =
667 static_cast<unsigned short>((m_cdcgp->getT0(iterSignalMapTrg->first.first) -
668 iterSignalMapTrg->second.m_driftTime) * m_tdcBinWidthInv + 0.5);
669 const CDCHit* cdcHit = m_cdcHits4Trg.appendNew(tdcCount, adcCount, iterSignalMapTrg->first.first);
672 m_simHits[iterSignalMapTrg->second.m_simHitIndex]->
addRelationTo(cdcHit);
674 if (rels.
size() != 0) {
677 double weight = rels.
weight(0);
696 double CDCDigitizerModule::smearDriftLength(
const double driftLength,
const double dDdt)
701 if (m_useSimpleDigitization) {
702 if (gRandom->Uniform() <= m_fraction) {
704 resolution = m_resolution1;
707 resolution = m_resolution2;
710 const unsigned short leftRight = m_posFlag;
711 double alpha = m_cdcgp->getAlpha(m_posWire, m_momentum);
712 double theta = m_cdcgp->getTheta(m_momentum);
713 resolution = m_cdcgp->getSigma(driftLength, m_wireID.getICLayer(), leftRight, alpha, theta);
714 resolution *= m_totalFudgeFactor;
719 double diff = resolution - dDdt * m_tdcResol;
721 resolution = sqrt(diff * (resolution + dDdt * m_tdcResol));
726 #if defined(CDC_DEBUG)
728 cout <<
"CDCDigitizerModule::smearDriftLength" << endl;
729 cout <<
"tdcResol= " << m_tdcResol << endl;
730 cout <<
"dDdt,resolution= " << dDdt <<
" " << resolution << endl;
734 double newDL = gRandom->Gaus(driftLength + mean , resolution);
735 while (newDL <= 0.) newDL = gRandom->Gaus(driftLength + mean, resolution);
741 double CDCDigitizerModule::getdDdt(
const double driftL)
747 double dDdt = m_driftV;
749 if (!m_useSimpleDigitization) {
750 const unsigned short layer = m_wireID.getICLayer();
751 const unsigned short leftRight = m_posFlag;
752 double alpha = m_cdcgp->getAlpha(m_posWire, m_momentum);
753 double theta = m_cdcgp->getTheta(m_momentum);
754 double t = m_cdcgp->getDriftTime(driftL, layer, leftRight, alpha, theta);
755 dDdt = m_cdcgp->getDriftV(t, layer, leftRight, alpha, theta);
757 #if defined(CDC_DEBUG)
759 cout <<
"CDCDigitizerModule::getdDdt" << endl;
760 cout <<
"**layer= " << layer << endl;
761 cout <<
"alpha= " << 180.*alpha / M_PI << std::endl;
764 for (
int i = 0; i < 1000; ++i) {
766 double d = m_cdcgp->getDriftLength(t, layer, lr, alpha, theta);
767 cout << t <<
" " << d << endl;
773 for (
int i = 0; i < 100; ++i) {
775 double d = m_cdcgp->getDriftLength(t, layer, lr, alpha, theta);
776 cout << t <<
" " << d << endl;
787 double CDCDigitizerModule::getDriftTime(
const double driftLength,
const bool addTof,
const bool addDelay)
798 if (m_useSimpleDigitization) {
799 driftT = (driftLength / Unit::cm) * m_driftVInv;
801 #if defined(CDC_DEBUG)
803 cout <<
"CDCDigitizerModule::getDriftTime" << endl;
804 cout <<
"driftvinv= " << m_driftVInv << endl;
807 const unsigned short layer = m_wireID.getICLayer();
808 const unsigned short leftRight = m_posFlag;
809 double alpha = m_cdcgp->getAlpha(m_posWire, m_momentum);
810 double theta = m_cdcgp->getTheta(m_momentum);
811 driftT = m_cdcgp->getDriftTime(driftLength, layer, leftRight, alpha, theta);
816 driftT += m_flightTime;
822 TVector3 backWirePos = m_cdcgp->wireBackwardPosition(m_wireID, set);
824 double propLength = (m_posWire - backWirePos).Mag();
828 if (m_gcp->getSenseWireZposMode() == 1) {
829 const unsigned short layer = m_wireID.getICLayer();
830 propLength += m_cdcgp->getBwdDeltaZ(layer);
834 if (m_useSimpleDigitization) {
835 driftT += (propLength / Unit::cm) * m_propSpeedInv;
837 #if defined(CDC_DEBUG)
838 cout <<
"pseedinv= " << m_propSpeedInv << endl;
841 const unsigned short layer = m_wireID.getICLayer();
842 driftT += (propLength / Unit::cm) * m_cdcgp->getPropSpeedInv(layer);
843 #if defined(CDC_DEBUG)
844 cout <<
"layer,pseedinv= " << layer <<
" " << m_cdcgp->getPropSpeedInv(layer) << endl;
854 unsigned short CDCDigitizerModule::getADCCount(
const WireID& wid,
double dEinGeV,
double dx,
double costh)
856 unsigned short adcCount = 0;
857 if (dEinGeV <= 0. || dx <= 0.)
return adcCount;
862 double conv = (100.0 / 3.2);
863 const unsigned short layer = wid.
getICLayer();
864 const unsigned short cell = wid.
getIWire();
865 double dEInkeV = dEinGeV / Unit::keV;
866 if (m_spaceChargeEffect) {
867 if (m_useDB4EDepToADC) conv = m_cdcgp->getEDepToADCConvFactor(layer, cell, dEInkeV, dx, costh);
869 if (m_useDB4EDepToADC) conv = m_cdcgp->getEDepToADCMainFactor(layer, cell);
874 adcCount =
static_cast<unsigned short>(std::round(conv * dEInkeV));
880 void CDCDigitizerModule::setFEElectronics()
882 const double& off = m_tdcThresholdOffset;
883 const double& gA = m_analogGain;
884 const double& gD = m_digitalGain;
885 const double& adcBW = m_adcBinWidth;
886 const double convF = gA / gD / adcBW;
887 const double el1TrgLatency = m_cdcgp->getMeanT0();
888 B2DEBUG(29,
"L1TRGLatency= " << el1TrgLatency);
889 const double c = 32. * m_tdcBinWidth;
891 if (!m_fEElectronicsFromDB) B2FATAL(
"No FEEElectronics dbobject!");
893 int mode = (fp.getBoardID() == -1) ? 1 : 0;
894 int iNBoards =
static_cast<int>(nBoards);
898 for (
int bdi = 1; bdi < iNBoards; ++bdi) {
899 m_uprEdgeOfTimeWindow[bdi] = el1TrgLatency - c * (fp.getTrgDelay() + 1);
900 if (m_uprEdgeOfTimeWindow[bdi] < 0.) B2FATAL(
"CDCDigitizer: Upper edge of time window < 0!");
901 m_lowEdgeOfTimeWindow[bdi] = m_uprEdgeOfTimeWindow[bdi] - c * (fp.getWidthOfTimeWindow() + 1);
902 if (m_lowEdgeOfTimeWindow[bdi] > 0.) B2FATAL(
"CDCDigitizer: Lower edge of time window > 0!");
903 m_adcThresh[bdi] = fp.getADCThresh();
904 m_tdcThresh[bdi] = convF * (off - fp.getTDCThreshInMV());
905 m_widthOfTimeWindow[bdi] = fp.getWidthOfTimeWindow() + 1;
911 for (
const auto& fpp : (*m_fEElectronicsFromDB)) {
912 int bdi = fpp.getBoardID();
913 if (mode == 0 && bdi == 0)
continue;
914 if (mode == 1 && bdi == -1)
continue;
915 if (bdi < 0 || bdi >= iNBoards) B2FATAL(
"CDCDigitizer:: Invalid no. of FEE board!");
916 m_uprEdgeOfTimeWindow[bdi] = el1TrgLatency - c * (fpp.getTrgDelay() + 1);
917 if (m_uprEdgeOfTimeWindow[bdi] < 0.) B2FATAL(
"CDCDigitizer: Upper edge of time window < 0!");
918 m_lowEdgeOfTimeWindow[bdi] = m_uprEdgeOfTimeWindow[bdi] - c * (fpp.getWidthOfTimeWindow() + 1);
919 if (m_lowEdgeOfTimeWindow[bdi] > 0.) B2FATAL(
"CDCDigitizer: Lower edge of time window > 0!");
920 m_adcThresh[bdi] = fpp.getADCThresh();
921 m_tdcThresh[bdi] = convF * (off - fpp.getTDCThreshInMV());
922 m_widthOfTimeWindow[bdi] = fpp.getWidthOfTimeWindow() + 1;
926 B2DEBUG(29,
"mode= " << mode);
927 for (
int bdi = 1; bdi < iNBoards; ++bdi) {
928 B2DEBUG(29, bdi <<
" " << m_lowEdgeOfTimeWindow[bdi] <<
" " << m_uprEdgeOfTimeWindow[bdi] <<
" " << m_adcThresh[bdi] <<
" " <<
934 void CDCDigitizerModule::setRunGain()
936 m_runGain = (*m_runGainFromDB)->getRunGain();
937 m_runGain *= m_overallGainFactor;
940 void CDCDigitizerModule::addXTalk()
942 map<WireID, XTalkInfo> xTalkMap;
943 map<WireID, XTalkInfo> xTalkMap1;
944 map<WireID, XTalkInfo>::iterator iterXTalkMap1;
947 int OriginalNoOfHits = m_cdcHits.getEntries();
948 B2DEBUG(m_debugLevel4XTalk,
"\n \n" <<
"#CDCHits " << OriginalNoOfHits);
949 for (
const auto& aHit : m_cdcHits) {
950 if (m_issue2ndHitWarning && aHit.is2ndHit()) {
951 B2WARNING(
"2nd TDC hit found, but not ready for it!");
955 short tdcCount = aHit.getTDCCount();
956 short adcCount = aHit.getADCCount();
957 short tot = aHit.getTOT();
958 short board = m_cdcgp->getBoardID(wid);
959 short channel = m_cdcgp->getChannelID(wid);
960 const vector<pair<short, asicChannel>> xTalks = (*m_xTalkFromDB)->getLibraryCrossTalk(channel, tdcCount, adcCount, tot);
962 int nXTalks = xTalks.size();
963 for (
int i = 0; i < nXTalks; ++i) {
964 const unsigned short tdcCount4XTalk = xTalks[i].second.TDC;
966 B2DEBUG(m_debugLevel4XTalk,
"\n" <<
" signal: " << channel <<
" " << tdcCount <<
" " << adcCount <<
" " << tot);
968 B2DEBUG(m_debugLevel4XTalk,
"xtalk: " << xTalks[i].first <<
" " << tdcCount4XTalk <<
" " << xTalks[i].second.ADC <<
" " <<
969 xTalks[i].second.TOT);
970 WireID widx = m_cdcgp->getWireID(board, xTalks[i].first);
971 if (!m_cdcgp->isBadWire(widx)) {
972 if (m_includeEarlyXTalks || (xTalks[i].second.TDC <= tdcCount)) {
973 const double t0 = m_cdcgp->getT0(widx);
974 const double ULOfTDC = (t0 - m_lowEdgeOfTimeWindow[board]) * m_tdcBinWidthInv;
975 const double LLOfTDC = (t0 - m_uprEdgeOfTimeWindow[board]) * m_tdcBinWidthInv;
976 if (LLOfTDC <= tdcCount4XTalk && tdcCount4XTalk <= ULOfTDC) {
977 const unsigned short status = 0;
978 xTalkMap.insert(make_pair(widx,
XTalkInfo(tdcCount4XTalk, xTalks[i].second.ADC, xTalks[i].second.TOT, status)));
988 B2DEBUG(m_debugLevel4XTalk,
"#xtalk hits: " << xTalkMap.size());
989 for (
const auto& aHit : xTalkMap) {
992 iterXTalkMap1 = xTalkMap1.find(wid);
993 unsigned short tdcCount = aHit.second.m_tdc;
994 unsigned short adcCount = aHit.second.m_adc;
995 unsigned short tot = aHit.second.m_tot;
996 unsigned short status = aHit.second.m_status;
998 if (iterXTalkMap1 == xTalkMap1.end()) {
999 xTalkMap1.insert(make_pair(wid,
XTalkInfo(tdcCount, adcCount, tot, status)));
1002 if (tdcCount < iterXTalkMap1->second.m_tdc) {
1003 iterXTalkMap1->second.m_tdc = tdcCount;
1004 B2DEBUG(m_debugLevel4XTalk,
"TDC-count of current xtalk: " << tdcCount);
1006 iterXTalkMap1->second.m_adc += adcCount;
1007 iterXTalkMap1->second.m_tot += tot;
1012 B2DEBUG(m_debugLevel4XTalk,
"#xtalk1 hits: " << xTalkMap1.size());
1013 for (
const auto& aX : xTalkMap1) {
1015 const unsigned short tdc4Bg = aX.second.m_tdc;
1016 const unsigned short adc4Bg = aX.second.m_adc;
1017 const unsigned short tot4Bg = aX.second.m_tot;
1018 const unsigned short status4Bg = aX.second.m_status;
1020 for (
int iHit = 0; iHit < OriginalNoOfHits; ++iHit) {
1021 CDCHit& aH = *(m_cdcHits[iHit]);
1022 if (aH.
getID() != aX.first.getEWire()) {
1028 const unsigned short tot4Sg = aH.
getTOT();
1034 if (tdc4Sg < tdc4Bg) {
1038 for (
int i = relSimHits.size() - 1; i >= 0; --i) {
1039 relSimHits.remove(i);
1042 for (
int i = relMCParticles.size() - 1; i >= 0; --i) {
1043 relMCParticles.remove(i);
1053 unsigned short s1 = tdc4Sg;
1054 unsigned short s2 = tdc4Bg;
1055 unsigned short w1 = tot4Sg;
1056 unsigned short w2 = tot4Bg;
1057 if (tdc4Sg < tdc4Bg) {
1065 const unsigned short e1 = s1 - w1;
1066 const unsigned short e2 = s2 - w2;
1069 double pulseW = w1 + w2;
1072 }
else if (e1 <= s2) {
1076 unsigned short board = m_cdcgp->getBoardID(aX.first);
1077 aH.
setTOT(std::min(std::round(pulseW / 32.),
static_cast<double>(m_widthOfTimeWindow[board])));
1086 m_cdcHits.appendNew(tdc4Bg, adc4Bg, aX.first, status4Bg, tot4Bg);
1087 B2DEBUG(m_debugLevel4XTalk,
"appended tdc,adc,tot,wid,status= " << tdc4Bg <<
" " << adc4Bg <<
" " << tot4Bg <<
" " << aX.first <<
1092 B2DEBUG(m_debugLevel4XTalk,
"original #hits, #hits= " << OriginalNoOfHits <<
" " << m_cdcHits.getEntries());