9#include <cdc/modules/cdcDigitizer/CDCDigitizerModule.h>
10#include <cdc/modules/cdcDigitizer/EDepInGas.h>
11#include <cdc/utilities/ClosestApproach.h>
13#include <framework/datastore/RelationArray.h>
15#include <framework/gearbox/Unit.h>
16#include <framework/logging/Logger.h>
44 "Name of output array for trigger. Can contain several hits per wire, "
45 "if they correspond to different time windows of 32ns.",
46 string(
"CDCHits4Trg"));
50 "Name of relation between MCParticles and CDCSimHits used",
string(
""));
52 "Name of relation between the CDCSimHits and the CDCHits used",
string(
""));
54 "Optional name of relation between the first MCParticles and CDCHits used",
string(
"FirstMatchedParticles"));
56 "Optional name of relation between all MCParticles and CDCHits used",
string(
"AllMatchedParticles"));
61 "If true, a simple x-t with a constant velocity is used for the drift-length to -time conversion",
false);
64 addParam(
"Fraction",
m_fraction,
"Fraction of first Gaussian used to smear drift length in cm", 1.0);
65 addParam(
"Mean1",
m_mean1,
"Mean value of first Gaussian used to smear drift length in cm", 0.0000);
66 addParam(
"Resolution1",
m_resolution1,
"Resolution of first Gaussian used to smear drift length in cm", 0.0130);
67 addParam(
"Mean2",
m_mean2,
"Mean value of second Gaussian used to smear drift length in cm", 0.0000);
68 addParam(
"Resolution2",
m_resolution2,
"Resolution of second Gaussian used to smear drift length in cm", 0.0000);
72 "If false, drift length will not be smeared.",
true);
75 "If true, some hits will be removed for better charge asymmetry simulation.",
false);
78 "Magnitude (w) of trigger timing jitter (ns). The trigger timing is randuminzed uniformly in a time window of [-w/2, +w/2].",
81 addParam(
"AddTimeWalk",
m_addTimeWalk,
"A switch for time-walk (pulse-heght dep. delay); true: on; false: off",
true);
83 "A switch used to control adding propagation delay in the wire into the final drift time or not; this is for signal hits.",
true);
85 "The same switch but for beam bg. hits.",
true);
87 "A switch used to control adding time of flight into the final drift time or not; this is for signal hits.",
true);
89 "The same switch but for beam bg. hits.",
true);
92 "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.",
96 "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.",
103 "TDC threshold (dE in eV) for Layers#8-56. The value corresponds to He-C2H6 gas", 250.);
105 "Same as TDCThreshold4Outer but for Layers#0-7,", 150.);
107 "Mode for extracting energy deposit in gas from energy deposit in gas+wire; =0: scaling using electron density; 1: scaling using most probab. energy deposit; 2: similar to 2 but slightly different; 3: extraction based on probability; 4: regeneration following probability",
112 "Threshold for ADC-count (in unit of count). ADC-count < threshold is treated as count=0.", 2);
113 addParam(
"tMin",
m_tMin,
"Lower edge of time window in ns; valid only for UseDB4FEE=false", -100.);
114 addParam(
"tMaxOuter",
m_tMaxOuter,
"Upper edge of time window in ns for the normal-cell layers; valid only for UseDB4FEE=false",
116 addParam(
"tMaxInner",
m_tMaxInner,
"Upper edge of time window in ns for the small-cell layers; valid only for UseDB4FEE=false",
124 addParam(
"UseDB4FEE",
m_useDB4FEE,
"Fetch and use FEE params. from database or not",
true);
131 addParam(
"Randomization",
m_randomization,
"Randomize timing with other sub-detectors; valid only for Synchronization=false",
136 "L1 trigger timing offset in count, [0,7] in a trigger bin. The default value is from exp14, while the value from exp12 is 2. This run dependence may be taken into account later if needed",
148 "Additional fudge factor for space resol. (common to all cells)", 1.);
151 "Degree of space charge effect on timing threshold; specify the range [0,1]; =1: full effect on threshold; =0: no effect",
154 addParam(
"AddXTalk",
m_addXTalk,
"A switch for crosstalk; true: on; false: off",
true);
165 "Effective energy (keV) needed for one electron production for gas gain smearing; average for alpha- and beta-sources.",
174 "Switch to store all MCRelations for the first three SimHits instead of only the first",
false);
176#if defined(CDC_DEBUG)
178 cout <<
"CDCDigitizer constructor" << endl;
221 if ((*m_fEElectronicsFromDB).isValid()) {
225 B2FATAL(
"CDCDigitizer:: CDCFEElectronics not valid !");
243 if ((*m_runGainFromDB).isValid()) {
246 B2FATAL(
"CDCDedxRunGain invalid!");
250 if ((*m_gain0FromDB).isValid()) {
253 B2FATAL(
"CDCDedxScaleFactor invalid!");
257 if ((*m_wireGainFromDB).isValid()) {
261 B2FATAL(
"CDCDedxWireGain invalid!");
267 if ((*m_xTalkFromDB).isValid()) {
269 B2FATAL(
"CDCCrossTalkLibrary invalid!");
274 if ((*m_corrToThresholdFromDB).isValid()) {
276 B2FATAL(
"CDCCorrToThresholds invalid!");
281 B2FATAL(
"CDCAlphaScaleFactorForAsymmetry invalid!");
285#if defined(CDC_DEBUG)
287 cout <<
"CDCDigitizer initialize" << endl;
291 cout <<
"m_driftV= " <<
m_driftV << endl;
304 ushort firstLayerOffset =
m_cdcgp->getOffsetOfFirstLayer();
305 if (
m_cdcgp->getEDepToADCMainFactor(firstLayerOffset, 0) == 0.) {
306 B2FATAL(
"CDCEDepToADCConversion payloads are unavailable!");
338 map<WireID, SignalInfo> signalMap;
341 map<WireID, unsigned short> adcMap;
347 map<pair<WireID, unsigned>,
SignalInfo> signalMapTrg;
348 map<pair<WireID, unsigned>,
SignalInfo>::iterator iterSignalMapTrg;
351 map<WireID, std::set<int>> particleMap;
352 map<WireID, std::set<int>>::iterator iterParticleMap;
362 B2DEBUG(
m_debugLevel,
"SimClockState unavailable so switched the mode from synchro to random.");
365 trigBin = gRandom->Integer(4);
367 if (trigBin < 0 || trigBin > 3) B2ERROR(
"Invalid trigger bin; must be an integer [0,3]!");
372 for (
unsigned short bd = 1; bd < c_nBoards; ++bd) {
375 B2DEBUG(
m_debugLevel, bd <<
" " << tMinInCount <<
" " << tMaxInCount);
386 B2DEBUG(
m_debugLevel,
"Number of CDCSimHits in the current event: " << nHits);
387 for (
int iHits = 0; iHits < nHits; ++iHits) {
393 if (
m_wireID.getISuperLayer() <
m_cdcgp->getOffsetOfFirstSuperLayer()) {
394 B2FATAL(
"SimHit with wireID " <<
m_wireID <<
" is in CDC SuperLayer: " <<
m_wireID.getISuperLayer() <<
" which should not happen.");
419 if ((bwpAlign - bwp).Mag() == 0. && (fwpAlign - fwp).Mag() == 0.)
m_align =
false;
429 double bckYSag = bwp.
Y();
430 double forYSag = fwp.
Y();
435 const int layerID =
m_wireID.getICLayer();
436 const int wireID =
m_wireID.getIWire();
437 m_cdcgp->getWireSagEffect(set, layerID, wireID, zpos, bckYSag, forYSag);
454 double deltaTime = 0.;
464 double dDdt =
getdDdt(hitDriftLength);
476 double hitDriftTime =
getDriftTime(hitDriftLength, addTof, addDelay);
484 hitDriftTime += trigTiming;
494 if (hitDriftTime < tMin || hitDriftTime > tMax)
continue;
500 if (
m_cdcgp->getMaterialDefinitionMode() != 2) {
505 double convFactorForThreshold = 1;
507 unsigned short adcCount = 0;
511 if (adcCount < adcTh) adcCount = 0;
513 if (iterADCMap == adcMap.end()) {
514 adcMap.insert(make_pair(
m_wireID, adcCount));
517 iterADCMap->second += adcCount;
523 double dEThreshold = 0.;
530 dEThreshold *= (*m_corrToThresholdFromDB)->getParam(
m_wireID.getICLayer());
531 B2DEBUG(
m_debugLevel,
"hitdE,dEThreshold,driftLength " << hitdE <<
" " << dEThreshold <<
" " << hitDriftLength);
533 if (hitdE < dEThreshold) {
534 B2DEBUG(
m_debugLevel,
"Below Ethreshold: " << hitdE <<
" " << dEThreshold);
539 unsigned short trigWindow = floor((hitDriftTime - tMin) *
m_tdcBinWidthInv / 32);
540 iterSignalMapTrg = signalMapTrg.find(make_pair(
m_wireID, trigWindow));
541 if (iterSignalMapTrg == signalMapTrg.end()) {
544 signalMapTrg.insert(make_pair(make_pair(
m_wireID, trigWindow),
547 if (hitDriftTime < iterSignalMapTrg->second.m_driftTime) {
548 iterSignalMapTrg->second.m_driftTime = hitDriftTime;
549 iterSignalMapTrg->second.m_simHitIndex = iHits;
552 iterSignalMapTrg->second.m_charge += adcCount;
566 if (eff < gRandom->Uniform())
continue;
570 const double a = bwpAlign.
X();
571 const double b = bwpAlign.
Y();
572 const double c = bwpAlign.
Z();
573 const B2Vector3D fmbAlign = fwpAlign - bwpAlign;
574 const double lmn = 1. / fmbAlign.
Mag();
575 const double l = fmbAlign.
X() * lmn;
576 const double m = fmbAlign.
Y() * lmn;
577 const double n = fmbAlign.
Z() * lmn;
582 double sub = l * dx + m * dy + n * dz;
583 const double driftLFromIn =
sqrt(dx * dx + dy * dy + dz * dz - sub * sub);
588 sub = l * dx + m * dy + n * dz;
589 const double driftLFromOut =
sqrt(dx * dx + dy * dy + dz * dz - sub * sub);
591 const double maxDriftL = std::max(driftLFromIn, driftLFromOut);
593 B2DEBUG(
m_debugLevel,
"driftLFromIn= " << driftLFromIn <<
" driftLFromOut= " << driftLFromOut <<
" minDriftL= " << minDriftL <<
598 iterSignalMap = signalMap.find(
m_wireID);
601 iterParticleMap = particleMap.find(
m_wireID);
605 if (rels.
size() != 0) {
608 mcIndex = int(mcparticle->
getIndex());
613 if (iterParticleMap == particleMap.end()) {
614 std::set<int> vecmc = {mcIndex};
615 particleMap.insert(make_pair(
m_wireID, vecmc));
617 iterParticleMap->second.insert(mcIndex);
622 if (iterSignalMap == signalMap.end()) {
625 signalMap.insert(make_pair(
m_wireID,
SignalInfo(iHits, hitDriftTime, adcCount, maxDriftL, minDriftL)));
629 if (hitDriftTime < iterSignalMap->second.m_driftTime) {
630 iterSignalMap->second.m_driftTime3 = iterSignalMap->second.m_driftTime2;
631 iterSignalMap->second.m_simHitIndex3 = iterSignalMap->second.m_simHitIndex2;
632 iterSignalMap->second.m_driftTime2 = iterSignalMap->second.m_driftTime;
633 iterSignalMap->second.m_simHitIndex2 = iterSignalMap->second.m_simHitIndex;
634 iterSignalMap->second.m_driftTime = hitDriftTime;
635 iterSignalMap->second.m_simHitIndex = iHits;
636 B2DEBUG(
m_debugLevel,
"hitDriftTime of current Signal: " << hitDriftTime <<
", hitDriftLength: " << hitDriftLength);
637 }
else if (hitDriftTime < iterSignalMap->second.m_driftTime2) {
638 iterSignalMap->second.m_driftTime3 = iterSignalMap->second.m_driftTime2;
639 iterSignalMap->second.m_simHitIndex3 = iterSignalMap->second.m_simHitIndex2;
640 iterSignalMap->second.m_driftTime2 = hitDriftTime;
641 iterSignalMap->second.m_simHitIndex2 = iHits;
642 }
else if (hitDriftTime < iterSignalMap->second.m_driftTime3) {
643 iterSignalMap->second.m_driftTime3 = hitDriftTime;
644 iterSignalMap->second.m_simHitIndex3 = iHits;
648 iterSignalMap->second.m_charge += adcCount;
651 if (iterSignalMap->second.m_maxDriftL < maxDriftL) iterSignalMap->second.m_maxDriftL = maxDriftL;
652 if (iterSignalMap->second.m_minDriftL > minDriftL) iterSignalMap->second.m_minDriftL = minDriftL;
653 B2DEBUG(
m_debugLevel,
"maxDriftL in struct= " << iterSignalMap->second.m_maxDriftL <<
"minDriftL in struct= " <<
654 iterSignalMap->second.m_minDriftL);
662 unsigned int iCDCHits = 0;
666 for (iterSignalMap = signalMap.begin(); iterSignalMap != signalMap.end(); ++iterSignalMap) {
671 iterADCMap = adcMap.find(iterSignalMap->first);
672 unsigned short adcCount = iterADCMap != adcMap.end() ? iterADCMap->second : 0;
685 B2DEBUG(
m_debugLevel,
"timewalk= " <<
m_cdcgp->getTimeWalk(iterSignalMap->first, adcCount));
686 iterSignalMap->second.m_driftTime +=
m_cdcgp->getTimeWalk(iterSignalMap->first, adcCount);
691 iterSignalMap->second.m_driftTime < 0.) {
697 int iHits = iterSignalMap->second.m_simHitIndex;
705 double random = gRandom->Uniform();
706 if ((Scale < 1) && (alpha > 0)) {
707 if (random > Scale) continue ;
709 if ((Scale > 1) && (alpha < 0)) {
710 if (random > 1. / Scale) continue ;
715 unsigned short tdcCount =
static_cast<unsigned short>((
getPositiveT0(iterSignalMap->first) - iterSignalMap->second.m_driftTime) *
719 double deltaDL = iterSignalMap->second.m_maxDriftL - iterSignalMap->second.m_minDriftL;
724 const unsigned short boardID =
m_cdcgp->getBoardID(iterSignalMap->first);
725 unsigned short tot = std::min(std::round(5.92749 * deltaDL + 2.59706),
static_cast<double>(
m_widthOfTimeWindowInCount[boardID]));
727 tot = std::min(
static_cast<int>(tot),
static_cast<int>(adcCount /
m_adcThresh[boardID]));
730 CDCHit* firstHit =
m_cdcHits.appendNew(tdcCount, adcCount, iterSignalMap->first, 0, tot);
733 cdcSimHitsToCDCHits.
add(iterSignalMap->second.m_simHitIndex, iCDCHits);
737 if (rels.
size() != 0) {
740 double weight = rels.
weight(0);
746 iterParticleMap = particleMap.find(iterSignalMap->first);
747 if (iterParticleMap != particleMap.end()) {
748 std::set<int> vv = iterParticleMap->second;
749 for (std::set<int>::iterator it = vv.begin(); it != vv.end(); ++it) {
761 if (iterSignalMap->second.m_simHitIndex >= 0) {
763 if (rels1.
size() != 0) {
766 double weight = rels1.
weight(0);
772 if (iterSignalMap->second.m_simHitIndex2 >= 0) {
774 if (rels2.
size() != 0) {
777 double weight = rels2.
weight(0);
783 if (iterSignalMap->second.m_simHitIndex3 >= 0) {
785 if (rels3.
size() != 0) {
788 double weight = rels3.
weight(0);
797 if (
m_output2ndHit && iterSignalMap->second.m_simHitIndex2 >= 0) {
798 unsigned short tdcCount2 =
static_cast<unsigned short>((
getPositiveT0(iterSignalMap->first) - iterSignalMap->second.m_driftTime2) *
800 if (tdcCount2 != tdcCount) {
801 CDCHit* secondHit =
m_cdcHits.appendNew(tdcCount2, adcCount, iterSignalMap->first, 0, tot);
814 cdcSimHitsToCDCHits.
add(iterSignalMap->second.m_simHitIndex2, iCDCHits);
819 if (rels.
size() != 0) {
822 double weight = rels.
weight(0);
827 if (iterSignalMap->second.m_simHitIndex3 >= 0) {
828 unsigned short tdcCount3 =
static_cast<unsigned short>((
getPositiveT0(iterSignalMap->first) - iterSignalMap->second.m_driftTime3) *
831 if (tdcCount3 != tdcCount) {
832 CDCHit* secondHit =
m_cdcHits.appendNew(tdcCount3, adcCount, iterSignalMap->first, 0, tot);
841 cdcSimHitsToCDCHits.
add(iterSignalMap->second.m_simHitIndex3, iCDCHits);
846 if (rels.
size() != 0) {
849 double weight = rels.
weight(0);
870 for (iterSignalMapTrg = signalMapTrg.begin(); iterSignalMapTrg != signalMapTrg.end(); ++iterSignalMapTrg) {
879 unsigned short adcCount = iterSignalMapTrg->second.m_charge;
880 unsigned short tdcCount =
881 static_cast<unsigned short>((
getPositiveT0(iterSignalMapTrg->first.first) -
883 const CDCHit* cdcHit =
m_cdcHits4Trg.appendNew(tdcCount, adcCount, iterSignalMapTrg->first.first);
886 m_simHits[iterSignalMapTrg->second.m_simHitIndex]->addRelationTo(cdcHit);
888 if (rels.
size() != 0) {
891 double weight = rels.
weight(0);
924 const unsigned short leftRight =
m_posFlag;
927 resolution =
m_cdcgp->getSigma(driftLength,
m_wireID.getICLayer(), leftRight, alpha, theta);
940#if defined(CDC_DEBUG)
942 cout <<
"CDCDigitizerModule::smearDriftLength" << endl;
944 cout <<
"dDdt,resolution= " << dDdt <<
" " << resolution << endl;
948 double newDL = gRandom->Gaus(driftLength + mean, resolution);
949 while (newDL <= 0.) newDL = gRandom->Gaus(driftLength + mean, resolution);
964 const unsigned short layer =
m_wireID.getICLayer();
965 const unsigned short leftRight =
m_posFlag;
968 double t =
m_cdcgp->getDriftTime(driftL, layer, leftRight, alpha, theta);
969 dDdt =
m_cdcgp->getDriftV(t, layer, leftRight, alpha, theta);
971#if defined(CDC_DEBUG)
973 cout <<
"CDCDigitizerModule::getdDdt" << endl;
974 cout <<
"**layer= " << layer << endl;
975 cout <<
"alpha= " << 180.*alpha / M_PI << std::endl;
978 for (
int i = 0; i < 1000; ++i) {
980 double d =
m_cdcgp->getDriftLength(t, layer, lr, alpha, theta);
981 cout << t <<
" " << d << endl;
987 for (
int i = 0; i < 100; ++i) {
989 double d =
m_cdcgp->getDriftLength(t, layer, lr, alpha, theta);
990 cout << t <<
" " << d << endl;
1015#if defined(CDC_DEBUG)
1016 cout <<
" " << endl;
1017 cout <<
"CDCDigitizerModule::getDriftTime" << endl;
1021 const unsigned short layer =
m_wireID.getICLayer();
1022 const unsigned short leftRight =
m_posFlag;
1025 driftT =
m_cdcgp->getDriftTime(driftLength, layer, leftRight, alpha, theta);
1038 double propLength = (
m_posWire - backWirePos).Mag();
1042 if (
m_gcp->getSenseWireZposMode() == 1) {
1043 const unsigned short layer =
m_wireID.getICLayer();
1044 propLength +=
m_cdcgp->getBwdDeltaZ(layer);
1051#if defined(CDC_DEBUG)
1055 const unsigned short layer =
m_wireID.getICLayer();
1057#if defined(CDC_DEBUG)
1058 cout <<
"layer,pseedinv= " << layer <<
" " <<
m_cdcgp->getPropSpeedInv(layer) << endl;
1068 unsigned short& adcCount,
double& convFactorForThreshold)
1070 static double conv00 = (100.0 / 3.2);
1071 convFactorForThreshold = conv00;
1073 if (dEinGeV <= 0. || dx <= 0.)
return;
1075 const unsigned short layer = wid.
getICLayer();
1076 const unsigned short cell = wid.
getIWire();
1079 double conv = conv00;
1082 conv =
m_cdcgp->getEDepToADCConvFactor(layer, cell, dEInkeV, dx, costh);
1083 double conv0 =
m_cdcgp->getEDepToADCMainFactor(layer, cell, costh);
1088 convFactorForThreshold = conv;
1091 if (convFactorForThreshold > 0.) {
1094 convFactorForThreshold = conv00;
1100 if (20 <= nElectrons) {
1101 relGain = std::max(0., gRandom->Gaus(1.,
sqrt(1. / (nElectrons * (1. +
m_thetaOfPolya)))));
1102 }
else if (1 <= nElectrons) {
1103 for (
int i = 1; i <= nElectrons; ++i) {
1106 relGain /= nElectrons;
1114 conv *= max(0., gRandom->Gaus(1.,
m_cdcgp->getEDepToADCSigma(layer, cell)));
1120 adcCount =
static_cast<unsigned short>(std::round(conv * dEInkeV));
1133 gRandom->RndmArray(2, urndm);
1134 x = xmax * urndm[0];
1137 y = ymax * urndm[1];
1150 const double convF = gA / gD / adcBW;
1151 const double el1TrgLatency =
m_cdcgp->getMeanT0();
1152 B2DEBUG(
m_debugLevel,
"L1TRGLatency= " << el1TrgLatency);
1157 int mode = (fp.getBoardID() == -1) ? 1 : 0;
1158 int iNBoards =
static_cast<int>(c_nBoards);
1162 for (
int bdi = 1; bdi < iNBoards; ++bdi) {
1168 m_tdcThresh[bdi] = convF * (off - fp.getTDCThreshInMV());
1177 int bdi = fpp.getBoardID();
1178 if (mode == 0 && bdi == 0)
continue;
1179 if (mode == 1 && bdi == -1)
continue;
1180 if (bdi < 0 || bdi >= iNBoards) B2FATAL(
"CDCDigitizer:: Invalid no. of FEE board!");
1186 m_tdcThresh[bdi] = convF * (off - fpp.getTDCThreshInMV());
1193 for (
int bdi = 1; bdi < iNBoards; ++bdi) {
1206 const int nLyrs = c_maxNSenseLayers;
1208 int nGoodL[nLyrs] = {};
1209 float wgL[nLyrs] = {};
1210 int nGoodSL[c_nSuperLayers] = {};
1211 float wgSL[c_nSuperLayers] = {};
1215 for (
int lyr = 0; lyr < nLyrs; ++lyr) {
1216 int nWs =
m_cdcgp->nWiresInLayer(lyr);
1217 for (
int w = 0; w < nWs; ++w) {
1219 float wg = (*m_wireGainFromDB)->getWireGain(iw);
1234 for (
int lyr = 0; lyr < nLyrs; ++lyr) {
1235 if (nGoodL[lyr] > 0) wgL[lyr] /= nGoodL[lyr];
1236 B2DEBUG(
m_debugLevel,
"lyr,ngood,gain= " << lyr <<
" " << nGoodL[lyr] <<
" " << wgL[lyr]);
1239 for (
unsigned int sl = 0; sl < c_nSuperLayers; ++sl) {
1240 if (nGoodSL[sl] > 0) wgSL[sl] /= nGoodSL[sl];
1241 B2DEBUG(
m_debugLevel,
"slyr,ngood,gain= " << sl <<
" " << nGoodSL[sl] <<
" " << wgSL[sl]);
1249 B2FATAL(
"No good wires !");
1251 B2DEBUG(
m_debugLevel,
"ngoodAll,gain= " << nGoodAll <<
" " << wgAll);
1254 for (
int lyr = 0; lyr < nLyrs; ++lyr) {
1255 int nWs =
m_cdcgp->nWiresInLayer(lyr);
1256 for (
int w = 0; w < nWs; ++w) {
1269 for (
int lyr = 0; lyr < nLyrs; ++lyr) {
1270 int nWs =
m_cdcgp->nWiresInLayer(lyr);
1271 for (
int w = 0; w < nWs; ++w) {
1273 B2WARNING(
"Gain for lyr and wire " << lyr <<
" " << w <<
"not > 0. Strange! Replace it with " << wgAll <<
".");
1280 m_runGain = (*m_runGainFromDB)->getRunGain();
1281 double cgain = (*m_gain0FromDB)->getScaleFactor();
1284 for (
int lyr = 0; lyr < nLyrs; ++lyr) {
1285 int nWs =
m_cdcgp->nWiresInLayer(lyr);
1286 for (
int w = 0; w < nWs; ++w) {
1296 map<WireID, XTalkInfo> xTalkMap;
1297 map<WireID, XTalkInfo> xTalkMap1;
1301 int OriginalNoOfHits =
m_cdcHits.getEntries();
1305 B2WARNING(
"2nd TDC hit found, but not ready for it!");
1307 WireID wid(aHit.getID());
1309 short tdcCount = aHit.getTDCCount();
1310 short adcCount = aHit.getADCCount();
1311 short tot = aHit.getTOT();
1312 short board =
m_cdcgp->getBoardID(wid);
1313 short channel =
m_cdcgp->getChannelID(wid);
1314 const vector<pair<short, asicChannel>> xTalks = (*m_xTalkFromDB)->getLibraryCrossTalk(channel, tdcCount, adcCount, tot);
1316 int nXTalks = xTalks.size();
1317 for (
int i = 0; i < nXTalks; ++i) {
1318 const unsigned short tdcCount4XTalk = xTalks[i].second.TDC;
1320 B2DEBUG(
m_debugLevel4XTalk,
"\n" <<
" signal: " << channel <<
" " << tdcCount <<
" " << adcCount <<
" " << tot);
1322 B2DEBUG(
m_debugLevel4XTalk,
"xtalk: " << xTalks[i].first <<
" " << tdcCount4XTalk <<
" " << xTalks[i].second.ADC <<
" " <<
1323 xTalks[i].second.TOT);
1325 if (!
m_cdcgp->isBadWire(widx)) {
1330 if (LLOfTDC <= tdcCount4XTalk && tdcCount4XTalk <= ULOfTDC) {
1331 const unsigned short status = 0;
1332 xTalkMap.insert(make_pair(widx,
XTalkInfo(tdcCount4XTalk, xTalks[i].second.ADC, xTalks[i].second.TOT, status)));
1343 for (
const auto& aHit : xTalkMap) {
1346 iterXTalkMap1 = xTalkMap1.find(wid);
1347 unsigned short tdcCount = aHit.second.m_tdc;
1348 unsigned short adcCount = aHit.second.m_adc;
1349 unsigned short tot = aHit.second.m_tot;
1350 unsigned short status = aHit.second.m_status;
1352 if (iterXTalkMap1 == xTalkMap1.end()) {
1353 xTalkMap1.insert(make_pair(wid,
XTalkInfo(tdcCount, adcCount, tot, status)));
1356 if (tdcCount < iterXTalkMap1->second.m_tdc) {
1357 iterXTalkMap1->second.m_tdc = tdcCount;
1360 iterXTalkMap1->second.m_adc += adcCount;
1361 iterXTalkMap1->second.m_tot += tot;
1367 for (
const auto& aX : xTalkMap1) {
1369 const unsigned short tdc4Bg = aX.second.m_tdc;
1370 const unsigned short adc4Bg = aX.second.m_adc;
1371 const unsigned short tot4Bg = aX.second.m_tot;
1372 const unsigned short status4Bg = aX.second.m_status;
1374 for (
int iHit = 0; iHit < OriginalNoOfHits; ++iHit) {
1376 if (aH.
getID() != aX.first.getEWire()) {
1382 const unsigned short tot4Sg = aH.
getTOT();
1388 if (tdc4Sg < tdc4Bg) {
1392 for (
int i = relSimHits.size() - 1; i >= 0; --i) {
1393 relSimHits.remove(i);
1396 for (
int i = relMCParticles.size() - 1; i >= 0; --i) {
1397 relMCParticles.remove(i);
1407 unsigned short s1 = tdc4Sg;
1408 unsigned short s2 = tdc4Bg;
1409 unsigned short w1 = tot4Sg;
1410 unsigned short w2 = tot4Bg;
1411 if (tdc4Sg < tdc4Bg) {
1419 const unsigned short e1 = s1 - w1;
1420 const unsigned short e2 = s2 - w2;
1423 double pulseW = w1 + w2;
1426 }
else if (e1 <= s2) {
1430 unsigned short board =
m_cdcgp->getBoardID(aX.first);
1440 m_cdcHits.appendNew(tdc4Bg, adc4Bg, aX.first, status4Bg, tot4Bg);
1441 B2DEBUG(
m_debugLevel4XTalk,
"appended tdc,adc,tot,wid,status= " << tdc4Bg <<
" " << adc4Bg <<
" " << tot4Bg <<
" " << aX.first <<
1452 double t0 =
m_cdcgp->getT0(wid);
DataType Z() const
access variable Z (= .at(2) without boundary check)
DataType X() const
access variable X (= .at(0) without boundary check)
DataType Y() const
access variable Y (= .at(1) without boundary check)
DataType Mag() const
The magnitude (rho in spherical coordinate system).
void SetY(DataType y)
set Y/2nd-coordinate
bool m_alphaCorrection
Switch for alpha correction.
double m_tdcThreshold4Inner
TDC threshold for inner layers in unit of eV.
double m_driftV
Nominal drift velocity (in cm/ns)
double Polya(double xmax=10)
Generate random number according to Polya distribution.
unsigned short m_boardID
FEE board ID.
int m_eDepInGasMode
Mode for extracting dE(gas) from dE(gas+wire)
double m_tMaxInner
Upper edge of time window in ns for the inner layers.
double getdDdt(double driftLength)
The method to get dD/dt.
bool m_outputNegativeDriftTime
A switch to output negative drift time to CDCHit.
unsigned short m_widthOfTimeWindowInCount[c_nBoards]
Width of time window.
bool m_extraADCSmearing
Switch for extra ADC smearing.
int m_debugLevel
Debug level.
bool m_includeEarlyXTalks
Flag to switch on/off xtalks earlier than the hit.
void makeSignalsAfterShapers(const WireID &wid, double edep, double dx, double costh, unsigned short &adcCount, double &convFactorForThreshold)
Function to write ADC-count and conversion factor for threshold.
double m_propSpeedInv
Inv.
int m_trgTimingOffsetInCount
Trigger timing offset in unit of count.
bool m_addInWirePropagationDelay4Bg
A switch used to control adding propagation delay into the total drift time or not for beam bg.
double getSemiTotalGain(int clayer, int cell) const
Return semi-total gain of the specified wire.
bool m_useDB4RunGain
Fetch run gain from DB.
DBObjPtr< CDCCrossTalkLibrary > * m_xTalkFromDB
Pointer to cross-talk from DB.
double m_analogGain
analog gain (V/pC)
double m_tdcBinWidthInv
m_tdcBinWidth^-1 (in ns^-1)
float m_uprEdgeOfTimeWindow[c_nBoards]
Upper edge of time-window.
double m_tMin
Lower edge of time window in ns.
bool m_output2ndHit
A switch to output 2nd hit.
int m_adcThreshold
Threshold for ADC in unit of count.
double m_resolution2
Resolution of the second Gaussian used to smear drift length.
bool m_addTimeOfFlight
A switch used to control adding time of flight into the total drift time or not.
double m_trigTimeJitter
Magnitude of trigger timing jitter (ns).
void initialize() override
Initialize variables, print info, and start CPU clock.
double smearDriftLength(double driftLength, double dDdt)
Method used to smear the drift length.
std::string m_MCParticlesToSimHitsName
Relation for origin of incoming SimHits.
void addXTalk()
Add crosstalk.
double m_digitalGain
digital gain (V/pC)
double m_totalFudgeFactor
total fudge factor for space resol.
std::string m_OptionalAllMCParticlesToHitsName
Relation name for optional matching of all MCParticles.
bool m_addInWirePropagationDelay
A switch used to control adding propagation delay into the total drift time or not.
bool m_align
A switch to control alignment.
void event() override
Actual digitization of all hits in the CDC.
unsigned short m_trgDelayInCount[c_nBoards]
Trigger delay in frontend electronics in count.
double m_globalTime
global time of this hit
bool m_doSmearing
A switch to control drift length smearing.
std::string m_inputCDCSimHitsName
Input array name.
DBObjPtr< CDCDedxWireGain > * m_wireGainFromDB
Pointer to wire gain from DB.
B2Vector3D m_posWire
wire position of this hit
bool m_addTimeOfFlight4Bg
A switch used to control adding time of flight into the total drift time or not for beam bg.
bool m_addXTalk
Flag to switch on/off crosstalk.
bool m_addTimeWalk
A switch used to control adding time-walk delay into the total drift time or not.
bool m_treatNegT0WiresAsGood
A switch for negative-t0 wires.
double m_addFudgeFactorForSigma
additional fudge factor for space resol.
bool m_matchAllMCParticles
A switch to match all particles to a hit, regardless whether they produced a hit or not.
bool m_useDB4EDepToADC
Fetch edep-to-ADC conversion params.
std::string m_SimHitsTOCDCHitsName
Relation for outgoing CDCHits.
int m_offsetForTriggerBin
Input to getCDCTriggerBin(offset)
void setSemiTotalGain()
Set semi-total gain (from DB)
double m_driftLength
drift length of this hit
double m_tdcResol
TDC resolution (in ns)
int m_tSimMode
Timing simulation mode.
bool m_useDB4FEE
Fetch FEE params from DB.
double m_adcBinWidth
ADC bin width (mV)
double m_runGain
run gain.
double m_thetaOfPolya
theta of Polya function for gas gain smearing
bool m_matchFirstMCParticles
A switch to match first three MCParticles, not just the one with smallest drift time.
bool m_spaceChargeEffect
Space charge effect.
double m_driftVInv
m_driftV^-1 (in ns/cm)
StoreArray< CDCSimHit > m_simHits
CDCSimHit array.
DBObjPtr< CDCCorrToThresholds > * m_corrToThresholdFromDB
Pointer to threshold correction from DB.
std::string m_OptionalFirstMCParticlesToHitsName
Relation name for optional matching of up to first three MCParticles.
StoreArray< CDCHit > m_cdcHits4Trg
CDCHit4trg array.
float m_semiTotalGain[c_maxNSenseLayers][c_maxNDriftCells]
total gain per wire
DBObjPtr< CDCDedxRunGain > * m_runGainFromDB
Pointer to run gain from DB.
double m_degOfSPEOnThreshold
Degree of space charge effect on timing threshold.
CDCDigitizerModule()
Constructor.
bool m_randomization
Flag to switch on/off timing randomization.
CDCSimHit * m_aCDCSimHit
Pointer to CDCSimHit.
WireID m_wireID
WireID of this hit.
DBObjPtr< CDCDedxScaleFactor > * m_gain0FromDB
Pointer to overall gain factor from DB.
double m_effWForGasGainSmearing
Effective energy (keV) for one electron prod.
bool m_issue2ndHitWarning
Flag to switch on/off a warning on the 2nd TDC hit.
double getPositiveT0(const WireID &)
Modify t0 for negative-t0 case.
double m_tdcThreshold4Outer
TDC threshold for outer layers in unit of eV.
double m_tdcThresholdOffset
Offset for TDC(digital) threshold (mV)
double m_fraction
Fraction of the first Gaussian used to smear drift length.
int m_debugLevel4XTalk
Debug level for crosstalk.
int m_shiftOfTimeWindowIn32Count
Shift of time window for synchronization in 32count.
unsigned short m_posFlag
left or right flag of this hit
std::string m_outputCDCHitsName4Trg
Output array name for trigger.
StoreArray< CDCHit > m_cdcHits
CDCHit array.
double m_mean1
Mean value of the first Gaussian used to smear drift length.
bool m_synchronization
Flag to switch on/off timing synchronization.
float m_lowEdgeOfTimeWindow[c_nBoards]
Lower edge of time-window.
double getDriftTime(double driftLength, bool addTof, bool addDelay)
The method to get drift time based on drift length.
std::string m_outputCDCHitsName
Output array name.
float m_tdcThresh[c_nBoards]
Threshold for timing-signal.
B2Vector3D m_momentum
3-momentum of this hit
StoreArray< MCParticle > m_mcParticles
Set edep-to-ADC conversion params.
bool m_gasGainSmearing
Switch for gas gain smearing.
double m_tMaxOuter
Upper edge of time window in ns for the outer layers.
DBArray< CDCFEElectronics > * m_fEElectronicsFromDB
Pointer to FE electronics params.
double m_resolution1
Resolution of the first Gaussian used to smear drift length.
B2Vector3D m_posTrack
track position of this hit
bool m_useSimpleDigitization
Use float Gaussian Smearing instead of proper digitization.
CDC::CDCGeometryPar * m_cdcgp
Cached Pointer to CDCGeometryPar.
double m_overallGainFactor
Overall gain factor.
void setFEElectronics()
Set FEE parameters (from DB)
double m_mean2
Mean value of the second Gaussian used to smear drift length.
OptionalDBObjPtr< CDCAlphaScaleFactorForAsymmetry > m_alphaScaleFactorsFromDB
the ratio of data to MC, for eff(alpha<0)/eff(alpha>0).
double m_tdcBinWidth
Width of a TDC bin (in ns)
CDC::CDCGeoControlPar * m_gcp
Cached pointer to CDCGeoControlPar.
bool m_correctForWireSag
A switch to control wire sag.
StoreObjPtr< SimClockState > m_simClockState
generated hardware clock state
double m_flightTime
flight time of this hit
float m_adcThresh[c_nBoards]
Threshold for FADC.
Database object for Fron-endt electronics params.
Class containing the result of the unpacker in raw data and the result of the digitizer in simulation...
void setTDCCount(short tdcCount)
Setter for TDC count.
void set2ndHitFlag()
Setter for 2nd hit flag.
short getTDCCount() const
Getter for TDC count.
void setTOT(unsigned short tot)
Setter for TOT.
unsigned short getID() const
Getter for encoded wire number.
unsigned short getStatus() const
Getter for CDCHit status.
void setADCCount(unsigned short adcCount)
Setter for ADC count.
unsigned short getADCCount() const
Getter for integrated charge.
unsigned short getTOT() const
Getter for TOT.
void setOtherHitIndices(CDCHit *otherHit)
Setter for the other hit indices.
void setStatus(unsigned short status)
Setter for CDCHit status.
static CDCGeoControlPar & getInstance()
Static method to get a reference to the CDCGeoControlPar instance.
The Class for CDC Geometry Parameters.
double getNominalPropSpeed() const
Return the nominal propagation speed of the sense wire (default: 27.25 cm/nsec).
EWirePosition
Wire position set.
double getTdcBinWidth() const
Return TDC bin width (nsec).
double getNominalDriftV() const
Return the nominal drift velocity of He-ethane gas (default: 4.0x10^-3 cm/nsec).
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
The Class for Energy deposit in the gas.
static EDepInGas & getInstance()
Static method to get a reference to the EDepInGas instance.
double getEDepInGas(int mode, int pdg, double p, double dx, double e3) const
Return the energy deosite in the gas.
Class for accessing arrays of objects in the database.
Class for accessing objects in the database.
@ c_WriteOut
Object/array should be saved by output modules.
@ c_Event
Different object in each event, all objects/arrays are invalidated after event() function has been ca...
A Class to store the Monte Carlo particle information.
int getIndex() const
Get 1-based index of the particle in the corresponding MCParticle list.
void setDescription(const std::string &description)
Sets the description of the module.
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Low-level class to create/modify relations between StoreArrays.
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
Class for type safe access to objects that are referred to in relations.
size_t size() const
Get number of relations.
float weight(int index) const
Get weight with index.
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
RelationVector< FROM > getRelationsFrom(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from another store array to this object.
static const double keV
[kiloelectronvolt]
static const double eV
[electronvolt]
static const double cm
Standard units with the value = 1.
Class to identify a wire inside the CDC.
unsigned short getICLayer() const
Getter for continuous layer numbering.
unsigned short getIWire() const
Getter for wire within the layer.
unsigned short getISuperLayer() const
Getter for Super-Layer.
void addParam(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
B2Vector3< double > B2Vector3D
typedef for common usage with double
double sqrt(double a)
sqrt for double
double ClosestApproach(const B2Vector3D &bwp, const B2Vector3D &fwp, const B2Vector3D &posIn, const B2Vector3D &posOut, B2Vector3D &hitPosition, B2Vector3D &wirePosition)
Returns a closest distance between a track and a wire.
Abstract base class for different kinds of events.
Structure for saving the signal information.
Structure for saving the x-talk information.