Belle II Software development
CDCDigitizerModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <cdc/modules/cdcDigitizer/CDCDigitizerModule.h>
10#include <cdc/modules/cdcDigitizer/EDepInGas.h>
11#include <cdc/utilities/ClosestApproach.h>
12
13#include <framework/datastore/RelationArray.h>
14//#include <framework/datastore/RelationIndex.h>
15#include <framework/gearbox/Unit.h>
16#include <framework/logging/Logger.h>
17
18#include <TRandom.h>
19#include <map>
20#include <utility>
21
22using namespace std;
23using namespace Belle2;
24using namespace CDC;
25
26// register module
27REG_MODULE(CDCDigitizer);
32 m_tdcResol(0.9825), m_driftV(4.0e-3),
33 m_driftVInv(250.0), m_propSpeedInv(27.25), m_align(true)
34{
35 // Set description
36 setDescription("Creates CDCHits from CDCSimHits.");
38
39 // Add parameters
40 // I/O
41 addParam("InputCDCSimHitsName", m_inputCDCSimHitsName, "Name of input array. Should consist of CDCSimHits.", string(""));
42 addParam("OutputCDCHitsName", m_outputCDCHitsName, "Name of output array. Will consist of CDCHits.", string(""));
43 addParam("OutputCDCHitsName4Trg", m_outputCDCHitsName4Trg,
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"));
47
48 //Relations
49 addParam("MCParticlesToCDCSimHitsName", m_MCParticlesToSimHitsName,
50 "Name of relation between MCParticles and CDCSimHits used", string(""));
51 addParam("CDCSimHistToCDCHitsName", m_SimHitsTOCDCHitsName,
52 "Name of relation between the CDCSimHits and the CDCHits used", string(""));
53 addParam("OptionalFirstMCParticlesToHitsName", m_OptionalFirstMCParticlesToHitsName,
54 "Optional name of relation between the first MCParticles and CDCHits used", string("FirstMatchedParticles"));
55 addParam("OptionalAllMCParticlesToHitsName", m_OptionalAllMCParticlesToHitsName,
56 "Optional name of relation between all MCParticles and CDCHits used", string("AllMatchedParticles"));
57
58
59 //Parameters for Digitization
60 addParam("UseSimpleDigitization", m_useSimpleDigitization,
61 "If true, a simple x-t with a constant velocity is used for the drift-length to -time conversion", false);
62
63 //float Gauss Parameters
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);
69
70 //Switch to control smearing
71 addParam("DoSmearing", m_doSmearing,
72 "If false, drift length will not be smeared.", true);
73
74 addParam("DoAlphaCorrection", m_alphaCorrection,
75 "If true, some hits will be removed for better charge asymmetry simulation.", false);
76
77 addParam("TrigTimeJitter", m_trigTimeJitter,
78 "Magnitude (w) of trigger timing jitter (ns). The trigger timing is randuminzed uniformly in a time window of [-w/2, +w/2].",
79 0.);
80 //Switches to control time information handling
81 addParam("AddTimeWalk", m_addTimeWalk, "A switch for time-walk (pulse-heght dep. delay); true: on; false: off", true);
82 addParam("AddInWirePropagationDelay", m_addInWirePropagationDelay,
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);
84 addParam("AddInWirePropagationDelay4Bg", m_addInWirePropagationDelay4Bg,
85 "The same switch but for beam bg. hits.", true);
86 addParam("AddTimeOfFlight", m_addTimeOfFlight,
87 "A switch used to control adding time of flight into the final drift time or not; this is for signal hits.", true);
88 addParam("AddTimeOfFlight4Bg", m_addTimeOfFlight4Bg,
89 "The same switch but for beam bg. hits.", true);
90 addParam("OutputNegativeDriftTime", m_outputNegativeDriftTime, "Output hits with negative drift time", true);
91 addParam("Output2ndHit", m_output2ndHit,
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.",
93 false);
94 //Switch to control sense wire sag
95 addParam("CorrectForWireSag", m_correctForWireSag,
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.",
97 true);
98 //Switch for negative-t0 wires
99 addParam("TreatNegT0WiresAsGood", m_treatNegT0WiresAsGood, "Treat wires with negative t0 (calibrated) as good wire4s.", true);
100
101 //Threshold
102 addParam("TDCThreshold4Outer", m_tdcThreshold4Outer,
103 "TDC threshold (dE in eV) for Layers#8-56. The value corresponds to He-C2H6 gas", 250.);
104 addParam("TDCThreshold4Inner", m_tdcThreshold4Inner,
105 "Same as TDCThreshold4Outer but for Layers#0-7,", 150.);
106 addParam("EDepInGasMode", m_eDepInGasMode,
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",
108 0);
109
110 //ADC Threshold
111 addParam("ADCThreshold", m_adcThreshold,
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",
115 500.);
116 addParam("tMaxInner", m_tMaxInner, "Upper edge of time window in ns for the small-cell layers; valid only for UseDB4FEE=false",
117 300.);
118 // The following doesn't make any sense. The only reasonable steerable would be a switch to decide if the jitter shall be
119 // activated. Then there has to be event by event jitter.
120 /* addParam("EventTime", m_eventTime,
121 "It is a timing of event, which includes a time jitter due to the trigger system, set in ns", float(0.0));*/
122
123 //Switch for database
124 addParam("UseDB4FEE", m_useDB4FEE, "Fetch and use FEE params. from database or not", true);
125 addParam("UseDB4EDepToADC", m_useDB4EDepToADC, "Uuse edep-to-ADC conversion params. from database or not", true);
126 addParam("UseDB4RunGain", m_useDB4RunGain, "Fetch and use run gain from database or not", true);
127 addParam("OverallGainFactor", m_overallGainFactor, "Overall gain factor for adjustment", 1.0);
128
129 //Switch for synchronization
130 addParam("Synchronization", m_synchronization, "Synchronize timing with other sub-detectors", m_synchronization);
131 addParam("Randomization", m_randomization, "Randomize timing with other sub-detectors; valid only for Synchronization=false",
133 addParam("OffsetForGetTriggerBin", m_offsetForTriggerBin, "Input to getCDCTriggerBin(offset), either of 0,1,2 or 3",
135 addParam("TrgTimingOffsetInCount", m_trgTimingOffsetInCount,
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",
138 addParam("ShiftOfTimeWindowIn32Count", m_shiftOfTimeWindowIn32Count,
139 "Shift of time window in 32count for synchronization (L1 timing=0)", m_shiftOfTimeWindowIn32Count);
140
141 //Some FEE params.
142 addParam("TDCThresholdOffset", m_tdcThresholdOffset, "Offset for TDC (digital) threshold (mV)", 3828.);
143 addParam("AnalogGain", m_analogGain, "Analog gain (V/pC)", 1.09);
144 addParam("DigitalGain", m_digitalGain, "Digital gain (V/pC)", 7.);
145 addParam("ADCBinWidth", m_adcBinWidth, "ADC bin width (mV)", 2.);
146
147 addParam("AddFudgeFactorForSigma", m_addFudgeFactorForSigma,
148 "Additional fudge factor for space resol. (common to all cells)", 1.);
149 addParam("SpaceChargeEffect", m_spaceChargeEffect, "Switch for space charge effect", true);
150 addParam("DegOfSPEOnThreshold", m_degOfSPEOnThreshold,
151 "Degree of space charge effect on timing threshold; specify the range [0,1]; =1: full effect on threshold; =0: no effect",
153
154 addParam("AddXTalk", m_addXTalk, "A switch for crosstalk; true: on; false: off", true);
155 addParam("Issue2ndHitWarning", m_issue2ndHitWarning, "=true: issue a warning when a 2nd TDC hit is found.", true);
156 addParam("IncludeEarlyXTalks", m_includeEarlyXTalks, "=true: include earlier x-talks as well than the signal hit in question.",
157 true);
158 addParam("DebugLevel", m_debugLevel, "Debug level; 20-29 are usable.", 20);
159 addParam("DebugLevel4XTalk", m_debugLevel4XTalk, "Debug level for crosstalk; 20-29 are usable.", 21);
160
161 //Gain smearing
162 addParam("GasGainSmearing", m_gasGainSmearing, "Switch for gas gain smearing for ADC simulation; true: on; false: off",
164 addParam("EffWForGasGainSmearing", m_effWForGasGainSmearing,
165 "Effective energy (keV) needed for one electron production for gas gain smearing; average for alpha- and beta-sources.",
167 addParam("ThetaOfPolyaFunction", m_thetaOfPolya, "Theta of Polya function for gas gain smearing", m_thetaOfPolya);
168 addParam("ExtraADCSmearing", m_extraADCSmearing, "Switch for extra ADC smearing; true: on; false: off", m_extraADCSmearing);
169 // addParam("SigmaForExtraADCSmearing", m_sigmaForExtraADCSmearing, "Gaussian sigma for extra ADC smearing; specify range [0,1]", m_sigmaForExtraADCSmearing);
170
171 // Switch for optional relations
172 addParam("MatchAllMCParticles", m_matchAllMCParticles, "Switch to store all MCRelations that produced a SimHit", false);
173 addParam("MatchFirstMCParticles", m_matchFirstMCParticles,
174 "Switch to store all MCRelations for the first three SimHits instead of only the first", false);
175
176#if defined(CDC_DEBUG)
177 cout << " " << endl;
178 cout << "CDCDigitizer constructor" << endl;
179#endif
180}
181
183{
185 m_simClockState.isOptional();
186
187 // Register the arrays in the DataStore, that are to be added in this module.
188 m_cdcHits.registerInDataStore(m_outputCDCHitsName);
189 m_simHits.registerRelationTo(m_cdcHits);
190 m_mcParticles.registerRelationTo(m_cdcHits);
193
194 // Arrays for trigger.
195 m_cdcHits4Trg.registerInDataStore(m_outputCDCHitsName4Trg);
196 m_simHits.registerRelationTo(m_cdcHits4Trg);
197 m_mcParticles.registerRelationTo(m_cdcHits4Trg);
198
200 CDCGeometryPar& cdcgp = *m_cdcgp;
204 m_driftV = cdcgp.getNominalDriftV();
205 m_driftVInv = 1. / m_driftV;
206 m_propSpeedInv = 1. / cdcgp.getNominalPropSpeed();
208 m_totalFudgeFactor = m_cdcgp->getFudgeFactorForSigma(2);
210 B2DEBUG(m_debugLevel, "totalFugeF in Digi= " << m_totalFudgeFactor);
211 /*
212 m_fraction = 1.0;
213 m_resolution1 = cdcgp.getNominalSpaceResol();
214 m_resolution2 = 0.;
215 m_mean1 = 0.;
216 m_mean2 = 0.;
217 */
218
219 if (m_useDB4FEE) {
221 if ((*m_fEElectronicsFromDB).isValid()) {
222 (*m_fEElectronicsFromDB).addCallback(this, &CDCDigitizerModule::setFEElectronics);
224 } else {
225 B2FATAL("CDCDigitizer:: CDCFEElectronics not valid !");
226 }
227 }
228
229 /*
230 if (m_useDB4EDepToADC) {
231 m_eDepToADCConversionsFromDB = new DBObjPtr<CDCEDepToADCConversions>;
232 if ((*m_eDepToADCConversionsFromDB).isValid()) {
233 (*m_eDepToADCConversionsFromDB).addCallback(this, &CDCDigitizerModule::setEDepToADCConversions);
234 setEDepToADCConversions();
235 } else {
236 B2FATAL("CDCDigitizer:: CDCEDepToADCConversions not valid !");
237 }
238 }
239 */
240
241 if (m_useDB4RunGain) {
243 if ((*m_runGainFromDB).isValid()) {
244 (*m_runGainFromDB).addCallback(this, &CDCDigitizerModule::setSemiTotalGain);
245 } else {
246 B2FATAL("CDCDedxRunGain invalid!");
247 }
248
250 if ((*m_gain0FromDB).isValid()) {
251 (*m_gain0FromDB).addCallback(this, &CDCDigitizerModule::setSemiTotalGain);
252 } else {
253 B2FATAL("CDCDedxScaleFactor invalid!");
254 }
255
257 if ((*m_wireGainFromDB).isValid()) {
258 (*m_wireGainFromDB).addCallback(this, &CDCDigitizerModule::setSemiTotalGain);
260 } else {
261 B2FATAL("CDCDedxWireGain invalid!");
262 }
263 }
264
265 if (m_addXTalk) {
267 if ((*m_xTalkFromDB).isValid()) {
268 } else {
269 B2FATAL("CDCCrossTalkLibrary invalid!");
270 }
271 }
272
274 if ((*m_corrToThresholdFromDB).isValid()) {
275 } else {
276 B2FATAL("CDCCorrToThresholds invalid!");
277 }
278
279 if (m_alphaCorrection) {
280 if (!m_alphaScaleFactorsFromDB.isValid()) {
281 B2FATAL("CDCAlphaScaleFactorForAsymmetry invalid!");
282 }
283 }
284
285#if defined(CDC_DEBUG)
286 cout << " " << endl;
287 cout << "CDCDigitizer initialize" << endl;
288 // cout << "m_tdcOffset= " << m_tdcOffset << endl;
289 cout << "m_tdcBinWidth= " << m_tdcBinWidth << endl;
290 cout << "m_tdcResol= " << m_tdcResol << endl;
291 cout << "m_driftV= " << m_driftV << endl;
292 cout << "m_driftVInv= " << m_driftVInv << endl;
293 cout << "m_propSpeedInv= " << m_propSpeedInv << endl;
294 /*
295 cout << "m_fraction= " << m_fraction << endl;
296 cout << "m_resolution1= " << m_resolution1 << endl;
297 cout << "m_resolution2= " << m_resolution2 << endl;
298 cout << "m_mean1= " << m_mean1 << endl;
299 cout << "m_mean2= " << m_mean2 << endl;
300 */
301#endif
302
303 if (m_useDB4EDepToADC) {
304 ushort firstLayerOffset = m_cdcgp->getOffsetOfFirstLayer();
305 if (m_cdcgp->getEDepToADCMainFactor(firstLayerOffset, 0) == 0.) {
306 B2FATAL("CDCEDepToADCConversion payloads are unavailable!");
307 }
308 }
309
310 // Set timing sim. mode
311 if (m_useDB4FEE) {
312 if (m_synchronization) { // synchronization
313 m_tSimMode = 0;
314 } else {
315 if (m_randomization) { // radomization
316 m_tSimMode = 1;
317 } else {
318 m_tSimMode = 2; // old sim.
319 }
320 }
321 } else {
322 m_tSimMode = 3; // old sim. w/o relying on fee db
323 }
324 B2DEBUG(m_debugLevel, "timing sim. mode= " << m_tSimMode);
325 if (m_tSimMode < 0 || m_tSimMode > 3) B2FATAL("invalid timing sim. mode= " << m_tSimMode);
326}
327
329{
330 // Get SimHit array, MCParticle array, and relation between the two.
331 RelationArray mcParticlesToCDCSimHits(m_mcParticles, m_simHits); //RelationArray created by CDC SensitiveDetector
332
333
334 //--- Start Digitization --------------------------------------------------------------------------------------------
335 // Merge the hits in the same cell and save them into CDC signal map.
336
337 // Define signal map
338 map<WireID, SignalInfo> signalMap;
340 // Define adc map
341 map<WireID, unsigned short> adcMap;
343 // map<WireID, double> adcMap;
344 // map<WireID, double>::iterator iterADCMap;
345
346 // signal map for trigger
347 map<pair<WireID, unsigned>, SignalInfo> signalMapTrg;
348 map<pair<WireID, unsigned>, SignalInfo>::iterator iterSignalMapTrg;
349
350 // signal map for all MCParticles: Wire <-> MCArrayIndex
351 map<WireID, std::set<int>> particleMap;
352 map<WireID, std::set<int>>::iterator iterParticleMap;
353
354
355 // Set time window per event
356 if (m_tSimMode == 0 || m_tSimMode == 1) {
357 int trigBin = 0;
358 if (m_simClockState.isValid()) {
359 trigBin = m_simClockState->getCDCTriggerBin(m_offsetForTriggerBin);
360 } else {
361 if (m_tSimMode == 0) {
362 B2DEBUG(m_debugLevel, "SimClockState unavailable so switched the mode from synchro to random.");
363 m_tSimMode = 1;
364 }
365 trigBin = gRandom->Integer(4);
366 }
367 if (trigBin < 0 || trigBin > 3) B2ERROR("Invalid trigger bin; must be an integer [0,3]!");
368 unsigned short offs = 8 * trigBin + m_trgTimingOffsetInCount;
369 B2DEBUG(m_debugLevel, "tSimMode,trigBin,offs= " << m_tSimMode << " " << trigBin << " " << offs);
370
371 //TODO: simplify the following 7 lines and setFEElectronics()
372 for (unsigned short bd = 1; bd < c_nBoards; ++bd) {
373 const short tMaxInCount = 32 * (m_shiftOfTimeWindowIn32Count - m_trgDelayInCount[bd]) - offs;
374 const short tMinInCount = tMaxInCount - 32 * m_widthOfTimeWindowInCount[bd];
375 B2DEBUG(m_debugLevel, bd << " " << tMinInCount << " " << tMaxInCount);
376 m_uprEdgeOfTimeWindow[bd] = m_tdcBinWidth * tMaxInCount;
377 m_lowEdgeOfTimeWindow[bd] = m_tdcBinWidth * tMinInCount;
378 }
379 }
380
381 // Set trigger timing jitter for this event
382 double trigTiming = m_trigTimeJitter == 0. ? 0. : m_trigTimeJitter * (gRandom->Uniform() - 0.5);
383 // std::cout << "trigTiming= " << trigTiming << std::endl;
384 // Loop over all hits
385 int nHits = m_simHits.getEntries();
386 B2DEBUG(m_debugLevel, "Number of CDCSimHits in the current event: " << nHits);
387 for (int iHits = 0; iHits < nHits; ++iHits) {
388 // Get a hit
389 m_aCDCSimHit = m_simHits[iHits];
390
391 // Hit geom. info
392 m_wireID = m_aCDCSimHit->getWireID();
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.");
395 }
396 // B2DEBUG(29, "Encoded wire number of current CDCSimHit: " << m_wireID);
397
398 m_posFlag = m_aCDCSimHit->getLeftRightPassageRaw();
399 m_boardID = m_cdcgp->getBoardID(m_wireID);
400 // B2DEBUG(29, "m_boardID= " << m_boardID);
401 m_posWire = m_aCDCSimHit->getPosWire();
402 m_posTrack = m_aCDCSimHit->getPosTrack();
403 m_momentum = m_aCDCSimHit->getMomentum();
404 m_flightTime = m_aCDCSimHit->getFlightTime();
405 m_globalTime = m_aCDCSimHit->getGlobalTime();
406 m_driftLength = m_aCDCSimHit->getDriftLength() * Unit::cm;
407
408 //include alignment effects
409 //basically align flag should be always on since on/off is controlled by the input alignment.xml file itself.
410 m_align = true;
411
412 B2Vector3D bwpAlign = m_cdcgp->wireBackwardPosition(m_wireID, CDCGeometryPar::c_Aligned);
413 B2Vector3D fwpAlign = m_cdcgp->wireForwardPosition(m_wireID, CDCGeometryPar::c_Aligned);
414
415 B2Vector3D bwp = m_cdcgp->wireBackwardPosition(m_wireID);
416 B2Vector3D fwp = m_cdcgp->wireForwardPosition(m_wireID);
417
418 //skip correction for wire-position alignment if unnecessary
419 if ((bwpAlign - bwp).Mag() == 0. && (fwpAlign - fwp).Mag() == 0.) m_align = false;
420 // std::cout << "a m_align= " << m_align << std::endl;
421
423
424 bwp = bwpAlign;
425 fwp = fwpAlign;
426
428 double zpos = m_posWire.Z();
429 double bckYSag = bwp.Y();
430 double forYSag = fwp.Y();
431
432 // CDCGeometryPar::EWirePosition set = m_align ?
433 // CDCGeometryPar::c_Aligned : CDCGeometryPar::c_Base;
434 CDCGeometryPar::EWirePosition set = CDCGeometryPar::c_Aligned;
435 const int layerID = m_wireID.getICLayer();
436 const int wireID = m_wireID.getIWire();
437 m_cdcgp->getWireSagEffect(set, layerID, wireID, zpos, bckYSag, forYSag);
438 bwp.SetY(bckYSag);
439 fwp.SetY(forYSag);
440 }
441
442 const B2Vector3D L = 5. * m_momentum.Unit(); //(cm) tentative
443 B2Vector3D posIn = m_posTrack - L;
444 B2Vector3D posOut = m_posTrack + L;
445 B2Vector3D posTrack = m_posTrack;
446 B2Vector3D posWire = m_posWire;
447
448 // m_driftLength = m_cdcgp->ClosestApproach(bwp, fwp, posIn, posOut, posTrack, posWire);
449 m_driftLength = ClosestApproach(bwp, fwp, posIn, posOut, posTrack, posWire);
450 // std::cout << "base-dl, sag-dl, diff= " << m_aCDCSimHit->getDriftLength() <<" "<< m_driftLength <<" "<< m_driftLength - m_aCDCSimHit->getDriftLength() << std::endl;
451 m_posTrack = posTrack;
452 m_posWire = posWire;
453
454 double deltaTime = 0.; //tentative (probably ok...)
455 // double deltaTime = (posTrack - m_posTrack).Mag() / speed;
456 m_flightTime += deltaTime;
457 m_globalTime += deltaTime;
458 m_posFlag = m_cdcgp->getNewLeftRightRaw(m_posWire, m_posTrack, m_momentum);
459 }
460
461 // Calculate measurement time.
462 // Smear drift length
463 double hitDriftLength = m_driftLength;
464 double dDdt = getdDdt(hitDriftLength);
465 if (m_doSmearing) {
466 hitDriftLength = smearDriftLength(hitDriftLength, dDdt);
467 }
468
469 //set flags
470 bool addTof = m_addTimeOfFlight4Bg;
471 bool addDelay = m_addInWirePropagationDelay4Bg;
472 if (m_aCDCSimHit->getBackgroundTag() == 0) {
473 addTof = m_addTimeOfFlight;
475 }
476 double hitDriftTime = getDriftTime(hitDriftLength, addTof, addDelay);
477
478 //add randamized event time for a beam bg. hit
479 if (m_aCDCSimHit->getBackgroundTag() != 0) {
480 hitDriftTime += m_globalTime - m_flightTime;
481 }
482
483 //add trigger timing jitter
484 hitDriftTime += trigTiming;
485
486 //apply time window cut
487 double tMin = m_tMin;
488 double tMax = m_tMaxOuter;
489 if (m_wireID.getISuperLayer() == 0) tMax = m_tMaxInner;
490 if (m_tSimMode <= 2) {
493 }
494 if (hitDriftTime < tMin || hitDriftTime > tMax) continue;
495
496 //Sum ADC count
497 const double stepLength = m_aCDCSimHit->getStepLength() * Unit::cm;
498 const double costh = m_momentum.Z() / m_momentum.Mag();
499 double hitdE = m_aCDCSimHit->getEnergyDep();
500 if (m_cdcgp->getMaterialDefinitionMode() != 2) { // for non wire-by-wire mode
501 static EDepInGas& edpg = EDepInGas::getInstance();
502 hitdE = edpg.getEDepInGas(m_eDepInGasMode, m_aCDCSimHit->getPDGCode(), m_momentum.Mag(), stepLength, hitdE);
503 }
504
505 double convFactorForThreshold = 1;
506 //TODO: modify the following function so that it can output timing signal in Volt in future
507 unsigned short adcCount = 0;
508 makeSignalsAfterShapers(m_wireID, hitdE, stepLength, costh, adcCount, convFactorForThreshold);
509 const unsigned short adcTh = m_useDB4FEE ? m_adcThresh[m_boardID] : m_adcThreshold;
510 // B2DEBUG(29, "adcTh,adcCount,convFactorForThreshold= " << adcTh <<" "<< adcCount <<" "<< convFactorForThreshold);
511 if (adcCount < adcTh) adcCount = 0;
512 iterADCMap = adcMap.find(m_wireID);
513 if (iterADCMap == adcMap.end()) {
514 adcMap.insert(make_pair(m_wireID, adcCount));
515 // adcMap.insert(make_pair(m_wireID, hitdE));
516 } else {
517 iterADCMap->second += adcCount;
518 // iterADCMap->second += hitdE;
519 }
520
521 //Apply energy threshold
522 // If hitdE < dEThreshold, the hit is ignored
523 double dEThreshold = 0.;
525 dEThreshold = m_tdcThresh[m_boardID] / convFactorForThreshold * Unit::keV;
526 } else {
527 dEThreshold = (m_wireID.getISuperLayer() == 0) ? m_tdcThreshold4Inner : m_tdcThreshold4Outer;
528 dEThreshold *= Unit::eV;
529 }
530 dEThreshold *= (*m_corrToThresholdFromDB)->getParam(m_wireID.getICLayer());
531 B2DEBUG(m_debugLevel, "hitdE,dEThreshold,driftLength " << hitdE << " " << dEThreshold << " " << hitDriftLength);
532
533 if (hitdE < dEThreshold) {
534 B2DEBUG(m_debugLevel, "Below Ethreshold: " << hitdE << " " << dEThreshold);
535 continue;
536 }
537
538 // add one hit per trigger time window to the trigger signal map
539 unsigned short trigWindow = floor((hitDriftTime - tMin) * m_tdcBinWidthInv / 32);
540 iterSignalMapTrg = signalMapTrg.find(make_pair(m_wireID, trigWindow));
541 if (iterSignalMapTrg == signalMapTrg.end()) {
542 // signalMapTrg.insert(make_pair(make_pair(m_wireID, trigWindow),
543 // SignalInfo(iHits, hitDriftTime, hitdE)));
544 signalMapTrg.insert(make_pair(make_pair(m_wireID, trigWindow),
545 SignalInfo(iHits, hitDriftTime, adcCount)));
546 } else {
547 if (hitDriftTime < iterSignalMapTrg->second.m_driftTime) {
548 iterSignalMapTrg->second.m_driftTime = hitDriftTime;
549 iterSignalMapTrg->second.m_simHitIndex = iHits;
550 }
551 // iterSignalMapTrg->second.m_charge += hitdE;
552 iterSignalMapTrg->second.m_charge += adcCount;
553 }
554
555 // Reject totally-dead wire; to be replaced by isDeadWire() in future
556 // N.B. The following lines for badwire must be after the above lines for trigger because badwires are different between trigger and tracking.
557 // Badwires for trigger are taken into account separately in the tsim module
558 if (m_cdcgp->isBadWire(m_wireID)) {
559 // std::cout<<"badwire= " << m_wireID.getICLayer() <<" "<< m_wireID.getIWire() << std::endl;
560 continue;
561 }
562 // Reject partly-dead wire as well
563 double eff = 1.;
564 if (m_cdcgp->isDeadWire(m_wireID, eff)) {
565 // std::cout << "wid,eff= " << m_wireID << " " << eff << std::endl;
566 if (eff < gRandom->Uniform()) continue;
567 }
568
569 // For TOT simulation, calculate drift length from In to the wire, and Out to the wire. The calculation is apprximate ignoring wire sag (this would be ok because TOT simulation is not required to be so accurate).
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;
578
579 double dx = m_aCDCSimHit->getPosIn().X() - a;
580 double dy = m_aCDCSimHit->getPosIn().Y() - b;
581 double dz = m_aCDCSimHit->getPosIn().Z() - c;
582 double sub = l * dx + m * dy + n * dz;
583 const double driftLFromIn = sqrt(dx * dx + dy * dy + dz * dz - sub * sub);
584
585 dx = m_aCDCSimHit->getPosOut().X() - a;
586 dy = m_aCDCSimHit->getPosOut().Y() - b;
587 dz = m_aCDCSimHit->getPosOut().Z() - c;
588 sub = l * dx + m * dy + n * dz;
589 const double driftLFromOut = sqrt(dx * dx + dy * dy + dz * dz - sub * sub);
590
591 const double maxDriftL = std::max(driftLFromIn, driftLFromOut);
592 const double minDriftL = m_driftLength;
593 B2DEBUG(m_debugLevel, "driftLFromIn= " << driftLFromIn << " driftLFromOut= " << driftLFromOut << " minDriftL= " << minDriftL <<
594 " maxDriftL= "
595 <<
596 maxDriftL << "m_driftLength= " << m_driftLength);
597
598 iterSignalMap = signalMap.find(m_wireID);
599
601 iterParticleMap = particleMap.find(m_wireID);
602 RelationVector<MCParticle> rels = m_aCDCSimHit->getRelationsFrom<MCParticle>();
603
604 int mcIndex = -1;
605 if (rels.size() != 0) {
606 if (rels.weight(0) > 0) {
607 const MCParticle* mcparticle = rels[0];
608 mcIndex = int(mcparticle->getIndex());
609 }
610 }
611
612 if (mcIndex >= 0) {
613 if (iterParticleMap == particleMap.end()) {
614 std::set<int> vecmc = {mcIndex};
615 particleMap.insert(make_pair(m_wireID, vecmc));
616 } else {
617 iterParticleMap->second.insert(mcIndex);
618 }
619 }
620 }
621
622 if (iterSignalMap == signalMap.end()) {
623 // new entry
624 // signalMap.insert(make_pair(m_wireID, SignalInfo(iHits, hitDriftTime, hitdE)));
625 signalMap.insert(make_pair(m_wireID, SignalInfo(iHits, hitDriftTime, adcCount, maxDriftL, minDriftL)));
626 B2DEBUG(m_debugLevel, "Creating new Signal with encoded wire number: " << m_wireID);
627 } else {
628 // ... smallest drift time has to be checked, ...
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;
645 }
646 // ... total charge has to be updated.
647 // iterSignalMap->second.m_charge += hitdE;
648 iterSignalMap->second.m_charge += adcCount;
649
650 // set max and min driftLs
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);
655 }
656
657 } // end loop over SimHits.
658
659 //--- Now Store the results into CDCHits and
660 // create corresponding relations between SimHits and CDCHits.
661
662 unsigned int iCDCHits = 0;
663 RelationArray cdcSimHitsToCDCHits(m_simHits, m_cdcHits); //SimHit<->CDCHit
664 RelationArray mcParticlesToCDCHits(m_mcParticles, m_cdcHits); //MCParticle<->CDCHit
665
666 for (iterSignalMap = signalMap.begin(); iterSignalMap != signalMap.end(); ++iterSignalMap) {
667
668 //add time-walk (here for simplicity)
669 // unsigned short adcCount = getADCCount(iterSignalMap->second.m_charge);
670 // unsigned short adcCount = iterSignalMap->second.m_charge;
671 iterADCMap = adcMap.find(iterSignalMap->first);
672 unsigned short adcCount = iterADCMap != adcMap.end() ? iterADCMap->second : 0;
673 /*
674 unsigned short adcCount = 0;
675 if (iterADCMap != adcMap.end()) {
676 adcCount = getADCCount(iterSignalMap->first, iterADCMap->second, 1., 0.);
677 unsigned short boardID = m_cdcgp->getBoardID(iterSignalMap->first);
678 // B2DEBUG(29, "boardID= " << boardID);
679 const unsigned short adcTh = m_useDB4FEE ? m_adcThresh[boardID] : m_adcThreshold;
680 if (adcCount < adcTh) adcCount = 0;
681 }
682 */
683
684 if (m_addTimeWalk) {
685 B2DEBUG(m_debugLevel, "timewalk= " << m_cdcgp->getTimeWalk(iterSignalMap->first, adcCount));
686 iterSignalMap->second.m_driftTime += m_cdcgp->getTimeWalk(iterSignalMap->first, adcCount);
687 }
688
689 //remove negative drift time (TDC) upon request
691 iterSignalMap->second.m_driftTime < 0.) {
692 continue;
693 }
694
695 //apply correction on alpha, to calibrate the charge asymmetry at hit-level
696 if (m_alphaCorrection) {
697 int iHits = iterSignalMap->second.m_simHitIndex;
698 m_aCDCSimHit = m_simHits[iHits];
699 m_posWire = m_aCDCSimHit->getPosWire();
700 m_momentum = m_aCDCSimHit->getMomentum();
701 int iLayer = m_aCDCSimHit->getWireID().getICLayer();
702 double alpha = m_cdcgp->getAlpha(m_posWire, m_momentum);
703
704 double Scale = m_alphaScaleFactorsFromDB->getScaleFactor(iLayer, alpha);
705 double random = gRandom->Uniform();
706 if ((Scale < 1) && (alpha > 0)) {
707 if (random > Scale) continue ; // remove this hit
708 }
709 if ((Scale > 1) && (alpha < 0)) {
710 if (random > 1. / Scale) continue ; // remove this hit
711 }
712 }
713
714 //N.B. No bias (+ or -0.5 count) is introduced on average in digitization by the real TDC (info. from KEK electronics division). So round off (t0 - drifttime) below.
715 unsigned short tdcCount = static_cast<unsigned short>((getPositiveT0(iterSignalMap->first) - iterSignalMap->second.m_driftTime) *
716 m_tdcBinWidthInv + 0.5);
717
718 //calculate tot; hard-coded currently
719 double deltaDL = iterSignalMap->second.m_maxDriftL - iterSignalMap->second.m_minDriftL;
720 if (deltaDL < 0.) {
721 B2DEBUG(m_debugLevel, "negative deltaDL= " << deltaDL);
722 deltaDL = 0.;
723 }
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]));
726 if (m_adcThresh[boardID] > 0) {
727 tot = std::min(static_cast<int>(tot), static_cast<int>(adcCount / m_adcThresh[boardID]));
728 }
729
730 CDCHit* firstHit = m_cdcHits.appendNew(tdcCount, adcCount, iterSignalMap->first, 0, tot);
731 // std::cout <<"firsthit?= " << firstHit->is2ndHit() << std::endl;
732 //set a relation: CDCSimHit -> CDCHit
733 cdcSimHitsToCDCHits.add(iterSignalMap->second.m_simHitIndex, iCDCHits);
734
735 //set a relation: MCParticle -> CDCHit
736 RelationVector<MCParticle> rels = m_simHits[iterSignalMap->second.m_simHitIndex]->getRelationsFrom<MCParticle>();
737 if (rels.size() != 0) {
738 //assumption: only one MCParticle
739 const MCParticle* mcparticle = rels[0];
740 double weight = rels.weight(0);
741 mcparticle->addRelationTo(firstHit, weight);
742 }
743
744 // Set relations to all particles that created a SimHit
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) {
750 // set all relations
751 int idx = *it;
752 MCParticle* part = m_mcParticles[idx - 1];
754 }
755 }
756 }
757
758 //set all relations to first hit if requested but dont create additional hits!
759 // relation 1
760 if (m_matchFirstMCParticles > 0) {
761 if (iterSignalMap->second.m_simHitIndex >= 0) {
762 RelationVector<MCParticle> rels1 = m_simHits[iterSignalMap->second.m_simHitIndex]->getRelationsFrom<MCParticle>();
763 if (rels1.size() != 0) {
764 //assumption: only one MCParticle
765 const MCParticle* mcparticle = rels1[0];
766 double weight = rels1.weight(0);
767 mcparticle->addRelationTo(firstHit, weight, m_OptionalFirstMCParticlesToHitsName);
768 }
769 }
770
771 // relation 2
772 if (iterSignalMap->second.m_simHitIndex2 >= 0) {
773 RelationVector<MCParticle> rels2 = m_simHits[iterSignalMap->second.m_simHitIndex2]->getRelationsFrom<MCParticle>();
774 if (rels2.size() != 0) {
775 //assumption: only one MCParticle
776 const MCParticle* mcparticle = rels2[0];
777 double weight = rels2.weight(0);
778 mcparticle->addRelationTo(firstHit, weight, m_OptionalFirstMCParticlesToHitsName);
779 }
780 }
781
782 // relation 3
783 if (iterSignalMap->second.m_simHitIndex3 >= 0) {
784 RelationVector<MCParticle> rels3 = m_simHits[iterSignalMap->second.m_simHitIndex3]->getRelationsFrom<MCParticle>();
785 if (rels3.size() != 0) {
786 //assumption: only one MCParticle
787 const MCParticle* mcparticle = rels3[0];
788 double weight = rels3.weight(0);
789 mcparticle->addRelationTo(firstHit, weight, m_OptionalFirstMCParticlesToHitsName);
790 }
791 }
792
793
794 }
795
796 //Set 2nd-hit related things if it exists
797 if (m_output2ndHit && iterSignalMap->second.m_simHitIndex2 >= 0) {
798 unsigned short tdcCount2 = static_cast<unsigned short>((getPositiveT0(iterSignalMap->first) - iterSignalMap->second.m_driftTime2) *
799 m_tdcBinWidthInv + 0.5);
800 if (tdcCount2 != tdcCount) {
801 CDCHit* secondHit = m_cdcHits.appendNew(tdcCount2, adcCount, iterSignalMap->first, 0, tot);
802 secondHit->set2ndHitFlag();
803 secondHit->setOtherHitIndices(firstHit);
804 // std::cout <<"2ndhit?= " << secondHit->is2ndHit() << std::endl;
805 // std::cout <<"1st-otherhitindex= " << firstHit->getOtherHitIndex() << std::endl;
806 // std::cout <<"2nd-otherhitindex= " << secondHit->getOtherHitIndex() << std::endl;
807 // secondHit->setOtherHitIndex(firstHit->getArrayIndex());
808 // firstHit->setOtherHitIndex(secondHit->getArrayIndex());
809 // std::cout <<"1st-otherhitindex= " << firstHit->getOtherHitIndex() << std::endl;
810 // std::cout <<"2nd-otherhitindex= " << secondHit->getOtherHitIndex() << std::endl;
811
812 //set a relation: CDCSimHit -> CDCHit
813 ++iCDCHits;
814 cdcSimHitsToCDCHits.add(iterSignalMap->second.m_simHitIndex2, iCDCHits);
815 // std::cout << "settdc2 " << firstHit->getTDCCount() << " " << secondHit->getTDCCount() << std::endl;
816
817 //set a relation: MCParticle -> CDCHit
818 rels = m_simHits[iterSignalMap->second.m_simHitIndex2]->getRelationsFrom<MCParticle>();
819 if (rels.size() != 0) {
820 //assumption: only one MCParticle
821 const MCParticle* mcparticle = rels[0];
822 double weight = rels.weight(0);
823 mcparticle->addRelationTo(secondHit, weight);
824 }
825 } else { //Check the 3rd hit when tdcCount = tdcCount2
826 // std::cout << "tdcCount1=2" << std::endl;
827 if (iterSignalMap->second.m_simHitIndex3 >= 0) {
828 unsigned short tdcCount3 = static_cast<unsigned short>((getPositiveT0(iterSignalMap->first) - iterSignalMap->second.m_driftTime3) *
829 m_tdcBinWidthInv + 0.5);
830 // std::cout << "tdcCount3= " << tdcCount3 << " " << tdcCount << std::endl;
831 if (tdcCount3 != tdcCount) {
832 CDCHit* secondHit = m_cdcHits.appendNew(tdcCount3, adcCount, iterSignalMap->first, 0, tot);
833 secondHit->set2ndHitFlag();
834 secondHit->setOtherHitIndices(firstHit);
835 // secondHit->setOtherHitIndex(firstHit->getArrayIndex());
836 // firstHit->setOtherHitIndex(secondHit->getArrayIndex());
837 // std::cout <<"2ndhit?= " << secondHit->is2ndHit() << std::endl;
838
839 //set a relation: CDCSimHit -> CDCHit
840 ++iCDCHits;
841 cdcSimHitsToCDCHits.add(iterSignalMap->second.m_simHitIndex3, iCDCHits);
842 // std::cout << "settdc3 " << firstHit->getTDCCount() << " " << secondHit->getTDCCount() << std::endl;
843
844 //set a relation: MCParticle -> CDCHit
845 rels = m_simHits[iterSignalMap->second.m_simHitIndex3]->getRelationsFrom<MCParticle>();
846 if (rels.size() != 0) {
847 //assumption: only one MCParticle
848 const MCParticle* mcparticle = rels[0];
849 double weight = rels.weight(0);
850 mcparticle->addRelationTo(secondHit, weight);
851 }
852 }
853 }
854 } //end of checking tdcCount 1=2 ?
855 } //end of 2nd hit setting
856
857 // std::cout <<"t0= " << m_cdcgp->getT0(iterSignalMap->first) << std::endl;
858 /* unsigned short tdcInCommonStop = static_cast<unsigned short>((m_tdcOffset - iterSignalMap->second.m_driftTime) * m_tdcBinWidthInv);
859 float driftTimeFromTDC = static_cast<float>(m_tdcOffset - (tdcInCommonStop + 0.5)) * m_tdcBinWidth;
860 std::cout <<"driftT bf digitization, TDC in common stop, digitized driftT = " << iterSignalMap->second.m_driftTime <<" "<< tdcInCommonStop <<" "<< driftTimeFromTDC << std::endl;
861 */
862 ++iCDCHits;
863 }
864
865 //Add crosstalk
866 if (m_addXTalk) addXTalk();
867
868 // Store the results with trigger time window in a separate array
869 // with corresponding relations.
870 for (iterSignalMapTrg = signalMapTrg.begin(); iterSignalMapTrg != signalMapTrg.end(); ++iterSignalMapTrg) {
871 /*
872 unsigned short adcCount = getADCCount(iterSignalMapTrg->first.first, iterSignalMapTrg->second.m_charge, 1., 0.);
873 unsigned short boardID = m_cdcgp->getBoardID(iterSignalMapTrg->first.first);
874 // B2DEBUG(29, "boardID= " << boardID);
875 const unsigned short adcTh = m_useDB4FEE ? m_adcThresh[boardID] : m_adcThreshold;
876 if (adcCount < adcTh) adcCount = 0;
877 */
878 // unsigned short adcCount = getADCCount(iterSignalMapTrg->second.m_charge);
879 unsigned short adcCount = iterSignalMapTrg->second.m_charge;
880 unsigned short tdcCount =
881 static_cast<unsigned short>((getPositiveT0(iterSignalMapTrg->first.first) -
882 iterSignalMapTrg->second.m_driftTime) * m_tdcBinWidthInv + 0.5);
883 const CDCHit* cdcHit = m_cdcHits4Trg.appendNew(tdcCount, adcCount, iterSignalMapTrg->first.first);
884
885 // relations
886 m_simHits[iterSignalMapTrg->second.m_simHitIndex]->addRelationTo(cdcHit);
887 RelationVector<MCParticle> rels = m_simHits[iterSignalMapTrg->second.m_simHitIndex]->getRelationsFrom<MCParticle>();
888 if (rels.size() != 0) {
889 //assumption: only one MCParticle
890 const MCParticle* mcparticle = rels[0];
891 double weight = rels.weight(0);
892 mcparticle->addRelationTo(cdcHit, weight);
893 }
894 }
895
896 /*
897 std::cout << " " << std::endl;
898 RelationIndex<MCParticle, CDCHit> mcp_to_hit(mcParticles, cdcHits);
899 if (!mcp_to_hit) B2FATAL("No MCParticle -> CDCHit relation founf!");
900 typedef RelationIndex<MCParticle, CDCHit>::Element RelationElement;
901 int ncdcHits = cdcHits.getEntries();
902 for (int j = 0; j < ncdcHits; ++j) {
903 for (const RelationElement& rel : mcp_to_hit.getElementsTo(cdcHits[j])) {
904 std::cout << j << " " << cdcHits[j]->is2ndHit() <<" "<< rel.from->getIndex() << " " << rel.weight << std::endl;
905 }
906 }
907 */
908}
909
910double CDCDigitizerModule::smearDriftLength(const double driftLength, const double dDdt)
911{
912 double mean = 0.;
913 double resolution;
914
916 if (gRandom->Uniform() <= m_fraction) {
917 mean = m_mean1;
918 resolution = m_resolution1;
919 } else {
920 mean = m_mean2;
921 resolution = m_resolution2;
922 }
923 } else {
924 const unsigned short leftRight = m_posFlag;
925 double alpha = m_cdcgp->getAlpha(m_posWire, m_momentum);
926 double theta = m_cdcgp->getTheta(m_momentum);
927 resolution = m_cdcgp->getSigma(driftLength, m_wireID.getICLayer(), leftRight, alpha, theta);
928 resolution *= m_totalFudgeFactor;
929 }
930
931 //subtract resol. due to digitization, which'll be added later in the digitization
932
933 double diff = resolution - dDdt * m_tdcResol;
934 if (diff > 0.) {
935 resolution = sqrt(diff * (resolution + dDdt * m_tdcResol));
936 } else {
937 resolution = 0.;
938 }
939
940#if defined(CDC_DEBUG)
941 cout << " " << endl;
942 cout << "CDCDigitizerModule::smearDriftLength" << endl;
943 cout << "tdcResol= " << m_tdcResol << endl;
944 cout << "dDdt,resolution= " << dDdt << " " << resolution << endl;
945#endif
946
947 // Smear drift length
948 double newDL = gRandom->Gaus(driftLength + mean, resolution);
949 while (newDL <= 0.) newDL = gRandom->Gaus(driftLength + mean, resolution);
950 // cout << "totalFugeF in Digi= " << m_totalFudgeFactor << endl;
951 return newDL;
952}
953
954
955double CDCDigitizerModule::getdDdt(const double driftL)
956{
957 //---------------------------------------------------------------------------------
958 // Calculates the 1'st derivative: dD/dt, where D: drift length before smearing; t: drift time
959 //---------------------------------------------------------------------------------
960
961 double dDdt = m_driftV;
962
964 const unsigned short layer = m_wireID.getICLayer();
965 const unsigned short leftRight = m_posFlag;
966 double alpha = m_cdcgp->getAlpha(m_posWire, m_momentum);
967 double theta = m_cdcgp->getTheta(m_momentum);
968 double t = m_cdcgp->getDriftTime(driftL, layer, leftRight, alpha, theta);
969 dDdt = m_cdcgp->getDriftV(t, layer, leftRight, alpha, theta);
970
971#if defined(CDC_DEBUG)
972 cout << " " << endl;
973 cout << "CDCDigitizerModule::getdDdt" << endl;
974 cout << "**layer= " << layer << endl;
975 cout << "alpha= " << 180.*alpha / M_PI << std::endl;
976 if (layer == 55) {
977 int lr = 0;
978 for (int i = 0; i < 1000; ++i) {
979 t = 1.0 * i;
980 double d = m_cdcgp->getDriftLength(t, layer, lr, alpha, theta);
981 cout << t << " " << d << endl;
982 }
983
984 cout << " " << endl;
985
986 lr = 1;
987 for (int i = 0; i < 100; ++i) {
988 t = 5 * i;
989 double d = m_cdcgp->getDriftLength(t, layer, lr, alpha, theta);
990 cout << t << " " << d << endl;
991 }
992 exit(-1);
993 }
994#endif
995 }
996
997 return dDdt;
998}
999
1000
1001double CDCDigitizerModule::getDriftTime(const double driftLength, const bool addTof, const bool addDelay)
1002{
1003 //---------------------------------------------------------------------------------
1004 // Method returning electron drift time (parameters: position in cm)
1005 // T(drift) = TOF + T(true drift time) + T(propagation delay in wire) - T(event),
1006 // T(event) is a timing of event, which includes a time jitter due to
1007 // the trigger system.
1008 //---------------------------------------------------------------------------------
1009
1010 double driftT = 0.;
1011
1013 driftT = (driftLength / Unit::cm) * m_driftVInv;
1014
1015#if defined(CDC_DEBUG)
1016 cout << " " << endl;
1017 cout << "CDCDigitizerModule::getDriftTime" << endl;
1018 cout << "driftvinv= " << m_driftVInv << endl;
1019#endif
1020 } else {
1021 const unsigned short layer = m_wireID.getICLayer();
1022 const unsigned short leftRight = m_posFlag;
1023 double alpha = m_cdcgp->getAlpha(m_posWire, m_momentum);
1024 double theta = m_cdcgp->getTheta(m_momentum);
1025 driftT = m_cdcgp->getDriftTime(driftLength, layer, leftRight, alpha, theta);
1026 // std::cout <<"alpha,theta,driftT= " << alpha <<" "<< theta <<" "<< driftT << std::endl;
1027 }
1028
1029 if (addTof) {
1030 driftT += m_flightTime; // in ns
1031 }
1032
1033 if (addDelay) {
1034 //calculate signal propagation length in the wire
1035 CDCGeometryPar::EWirePosition set = m_align ? CDCGeometryPar::c_Aligned : CDCGeometryPar::c_Base;
1036 B2Vector3D backWirePos = m_cdcgp->wireBackwardPosition(m_wireID, set);
1037
1038 double propLength = (m_posWire - backWirePos).Mag();
1039 // if (m_cdcgp->getSenseWireZposMode() == 1) {
1040 //TODO: replace the following with cached reference
1041 // std::cout << m_gcp->getInstance().getSenseWireZposMode() << std::endl;
1042 if (m_gcp->getSenseWireZposMode() == 1) {
1043 const unsigned short layer = m_wireID.getICLayer();
1044 propLength += m_cdcgp->getBwdDeltaZ(layer);
1045 }
1046 // B2DEBUG(29, "Propagation in wire length: " << propLength);
1047
1049 driftT += (propLength / Unit::cm) * m_propSpeedInv;
1050
1051#if defined(CDC_DEBUG)
1052 cout << "pseedinv= " << m_propSpeedInv << endl;
1053#endif
1054 } else {
1055 const unsigned short layer = m_wireID.getICLayer();
1056 driftT += (propLength / Unit::cm) * m_cdcgp->getPropSpeedInv(layer);
1057#if defined(CDC_DEBUG)
1058 cout << "layer,pseedinv= " << layer << " " << m_cdcgp->getPropSpeedInv(layer) << endl;
1059#endif
1060 }
1061 }
1062
1063 return driftT;
1064}
1065
1066
1067void CDCDigitizerModule::makeSignalsAfterShapers(const WireID& wid, double dEinGeV, double dx, double costh,
1068 unsigned short& adcCount, double& convFactorForThreshold)
1069{
1070 static double conv00 = (100.0 / 3.2); //keV -> coun (original from some test beam results)
1071 convFactorForThreshold = conv00;
1072 adcCount = 0;
1073 if (dEinGeV <= 0. || dx <= 0.) return;
1074
1075 const unsigned short layer = wid.getICLayer();
1076 const unsigned short cell = wid.getIWire();
1077 double dEInkeV = dEinGeV / Unit::keV;
1078
1079 double conv = conv00;
1080 if (m_spaceChargeEffect) {
1081 if (m_useDB4EDepToADC) {
1082 conv = m_cdcgp->getEDepToADCConvFactor(layer, cell, dEInkeV, dx, costh);
1083 double conv0 = m_cdcgp->getEDepToADCMainFactor(layer, cell, costh);
1084 convFactorForThreshold = (conv0 + m_degOfSPEOnThreshold * (conv - conv0));
1085 }
1086 } else {
1087 if (m_useDB4EDepToADC) conv = m_cdcgp->getEDepToADCMainFactor(layer, cell, costh);
1088 convFactorForThreshold = conv;
1089 }
1090
1091 if (convFactorForThreshold > 0.) {
1092 convFactorForThreshold *= getSemiTotalGain(layer, cell);
1093 } else {
1094 convFactorForThreshold = conv00;
1095 }
1096
1097 if (m_gasGainSmearing) {
1098 const int nElectrons = std::round(dEInkeV / m_effWForGasGainSmearing);
1099 double relGain = 0;
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) {
1104 relGain += Polya();
1105 }
1106 relGain /= nElectrons;
1107 } else {
1108 relGain = 1;
1109 }
1110 conv *= relGain;
1111 }
1112
1113 if (m_extraADCSmearing) {
1114 conv *= max(0., gRandom->Gaus(1., m_cdcgp->getEDepToADCSigma(layer, cell)));
1115 }
1116
1117 conv *= getSemiTotalGain(layer, cell);
1118
1119 //The ADCcount is obtained by rounding-up (measured voltage)/bin in real ADC. This is true both for pedestal and signal voltages, so the pedestal-subtracted ADCcount (simulated here) is rounded.
1120 adcCount = static_cast<unsigned short>(std::round(conv * dEInkeV));
1121 return;
1122}
1123
1124
1125double CDCDigitizerModule::Polya(double xmax)
1126{
1127 double x = 0;
1128 double y = 1;
1129 double fx = 0;
1130 double urndm[2];
1131 static double ymax = pow(m_thetaOfPolya, m_thetaOfPolya) * exp(-m_thetaOfPolya);
1132 while (y > fx) {
1133 gRandom->RndmArray(2, urndm);
1134 x = xmax * urndm[0];
1135 double a = (1 + m_thetaOfPolya) * x;
1136 fx = pow(a, m_thetaOfPolya) * exp(-a);
1137 y = ymax * urndm[1];
1138 }
1139 return x;
1140}
1141
1142
1143// Set FEE parameters (from DB)
1145{
1146 const double& off = m_tdcThresholdOffset;
1147 const double& gA = m_analogGain;
1148 const double& gD = m_digitalGain;
1149 const double& adcBW = m_adcBinWidth;
1150 const double convF = gA / gD / adcBW;
1151 const double el1TrgLatency = m_cdcgp->getMeanT0(); // ns
1152 B2DEBUG(m_debugLevel, "L1TRGLatency= " << el1TrgLatency);
1153 const double c = 32. * m_tdcBinWidth;
1154
1155 if (!m_fEElectronicsFromDB) B2FATAL("No FEEElectronics dbobject!");
1156 const CDCFEElectronics& fp = *((*m_fEElectronicsFromDB)[0]);
1157 int mode = (fp.getBoardID() == -1) ? 1 : 0;
1158 int iNBoards = static_cast<int>(c_nBoards);
1159
1160 //set typical values for all channels first if mode=1
1161 if (mode == 1) {
1162 for (int bdi = 1; bdi < iNBoards; ++bdi) {
1163 m_uprEdgeOfTimeWindow[bdi] = el1TrgLatency - c * (fp.getTrgDelay() + 1);
1164 if (m_uprEdgeOfTimeWindow[bdi] < 0.) B2FATAL("CDCDigitizer: Upper edge of time window < 0!");
1165 m_lowEdgeOfTimeWindow[bdi] = m_uprEdgeOfTimeWindow[bdi] - c * (fp.getWidthOfTimeWindow() + 1);
1166 if (m_lowEdgeOfTimeWindow[bdi] > 0.) B2FATAL("CDCDigitizer: Lower edge of time window > 0!");
1167 m_adcThresh[bdi] = fp.getADCThresh();
1168 m_tdcThresh[bdi] = convF * (off - fp.getTDCThreshInMV());
1169 m_widthOfTimeWindowInCount[bdi] = fp.getWidthOfTimeWindow() + 1;
1170 m_trgDelayInCount [bdi] = fp.getTrgDelay();
1171 }
1172 }
1173
1174 //overwrite values for specific channels if mode=1
1175 //set typical values for all channels if mode=0
1176 for (const auto& fpp : (*m_fEElectronicsFromDB)) {
1177 int bdi = fpp.getBoardID();
1178 if (mode == 0 && bdi == 0) continue; //bdi=0 is dummy (not used)
1179 if (mode == 1 && bdi == -1) continue; //skip typical case
1180 if (bdi < 0 || bdi >= iNBoards) B2FATAL("CDCDigitizer:: Invalid no. of FEE board!");
1181 m_uprEdgeOfTimeWindow[bdi] = el1TrgLatency - c * (fpp.getTrgDelay() + 1);
1182 if (m_uprEdgeOfTimeWindow[bdi] < 0.) B2FATAL("CDCDigitizer: Upper edge of time window < 0!");
1183 m_lowEdgeOfTimeWindow[bdi] = m_uprEdgeOfTimeWindow[bdi] - c * (fpp.getWidthOfTimeWindow() + 1);
1184 if (m_lowEdgeOfTimeWindow[bdi] > 0.) B2FATAL("CDCDigitizer: Lower edge of time window > 0!");
1185 m_adcThresh[bdi] = fpp.getADCThresh();
1186 m_tdcThresh[bdi] = convF * (off - fpp.getTDCThreshInMV());
1187 m_widthOfTimeWindowInCount[bdi] = fpp.getWidthOfTimeWindow() + 1;
1188 m_trgDelayInCount [bdi] = fpp.getTrgDelay();
1189 }
1190
1191 //debug
1192 B2DEBUG(m_debugLevel, "mode= " << mode);
1193 for (int bdi = 1; bdi < iNBoards; ++bdi) {
1194 B2DEBUG(m_debugLevel, bdi << " " << m_lowEdgeOfTimeWindow[bdi] << " " << m_uprEdgeOfTimeWindow[bdi] << " " << m_adcThresh[bdi] <<
1195 " " <<
1196 m_tdcThresh[bdi]);
1197 }
1198}
1199
1200// Set Run-gain (from DB)
1202{
1203 B2DEBUG(m_debugLevel, " ");
1204
1205 //read individual wire gains
1206 const int nLyrs = c_maxNSenseLayers;
1207 B2DEBUG(m_debugLevel, "nLyrs= " << nLyrs);
1208 int nGoodL[nLyrs] = {};
1209 float wgL[nLyrs] = {};
1210 int nGoodSL[c_nSuperLayers] = {};
1211 float wgSL[c_nSuperLayers] = {};
1212 int nGoodAll = 0;
1213 float wgAll = 0;
1214 int iw = -1;
1215 for (int lyr = 0; lyr < nLyrs; ++lyr) {
1216 int nWs = m_cdcgp->nWiresInLayer(lyr);
1217 for (int w = 0; w < nWs; ++w) {
1218 ++iw;
1219 float wg = (*m_wireGainFromDB)->getWireGain(iw);
1220 m_semiTotalGain[lyr][w] = wg;
1221 if (wg > 0) {
1222 ++nGoodL[lyr];
1223 wgL[lyr] += wg;
1224 WireID wid(lyr, w);
1225 ++nGoodSL[wid.getISuperLayer()];
1226 wgSL[wid.getISuperLayer()] += wg;
1227 ++nGoodAll;
1228 wgAll += wg;
1229 }
1230 }
1231 }
1232
1233 //calculate mean gain per layer
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]);
1237 }
1238 //calculate mean gain per superlayer
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]);
1242 }
1243
1244
1245 //calculate mean gain over all wires
1246 if (nGoodAll > 0) {
1247 wgAll /= nGoodAll;
1248 } else {
1249 B2FATAL("No good wires !");
1250 }
1251 B2DEBUG(m_debugLevel, "ngoodAll,gain= " << nGoodAll << " " << wgAll);
1252
1253 //set gain also for bad/dead wires (bad/dead in terms of dE/dx pid)
1254 for (int lyr = 0; lyr < nLyrs; ++lyr) {
1255 int nWs = m_cdcgp->nWiresInLayer(lyr);
1256 for (int w = 0; w < nWs; ++w) {
1257 if (m_semiTotalGain[lyr][w] <= 0) {
1258 if (wgL[lyr] > 0) {
1259 m_semiTotalGain[lyr][w] = wgL[lyr];
1260 } else {
1261 WireID wid(lyr, w);
1262 m_semiTotalGain[lyr][w] = wgSL[wid.getISuperLayer()];
1263 }
1264 }
1265 }
1266 }
1267
1268 //check if all gains > 0
1269 for (int lyr = 0; lyr < nLyrs; ++lyr) {
1270 int nWs = m_cdcgp->nWiresInLayer(lyr);
1271 for (int w = 0; w < nWs; ++w) {
1272 if (m_semiTotalGain[lyr][w] <= 0) {
1273 B2WARNING("Gain for lyr and wire " << lyr << " " << w << "not > 0. Strange! Replace it with " << wgAll << ".");
1274 m_semiTotalGain[lyr][w] = wgAll;
1275 }
1276 }
1277 }
1278
1279//multiply common factor for all wires
1280 m_runGain = (*m_runGainFromDB)->getRunGain();
1281 double cgain = (*m_gain0FromDB)->getScaleFactor();
1282 B2DEBUG(m_debugLevel, "runGain, sf= " << m_runGain << " " << cgain);
1283 cgain *= m_runGain * m_overallGainFactor;
1284 for (int lyr = 0; lyr < nLyrs; ++lyr) {
1285 int nWs = m_cdcgp->nWiresInLayer(lyr);
1286 for (int w = 0; w < nWs; ++w) {
1287 m_semiTotalGain[lyr][w] *= cgain;
1288 B2DEBUG(m_debugLevel, "lyr,wire,gain= " << lyr << " " << w << " " << m_semiTotalGain[lyr][w]);
1289 }
1290 }
1291}
1292
1293
1295{
1296 map<WireID, XTalkInfo> xTalkMap;
1297 map<WireID, XTalkInfo> xTalkMap1;
1299
1300 // Loop over all cdc hits to create a xtalk map
1301 int OriginalNoOfHits = m_cdcHits.getEntries();
1302 B2DEBUG(m_debugLevel4XTalk, "\n \n" << "#CDCHits " << OriginalNoOfHits);
1303 for (const auto& aHit : m_cdcHits) {
1304 if (m_issue2ndHitWarning && aHit.is2ndHit()) {
1305 B2WARNING("2nd TDC hit found, but not ready for it!");
1306 }
1307 WireID wid(aHit.getID());
1308 // B2DEBUG(m_debugLevel4XTalk, "Encoded wireid of current CDCHit: " << wid);
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);
1315
1316 int nXTalks = xTalks.size();
1317 for (int i = 0; i < nXTalks; ++i) {
1318 const unsigned short tdcCount4XTalk = xTalks[i].second.TDC;
1319 if (i == 0) {
1320 B2DEBUG(m_debugLevel4XTalk, "\n" << " signal: " << channel << " " << tdcCount << " " << adcCount << " " << tot);
1321 }
1322 B2DEBUG(m_debugLevel4XTalk, "xtalk: " << xTalks[i].first << " " << tdcCount4XTalk << " " << xTalks[i].second.ADC << " " <<
1323 xTalks[i].second.TOT);
1324 WireID widx = m_cdcgp->getWireID(board, xTalks[i].first);
1325 if (!m_cdcgp->isBadWire(widx)) { // for non-bad wire
1326 if (m_includeEarlyXTalks || (xTalks[i].second.TDC <= tdcCount)) {
1327 const double t0 = getPositiveT0(widx);
1328 const double ULOfTDC = (t0 - m_lowEdgeOfTimeWindow[board]) * m_tdcBinWidthInv;
1329 const double LLOfTDC = (t0 - m_uprEdgeOfTimeWindow[board]) * m_tdcBinWidthInv;
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)));
1333 }
1334 }
1335 // } else {
1336 // cout<<"badwire= " << widx.getICLayer() <<" "<< widx.getIWire() << endl;
1337 }
1338 } //end of xtalk loop
1339 } //end of cdc hit loop
1340
1341 //Loop over all xtalk hits to creat a new xtalk map with only the fastest hits kept (approx.)
1342 B2DEBUG(m_debugLevel4XTalk, "#xtalk hits: " << xTalkMap.size());
1343 for (const auto& aHit : xTalkMap) {
1344 WireID wid = aHit.first;
1345
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;
1351
1352 if (iterXTalkMap1 == xTalkMap1.end()) { // new entry
1353 xTalkMap1.insert(make_pair(wid, XTalkInfo(tdcCount, adcCount, tot, status)));
1354 // B2DEBUG(m_debugLevel4XTalk, "Creating a new xtalk hit with encoded wire no.: " << wid);
1355 } else { // not new; check if fastest
1356 if (tdcCount < iterXTalkMap1->second.m_tdc) {
1357 iterXTalkMap1->second.m_tdc = tdcCount;
1358 B2DEBUG(m_debugLevel4XTalk, "TDC-count of current xtalk: " << tdcCount);
1359 }
1360 iterXTalkMap1->second.m_adc += adcCount;
1361 iterXTalkMap1->second.m_tot += tot; // approx.
1362 }
1363 } // end of xtalk loop
1364
1365 //add xtalk in the same way as the beam bg. overlay
1366 B2DEBUG(m_debugLevel4XTalk, "#xtalk1 hits: " << xTalkMap1.size());
1367 for (const auto& aX : xTalkMap1) {
1368 bool append = true;
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;
1373
1374 for (int iHit = 0; iHit < OriginalNoOfHits; ++iHit) {
1375 CDCHit& aH = *(m_cdcHits[iHit]);
1376 if (aH.getID() != aX.first.getEWire()) { //wire id unmatched
1377 continue;
1378 } else { //wire id matched
1379 append = false;
1380 const unsigned short tdc4Sg = aH.getTDCCount();
1381 const unsigned short adc4Sg = aH.getADCCount();
1382 const unsigned short tot4Sg = aH.getTOT();
1383 // B2DEBUG(m_debuglevel4XTalk, "Sg tdc,adc,tot= " << tdc4Sg << " " << adc4Sg << " " << tot4Sg);
1384 // B2DEBUG(m_debugLevel4XTalk, "Bg tdc,adc,tot= " << tdc4Bg << " " << adc4Bg << " " << tot4Bg);
1385
1386 // If the BG hit is faster than the true hit, the TDC count is replaced, and
1387 // the relations are removed. ADC counts are summed up.
1388 if (tdc4Sg < tdc4Bg) {
1389 aH.setTDCCount(tdc4Bg);
1390 aH.setStatus(status4Bg);
1391 auto relSimHits = aH.getRelationsFrom<CDCSimHit>();
1392 for (int i = relSimHits.size() - 1; i >= 0; --i) {
1393 relSimHits.remove(i);
1394 }
1395 auto relMCParticles = aH.getRelationsFrom<MCParticle>();
1396 for (int i = relMCParticles.size() - 1; i >= 0; --i) {
1397 relMCParticles.remove(i);
1398 }
1399 }
1400
1401 aH.setADCCount(adc4Sg + adc4Bg);
1402
1403 //Set TOT for signal+background case. It is assumed that the start timing
1404 //of a pulse (input to ADC) is given by the TDC-count. This is an
1405 //approximation because analog (for ADC) and digital (for TDC) parts are
1406 //different in the front-end electronics.
1407 unsigned short s1 = tdc4Sg; //start time of 1st pulse
1408 unsigned short s2 = tdc4Bg; //start time of 2nd pulse
1409 unsigned short w1 = tot4Sg; //its width
1410 unsigned short w2 = tot4Bg; //its width
1411 if (tdc4Sg < tdc4Bg) {
1412 s1 = tdc4Bg;
1413 w1 = tot4Bg;
1414 s2 = tdc4Sg;
1415 w2 = tot4Sg;
1416 }
1417 w1 *= 32;
1418 w2 *= 32;
1419 const unsigned short e1 = s1 - w1; //end time of 1st pulse
1420 const unsigned short e2 = s2 - w2; //end time of 2nd pulse
1421 // B2DEBUG(m_debuglevel4Xtalk, "s1,e1,w1,s2,e2,w2= " << s1 << " " << e1 << " " << w1 << " " << s2 << " " << e2 << " " << w2);
1422
1423 double pulseW = w1 + w2;
1424 if (e1 <= e2) {
1425 pulseW = w1;
1426 } else if (e1 <= s2) {
1427 pulseW = s1 - e2;
1428 }
1429
1430 unsigned short board = m_cdcgp->getBoardID(aX.first);
1431 aH.setTOT(std::min(std::round(pulseW / 32.), static_cast<double>(m_widthOfTimeWindowInCount[board])));
1432 B2DEBUG(m_debugLevel4XTalk, "replaced tdc,adc,tot,wid,status= " << aH.getTDCCount() << " " << aH.getADCCount() << " " << aH.getTOT()
1433 <<
1434 " " << aH.getID() << " " << aH.getStatus());
1435 break;
1436 }
1437 } //end of cdc hit loop
1438
1439 if (append) {
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 <<
1442 " " <<
1443 status4Bg);
1444 }
1445 } //end of x-talk loop
1446 B2DEBUG(m_debugLevel4XTalk, "original #hits, #hits= " << OriginalNoOfHits << " " << m_cdcHits.getEntries());
1447}
1448
1449
1451{
1452 double t0 = m_cdcgp->getT0(wid);
1453 if (t0 <= 0 && m_treatNegT0WiresAsGood) t0 = m_cdcgp->getMeanT0();
1454 // B2DEBUG(m_debugLevel, m_cdcgp->getT0(wid) <<" "<< m_cdcgp->getMeanT0() <<" "<< t0);
1455 return t0;
1456}
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition B2Vector3.h:435
DataType X() const
access variable X (= .at(0) without boundary check)
Definition B2Vector3.h:431
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition B2Vector3.h:433
DataType Mag() const
The magnitude (rho in spherical coordinate system).
Definition B2Vector3.h:159
void SetY(DataType y)
set Y/2nd-coordinate
Definition B2Vector3.h:459
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.
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.
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.
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_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.
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...
Definition CDCHit.h:40
void setTDCCount(short tdcCount)
Setter for TDC count.
Definition CDCHit.h:128
void set2ndHitFlag()
Setter for 2nd hit flag.
Definition CDCHit.h:113
short getTDCCount() const
Getter for TDC count.
Definition CDCHit.h:219
void setTOT(unsigned short tot)
Setter for TOT.
Definition CDCHit.h:160
unsigned short getID() const
Getter for encoded wire number.
Definition CDCHit.h:193
unsigned short getStatus() const
Getter for CDCHit status.
Definition CDCHit.h:199
void setADCCount(unsigned short adcCount)
Setter for ADC count.
Definition CDCHit.h:135
unsigned short getADCCount() const
Getter for integrated charge.
Definition CDCHit.h:230
unsigned short getTOT() const
Getter for TOT.
Definition CDCHit.h:248
void setOtherHitIndices(CDCHit *otherHit)
Setter for the other hit indices.
Definition CDCHit.h:147
void setStatus(unsigned short status)
Setter for CDCHit status.
Definition CDCHit.h:106
Example Detector.
Definition CDCSimHit.h:21
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.
Definition EDepInGas.h:20
static EDepInGas & getInstance()
Static method to get a reference to the EDepInGas instance.
Definition EDepInGas.cc:21
double getEDepInGas(int mode, int pdg, double p, double dx, double e3) const
Return the energy deosite in the gas.
Definition EDepInGas.cc:135
Class for accessing arrays of objects in the database.
Definition DBArray.h:26
Class for accessing objects in the database.
Definition DBObjPtr.h:21
@ c_WriteOut
Object/array should be saved by output modules.
Definition DataStore.h:70
@ c_Event
Different object in each event, all objects/arrays are invalidated after event() function has been ca...
Definition DataStore.h:59
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
int getIndex() const
Get 1-based index of the particle in the corresponding MCParticle list.
Definition MCParticle.h:219
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition Module.cc:208
Module()
Constructor.
Definition Module.cc:30
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition Module.h:80
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]
Definition Unit.h:113
static const double eV
[electronvolt]
Definition Unit.h:112
static const double cm
Standard units with the value = 1.
Definition Unit.h:47
Class to identify a wire inside the CDC.
Definition WireID.h:34
unsigned short getICLayer() const
Getter for continuous layer numbering.
Definition WireID.cc:24
unsigned short getIWire() const
Getter for wire within the layer.
Definition WireID.h:145
unsigned short getISuperLayer() const
Getter for Super-Layer.
Definition WireID.h:130
STL iterator class.
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition B2Vector3.h:516
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
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.
STL namespace.
Structure for saving the signal information.
Structure for saving the x-talk information.