Belle II Software  release-05-01-25
CDCDigitizerModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2012 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Guofu Cao, Martin Heck, CDC group *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <cdc/modules/cdcDigitizer/CDCDigitizerModule.h>
12 #include <cdc/utilities/ClosestApproach.h>
13 
14 #include <framework/datastore/RelationArray.h>
15 //#include <framework/datastore/RelationIndex.h>
16 #include <framework/gearbox/Unit.h>
17 #include <framework/logging/Logger.h>
18 
19 #include <TRandom.h>
20 #include <map>
21 #include <utility>
22 
23 using namespace std;
24 using namespace Belle2;
25 using namespace CDC;
26 
27 
28 // register module
29 REG_MODULE(CDCDigitizer)
31  m_cdcgp(), m_gcp(), m_aCDCSimHit(), m_posFlag(0),
32  m_driftLength(0.0), m_flightTime(0.0), m_globalTime(0.0),
33  m_tdcBinWidth(1.0), m_tdcBinWidthInv(1.0),
34  m_tdcResol(0.9825), m_driftV(4.0e-3),
35  m_driftVInv(250.0), m_propSpeedInv(27.25), m_align(true)
36 {
37  // Set description
38  setDescription("Creates CDCHits from CDCSimHits.");
39  setPropertyFlags(c_ParallelProcessingCertified);
40 
41  // Add parameters
42  // I/O
43  addParam("InputCDCSimHitsName", m_inputCDCSimHitsName, "Name of input array. Should consist of CDCSimHits.", string(""));
44  addParam("OutputCDCHitsName", m_outputCDCHitsName, "Name of output array. Will consist of CDCHits.", string(""));
45  addParam("OutputCDCHitsName4Trg", m_outputCDCHitsName4Trg,
46  "Name of output array for trigger. Can contain several hits per wire, "
47  "if they correspond to different time windows of 32ns.",
48  string("CDCHits4Trg"));
49 
50  //Relations
51  addParam("MCParticlesToCDCSimHitsName", m_MCParticlesToSimHitsName,
52  "Name of relation between MCParticles and CDCSimHits used", string(""));
53  addParam("CDCSimHistToCDCHitsName", m_SimHitsTOCDCHitsName,
54  "Name of relation between the CDCSimHits and the CDCHits used", string(""));
55 
56 
57  //Parameters for Digitization
58  addParam("UseSimpleDigitization", m_useSimpleDigitization,
59  "If true, a simple x-t with a constant velocity is used for the drift-length to -time conversion", false);
60 
61  //float Gauss Parameters
62  addParam("Fraction", m_fraction, "Fraction of first Gaussian used to smear drift length in cm", 1.0);
63  addParam("Mean1", m_mean1, "Mean value of first Gaussian used to smear drift length in cm", 0.0000);
64  addParam("Resolution1", m_resolution1, "Resolution of first Gaussian used to smear drift length in cm", 0.0130);
65  addParam("Mean2", m_mean2, "Mean value of second Gaussian used to smear drift length in cm", 0.0000);
66  addParam("Resolution2", m_resolution2, "Resolution of second Gaussian used to smear drift length in cm", 0.0000);
67 
68  //Switch to control smearing
69  addParam("DoSmearing", m_doSmearing,
70  "If false, drift length will not be smeared.", true);
71 
72  addParam("TrigTimeJitter", m_trigTimeJitter,
73  "Magnitude (w) of trigger timing jitter (ns). The trigger timing is randuminzed uniformly in a time window of [-w/2, +w/2].",
74  0.);
75  //Switches to control time information handling
76  addParam("AddTimeWalk", m_addTimeWalk, "A switch for time-walk (pulse-heght dep. delay); true: on; false: off", true);
77  addParam("AddInWirePropagationDelay", m_addInWirePropagationDelay,
78  "A switch used to control adding propagation delay in the wire into the final drift time or not; this is for signal hits.", true);
79  addParam("AddInWirePropagationDelay4Bg", m_addInWirePropagationDelay4Bg,
80  "The same switch but for beam bg. hits.", true);
81  addParam("AddTimeOfFlight", m_addTimeOfFlight,
82  "A switch used to control adding time of flight into the final drift time or not; this is for signal hits.", true);
83  addParam("AddTimeOfFlight4Bg", m_addTimeOfFlight4Bg,
84  "The same switch but for beam bg. hits.", true);
85  addParam("OutputNegativeDriftTime", m_outputNegativeDriftTime, "Output hits with negative drift time", true);
86  addParam("Output2ndHit", m_output2ndHit,
87  "Output the 2nd hit if exists in the time window. Note that it is not well-simulated at all, partly because no cross-talk betw. channels is simulated.",
88  false);
89  //Switch to control sense wire sag
90  addParam("CorrectForWireSag", m_correctForWireSag,
91  "A switch for sense wire sag effect; true: drift-time is calculated with the sag taken into account; false: not. Here, sag means the perturbative part which corresponds to alignment in case of wire-position. The main part (corresponding to design+displacement in wire-position) is taken into account in FullSim; you can control it via CDCJobCntlParModifier.",
92  true);
93 
94  //Threshold
95  addParam("TDCThreshold4Outer", m_tdcThreshold4Outer,
96  "TDC threshold (dE in eV) for Layers#8-56. The value corresponds to He-C2H6 gas; for the gas+wire (MaterialDefinitionMode=0) case, (this value)*f will be used, where f is specified by GasToGasWire",
97  25.0);
98  addParam("TDCThreshold4Inner", m_tdcThreshold4Inner,
99  "Same as TDCThreshold4Outer but for Layers#0-7,", 25.0);
100  addParam("GasToGasWire", m_gasToGasWire,
101  "(Approximate) ratio of dE in He/C2H6-gas to dE in gas+wire, where dE is energy deposit.", 1. / 1.478);
102 
103  //ADC Threshold
104  addParam("ADCThreshold", m_adcThreshold,
105  "Threshold for ADC-count (in unit of count). ADC-count < threshold is treated as count=0.", 2);
106  addParam("tMin", m_tMin, "Lower edge of time window in ns", -100.);
107  addParam("tMaxOuter", m_tMaxOuter, "Upper edge of time window in ns for the normal-cell layers", 500.);
108  addParam("tMaxInner", m_tMaxInner, "Upper edge of time window in ns for the small-cell layers", 300.);
109  // The following doesn't make any sense. The only reasonable steerable would be a switch to decide if the jitter shall be
110  // activated. Then there has to be event by event jitter.
111  /* addParam("EventTime", m_eventTime,
112  "It is a timing of event, which includes a time jitter due to the trigger system, set in ns", float(0.0));*/
113 
114  //Switch for database
115  addParam("UseDB4FEE", m_useDB4FEE, "Fetch and use FEE params. from database or not", true);
116  addParam("UseDB4EDepToADC", m_useDB4EDepToADC, "Uuse edep-to-ADC conversion params. from database or not", true);
117  addParam("UseDB4RunGain", m_useDB4RunGain, "Fetch and use run gain from database or not", true);
118  addParam("OverallGainFactor", m_overallGainFactor, "Overall gain factor for adjustment", 1.0);
119 
120  //Some FEE params.
121  addParam("TDCThresholdOffset", m_tdcThresholdOffset, "Offset for TDC (digital) threshold (mV)", 3820.);
122  addParam("AnalogGain", m_analogGain, "Analog gain (V/pC)", 1.1);
123  addParam("DigitalGain", m_digitalGain, "Digital gain (V/pC)", 15.);
124  addParam("ADCBinWidth", m_adcBinWidth, "ADC bin width (mV)", 2.);
125 
126  addParam("AddFudgeFactorForSigma", m_addFudgeFactorForSigma,
127  "Additional fudge factor for space resol. (common to all cells)", 1.);
128  addParam("SpaceChargeEffect", m_spaceChargeEffect, "Switch for space charge effect", true);
129 
130  addParam("AddXTalk", m_addXTalk, "A switch for crosstalk; true: on; false: off", true);
131  addParam("Issue2ndHitWarning", m_issue2ndHitWarning, "=true: issue a warning when a 2nd TDC hit is found.", true);
132  addParam("IncludeEarlyXTalks", m_includeEarlyXTalks, "=true: include earlier x-talks as well than the signal hit in question.",
133  true);
134  addParam("DebugLevel4XTalk", m_debugLevel4XTalk, "Debug level fr crosstalk; 20-29 are usable.", 20);
135 
136 
137 #if defined(CDC_DEBUG)
138  cout << " " << endl;
139  cout << "CDCDigitizer constructor" << endl;
140 #endif
141 }
142 
143 void CDCDigitizerModule::initialize()
144 {
145  m_simHits.isRequired(m_inputCDCSimHitsName);
146 
147  // Register the arrays in the DataStore, that are to be added in this module.
148  m_cdcHits.registerInDataStore(m_outputCDCHitsName);
149  m_simHits.registerRelationTo(m_cdcHits);
150  m_mcParticles.registerRelationTo(m_cdcHits);
151  // Arrays for trigger.
152  m_cdcHits4Trg.registerInDataStore(m_outputCDCHitsName4Trg);
153  m_simHits.registerRelationTo(m_cdcHits4Trg);
154  m_mcParticles.registerRelationTo(m_cdcHits4Trg);
155 
156  m_cdcgp = &(CDCGeometryPar::Instance());
157  CDCGeometryPar& cdcgp = *m_cdcgp;
158  m_tdcBinWidth = cdcgp.getTdcBinWidth();
159  m_tdcBinWidthInv = 1. / m_tdcBinWidth;
160  m_tdcResol = m_tdcBinWidth / sqrt(12.);
161  m_driftV = cdcgp.getNominalDriftV();
162  m_driftVInv = 1. / m_driftV;
163  m_propSpeedInv = 1. / cdcgp.getNominalPropSpeed();
164  m_scaleFac = 1.;
165  if (m_cdcgp->getMaterialDefinitionMode() == 0) { //gas+wire mode
166  m_scaleFac = m_gasToGasWire;
167  }
168  m_scaleFac *= Unit::GeV;
169  m_gcp = &(CDCGeoControlPar::getInstance());
170  m_totalFudgeFactor = m_cdcgp->getFudgeFactorForSigma(2);
171  m_totalFudgeFactor *= m_addFudgeFactorForSigma;
172  B2DEBUG(29, "totalFugeF in Digi= " << m_totalFudgeFactor);
173  /*
174  m_fraction = 1.0;
175  m_resolution1 = cdcgp.getNominalSpaceResol();
176  m_resolution2 = 0.;
177  m_mean1 = 0.;
178  m_mean2 = 0.;
179  */
180 
181  if (m_useDB4FEE) {
182  m_fEElectronicsFromDB = new DBArray<CDCFEElectronics>;
183  if ((*m_fEElectronicsFromDB).isValid()) {
184  (*m_fEElectronicsFromDB).addCallback(this, &CDCDigitizerModule::setFEElectronics);
185  setFEElectronics();
186  } else {
187  B2FATAL("CDCDigitizer:: CDCFEElectronics not valid !");
188  }
189  }
190 
191  /*
192  if (m_useDB4EDepToADC) {
193  m_eDepToADCConversionsFromDB = new DBObjPtr<CDCEDepToADCConversions>;
194  if ((*m_eDepToADCConversionsFromDB).isValid()) {
195  (*m_eDepToADCConversionsFromDB).addCallback(this, &CDCDigitizerModule::setEDepToADCConversions);
196  setEDepToADCConversions();
197  } else {
198  B2FATAL("CDCDigitizer:: CDCEDepToADCConversions not valid !");
199  }
200  }
201  */
202 
203  if (m_useDB4RunGain) {
204  m_runGainFromDB = new DBObjPtr<CDCDedxRunGain>;
205  if ((*m_runGainFromDB).isValid()) {
206  (*m_runGainFromDB).addCallback(this, &CDCDigitizerModule::setRunGain);
207  setRunGain();
208  } else {
209  B2FATAL("CDCDedxRunGain invalid!");
210  }
211  }
212  B2DEBUG(29, "run-gain = " << m_runGain);
213 
214  if (m_addXTalk) {
215  m_xTalkFromDB = new DBObjPtr<CDCCrossTalkLibrary>;
216  if ((*m_xTalkFromDB).isValid()) {
217  } else {
218  B2FATAL("CDCCrossTalkLibrary invalid!");
219  }
220  }
221 
222 #if defined(CDC_DEBUG)
223  cout << " " << endl;
224  cout << "CDCDigitizer initialize" << endl;
225  // cout << "m_tdcOffset= " << m_tdcOffset << endl;
226  cout << "m_tdcBinWidth= " << m_tdcBinWidth << endl;
227  cout << "m_tdcResol= " << m_tdcResol << endl;
228  cout << "m_driftV= " << m_driftV << endl;
229  cout << "m_driftVInv= " << m_driftVInv << endl;
230  cout << "m_propSpeedInv= " << m_propSpeedInv << endl;
231  /*
232  cout << "m_fraction= " << m_fraction << endl;
233  cout << "m_resolution1= " << m_resolution1 << endl;
234  cout << "m_resolution2= " << m_resolution2 << endl;
235  cout << "m_mean1= " << m_mean1 << endl;
236  cout << "m_mean2= " << m_mean2 << endl;
237  */
238 #endif
239 
240  if (m_useDB4EDepToADC) {
241  if (m_cdcgp->getEDepToADCMainFactor(0, 0) == 0.) {
242  B2FATAL("CDCEDepToADCConversion payloads are unavailable!");
243  }
244  }
245 
246 }
247 
248 void CDCDigitizerModule::event()
249 {
250  // Get SimHit array, MCParticle array, and relation between the two.
251  RelationArray mcParticlesToCDCSimHits(m_mcParticles, m_simHits); //RelationArray created by CDC SensitiveDetector
252 
253 
254  //--- Start Digitization --------------------------------------------------------------------------------------------
255  // Merge the hits in the same cell and save them into CDC signal map.
256 
257  // Define signal map
258  map<WireID, SignalInfo> signalMap;
259  map<WireID, SignalInfo>::iterator iterSignalMap;
260  // Define adc map
261  map<WireID, unsigned short> adcMap;
262  map<WireID, unsigned short>::iterator iterADCMap;
263  // map<WireID, double> adcMap;
264  // map<WireID, double>::iterator iterADCMap;
265 
266  // signal map for trigger
267  map<pair<WireID, unsigned>, SignalInfo> signalMapTrg;
268  map<pair<WireID, unsigned>, SignalInfo>::iterator iterSignalMapTrg;
269 
270  // Set trigger timing jitter for this event
271  double trigTiming = m_trigTimeJitter == 0. ? 0. : m_trigTimeJitter * (gRandom->Uniform() - 0.5);
272  // std::cout << "trigTiming= " << trigTiming << std::endl;
273  // Loop over all hits
274  int nHits = m_simHits.getEntries();
275  B2DEBUG(29, "Number of CDCSimHits in the current event: " << nHits);
276  for (int iHits = 0; iHits < nHits; ++iHits) {
277  // Get a hit
278  m_aCDCSimHit = m_simHits[iHits];
279 
280  // Hit geom. info
281  m_wireID = m_aCDCSimHit->getWireID();
282  // B2DEBUG(29, "Encoded wire number of current CDCSimHit: " << m_wireID);
283 
284  m_posFlag = m_aCDCSimHit->getLeftRightPassageRaw();
285  m_boardID = m_cdcgp->getBoardID(m_wireID);
286  // B2DEBUG(29, "m_boardID= " << m_boardID);
287  m_posWire = m_aCDCSimHit->getPosWire();
288  m_posTrack = m_aCDCSimHit->getPosTrack();
289  m_momentum = m_aCDCSimHit->getMomentum();
290  m_flightTime = m_aCDCSimHit->getFlightTime();
291  m_globalTime = m_aCDCSimHit->getGlobalTime();
292  m_driftLength = m_aCDCSimHit->getDriftLength() * Unit::cm;
293 
294  //include alignment effects
295  //basically align flag should be always on since on/off is controlled by the input alignment.xml file itself.
296  m_align = true;
297 
298  TVector3 bwpAlign = m_cdcgp->wireBackwardPosition(m_wireID, CDCGeometryPar::c_Aligned);
299  TVector3 fwpAlign = m_cdcgp->wireForwardPosition(m_wireID, CDCGeometryPar::c_Aligned);
300 
301  TVector3 bwp = m_cdcgp->wireBackwardPosition(m_wireID);
302  TVector3 fwp = m_cdcgp->wireForwardPosition(m_wireID);
303 
304  //skip correction for wire-position alignment if unnecessary
305  if ((bwpAlign - bwp).Mag() == 0. && (fwpAlign - fwp).Mag() == 0.) m_align = false;
306  // std::cout << "a m_align= " << m_align << std::endl;
307 
308  if (m_align || m_correctForWireSag) {
309 
310  bwp = bwpAlign;
311  fwp = fwpAlign;
312 
313  if (m_correctForWireSag) {
314  double zpos = m_posWire.z();
315  double bckYSag = bwp.y();
316  double forYSag = fwp.y();
317 
318  // CDCGeometryPar::EWirePosition set = m_align ?
319  // CDCGeometryPar::c_Aligned : CDCGeometryPar::c_Base;
320  CDCGeometryPar::EWirePosition set = CDCGeometryPar::c_Aligned;
321  const int layerID = m_wireID.getICLayer();
322  const int wireID = m_wireID.getIWire();
323  m_cdcgp->getWireSagEffect(set, layerID, wireID, zpos, bckYSag, forYSag);
324  bwp.SetY(bckYSag);
325  fwp.SetY(forYSag);
326  }
327 
328  const TVector3 L = 5. * m_momentum.Unit(); //(cm) tentative
329  TVector3 posIn = m_posTrack - L;
330  TVector3 posOut = m_posTrack + L;
331  TVector3 posTrack = m_posTrack;
332  TVector3 posWire = m_posWire;
333 
334  // m_driftLength = m_cdcgp->ClosestApproach(bwp, fwp, posIn, posOut, posTrack, posWire);
335  m_driftLength = ClosestApproach(bwp, fwp, posIn, posOut, posTrack, posWire);
336  // std::cout << "base-dl, sag-dl, diff= " << m_aCDCSimHit->getDriftLength() <<" "<< m_driftLength <<" "<< m_driftLength - m_aCDCSimHit->getDriftLength() << std::endl;
337  m_posTrack = posTrack;
338  m_posWire = posWire;
339 
340  double deltaTime = 0.; //tentative (probably ok...)
341  // double deltaTime = (posTrack - m_posTrack).Mag() / speed;
342  m_flightTime += deltaTime;
343  m_globalTime += deltaTime;
344  m_posFlag = m_cdcgp->getNewLeftRightRaw(m_posWire, m_posTrack, m_momentum);
345  }
346 
347  // Calculate measurement time.
348  // Smear drift length
349  double hitDriftLength = m_driftLength;
350  double dDdt = getdDdt(hitDriftLength);
351  if (m_doSmearing) {
352  hitDriftLength = smearDriftLength(hitDriftLength, dDdt);
353  }
354 
355  //set flags
356  bool addTof = m_addTimeOfFlight4Bg;
357  bool addDelay = m_addInWirePropagationDelay4Bg;
358  if (m_aCDCSimHit->getBackgroundTag() == 0) {
359  addTof = m_addTimeOfFlight;
360  addDelay = m_addInWirePropagationDelay;
361  }
362  double hitDriftTime = getDriftTime(hitDriftLength, addTof, addDelay);
363 
364  //add randamized event time for a beam bg. hit
365  if (m_aCDCSimHit->getBackgroundTag() != 0) {
366  hitDriftTime += m_globalTime - m_flightTime;
367  }
368 
369  //add trigger timing jitter
370  hitDriftTime += trigTiming;
371 
372  //apply time window cut
373  double tMin = m_tMin;
374  double tMax = m_tMaxOuter;
375  if (m_wireID.getISuperLayer() == 0) tMax = m_tMaxInner;
376  if (m_useDB4FEE) {
377  tMin = m_lowEdgeOfTimeWindow[m_boardID];
378  tMax = m_uprEdgeOfTimeWindow[m_boardID];
379  }
380  if (hitDriftTime < tMin || hitDriftTime > tMax) continue;
381 
382  //Sum ADC count
383  const double stepLength = m_aCDCSimHit->getStepLength() * Unit::cm;
384  const double costh = m_momentum.z() / m_momentum.Mag();
385  const double hitdE = m_scaleFac * m_aCDCSimHit->getEnergyDep();
386  // B2DEBUG(29, "m_scaleFac,UnitGeV= " << m_scaleFac <<" "<< Unit::GeV);
387  unsigned short layerID = m_wireID.getICLayer();
388  unsigned short cellID = m_wireID.getIWire();
389  // unsigned short adcCount = getADCCount(layerID, cellID, hitdE, stepLength, costh);
390  unsigned short adcCount = getADCCount(m_wireID, hitdE, stepLength, costh);
391  const unsigned short adcTh = m_useDB4FEE ? m_adcThresh[m_boardID] : m_adcThreshold;
392  // B2DEBUG(29, "adcTh,adcCount= " << adcTh <<" "<< adcCount);
393  if (adcCount < adcTh) adcCount = 0;
394  iterADCMap = adcMap.find(m_wireID);
395  if (iterADCMap == adcMap.end()) {
396  adcMap.insert(make_pair(m_wireID, adcCount));
397  // adcMap.insert(make_pair(m_wireID, hitdE));
398  } else {
399  iterADCMap->second += adcCount;
400  // iterADCMap->second += hitdE;
401  }
402 
403  //Apply energy threshold
404  // If hitdE < dEThreshold, the hit is ignored
405  // M. Uchida 2012.08.31
406  double dEThreshold = 0.;
407  if (m_useDB4FEE && m_useDB4EDepToADC) {
408  dEThreshold = (m_tdcThresh[m_boardID] / m_cdcgp->getEDepToADCMainFactor(layerID, cellID)) * Unit::keV;
409  } else {
410  dEThreshold = (m_wireID.getISuperLayer() == 0) ? m_tdcThreshold4Inner : m_tdcThreshold4Outer;
411  dEThreshold *= Unit::eV;
412  }
413  B2DEBUG(29, "hitdE,dEThreshold,driftLength " << hitdE << " " << dEThreshold << " " << hitDriftLength);
414 
415  if (hitdE < dEThreshold) {
416  B2DEBUG(29, "Below Ethreshold: " << hitdE << " " << dEThreshold);
417  continue;
418  }
419 
420  // add one hit per trigger time window to the trigger signal map
421  unsigned short trigWindow = floor((hitDriftTime - m_tMin) * m_tdcBinWidthInv / 32);
422  iterSignalMapTrg = signalMapTrg.find(make_pair(m_wireID, trigWindow));
423  if (iterSignalMapTrg == signalMapTrg.end()) {
424  // signalMapTrg.insert(make_pair(make_pair(m_wireID, trigWindow),
425  // SignalInfo(iHits, hitDriftTime, hitdE)));
426  signalMapTrg.insert(make_pair(make_pair(m_wireID, trigWindow),
427  SignalInfo(iHits, hitDriftTime, adcCount)));
428  } else {
429  if (hitDriftTime < iterSignalMapTrg->second.m_driftTime) {
430  iterSignalMapTrg->second.m_driftTime = hitDriftTime;
431  iterSignalMapTrg->second.m_simHitIndex = iHits;
432  }
433  // iterSignalMapTrg->second.m_charge += hitdE;
434  iterSignalMapTrg->second.m_charge += adcCount;
435  }
436 
437  // Reject totally-dead wire; to be replaced by isDeadWire() in future
438  if (m_cdcgp->isBadWire(m_wireID)) {
439  // std::cout<<"badwire= " << m_wireID.getICLayer() <<" "<< m_wireID.getIWire() << std::endl;
440  continue;
441  }
442  // Reject partly-dead wire as well
443  double eff = 1.;
444  if (m_cdcgp->isDeadWire(m_wireID, eff)) {
445  // std::cout << "wid,eff= " << m_wireID << " " << eff << std::endl;
446  if (eff < gRandom->Uniform()) continue;
447  }
448 
449  // 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).
450  const double a = bwpAlign.X();
451  const double b = bwpAlign.Y();
452  const double c = bwpAlign.Z();
453  const TVector3 fmbAlign = fwpAlign - bwpAlign;
454  const double lmn = 1. / fmbAlign.Mag();
455  const double l = fmbAlign.X() * lmn;
456  const double m = fmbAlign.Y() * lmn;
457  const double n = fmbAlign.Z() * lmn;
458 
459  double dx = m_aCDCSimHit->getPosIn().X() - a;
460  double dy = m_aCDCSimHit->getPosIn().Y() - b;
461  double dz = m_aCDCSimHit->getPosIn().Z() - c;
462  double sub = l * dx + m * dy + n * dz;
463  const double driftLFromIn = sqrt(dx * dx + dy * dy + dz * dz - sub * sub);
464 
465  dx = m_aCDCSimHit->getPosOut().X() - a;
466  dy = m_aCDCSimHit->getPosOut().Y() - b;
467  dz = m_aCDCSimHit->getPosOut().Z() - c;
468  sub = l * dx + m * dy + n * dz;
469  const double driftLFromOut = sqrt(dx * dx + dy * dy + dz * dz - sub * sub);
470 
471  const double maxDriftL = std::max(driftLFromIn, driftLFromOut);
472  const double minDriftL = m_driftLength;
473  B2DEBUG(29, "driftLFromIn= " << driftLFromIn << " driftLFromOut= " << driftLFromOut << " minDriftL= " << minDriftL << " maxDriftL= "
474  <<
475  maxDriftL << "m_driftLength= " << m_driftLength);
476 
477  iterSignalMap = signalMap.find(m_wireID);
478 
479  if (iterSignalMap == signalMap.end()) {
480  // new entry
481  // signalMap.insert(make_pair(m_wireID, SignalInfo(iHits, hitDriftTime, hitdE)));
482  signalMap.insert(make_pair(m_wireID, SignalInfo(iHits, hitDriftTime, adcCount, maxDriftL, minDriftL)));
483  B2DEBUG(29, "Creating new Signal with encoded wire number: " << m_wireID);
484  } else {
485  // ... smallest drift time has to be checked, ...
486  if (hitDriftTime < iterSignalMap->second.m_driftTime) {
487  iterSignalMap->second.m_driftTime3 = iterSignalMap->second.m_driftTime2;
488  iterSignalMap->second.m_simHitIndex3 = iterSignalMap->second.m_simHitIndex2;
489  iterSignalMap->second.m_driftTime2 = iterSignalMap->second.m_driftTime;
490  iterSignalMap->second.m_simHitIndex2 = iterSignalMap->second.m_simHitIndex;
491  iterSignalMap->second.m_driftTime = hitDriftTime;
492  iterSignalMap->second.m_simHitIndex = iHits;
493  B2DEBUG(29, "hitDriftTime of current Signal: " << hitDriftTime << ", hitDriftLength: " << hitDriftLength);
494  } else if (hitDriftTime < iterSignalMap->second.m_driftTime2) {
495  iterSignalMap->second.m_driftTime3 = iterSignalMap->second.m_driftTime2;
496  iterSignalMap->second.m_simHitIndex3 = iterSignalMap->second.m_simHitIndex2;
497  iterSignalMap->second.m_driftTime2 = hitDriftTime;
498  iterSignalMap->second.m_simHitIndex2 = iHits;
499  } else if (hitDriftTime < iterSignalMap->second.m_driftTime3) {
500  iterSignalMap->second.m_driftTime3 = hitDriftTime;
501  iterSignalMap->second.m_simHitIndex3 = iHits;
502  }
503  // ... total charge has to be updated.
504  // iterSignalMap->second.m_charge += hitdE;
505  iterSignalMap->second.m_charge += adcCount;
506 
507  // set max and min driftLs
508  if (iterSignalMap->second.m_maxDriftL < maxDriftL) iterSignalMap->second.m_maxDriftL = maxDriftL;
509  if (iterSignalMap->second.m_minDriftL > minDriftL) iterSignalMap->second.m_minDriftL = minDriftL;
510  B2DEBUG(29, "maxDriftL in struct= " << iterSignalMap->second.m_maxDriftL << "minDriftL in struct= " <<
511  iterSignalMap->second.m_minDriftL);
512  }
513 
514  } // end loop over SimHits.
515 
516  //--- Now Store the results into CDCHits and
517  // create corresponding relations between SimHits and CDCHits.
518 
519  unsigned int iCDCHits = 0;
520  RelationArray cdcSimHitsToCDCHits(m_simHits, m_cdcHits); //SimHit<->CDCHit
521  RelationArray mcParticlesToCDCHits(m_mcParticles, m_cdcHits); //MCParticle<->CDCHit
522 
523  for (iterSignalMap = signalMap.begin(); iterSignalMap != signalMap.end(); ++iterSignalMap) {
524 
525  //add time-walk (here for simplicity)
526  // unsigned short adcCount = getADCCount(iterSignalMap->second.m_charge);
527  // unsigned short adcCount = iterSignalMap->second.m_charge;
528  iterADCMap = adcMap.find(iterSignalMap->first);
529  unsigned short adcCount = iterADCMap != adcMap.end() ? iterADCMap->second : 0;
530  /*
531  unsigned short adcCount = 0;
532  if (iterADCMap != adcMap.end()) {
533  adcCount = getADCCount(iterSignalMap->first, iterADCMap->second, 1., 0.);
534  unsigned short boardID = m_cdcgp->getBoardID(iterSignalMap->first);
535  // B2DEBUG(29, "boardID= " << boardID);
536  const unsigned short adcTh = m_useDB4FEE ? m_adcThresh[boardID] : m_adcThreshold;
537  if (adcCount < adcTh) adcCount = 0;
538  }
539  */
540 
541  if (m_addTimeWalk) {
542  B2DEBUG(29, "timewalk= " << m_cdcgp->getTimeWalk(iterSignalMap->first, adcCount));
543  iterSignalMap->second.m_driftTime += m_cdcgp->getTimeWalk(iterSignalMap->first, adcCount);
544  }
545 
546  //remove negative drift time (TDC) upon request
547  if (!m_outputNegativeDriftTime &&
548  iterSignalMap->second.m_driftTime < 0.) {
549  continue;
550  }
551 
552  //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.
553  unsigned short tdcCount = static_cast<unsigned short>((m_cdcgp->getT0(iterSignalMap->first) - iterSignalMap->second.m_driftTime) *
554  m_tdcBinWidthInv + 0.5);
555 
556  //calculate tot; hard-coded currently
557  double deltaDL = iterSignalMap->second.m_maxDriftL - iterSignalMap->second.m_minDriftL;
558  if (deltaDL < 0.) {
559  B2DEBUG(29, "negative deltaDL= " << deltaDL);
560  deltaDL = 0.;
561  }
562  const unsigned short boardID = m_cdcgp->getBoardID(iterSignalMap->first);
563  unsigned short tot = std::min(std::round(5.92749 * deltaDL + 2.59706), static_cast<double>(m_widthOfTimeWindow[boardID]));
564  if (m_adcThresh[boardID] > 0) {
565  tot = std::min(static_cast<int>(tot), static_cast<int>(adcCount / m_adcThresh[boardID]));
566  }
567 
568  CDCHit* firstHit = m_cdcHits.appendNew(tdcCount, adcCount, iterSignalMap->first, 0, tot);
569  // std::cout <<"firsthit?= " << firstHit->is2ndHit() << std::endl;
570  //set a relation: CDCSimHit -> CDCHit
571  cdcSimHitsToCDCHits.add(iterSignalMap->second.m_simHitIndex, iCDCHits);
572 
573  //set a relation: MCParticle -> CDCHit
574  RelationVector<MCParticle> rels = m_simHits[iterSignalMap->second.m_simHitIndex]->getRelationsFrom<MCParticle>();
575  if (rels.size() != 0) {
576  //assumption: only one MCParticle
577  const MCParticle* mcparticle = rels[0];
578  double weight = rels.weight(0);
579  mcparticle->addRelationTo(firstHit, weight);
580  }
581 
582  //Set 2nd-hit related things if it exists
583  if (m_output2ndHit && iterSignalMap->second.m_simHitIndex2 >= 0) {
584  unsigned short tdcCount2 = static_cast<unsigned short>((m_cdcgp->getT0(iterSignalMap->first) - iterSignalMap->second.m_driftTime2) *
585  m_tdcBinWidthInv + 0.5);
586  if (tdcCount2 != tdcCount) {
587  CDCHit* secondHit = m_cdcHits.appendNew(tdcCount2, adcCount, iterSignalMap->first, 0, tot);
588  secondHit->set2ndHitFlag();
589  secondHit->setOtherHitIndices(firstHit);
590  // std::cout <<"2ndhit?= " << secondHit->is2ndHit() << std::endl;
591  // std::cout <<"1st-otherhitindex= " << firstHit->getOtherHitIndex() << std::endl;
592  // std::cout <<"2nd-otherhitindex= " << secondHit->getOtherHitIndex() << std::endl;
593  // secondHit->setOtherHitIndex(firstHit->getArrayIndex());
594  // firstHit->setOtherHitIndex(secondHit->getArrayIndex());
595  // std::cout <<"1st-otherhitindex= " << firstHit->getOtherHitIndex() << std::endl;
596  // std::cout <<"2nd-otherhitindex= " << secondHit->getOtherHitIndex() << std::endl;
597 
598  //set a relation: CDCSimHit -> CDCHit
599  ++iCDCHits;
600  cdcSimHitsToCDCHits.add(iterSignalMap->second.m_simHitIndex2, iCDCHits);
601  // std::cout << "settdc2 " << firstHit->getTDCCount() << " " << secondHit->getTDCCount() << std::endl;
602 
603  //set a relation: MCParticle -> CDCHit
604  rels = m_simHits[iterSignalMap->second.m_simHitIndex2]->getRelationsFrom<MCParticle>();
605  if (rels.size() != 0) {
606  //assumption: only one MCParticle
607  const MCParticle* mcparticle = rels[0];
608  double weight = rels.weight(0);
609  mcparticle->addRelationTo(secondHit, weight);
610  }
611  } else { //Check the 3rd hit when tdcCount = tdcCount2
612  // std::cout << "tdcCount1=2" << std::endl;
613  if (iterSignalMap->second.m_simHitIndex3 >= 0) {
614  unsigned short tdcCount3 = static_cast<unsigned short>((m_cdcgp->getT0(iterSignalMap->first) - iterSignalMap->second.m_driftTime3) *
615  m_tdcBinWidthInv + 0.5);
616  // std::cout << "tdcCount3= " << tdcCount3 << " " << tdcCount << std::endl;
617  if (tdcCount3 != tdcCount) {
618  CDCHit* secondHit = m_cdcHits.appendNew(tdcCount3, adcCount, iterSignalMap->first, 0, tot);
619  secondHit->set2ndHitFlag();
620  secondHit->setOtherHitIndices(firstHit);
621  // secondHit->setOtherHitIndex(firstHit->getArrayIndex());
622  // firstHit->setOtherHitIndex(secondHit->getArrayIndex());
623  // std::cout <<"2ndhit?= " << secondHit->is2ndHit() << std::endl;
624 
625  //set a relation: CDCSimHit -> CDCHit
626  ++iCDCHits;
627  cdcSimHitsToCDCHits.add(iterSignalMap->second.m_simHitIndex3, iCDCHits);
628  // std::cout << "settdc3 " << firstHit->getTDCCount() << " " << secondHit->getTDCCount() << std::endl;
629 
630  //set a relation: MCParticle -> CDCHit
631  rels = m_simHits[iterSignalMap->second.m_simHitIndex3]->getRelationsFrom<MCParticle>();
632  if (rels.size() != 0) {
633  //assumption: only one MCParticle
634  const MCParticle* mcparticle = rels[0];
635  double weight = rels.weight(0);
636  mcparticle->addRelationTo(secondHit, weight);
637  }
638  }
639  }
640  } //end of checking tdcCount 1=2 ?
641  } //end of 2nd hit setting
642 
643  // std::cout <<"t0= " << m_cdcgp->getT0(iterSignalMap->first) << std::endl;
644  /* unsigned short tdcInCommonStop = static_cast<unsigned short>((m_tdcOffset - iterSignalMap->second.m_driftTime) * m_tdcBinWidthInv);
645  float driftTimeFromTDC = static_cast<float>(m_tdcOffset - (tdcInCommonStop + 0.5)) * m_tdcBinWidth;
646  std::cout <<"driftT bf digitization, TDC in common stop, digitized driftT = " << iterSignalMap->second.m_driftTime <<" "<< tdcInCommonStop <<" "<< driftTimeFromTDC << std::endl;
647  */
648  ++iCDCHits;
649  }
650 
651  //Add crosstalk
652  if (m_addXTalk) addXTalk();
653 
654  // Store the results with trigger time window in a separate array
655  // with corresponding relations.
656  for (iterSignalMapTrg = signalMapTrg.begin(); iterSignalMapTrg != signalMapTrg.end(); ++iterSignalMapTrg) {
657  /*
658  unsigned short adcCount = getADCCount(iterSignalMapTrg->first.first, iterSignalMapTrg->second.m_charge, 1., 0.);
659  unsigned short boardID = m_cdcgp->getBoardID(iterSignalMapTrg->first.first);
660  // B2DEBUG(29, "boardID= " << boardID);
661  const unsigned short adcTh = m_useDB4FEE ? m_adcThresh[boardID] : m_adcThreshold;
662  if (adcCount < adcTh) adcCount = 0;
663  */
664  // unsigned short adcCount = getADCCount(iterSignalMapTrg->second.m_charge);
665  unsigned short adcCount = iterSignalMapTrg->second.m_charge;
666  unsigned short tdcCount =
667  static_cast<unsigned short>((m_cdcgp->getT0(iterSignalMapTrg->first.first) -
668  iterSignalMapTrg->second.m_driftTime) * m_tdcBinWidthInv + 0.5);
669  const CDCHit* cdcHit = m_cdcHits4Trg.appendNew(tdcCount, adcCount, iterSignalMapTrg->first.first);
670 
671  // relations
672  m_simHits[iterSignalMapTrg->second.m_simHitIndex]->addRelationTo(cdcHit);
673  RelationVector<MCParticle> rels = m_simHits[iterSignalMapTrg->second.m_simHitIndex]->getRelationsFrom<MCParticle>();
674  if (rels.size() != 0) {
675  //assumption: only one MCParticle
676  const MCParticle* mcparticle = rels[0];
677  double weight = rels.weight(0);
678  mcparticle->addRelationTo(cdcHit, weight);
679  }
680  }
681 
682  /*
683  std::cout << " " << std::endl;
684  RelationIndex<MCParticle, CDCHit> mcp_to_hit(mcParticles, cdcHits);
685  if (!mcp_to_hit) B2FATAL("No MCParticle -> CDCHit relation founf!");
686  typedef RelationIndex<MCParticle, CDCHit>::Element RelationElement;
687  int ncdcHits = cdcHits.getEntries();
688  for (int j = 0; j < ncdcHits; ++j) {
689  for (const RelationElement& rel : mcp_to_hit.getElementsTo(cdcHits[j])) {
690  std::cout << j << " " << cdcHits[j]->is2ndHit() <<" "<< rel.from->getIndex() << " " << rel.weight << std::endl;
691  }
692  }
693  */
694 }
695 
696 double CDCDigitizerModule::smearDriftLength(const double driftLength, const double dDdt)
697 {
698  double mean = 0.;
699  double resolution;
700 
701  if (m_useSimpleDigitization) {
702  if (gRandom->Uniform() <= m_fraction) {
703  mean = m_mean1;
704  resolution = m_resolution1;
705  } else {
706  mean = m_mean2;
707  resolution = m_resolution2;
708  }
709  } else {
710  const unsigned short leftRight = m_posFlag;
711  double alpha = m_cdcgp->getAlpha(m_posWire, m_momentum);
712  double theta = m_cdcgp->getTheta(m_momentum);
713  resolution = m_cdcgp->getSigma(driftLength, m_wireID.getICLayer(), leftRight, alpha, theta);
714  resolution *= m_totalFudgeFactor;
715  }
716 
717  //subtract resol. due to digitization, which'll be added later in the digitization
718 
719  double diff = resolution - dDdt * m_tdcResol;
720  if (diff > 0.) {
721  resolution = sqrt(diff * (resolution + dDdt * m_tdcResol));
722  } else {
723  resolution = 0.;
724  }
725 
726 #if defined(CDC_DEBUG)
727  cout << " " << endl;
728  cout << "CDCDigitizerModule::smearDriftLength" << endl;
729  cout << "tdcResol= " << m_tdcResol << endl;
730  cout << "dDdt,resolution= " << dDdt << " " << resolution << endl;
731 #endif
732 
733  // Smear drift length
734  double newDL = gRandom->Gaus(driftLength + mean , resolution);
735  while (newDL <= 0.) newDL = gRandom->Gaus(driftLength + mean, resolution);
736  // cout << "totalFugeF in Digi= " << m_totalFudgeFactor << endl;
737  return newDL;
738 }
739 
740 
741 double CDCDigitizerModule::getdDdt(const double driftL)
742 {
743  //---------------------------------------------------------------------------------
744  // Calculates the 1'st derivative: dD/dt, where D: drift length before smearing; t: drift time
745  //---------------------------------------------------------------------------------
746 
747  double dDdt = m_driftV;
748 
749  if (!m_useSimpleDigitization) {
750  const unsigned short layer = m_wireID.getICLayer();
751  const unsigned short leftRight = m_posFlag;
752  double alpha = m_cdcgp->getAlpha(m_posWire, m_momentum);
753  double theta = m_cdcgp->getTheta(m_momentum);
754  double t = m_cdcgp->getDriftTime(driftL, layer, leftRight, alpha, theta);
755  dDdt = m_cdcgp->getDriftV(t, layer, leftRight, alpha, theta);
756 
757 #if defined(CDC_DEBUG)
758  cout << " " << endl;
759  cout << "CDCDigitizerModule::getdDdt" << endl;
760  cout << "**layer= " << layer << endl;
761  cout << "alpha= " << 180.*alpha / M_PI << std::endl;
762  if (layer == 55) {
763  int lr = 0;
764  for (int i = 0; i < 1000; ++i) {
765  t = 1.0 * i;
766  double d = m_cdcgp->getDriftLength(t, layer, lr, alpha, theta);
767  cout << t << " " << d << endl;
768  }
769 
770  cout << " " << endl;
771 
772  lr = 1;
773  for (int i = 0; i < 100; ++i) {
774  t = 5 * i;
775  double d = m_cdcgp->getDriftLength(t, layer, lr, alpha, theta);
776  cout << t << " " << d << endl;
777  }
778  exit(-1);
779  }
780 #endif
781  }
782 
783  return dDdt;
784 }
785 
786 
787 double CDCDigitizerModule::getDriftTime(const double driftLength, const bool addTof, const bool addDelay)
788 {
789  //---------------------------------------------------------------------------------
790  // Method returning electron drift time (parameters: position in cm)
791  // T(drift) = TOF + T(true drift time) + T(propagation delay in wire) - T(event),
792  // T(event) is a timing of event, which includes a time jitter due to
793  // the trigger system.
794  //---------------------------------------------------------------------------------
795 
796  double driftT = 0.;
797 
798  if (m_useSimpleDigitization) {
799  driftT = (driftLength / Unit::cm) * m_driftVInv;
800 
801 #if defined(CDC_DEBUG)
802  cout << " " << endl;
803  cout << "CDCDigitizerModule::getDriftTime" << endl;
804  cout << "driftvinv= " << m_driftVInv << endl;
805 #endif
806  } else {
807  const unsigned short layer = m_wireID.getICLayer();
808  const unsigned short leftRight = m_posFlag;
809  double alpha = m_cdcgp->getAlpha(m_posWire, m_momentum);
810  double theta = m_cdcgp->getTheta(m_momentum);
811  driftT = m_cdcgp->getDriftTime(driftLength, layer, leftRight, alpha, theta);
812  // std::cout <<"alpha,theta,driftT= " << alpha <<" "<< theta <<" "<< driftT << std::endl;
813  }
814 
815  if (addTof) {
816  driftT += m_flightTime; // in ns
817  }
818 
819  if (addDelay) {
820  //calculate signal propagation length in the wire
821  CDCGeometryPar::EWirePosition set = m_align ? CDCGeometryPar::c_Aligned : CDCGeometryPar::c_Base;
822  TVector3 backWirePos = m_cdcgp->wireBackwardPosition(m_wireID, set);
823 
824  double propLength = (m_posWire - backWirePos).Mag();
825  // if (m_cdcgp->getSenseWireZposMode() == 1) {
826  //TODO: replace the following with cached reference
827  // std::cout << m_gcp->getInstance().getSenseWireZposMode() << std::endl;
828  if (m_gcp->getSenseWireZposMode() == 1) {
829  const unsigned short layer = m_wireID.getICLayer();
830  propLength += m_cdcgp->getBwdDeltaZ(layer);
831  }
832  // B2DEBUG(29, "Propagation in wire length: " << propLength);
833 
834  if (m_useSimpleDigitization) {
835  driftT += (propLength / Unit::cm) * m_propSpeedInv;
836 
837 #if defined(CDC_DEBUG)
838  cout << "pseedinv= " << m_propSpeedInv << endl;
839 #endif
840  } else {
841  const unsigned short layer = m_wireID.getICLayer();
842  driftT += (propLength / Unit::cm) * m_cdcgp->getPropSpeedInv(layer);
843 #if defined(CDC_DEBUG)
844  cout << "layer,pseedinv= " << layer << " " << m_cdcgp->getPropSpeedInv(layer) << endl;
845 #endif
846  }
847  }
848 
849  return driftT;
850 }
851 
852 
853 //unsigned short CDCDigitizerModule::getADCCount(unsigned short layer, unsigned short cell, double edep, double dx, double costh)
854 unsigned short CDCDigitizerModule::getADCCount(const WireID& wid, double dEinGeV, double dx, double costh)
855 {
856  unsigned short adcCount = 0;
857  if (dEinGeV <= 0. || dx <= 0.) return adcCount;
858 
859  // The current value is taken from E. Nakano from some test beam results. This should be a somewhat realistic value, but doesn't need to be exact.
860  // Similar as geometry parameters are for the ideal geometry, not the real one.
861  // const double conversionEDepToADC = (100.0 / 3.2) * 1e6;
862  double conv = (100.0 / 3.2); //keV -> count
863  const unsigned short layer = wid.getICLayer();
864  const unsigned short cell = wid.getIWire();
865  double dEInkeV = dEinGeV / Unit::keV;
866  if (m_spaceChargeEffect) {
867  if (m_useDB4EDepToADC) conv = m_cdcgp->getEDepToADCConvFactor(layer, cell, dEInkeV, dx, costh);
868  } else {
869  if (m_useDB4EDepToADC) conv = m_cdcgp->getEDepToADCMainFactor(layer, cell);
870  }
871  conv *= m_runGain;
872 
873  //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.
874  adcCount = static_cast<unsigned short>(std::round(conv * dEInkeV));
875  return adcCount;
876 }
877 
878 
879 // Set FEE parameters (from DB)
880 void CDCDigitizerModule::setFEElectronics()
881 {
882  const double& off = m_tdcThresholdOffset;
883  const double& gA = m_analogGain;
884  const double& gD = m_digitalGain;
885  const double& adcBW = m_adcBinWidth;
886  const double convF = gA / gD / adcBW;
887  const double el1TrgLatency = m_cdcgp->getMeanT0(); // ns
888  B2DEBUG(29, "L1TRGLatency= " << el1TrgLatency);
889  const double c = 32. * m_tdcBinWidth;
890 
891  if (!m_fEElectronicsFromDB) B2FATAL("No FEEElectronics dbobject!");
892  const CDCFEElectronics& fp = *((*m_fEElectronicsFromDB)[0]);
893  int mode = (fp.getBoardID() == -1) ? 1 : 0;
894  int iNBoards = static_cast<int>(nBoards);
895 
896  //set typical values for all channels first if mode=1
897  if (mode == 1) {
898  for (int bdi = 1; bdi < iNBoards; ++bdi) {
899  m_uprEdgeOfTimeWindow[bdi] = el1TrgLatency - c * (fp.getTrgDelay() + 1);
900  if (m_uprEdgeOfTimeWindow[bdi] < 0.) B2FATAL("CDCDigitizer: Upper edge of time window < 0!");
901  m_lowEdgeOfTimeWindow[bdi] = m_uprEdgeOfTimeWindow[bdi] - c * (fp.getWidthOfTimeWindow() + 1);
902  if (m_lowEdgeOfTimeWindow[bdi] > 0.) B2FATAL("CDCDigitizer: Lower edge of time window > 0!");
903  m_adcThresh[bdi] = fp.getADCThresh();
904  m_tdcThresh[bdi] = convF * (off - fp.getTDCThreshInMV());
905  m_widthOfTimeWindow[bdi] = fp.getWidthOfTimeWindow() + 1;
906  }
907  }
908 
909  //ovewrite values for specific channels if mode=1
910  //set typical values for all channels if mode=0
911  for (const auto& fpp : (*m_fEElectronicsFromDB)) {
912  int bdi = fpp.getBoardID();
913  if (mode == 0 && bdi == 0) continue; //bdi=0 is dummy (not used)
914  if (mode == 1 && bdi == -1) continue; //skip typical case
915  if (bdi < 0 || bdi >= iNBoards) B2FATAL("CDCDigitizer:: Invalid no. of FEE board!");
916  m_uprEdgeOfTimeWindow[bdi] = el1TrgLatency - c * (fpp.getTrgDelay() + 1);
917  if (m_uprEdgeOfTimeWindow[bdi] < 0.) B2FATAL("CDCDigitizer: Upper edge of time window < 0!");
918  m_lowEdgeOfTimeWindow[bdi] = m_uprEdgeOfTimeWindow[bdi] - c * (fpp.getWidthOfTimeWindow() + 1);
919  if (m_lowEdgeOfTimeWindow[bdi] > 0.) B2FATAL("CDCDigitizer: Lower edge of time window > 0!");
920  m_adcThresh[bdi] = fpp.getADCThresh();
921  m_tdcThresh[bdi] = convF * (off - fpp.getTDCThreshInMV());
922  m_widthOfTimeWindow[bdi] = fpp.getWidthOfTimeWindow() + 1;
923  }
924 
925  //debug
926  B2DEBUG(29, "mode= " << mode);
927  for (int bdi = 1; bdi < iNBoards; ++bdi) {
928  B2DEBUG(29, bdi << " " << m_lowEdgeOfTimeWindow[bdi] << " " << m_uprEdgeOfTimeWindow[bdi] << " " << m_adcThresh[bdi] << " " <<
929  m_tdcThresh[bdi]);
930  }
931 }
932 
933 // Set Run-gain (from DB)
934 void CDCDigitizerModule::setRunGain()
935 {
936  m_runGain = (*m_runGainFromDB)->getRunGain();
937  m_runGain *= m_overallGainFactor;
938 }
939 
940 void CDCDigitizerModule::addXTalk()
941 {
942  map<WireID, XTalkInfo> xTalkMap;
943  map<WireID, XTalkInfo> xTalkMap1;
944  map<WireID, XTalkInfo>::iterator iterXTalkMap1;
945 
946  // Loop over all cdc hits to create a xtalk map
947  int OriginalNoOfHits = m_cdcHits.getEntries();
948  B2DEBUG(m_debugLevel4XTalk, "\n \n" << "#CDCHits " << OriginalNoOfHits);
949  for (const auto& aHit : m_cdcHits) {
950  if (m_issue2ndHitWarning && aHit.is2ndHit()) {
951  B2WARNING("2nd TDC hit found, but not ready for it!");
952  }
953  WireID wid(aHit.getID());
954  // B2DEBUG(m_debugLevel4XTalk, "Encoded wireid of current CDCHit: " << wid);
955  short tdcCount = aHit.getTDCCount();
956  short adcCount = aHit.getADCCount();
957  short tot = aHit.getTOT();
958  short board = m_cdcgp->getBoardID(wid);
959  short channel = m_cdcgp->getChannelID(wid);
960  const vector<pair<short, asicChannel>> xTalks = (*m_xTalkFromDB)->getLibraryCrossTalk(channel, tdcCount, adcCount, tot);
961 
962  int nXTalks = xTalks.size();
963  for (int i = 0; i < nXTalks; ++i) {
964  const unsigned short tdcCount4XTalk = xTalks[i].second.TDC;
965  if (i == 0) {
966  B2DEBUG(m_debugLevel4XTalk, "\n" << " signal: " << channel << " " << tdcCount << " " << adcCount << " " << tot);
967  }
968  B2DEBUG(m_debugLevel4XTalk, "xtalk: " << xTalks[i].first << " " << tdcCount4XTalk << " " << xTalks[i].second.ADC << " " <<
969  xTalks[i].second.TOT);
970  WireID widx = m_cdcgp->getWireID(board, xTalks[i].first);
971  if (!m_cdcgp->isBadWire(widx)) { // for non-bad wire
972  if (m_includeEarlyXTalks || (xTalks[i].second.TDC <= tdcCount)) {
973  const double t0 = m_cdcgp->getT0(widx);
974  const double ULOfTDC = (t0 - m_lowEdgeOfTimeWindow[board]) * m_tdcBinWidthInv;
975  const double LLOfTDC = (t0 - m_uprEdgeOfTimeWindow[board]) * m_tdcBinWidthInv;
976  if (LLOfTDC <= tdcCount4XTalk && tdcCount4XTalk <= ULOfTDC) {
977  const unsigned short status = 0;
978  xTalkMap.insert(make_pair(widx, XTalkInfo(tdcCount4XTalk, xTalks[i].second.ADC, xTalks[i].second.TOT, status)));
979  }
980  }
981  // } else {
982  // cout<<"badwire= " << widx.getICLayer() <<" "<< widx.getIWire() << endl;
983  }
984  } //end of xtalk loop
985  } //end of cdc hit loop
986 
987  //Loop over all xtalk hits to creat a new xtalk map with only the fastest hits kept (approx.)
988  B2DEBUG(m_debugLevel4XTalk, "#xtalk hits: " << xTalkMap.size());
989  for (const auto& aHit : xTalkMap) {
990  WireID wid = aHit.first;
991 
992  iterXTalkMap1 = xTalkMap1.find(wid);
993  unsigned short tdcCount = aHit.second.m_tdc;
994  unsigned short adcCount = aHit.second.m_adc;
995  unsigned short tot = aHit.second.m_tot;
996  unsigned short status = aHit.second.m_status;
997 
998  if (iterXTalkMap1 == xTalkMap1.end()) { // new entry
999  xTalkMap1.insert(make_pair(wid, XTalkInfo(tdcCount, adcCount, tot, status)));
1000  // B2DEBUG(m_debugLevel4XTalk, "Creating a new xtalk hit with encoded wire no.: " << wid);
1001  } else { // not new; check if fastest
1002  if (tdcCount < iterXTalkMap1->second.m_tdc) {
1003  iterXTalkMap1->second.m_tdc = tdcCount;
1004  B2DEBUG(m_debugLevel4XTalk, "TDC-count of current xtalk: " << tdcCount);
1005  }
1006  iterXTalkMap1->second.m_adc += adcCount;
1007  iterXTalkMap1->second.m_tot += tot; // approx.
1008  }
1009  } // end of xtalk loop
1010 
1011  //add xtalk in the same way as the beam bg. overlay
1012  B2DEBUG(m_debugLevel4XTalk, "#xtalk1 hits: " << xTalkMap1.size());
1013  for (const auto& aX : xTalkMap1) {
1014  bool append = true;
1015  const unsigned short tdc4Bg = aX.second.m_tdc;
1016  const unsigned short adc4Bg = aX.second.m_adc;
1017  const unsigned short tot4Bg = aX.second.m_tot;
1018  const unsigned short status4Bg = aX.second.m_status;
1019 
1020  for (int iHit = 0; iHit < OriginalNoOfHits; ++iHit) {
1021  CDCHit& aH = *(m_cdcHits[iHit]);
1022  if (aH.getID() != aX.first.getEWire()) { //wire id unmatched
1023  continue;
1024  } else { //wire id matched
1025  append = false;
1026  const unsigned short tdc4Sg = aH.getTDCCount();
1027  const unsigned short adc4Sg = aH.getADCCount();
1028  const unsigned short tot4Sg = aH.getTOT();
1029  // B2DEBUG(m_debuglevel4XTalk, "Sg tdc,adc,tot= " << tdc4Sg << " " << adc4Sg << " " << tot4Sg);
1030  // B2DEBUG(m_debugLevel4XTalk, "Bg tdc,adc,tot= " << tdc4Bg << " " << adc4Bg << " " << tot4Bg);
1031 
1032  // If the BG hit is faster than the true hit, the TDC count is replaced, and
1033  // the relations are removed. ADC counts are summed up.
1034  if (tdc4Sg < tdc4Bg) {
1035  aH.setTDCCount(tdc4Bg);
1036  aH.setStatus(status4Bg);
1037  auto relSimHits = aH.getRelationsFrom<CDCSimHit>();
1038  for (int i = relSimHits.size() - 1; i >= 0; --i) {
1039  relSimHits.remove(i);
1040  }
1041  auto relMCParticles = aH.getRelationsFrom<MCParticle>();
1042  for (int i = relMCParticles.size() - 1; i >= 0; --i) {
1043  relMCParticles.remove(i);
1044  }
1045  }
1046 
1047  aH.setADCCount(adc4Sg + adc4Bg);
1048 
1049  //Set TOT for signal+background case. It is assumed that the start timing
1050  //of a pulse (input to ADC) is given by the TDC-count. This is an
1051  //approximation becasue analog (for ADC) and digital (for TDC) parts are
1052  //different in the front-end electronics.
1053  unsigned short s1 = tdc4Sg; //start time of 1st pulse
1054  unsigned short s2 = tdc4Bg; //start time of 2nd pulse
1055  unsigned short w1 = tot4Sg; //its width
1056  unsigned short w2 = tot4Bg; //its width
1057  if (tdc4Sg < tdc4Bg) {
1058  s1 = tdc4Bg;
1059  w1 = tot4Bg;
1060  s2 = tdc4Sg;
1061  w2 = tot4Sg;
1062  }
1063  w1 *= 32;
1064  w2 *= 32;
1065  const unsigned short e1 = s1 - w1; //end time of 1st pulse
1066  const unsigned short e2 = s2 - w2; //end time of 2nd pulse
1067  // B2DEBUG(m_debuglevel4Xtalk, "s1,e1,w1,s2,e2,w2= " << s1 << " " << e1 << " " << w1 << " " << s2 << " " << e2 << " " << w2);
1068 
1069  double pulseW = w1 + w2;
1070  if (e1 <= e2) {
1071  pulseW = w1;
1072  } else if (e1 <= s2) {
1073  pulseW = s1 - e2;
1074  }
1075 
1076  unsigned short board = m_cdcgp->getBoardID(aX.first);
1077  aH.setTOT(std::min(std::round(pulseW / 32.), static_cast<double>(m_widthOfTimeWindow[board])));
1078  B2DEBUG(m_debugLevel4XTalk, "replaced tdc,adc,tot,wid,status= " << aH.getTDCCount() << " " << aH.getADCCount() << " " << aH.getTOT()
1079  <<
1080  " " << aH.getID() << " " << aH.getStatus());
1081  break;
1082  }
1083  } //end of cdc hit loop
1084 
1085  if (append) {
1086  m_cdcHits.appendNew(tdc4Bg, adc4Bg, aX.first, status4Bg, tot4Bg);
1087  B2DEBUG(m_debugLevel4XTalk, "appended tdc,adc,tot,wid,status= " << tdc4Bg << " " << adc4Bg << " " << tot4Bg << " " << aX.first <<
1088  " " <<
1089  status4Bg);
1090  }
1091  } //end of x-talk loop
1092  B2DEBUG(m_debugLevel4XTalk, "original #hits, #hits= " << OriginalNoOfHits << " " << m_cdcHits.getEntries());
1093 }
Belle2::RelationVector::size
size_t size() const
Get number of relations.
Definition: RelationVector.h:98
Belle2::CDCHit::setStatus
void setStatus(unsigned short status)
Setter for CDCHit status.
Definition: CDCHit.h:117
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:72
Belle2::WireID
Class to identify a wire inside the CDC.
Definition: WireID.h:44
Belle2::CDCHit::setOtherHitIndices
void setOtherHitIndices(CDCHit *otherHit)
Setter for the other hit indices.
Definition: CDCHit.h:158
Belle2::CDCHit::getID
unsigned short getID() const
Getter for encoded wire number.
Definition: CDCHit.h:204
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::CDCHit
Class containing the result of the unpacker in raw data and the result of the digitizer in simulation...
Definition: CDCHit.h:51
Belle2::CDCHit::getTDCCount
short getTDCCount() const
Getter for TDC count.
Definition: CDCHit.h:230
Belle2::DBArray
Class for accessing arrays of objects in the database.
Definition: DBArray.h:36
Belle2::CDC::CDCGeometryPar::getNominalPropSpeed
double getNominalPropSpeed() const
Return the nominal propagation speed of the sense wire (default: 27.25 cm/nsec).
Definition: CDCGeometryPar.h:752
Belle2::RelationsInterface::addRelationTo
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).
Definition: RelationsObject.h:144
Belle2::CDC::CDCGeometryPar::getTdcBinWidth
double getTdcBinWidth() const
Return TDC bin width (nsec).
Definition: CDCGeometryPar.h:732
Belle2::CDCDigitizerModule::SignalInfo
Structure for saving the signal information.
Definition: CDCDigitizerModule.h:237
Belle2::CDCSimHit
Example Detector.
Definition: CDCSimHit.h:33
Belle2::CDC::ClosestApproach
double ClosestApproach(const TVector3 &bwp, const TVector3 &fwp, const TVector3 &posIn, const TVector3 &posOut, TVector3 &hitPosition, TVector3 &wirePosition)
Returns a closest distance between a track and a wire.
Definition: ClosestApproach.cc:27
Belle2::CDCHit::getStatus
unsigned short getStatus() const
Getter for CDCHit status.
Definition: CDCHit.h:210
Belle2::CDCHit::getTOT
unsigned short getTOT() const
Getter for TOT.
Definition: CDCHit.h:259
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:31
Belle2::CDCHit::getADCCount
unsigned short getADCCount() const
Getter for integrated charge.
Definition: CDCHit.h:241
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::CDC::CDCGeometryPar
The Class for CDC Geometry Parameters.
Definition: CDCGeometryPar.h:75
Belle2::RelationVector
Class for type safe access to objects that are referred to in relations.
Definition: DataStore.h:38
Belle2::CDCHit::setTDCCount
void setTDCCount(short tdcCount)
Setter for TDC count.
Definition: CDCHit.h:139
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CDCHit::set2ndHitFlag
void set2ndHitFlag()
Setter for 2nd hit flag.
Definition: CDCHit.h:124
Belle2::CDCDigitizerModule
The Class for Detailed Digitization of CDC.
Definition: CDCDigitizerModule.h:60
Belle2::RelationsInterface::getRelationsFrom
RelationVector< FROM > getRelationsFrom(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from another store array to this object.
Definition: RelationsObject.h:214
Belle2::CDCHit::setTOT
void setTOT(unsigned short tot)
Setter for TOT.
Definition: CDCHit.h:171
Belle2::RelationArray::add
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
Definition: RelationArray.h:291
Belle2::CDC::CDCGeometryPar::getNominalDriftV
double getNominalDriftV() const
Return the nominal drift velocity of He-ethane gas (default: 4.0x10^-3 cm/nsec).
Definition: CDCGeometryPar.h:742
Belle2::CDCHit::setADCCount
void setADCCount(unsigned short adcCount)
Setter for ADC count.
Definition: CDCHit.h:146
Belle2::CDC::CDCGeometryPar::EWirePosition
EWirePosition
Wire position set.
Definition: CDCGeometryPar.h:80
Belle2::CDCDigitizerModule::XTalkInfo
Structure for saving the x-talk information.
Definition: CDCDigitizerModule.h:257
Belle2::RelationVector::weight
float weight(int index) const
Get weight with index.
Definition: RelationVector.h:120
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::WireID::getIWire
unsigned short getIWire() const
Getter for wire within the layer.
Definition: WireID.h:155
Belle2::WireID::getICLayer
unsigned short getICLayer() const
Getter for continuous layer numbering.
Definition: WireID.cc:26
Belle2::DBAccessorBase::addCallback
void addCallback(std::function< void(const std::string &)> callback, bool onDestruction=false)
Add a callback method.
Definition: DBAccessorBase.h:108
Belle2::CDCFEElectronics
Database object for Fron-endt electronics params.
Definition: CDCFEElectronics.h:30