All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CRTDetSimAlg.cc
Go to the documentation of this file.
1 #ifndef IC_CRTDETSIMALG_CC
2 #define IC_CRTDETSIMALG_CC
3 
4 #include "icaruscode/CRT/CRTUtils/CRTDetSimAlg.h"
5 
6 namespace icarus{
7  namespace crt {
8 
9  bool TimeOrderCRTData(std::pair<ChanData, AuxDetIDE> crtdat1,
10  std::pair<ChanData, AuxDetIDE> crtdat2) {
11  return ( crtdat1.first.ts < crtdat2.first.ts );
12  }//TimeOrderCRTData()
13 
14  //-------------------------------------------------------------------------------------------
15  //constructor
16  CRTDetSimAlg::CRTDetSimAlg(fhicl::ParameterSet const & p, CLHEP::HepRandomEngine& engine) :
17  fNsim_m(0), fNsim_d(0), fNsim_c(0), fNchandat_m(0), fNchandat_d(0), fNchandat_c(0),
18  fNmissthr_c(0), fNmissthr_d(0), fNmissthr_m(0), fNmiss_strcoin_c(0), fNdual_m(0),
19  fHasFilledTaggers(false), fRandEngine(engine)
20  {
21 
22  this->reconfigure(p);
23  fRegCounts.clear();
24  fRegions.clear();
25  fTaggers.clear();
26  fCrtutils = new CRTCommonUtils();
27  }
28  //----------------------------------------------------------------------
29  //getting parameter values from FHiCL (called by constructor)
30  void CRTDetSimAlg::reconfigure(fhicl::ParameterSet const & p) {
31  fVerbose = p.get<bool>("Verbose");
32  fUltraVerbose = p.get<bool>("UltraVerbose");
33  fGlobalT0Offset = p.get<double>("GlobalT0Offset");
34  fTDelayNorm = p.get<double>("TDelayNorm");
35  fTDelayShift = p.get<double>("TDelayShift");
36  fTDelaySigma = p.get<double>("TDelaySigma");
37  fTDelayOffset = p.get<double>("TDelayOffset");
38  fTDelayRMSGausNorm = p.get<double>("TDelayRMSGausNorm");
39  fTDelayRMSGausShift = p.get<double>("TDelayRMSGausShift");
40  fTDelayRMSGausSigma = p.get<double>("TDelayRMSGausSigma");
41  fTDelayRMSExpNorm = p.get<double>("TDelayRMSExpNorm");
42  fTDelayRMSExpShift = p.get<double>("TDelayRMSExpShift");
43  fTDelayRMSExpScale = p.get<double>("TDelayRMSExpScale");
44  fPropDelay = p.get<double>("PropDelay");
45  fPropDelayError = p.get<double>("PropDelayError");
46  fTResInterpolator = p.get<double>("TResInterpolator");
47  fUseEdep = p.get<bool>("UseEdep");
48  fQ0 = p.get<double>("Q0");
49  fQPed = p.get<double>("QPed");
50  fQSlope = p.get<double>("QSlope");
51  fQRMS = p.get<double>("QRMS");
52  fQMax = p.get<uint16_t>("QMax");
53  fQThresholdC = p.get<uint16_t>("QThresholdC");
54  fQThresholdM = p.get<uint16_t>("QThresholdM");
55  fQThresholdD = p.get<uint16_t>("QThresholdD");
56  fStripCoincidenceWindow = p.get<double>("StripCoincidenceWindow");
57  fApplyStripCoinC = p.get<bool>("ApplyStripCoincidenceC");
58  fApplyCoincidenceC = p.get<bool>("ApplyCoincidenceC");
59  fApplyCoincidenceM = p.get<bool>("ApplyCoincidenceM");
60  fApplyCoincidenceD = p.get<bool>("ApplyCoincidenceD");
61  fLayerCoincidenceWindowC = p.get<double>("LayerCoincidenceWindowC");
62  fLayerCoincidenceWindowM = p.get<double>("LayerCoincidenceWindowM");
63  fLayerCoincidenceWindowD = p.get<double>("LayerCoincidenceWindowD");
64  //fAbsLenEffC = p.get<double>("AbsLenEffC");
65  //fAbsLenEffM = p.get<double>("AbsLenEffM");
66  //fAbsLenEffD = p.get<double>("AbsLenEffD");
67  fDeadTime = p.get<double>("DeadTime");
68  fBiasTime = p.get<double>("BiasTime");
69  fUseBirks = p.get<bool>("UseBirks");
70  fKbirks = p.get<double>("Kbirks");
71  }//CRTDetSim::reconfigure()
72 
73  //-------------------------------------------------------------------------------------------------------------
74  // function called after loop over AuxDetSimChannels where FillTaggers was used to perform first detsim step.
75  // this function applies trigger logic, deadtime, and close-in-time signal biasing effects. it produces the
76  // "triggered data" products which make it into the event
77  vector<pair<CRTData,vector<sim::AuxDetIDE>>> CRTDetSimAlg::CreateData()
78  {
79  vector<pair<CRTData, vector<AuxDetIDE>>> dataCol;
80 
81  if(!fHasFilledTaggers)
82  throw cet::exception("CRTDetSimAlg") << "CreateData() called with empty taggers map!";
83  if(fTaggers.size()==0)
84  return dataCol;
85 
86  int eve=1;
87 
88  int ncombined_c=0, ncombined_m=0, ncombined_d=0; //channel signals close in time, biasing first signal entering into track and hold circuit
89  int nmiss_lock_c=0, nmiss_lock_d=0, nmiss_lock_m=0; //channel signals missed due to channel already locked, but before readout (dead time)
90  int nmiss_dead_c=0, nmiss_dead_d=0, nmiss_dead_m=0; //track losses from deadtime
91  int nmiss_opencoin_c = 0, nmiss_opencoin_d = 0; //track losses where no coincidence possible
92  int nmiss_coin_c = 0;
93  int nmiss_coin_d = 0;
94  int nmiss_coin_m = 0;
95  int event = 0; //FEB entry number (can have multiple per art event)
96  int nhit_m=0, nhit_c=0, nhit_d=0; //channels with signals contained in readout of FEB
97  int neve_m=0, neve_c=0, neve_d=0; //no. readouts of FEBs for each subsystem
98 
99  // Loop over all FEBs (key for taggers) with a hit and check coincidence requirement.
100  // For each FEB, find channel providing trigger and determine if
101  // other hits are in concidence with the trigger (keep)
102  // or if hits occur during R/O (dead time) (lost)
103  // or if hits are part of a different event (keep for now)
104  // First apply dead time correction, biasing effect if configured to do so.
105  // Front-end logic: For CERN or DC modules require at least one hit in each X-X layer.
106  if (fUltraVerbose) std::cout << '\n' << "about to loop over taggers (size " << fTaggers.size() << " )" << std::endl;
107 
108  for (auto trg : fTaggers) {
109  //if(trg.second.data.size()!=trg.second.ide.size())
110  // std::cout << "WARNING DATA AND INDEX VECTOR SIZE MISMATCH!" << std::endl;
111 
112  event = 0;
113  ChanData *chanTrigData(nullptr);
114  ChanData *chanTmpData(nullptr);
115  set<int> trackNHold = {}; //set of channels close in time to triggered readout above threshold
116  set<int> layerNHold = {}; //layers with channels above threshold (used for checking layer-layer coincidence)
117  if((int)trg.first==INT_MAX) std::cout << "WARNING: bad mac5 found!" << std::endl;
118  bool minosPairFound = false, istrig=false;
119  vector<ChanData> passingData; //data to be included in "readout" of FEB
120  vector<AuxDetIDE> passingIDE;
121  uint64_t ttrig=0.0, ttmp=0.0; //time stamps on trigger channel, channel considered as part of readout
122  AuxDetIDE idetmp;
123  size_t trigIndex = 0;
124  uint16_t adc[64];
125 
126  //check "open" coincidence (just check if coincdence possible w/hit in both layers)
127  if (trg.second.type=='c' && fApplyCoincidenceC && trg.second.layerid.size()<2) {
128  nmiss_opencoin_c++;
129  continue;
130  }
131  if (trg.second.type=='d' && fApplyCoincidenceD && trg.second.layerid.size()<2) {
132  nmiss_opencoin_d++;
133  continue;
134  }
135 
136  //time order ChannalData objects in this FEB by T0
137  std::sort((trg.second.data).begin(),(trg.second.data).end(),TimeOrderCRTData);
138 
139  if (fUltraVerbose) std::cout << "processing data for FEB " << (int)trg.first << " with "
140  << trg.second.data.size() << " entries..." << '\n'
141  << " type: " << trg.second.type << '\n'
142  << " region: " << trg.second.reg << '\n'
143  << " layerID: " << *trg.second.layerid.begin() << '\n' << std::endl;
144 
145  //outer (primary) loop over all data products for this FEB
146  for ( size_t i=0; i< trg.second.data.size(); i++ ) {
147 
148  //get data for earliest entry
149  if(i==0) {
150  chanTrigData = &(trg.second.data[0].first);
151  ttrig = chanTrigData->ts; //time stamp [ns]
152  trackNHold.insert(chanTrigData->channel);
153  layerNHold.insert(trg.second.chanlayer[chanTrigData->channel]);
154  passingData.push_back(*chanTrigData);
155  passingIDE.push_back(trg.second.data[0].second);
156  ttmp = ttrig;
157  idetmp = passingIDE.back();
158  if(trg.second.type!='m')
159  continue;
160  }
161 
162  else {
163  chanTmpData = &(trg.second.data[i].first);
164  ttmp = chanTmpData->ts;
165  idetmp = trg.second.data[i].second;
166  }
167 
168  //for C and D modules only and coin. enabled, if assumed trigger channel has no coincidence
169  // set trigger channel to tmp channel and try again
170  if ( layerNHold.size()==1 &&
171  ( (trg.second.type=='c' && fApplyCoincidenceC && ttmp-ttrig>fLayerCoincidenceWindowC) ||
172  (trg.second.type=='d' && fApplyCoincidenceD && ttmp-ttrig>fLayerCoincidenceWindowD )))
173  {
174  trigIndex++;
175  chanTrigData = &(trg.second.data[trigIndex].first);
176  i = trigIndex; //+1
177  ttrig = chanTrigData->ts;
178  trackNHold.clear();
179  layerNHold.clear();
180  passingData.clear();
181  passingIDE.clear();
182  trackNHold.insert(chanTrigData->channel);
183  layerNHold.insert(trg.second.chanlayer[chanTrigData->channel]);
184  passingData.push_back(*chanTrigData);
185  passingIDE.push_back(trg.second.data[trigIndex].second);
186  if(trg.second.type=='c') nmiss_coin_c++;
187  if(trg.second.type=='d') nmiss_coin_d++;
188  continue;
189  }
190 
191  if(fUltraVerbose)
192  std::cout << "checking for coincidence..." << std::endl;
193  //check if coincidence condtion met
194  //for c and d modules, just need time stamps within tagger obj
195  //for m modules, need to check coincidence with other tagger objs
196  if (trg.second.type=='m' && !minosPairFound && fApplyCoincidenceM) {
197  for (auto trg2 :fTaggers) {
198 
199  if( trg2.second.type!='m' || //is other mod 'm' type
200  trg.second.modid == trg2.second.modid || //other mod not same as this one
201  trg.second.reg != trg2.second.reg || //other mod is in same region
202  *trg2.second.layerid.begin() == *trg.second.layerid.begin()) //modules are in adjacent layers
203  continue;
204 
205  //find entry within coincidence window starting with this FEB's
206  //triggering channel in coincidence candidate's FEB
207  for ( size_t j=0; j< trg2.second.data.size(); j++ ) {
208  uint64_t t2tmp = trg2.second.data[j].first.ts; //in us
209  if ( util::absDiff(t2tmp,ttrig) < fLayerCoincidenceWindowM) {
210  minosPairFound = true;
211  break;
212  }
213  }
214  //we found a valid pair so move on to next step
215  if (minosPairFound)
216  break;
217  }//inner loop over febs (taggers)
218 
219  //if no coincidence pairs found, reinitialize and move to next FEB
220  if(!minosPairFound) {
221  if(fUltraVerbose)
222  std::cout << "MINOS pair NOT found! Skipping to next FEB..." << std::endl;
223  if(trg.second.data.size()==1) continue;
224  trigIndex++;
225  chanTrigData = &(trg.second.data[trigIndex].first);
226  i = trigIndex;//+1;
227  ttrig = chanTrigData->ts;
228  trackNHold.clear();
229  layerNHold.clear();
230  passingData.clear();
231  passingIDE.clear();
232  trackNHold.insert(chanTrigData->channel);
233  layerNHold.insert(trg.second.chanlayer[chanTrigData->channel]);
234  passingData.push_back(*chanTrigData);
235  passingIDE.push_back(trg.second.data[trigIndex].second);
236  nmiss_coin_m++;
237  continue;
238  }
239  else
240  istrig = true;
241  }//if minos module and no pair yet found
242 
243  if(fUltraVerbose)
244  std::cout << "done checking coinceidence...moving on to latency effects..." << std::endl;
245  if(!minosPairFound && trg.second.type=='m')
246  std::cout << "WARNING: !minosPairFound...should not see this message!" << std::endl;
247 
248  int adctmp = 0;
249  //currently assuming layer coincidence window is same as track and hold window (FIX ME!)
250  if (i>0 && ((trg.second.type=='c' && ttmp < ttrig + fLayerCoincidenceWindowC) ||
251  (trg.second.type=='d' && ttmp < ttrig + fLayerCoincidenceWindowD) ||
252  (trg.second.type=='m' && ttmp < ttrig + fLayerCoincidenceWindowM)) )
253  {
254 
255  //if channel not locked
256  if ((trackNHold.insert(chanTmpData->channel)).second) {
257 
258  //channel added to vector of ChannelData for current FEB readout
259  passingData.push_back(*chanTmpData);
260  passingIDE.push_back(idetmp);
261 
262  //if not m module, check to see if strip is first in time in adjacent layer w.r.t. trigger strip
263  if (layerNHold.insert(trg.second.chanlayer[chanTmpData->channel]).second
264  && trg.second.type != 'm')
265  {
266  //flagging strips which produce triggering condition (layer-layer coincidence)
267  istrig = true;
268  }
269  } //end if channel not locked
270 
271  //check for signal biasing
272  else if (ttmp < ttrig + fBiasTime) {
273  for(size_t dat=0; dat<passingData.size(); dat++) {
274  if(passingData[dat].channel!=chanTmpData->channel)
275  continue;
276  adctmp = passingData[dat].adc;
277  adctmp += chanTmpData->adc;
278  if(adctmp>fQMax) adctmp = fQMax;
279  passingData[dat].adc = adctmp;
280  passingIDE.push_back(idetmp);
281  break;
282 
283  }
284 
285  switch (trg.second.type) {
286  case 'c' : ncombined_c++; break;
287  case 'd' : ncombined_d++; break;
288  case 'm' : ncombined_m++; break;
289  }
290  } //channel is locked but hit close in time to bias pulse height
291 
292  else switch (trg.second.type) {
293  case 'c' : nmiss_lock_c++; break;
294  case 'd' : nmiss_lock_d++; break;
295  case 'm' : nmiss_lock_m++; break;
296  } //data is discarded (not writeable)
297 
298  }//if hits inside track and hold window
299 
300  else if ( i>0 && ttmp <= ttrig + fDeadTime ) {
301  switch (trg.second.type) {
302  case 'c' : nmiss_dead_c++; break;
303  case 'd' : nmiss_dead_d++; break;
304  case 'm' : nmiss_dead_m++; break;
305  }
306  } // hits occuring during digitization lost (dead time)
307 
308  //"read out" data for this event, first hit after dead time as next trigger channel
309  else if ( ttmp > ttrig + fDeadTime) {
310 
311  if(istrig)
312  {int regnum = fCrtutils->AuxDetRegionNameToNum(trg.second.reg);
313  if( (fRegions.insert(regnum)).second) fRegCounts[regnum] = 1;
314  else fRegCounts[regnum]++;
315 
316  if (fUltraVerbose) {
317  std::cout << "creating CRTData product just after deadtime" << '\n'
318  << " event: " << eve << '\n'
319  << " mac5: " << trg.first << '\n'
320  << " FEB entry: " << event << '\n'
321  << " trig time: " << ttrig << '\n'
322  << " trig channel: " << chanTrigData->channel << '\n'
323  << " passing data: " << std::endl;
324  for(size_t i=0; i<passingData.size(); i++) {
325  std::cout
326  << " index: " << i << '\n'
327  << " channel: " << passingData[i].channel << '\n'
328  << " t0: " << passingData[i].ts << '\n'
329  << " adc: " << passingData[i].adc << std::endl;
330  }
331  }
332 
333  if(passingData.size()>passingIDE.size())
334  std::cout << "data/IDE size mismatch!" << passingData.size()-passingIDE.size() << std::endl;
335  FillAdcArr(passingData,adc);
336  dataCol.push_back(std::make_pair(
337  FillCRTData(trg.first,event,ttrig,ttrig,adc),
338  passingIDE) );
339 
340  if (fUltraVerbose) std::cout << " ...success!" << std::endl;
341  event++;
342  if (trg.second.type=='c') {neve_c++; nhit_c+=passingData.size(); }
343  if (trg.second.type=='d') {neve_d++; nhit_d+=passingData.size(); }
344  if (trg.second.type=='m') {neve_m++; nhit_m+=passingData.size(); } }
345  trigIndex = i;
346  ttrig = ttmp;
347  chanTrigData = chanTmpData;
348  passingData.clear();
349  passingIDE.clear();
350  trackNHold.clear();
351  layerNHold.clear();
352  passingData.push_back(*chanTrigData);
353  passingIDE.push_back(idetmp);
354  trackNHold.insert(chanTrigData->channel);
355  layerNHold.insert(trg.second.chanlayer[chanTmpData->channel]);
356  minosPairFound = false;
357  istrig = false;
358  }
359 
360  if (!(ttmp > ttrig + fDeadTime) && i==trg.second.data.size()-1 && istrig) {
361 
362  if (fUltraVerbose) {
363  std::cout << "creating CRTData product at end of FEB events..." << '\n'
364  << " event: " << eve << '\n'
365  << " mac5: " << trg.first << '\n'
366  << " FEB entry: " << event << '\n'
367  << " trig time: " << ttrig << '\n'
368  << " trig channel: " << chanTrigData->channel << '\n'
369  << " passing data: " << std::endl;
370  for(size_t i=0; i<passingData.size(); i++) {
371  std::cout
372  << " index: " << i << '\n'
373  << " channel: " << passingData[i].channel << '\n'
374  << " t0: " << passingData[i].ts << '\n'
375  << " adc: " << passingData[i].adc << std::endl;
376  }
377  }
378 
379  int regnum = fCrtutils->AuxDetRegionNameToNum(trg.second.reg);
380  if( (fRegions.insert(regnum)).second) fRegCounts[regnum] = 1;
381  else fRegCounts[regnum]++;
382 
383  if(passingData.size()>passingIDE.size())
384  std::cout << "data/IDE size mismatch! " << passingData.size()-passingIDE.size() << std::endl;
385  FillAdcArr(passingData,adc);
386  dataCol.push_back( std::make_pair(
387  FillCRTData(trg.first,event,ttrig,ttrig,adc),
388  passingIDE) );
389  if (fUltraVerbose)
390  std::cout << " ...success!" << std::endl;
391  event++;
392  if (trg.second.type=='c') {neve_c++; nhit_c+=passingData.size(); }
393  if (trg.second.type=='d') {neve_d++; nhit_d+=passingData.size(); }
394  if (trg.second.type=='m') {neve_m++; nhit_m+=passingData.size(); }
395 
396  }//if last event and not already written
397  }//for data entries (hits)
398 
399  if(fUltraVerbose) std::cout << " outside loop over FEB data entries...moving on to next FEB..." << std::endl;
400 
401  } // for taggers
402 
403  if (fVerbose) {
404  std::cout << "|---------------------- CRTDetSim Summary ----------------------|" << '\n'
405  << " - - - EDeps from AuxDetIDE's - - -" << '\n'
406  << " CERN sim strip hits: " << fNsim_c << '\n'
407  << " MINOS sim strip hits: " << fNsim_m << '\n'
408  << " DC sim strip hits: " << fNsim_d << '\n' << '\n'
409  << " - - - Single Channel Threshold - - -" << '\n'
410  << " Pass:" << '\n'
411  << " CERN channel signals > thresh: " << fNchandat_c << '\n'
412  << " MINOS channel signals > thresh: " << fNchandat_m << '\n'
413  << " both ends > thresh: " << fNdual_m << '\n'
414  << " DC channel signals > thresh: " << fNchandat_d << '\n'
415  << " Lost:" << '\n'
416  << " CERN channel signals < thresh: " << fNmissthr_c << '\n'
417  << " MINOS channel signals < thresh: " << fNmissthr_m << '\n'
418  << " DC channel signals < thresh: " << fNmissthr_d << '\n' << '\n'
419  << " - - - System Specifc Coincidence Loss- - -" << '\n'
420  << " CERN fiber-fiber: " << fNmiss_strcoin_c << " (" << fNmiss_strcoin_c*2 << " channels signals lost)" << '\n'
421  << " CERN open layer-layer: " << nmiss_opencoin_c << '\n'
422  << " DC open layer-layer: " << nmiss_opencoin_d << '\n'
423  << " CERN layer-layer: " << nmiss_coin_c << " (" << 100.0*nmiss_coin_c/fNchandat_c << "%)" << '\n'
424  << " MINOS layer-layer: " << nmiss_coin_m << " (" << 100.0*nmiss_coin_m/fNchandat_m << "%)" << '\n'
425  << " DC layer-layer: " << nmiss_coin_d << " (" << 100.0*nmiss_coin_d/fNchandat_d << "%)" << '\n' << '\n'
426  << " - - - Front-End Electronics Effects - - -" << '\n'
427  << " CERN trackNHold loss: " << nmiss_lock_c << " (" << 100.0*nmiss_lock_c/fNchandat_c << "%)" << '\n'
428  << " MINOS trackNHold loss: " << nmiss_lock_m << " (" << 100.0*nmiss_lock_m/fNchandat_m << "%)" << '\n'
429  << " DC trackNHold loss: " << nmiss_lock_d << " (" << 100.0*nmiss_lock_d/fNchandat_d << "%)" << '\n'
430  << " CERN combined: " << ncombined_c << " (" << 100.0*ncombined_c/fNchandat_c << "%)" << '\n'
431  << " MINOS combined: " << ncombined_m << " (" << 100.0*ncombined_m/fNchandat_m << "%)" << '\n'
432  << " DC combined: " << ncombined_d << " (" << 100.0*ncombined_d/fNchandat_d << "%)" << '\n'
433  << " CERN deadTime loss: " << nmiss_dead_c << " (" << 100.0*nmiss_dead_c/fNchandat_c << "%)" << '\n'
434  << " MINOS deadTime loss: " << nmiss_dead_m << " (" << 100.0*nmiss_dead_m/fNchandat_m << "%)" << '\n'
435  << " DC deadTime loss: " << nmiss_dead_d << " (" << 100.0*nmiss_dead_d/fNchandat_d << "%)" << '\n' << '\n'
436  << " - - - Passing Hits Pushed To Event - - -" << '\n'
437  << " CERN strip hits remaining: " << nhit_c+ncombined_c << " (" << 100.0*(nhit_c+ncombined_c)/fNchandat_c << "%)" << '\n'
438  << " MINOS strip hits remaining: " << nhit_m+ncombined_m << " (" << 100.0*(nhit_m+ncombined_m)/fNchandat_m << "%)" << '\n'
439  << " DC strip hits remaining: " << nhit_d+ncombined_d << " (" << 100.0*(nhit_d+ncombined_d)/fNchandat_d << "%)" << '\n'
440  << " CERN channel signals generated: " << nhit_c << '\n'
441  << " MINOS channel signals generated: " << nhit_m << '\n'
442  << " DC channel signals generated: " << nhit_d << '\n'
443  << " CERN events: " << neve_c << '\n'
444  << " MINOS events: " << neve_m << '\n'
445  << " DC events: " << neve_d << '\n'
446  << " Total pushes to event: " << dataCol.size() << std::endl;
447 
448  if(dataCol.size()>0) {
449  std::map<int,int>::iterator it = fRegCounts.begin();
450  std::cout << '\n' << "FEB events per CRT region: " << '\n' << std::endl;
451 
452  do {
453  std::cout << fCrtutils->GetRegionNameFromNum((*it).first) << ": , events: " << (*it).second << '\n' << std::endl;
454  it++;
455  }
456  while ( it != fRegCounts.end() );
457  }//if any CRTData
458  } //if verbose
459 
460  return dataCol;
461 
462  }//end CreateData()
463 
464  //-----------------------------------------------------------------------------
465  // intented to be called within loop over AuxDetChannels and provided the
466  // AuxDetChannelID, AuxDetSensitiveChannelID, vector of AuxDetIDEs and
467  // the number of ides from the end of the vector to include in the detector sim.
468  // handles deposited energy -> light output at SiPM (including attenuation)
469  // -> PEs from the SiPM (or PMT in case of bottom CRT) with associated time stamps
470  void CRTDetSimAlg::FillTaggers(detinfo::DetectorClocksData const& clockData,
471  const uint32_t adid, const uint32_t adsid, const vector<sim::AuxDetIDE>& ides) {
472 
473  art::ServiceHandle<geo::Geometry> geoService;
474  detinfo::ElecClock trigClock = clockData.TriggerClock();
475  fHasFilledTaggers = true;
476 
477  const geo::AuxDetGeo& adGeo = geoService->AuxDet(adid); //pointer to module object
478  const geo::AuxDetSensitiveGeo& adsGeo = adGeo.SensitiveVolume(adsid); //pointer to strip object
479  const char auxDetType = fCrtutils->GetAuxDetType(adid); //CRT module type (c, d, or m)
480  const string region = fCrtutils->GetAuxDetRegion(adid); //CRT region
481 
482  int layid = INT_MAX; //set to 0 or 1 if layerid determined
483  uint8_t mac5=UINT8_MAX, mac5dual=UINT8_MAX; //front-end board ID, dual for MINOS modules (not cut)
484 
485  // Find the path to the strip geo node, to locate it in the hierarchy
486  set<string> volNames = { adsGeo.TotalVolume()->GetName() };
487  vector<vector<TGeoNode const*> > paths = geoService->FindAllVolumePaths(volNames);
488 
489  string path = "";
490  for (size_t inode=0; inode<paths.at(0).size(); inode++) {
491  path += paths.at(0).at(inode)->GetName();
492  if (inode < paths.at(0).size() - 1) {
493  path += "/";
494  }
495  }
496 
497  TGeoManager* manager = geoService->ROOTGeoManager();
498  manager->cd(path.c_str());
499  TGeoNode* nodeStrip = manager->GetCurrentNode();
500  TGeoNode* nodeInner = manager->GetMother(1);
501  TGeoNode* nodeModule = manager->GetMother(2);
502  double origin[3] = {0, 0, 0};
503  // Module position in parent (tagger) frame
504  double modulePosMother[3]; //position in CRT region volume
505  nodeModule->LocalToMaster(origin, modulePosMother);
506 
507  // strip position in module frame
508  double stripPosMother[3];
509  double stripPosModule[3];
510  nodeStrip->LocalToMaster(origin, stripPosMother);
511  nodeInner->LocalToMaster(stripPosMother,stripPosModule);
512 
513  // Determine layid
514  // for C modules, two diff. lay thick
515  // 1cm for top (y>0) and 1.5cm for bottom (y<0)
516  if (auxDetType == 'c' || auxDetType == 'd')
517  layid = (stripPosModule[1] > 0);
518 
519  if (auxDetType == 'm') {
520  //lateral stacks (6 total, 3 per side)
521  if ( region.find("West")!=std::string::npos || region.find("East")!=std::string::npos ) {
522  layid = ( modulePosMother[0]>0 );
523  }
524  //longitudinal walls
525  if ( region=="South" || region=="North" ) {
526  layid = ( modulePosMother[2]> 0 );
527  }
528  }
529 
530  if(layid==INT_MAX)
531  mf::LogError("CRT")
532  << "layid NOT SET!!!" << '\n'
533  << " ADType: " << auxDetType << '\n'
534  << " ADRegion: " << region;
535 
536  // ------------- loop over AuxDetIDEs ---------
537  // Simulate the CRT response for each hit in this strip
538  for (auto const& ide : ides ) {
539  // count true number of energy deposits for each strip
540  if (auxDetType=='c') fNsim_c++;
541  if (auxDetType=='d') fNsim_d++;
542  if (auxDetType=='m') fNsim_m++;
543 
544  // What is the distance from the hit (centroid of the entry
545  // and exit points) to the readout end?
546  double x = (ide.entryX + ide.exitX) / 2;
547  double y = (ide.entryY + ide.exitY) / 2;
548  double z = (ide.entryZ + ide.exitZ) / 2;
549  double world[3] = {x, y, z};
550  double svHitPosLocal[3];
551  double modHitPosLocal[3];
552  adsGeo.WorldToLocal(world, svHitPosLocal); //position in strip frame (origin at center)
553  adGeo.WorldToLocal(world, modHitPosLocal); //position in module frame (origin at center)
554 
555  //check hit point is contained within the strip according to geometry
556  if ( abs(svHitPosLocal[0])>adsGeo.HalfWidth1()+0.001 ||
557  abs(svHitPosLocal[1])>adsGeo.HalfHeight()+0.001 ||
558  abs(svHitPosLocal[2])>adsGeo.HalfLength()+0.001)
559  mf::LogError("CRT") << "HIT POINT OUTSIDE OF SENSITIVE VOLUME!" << '\n'
560  << " AD: " << adid << " , ADS: " << adsid << '\n'
561  << " Local position (x,y,z): ( " << svHitPosLocal[0]
562  << " , " << svHitPosLocal[1] << " , " << svHitPosLocal[2] << " )";
563 
564  //calculate Birks' correction factor
565  double dl = sqrt(pow(ide.entryX-ide.exitX,2)+pow(ide.entryY-ide.exitY,2)+pow(ide.entryZ-ide.exitZ,2));
566  double birksCorr = 1./(1+fKbirks*ide.energyDeposited/dl);
567 
568  // The expected number of PE, using a quadratic model for the distance
569  // dependence, and scaling linearly with deposited energy.
570  double qr = fUseEdep ? ide.energyDeposited / fQ0 : 1.0;
571  if (auxDetType == 'c'&& layid==0) qr *= 1.5; //c bottom layer strips 50% thicker
572 
573  //longitudinal distance (m) along the strip for fiber atten. calculation
574  //assuming SiPM is on +z end (also -z for m modules)
575  double distToReadout = abs( adsGeo.HalfLength() - svHitPosLocal[2])*0.01;
576  double distToReadout2 = abs(-adsGeo.HalfLength() - svHitPosLocal[2])*0.01;
577 
578  double npeExpected = GetLongAtten(distToReadout) * qr;
579  double npeExpected2 = GetLongAtten(distToReadout2) * qr;
580 
581  //Attenuation factor for transverse propegation in the bulk (c modules only)
582  double abs0=1.0, abs1=1.0;
583  if(auxDetType=='c'){
584  pair<double,double> tmp = GetTransAtten(svHitPosLocal[0]);
585  abs0 = tmp.first;
586  abs1 = tmp.second;
587  }
588 
589  //# photons arriving at SiPM
590  double npeExp0 = npeExpected * abs0;;
591  double npeExp1 = npeExpected * abs1;
592  double npeExp0Dual = npeExpected2 * abs0;
593 
594  //apply Birks' quenching if requested
595  if(fUseBirks){
596  npeExp0 *= birksCorr;
597  npeExp1 *= birksCorr;
598  npeExp0Dual *= birksCorr;
599  }
600 
601  //sanity check on simulated light output
602  if(npeExp0<0||npeExp1<0||npeExp0Dual<0)
603  mf::LogError("CRT") << "NEGATIVE PE!!!!!";
604 
605  // Observed PE (Poisson-fluctuated)
606  int npe0 = CLHEP::RandPoisson::shoot(&fRandEngine, npeExp0);
607  int npe1 = CLHEP::RandPoisson::shoot(&fRandEngine, npeExp1);
608  int npe0Dual = CLHEP::RandPoisson::shoot(&fRandEngine, npeExp0Dual);
609 
610  // Time relative to trigger [ns], accounting for propagation delay and 'walk'
611  // for the fixed-threshold discriminator
612  double tTrue = (ide.entryT + ide.exitT) / 2 + fGlobalT0Offset;
613  if(tTrue<0)
614  mf::LogError("CRTDetSim")
615  << "negative true time passed to time stamp generator!";
616  uint64_t t0 =
617  GetChannelTriggerTicks(trigClock, tTrue, npe0, distToReadout*100);
618  uint64_t t1 =
619  GetChannelTriggerTicks(trigClock, tTrue, npe1, distToReadout*100);
620  uint64_t t0Dual =
621  GetChannelTriggerTicks(trigClock, tTrue, npe0Dual, distToReadout2*100);
622 
623  // Time relative to PPS: Random for now! (FIXME)
624  //int ppsTicks =
625  // CLHEP::RandFlat::shootInt(&fRandEngine, trigClock.Frequency() * 1e6);
626 
627  // SiPM and ADC response: Npe to ADC counts
628  uint16_t q0 =
629  CLHEP::RandGauss::shoot(&fRandEngine, fQPed + fQSlope * npe0, fQRMS * sqrt(npe0));
630  uint16_t q1 =
631  CLHEP::RandGauss::shoot(&fRandEngine, fQPed + fQSlope * npe1, fQRMS * sqrt(npe1));
632  uint16_t q0Dual =
633  CLHEP::RandGauss::shoot(&fRandEngine, fQPed + fQSlope * npe0Dual, fQRMS * sqrt(npe0Dual));
634 
635  if(q0>fQMax) q0 = fQMax;
636  if(q1>fQMax) q1 = fQMax;
637  if(q0Dual>fQMax) q0Dual = fQMax;
638 
639  // Adjacent channels on a strip are numbered sequentially.
640  //
641  // In the AuxDetChannelMapAlg methods, channels are identified by an
642  // AuxDet name (retrievable given the hit AuxDet ID) which specifies a
643  // module, and a channel number from 0 to 32.
644 
645  int channel0ID=0, channel1ID=0;
646  pair<uint8_t,uint8_t> macs = fCrtutils->ADToMac(adid);
647  mac5 = macs.first;
648  int changroup = fCrtutils->ADToChanGroup(adid);
649 
650  switch (auxDetType){
651  case 'c' :
652  channel0ID = 2 * adsid + 0;
653  channel1ID = 2 * adsid + 1;
654  break;
655  case 'd' :
656  channel0ID = adsid;
657  break;
658  case 'm' :
659  channel0ID = adsid/2 + 10*(changroup-1);
660  if (fCrtutils->NFeb(adid)==2) {
661  mac5dual = macs.second;
662  }
663  break;
664  }
665 
666  if (mac5==UINT8_MAX) mf::LogError("CRT") << "mac addrs not set!";
667 
668  // Apply ADC threshold and strip-level coincidence (both fibers fire)
669  if (auxDetType=='c') {
670  if ((fApplyStripCoinC && q0>fQThresholdC && q1>fQThresholdC
671  && util::absDiff(t0,t1)<fStripCoincidenceWindow)||
672  (!fApplyStripCoinC && (q0>fQThresholdC || q1>fQThresholdC)) )
673  {
674  Tagger& tagger = fTaggers[mac5];
675  tagger.layerid.insert(layid);
676  tagger.chanlayer[channel0ID] = layid;
677  tagger.chanlayer[channel1ID] = layid;
678  tagger.reg = region;
679  tagger.type = 'c';
680  tagger.modid = adid;
681  if (q0>fQThresholdC) {
682  tagger.data.push_back(
683  std::make_pair(FillChanData(channel0ID,q0,t0),ide));
684  fNchandat_c++;
685  }
686  else fNmissthr_c++;
687  if (q1>fQThresholdC) {
688  tagger.data.push_back(
689  std::make_pair(FillChanData(channel1ID,q1,t1),ide));
690  fNchandat_c++;
691  }
692  else fNmissthr_c++;
693  //nchandat_c+=2;
694  }
695  }//if fiber-fiber coincidence
696 
697  if (auxDetType=='d' && q0 > fQThresholdD) {
698  Tagger& tagger = fTaggers[mac5];
699  tagger.layerid.insert(layid);
700  tagger.chanlayer[channel0ID] = layid;
701  tagger.reg = region;
702  tagger.type = 'd';
703  tagger.modid = adid;
704  tagger.data.push_back(
705  std::make_pair(FillChanData(channel0ID,q0,t0),ide));
706  fNchandat_d++;
707  }//if one strip above threshold
708 
709  if (auxDetType=='m') {
710  if(q0 > fQThresholdM) {
711  Tagger& tagger = fTaggers[mac5];
712  tagger.layerid.insert(layid);
713  tagger.chanlayer[channel0ID] = layid;
714  tagger.reg = region;
715  tagger.type = 'm';
716  tagger.modid = adid;
717  tagger.data.push_back(
718  std::make_pair(FillChanData(channel0ID,q0,t0),ide));
719  fNchandat_m++;
720  }
721  if(q0Dual > fQThresholdM && fCrtutils->NFeb(adid)==2) {
722  Tagger& tagger = fTaggers[mac5dual];
723  tagger.layerid.insert(layid);
724  tagger.chanlayer[channel0ID] = layid;
725  tagger.reg = region;
726  tagger.type = 'm';
727  tagger.modid = adid;
728  tagger.data.push_back(
729  std::make_pair(FillChanData(channel0ID,q0Dual,t0Dual),ide));
730  fNchandat_m++;
731  }
732  if(q0 > fQThresholdM && q0Dual > fQThresholdM) fNdual_m++;
733  }//if one strip above threshold at either end
734 
735  //counting losses
736  if (auxDetType == 'c') {
737  if ( fApplyStripCoinC && util::absDiff(t0,t1) >= fStripCoincidenceWindow ) fNmiss_strcoin_c++;
738  }
739  if (auxDetType == 'd' && q0 < fQThresholdD) fNmissthr_d++;
740  if (auxDetType == 'm') {
741  if( q0 < fQThresholdM) fNmissthr_m++;
742  if( q0Dual < fQThresholdM && fCrtutils->NFeb(adid)==2) fNmissthr_m++;
743  }
744 
745  //print detsim info (if enabled)
746  if (fUltraVerbose&&
747  (
748  (auxDetType=='c' && q0>fQThresholdC && q1>fQThresholdC && util::absDiff(t0, t1) < fStripCoincidenceWindow) ||
749  (auxDetType=='d' && q0>fQThresholdD ) ||
750  (auxDetType=='m' && (q0>fQThresholdM || q0Dual>fQThresholdM)) ))
751  std::cout << '\n'
752  << "CRT HIT VOL " << (adGeo.TotalVolume())->GetName() << "\n"
753  << "CRT HIT SENSITIVE VOL " << (adsGeo.TotalVolume())->GetName() << "\n"
754  << "CRT HIT AuxDetID " << adid << " / AuxDetSensitiveID " << adsid << "\n"
755  << "CRT module type: " << auxDetType << " , CRT region: " << region << '\n'
756  << "CRT channel: " << channel0ID << " , mac5: " << (int)mac5 << '\n'
757  << "CRT HIT POS (world coords) " << x << " " << y << " " << z << "\n"
758  << "CRT STRIP POS (module coords) " << svHitPosLocal[0] << " " << svHitPosLocal[1] << " " << svHitPosLocal[2] << "\n"
759  << "CRT MODULE POS (region coords) " << modHitPosLocal[0] << " " << modHitPosLocal[1] << " "<< modHitPosLocal[2] << " " << "\n"
760  << "CRT layer ID: " << layid << "\n"
761  << "CRT distToReadout: " << distToReadout << ", distToReadout2: " << distToReadout2 << ", qr = " << qr << '\n'
762  << "CRT abs0: " << abs0 << " , abs1: " << abs1 << '\n'
763  << "CRT npeExpected: " << npeExpected << " , npeExpected2: " << npeExpected2 << '\n'
764  //<< "CRT npeExp0: " << npeExp0 << " , npeExp1: " << npeExp1 << " , npeExp0Dual: " << npeExp0Dual << '\n'
765  << "CRT npeSiPM0: " << npe0 << " , npeSiPM1: " << npe1 << " , npeSiPM0Dual: " << npe0Dual << '\n'
766  << "CRT charge q0: " << q0 << ", q1: " << q1 << '\n'
767  //<< "CRT timing: tTrue: " << tTrue << ", t0: " << t0 << ", t1: " << t1 << ", dt: " << util::absDiff(t0,t1) << '\n'
768  << "CRT timing: tTrue: " << tTrue << ", t0: " << t0 << ", t1: " << t1 << ", |t0-t1|: " << util::absDiff(t0,t1) << '\n'
769  //<< " recoT-trueT = " << t0-tTrue << std::endl;
770  << " recoT0-trueT = " << t0-tTrue << ", recoT1-trueT = " << t1-tTrue << ", recoT0Dual-trueT = " << t0Dual-tTrue << ", recoT0-recoT0Dual = " << t0-t0Dual << std::endl;
771 
772  }//for AuxDetIDEs
773 
774  } //end FillTaggers
775 
776  //---------------------------------------------------------------
777  // for CERN modules only. Given position in local strip coordinates,
778  // calculate the transverse attentuation factor
779  std::pair<double,double> CRTDetSimAlg::GetTransAtten(const double pos) {
780  //attenuation coefficients for each fiber
781  double abs0=0.0, abs1=0.0;
782 
783  //coefficiencts for transverse attenuation in CERN modules
784  //from early beta-source tranverse scan data
785  double at0_c = 0.682976;
786  double at1_c = -0.0204477;
787  double at2_c = -0.000707564;
788  double at3_c = 0.000636617;
789  double at4_c = 0.000147957;
790  double at5_c = -3.89078e-05;
791 
792  double at0_r = 0.139941;
793  double at1_r = 0.168238;
794  double at2_r = -0.0198199;
795  double at3_r = 0.000781752;
796 
797  double at0_l = 8.78875;
798  double at1_l = 3.54602;
799  double at2_l = 0.595592;
800  double at3_l = 0.0449169;
801  double at4_l = 0.00127892;
802 
803  //hit between both fibers
804  if ( abs(pos) <= 5.5 ) {
805  abs0 = at5_c*pow(pos,5) + at4_c*pow(pos,4) + at3_c*pow(pos,3)
806  + at2_c*pow(pos,2) + at1_c*pos + at0_c;
807  abs1 = -1*at5_c*pow(pos,5) + at4_c*pow(pos,4) - at3_c*pow(pos,3)
808  + at2_c*pow(pos,2) - at1_c*pos + at0_c;
809  }
810 
811  //hit to right of both fibers
812  if ( pos > 5.5 ) {
813  abs0 = at3_r*pow(pos,3) + at2_r*pow(pos,2) + at1_r*pos + at0_r;
814  abs1 = at4_l*pow(pos,4) - at3_l*pow(pos,3)
815  + at2_l*pow(pos,2) - at1_l*pos + at0_l;
816  }
817 
818  //hit to left of both fibers
819  if ( pos < -5.5 ) {
820  abs0 = at4_l*pow(pos,4) + at3_l*pow(pos,3) \
821  + at2_l*pow(pos,2) + at1_l*pos + at0_l;
822  abs1 = -1*at3_r*pow(pos,3) + at2_r*pow(pos,2) - at1_r*pos + at0_r;
823  }
824 
825  return std::make_pair(abs0,abs1);
826 
827  }
828 
829  //---------------------------------------------------------------
830  double CRTDetSimAlg::GetLongAtten(const double dist) {
831 
832  double pes=0.0;
833 
834  //coefficients for quadratic fit to MINOS test data w/S14
835  //obtained for normally incident cosmic muons
836  //p2_m * x^2 + p1_m * x + p0_m, x is distance from readout [m]
837  double p0_m = 36.5425; //initial light yield (pe) before any attenuation in m scintillator
838  double p1_m = -6.3895;
839  double p2_m = 0.3742;
840 
841  pes = p2_m * pow(dist,2) + p1_m * dist + p0_m;
842 
843  return pes;
844 
845  }
846 
847  //function for simulating time response
848  // takes true hit time, LY(PE) observed, and intitudinal distance from readout
849  // uses 12 FHiCL configurable parameters
850  // returns simulated time in units of clock ticks
851  uint64_t CRTDetSimAlg::GetChannelTriggerTicks(detinfo::ElecClock& clock,
852  double t0, float npeMean, float r)
853  {
854  // Hit timing, with smearing and NPE dependence
855  double tDelayMean =
856  fTDelayNorm *
857  exp(-0.5 * pow((npeMean - fTDelayShift) / fTDelaySigma, 2)) +
858  fTDelayOffset;
859 
860  double tDelayRMS =
861  fTDelayRMSGausNorm *
862  exp(-pow(npeMean - fTDelayRMSGausShift, 2) / fTDelayRMSGausSigma) +
863  fTDelayRMSExpNorm *
864  exp(-(npeMean - fTDelayRMSExpShift) / fTDelayRMSExpScale);
865 
866  double tDelay = CLHEP::RandGauss::shoot(&fRandEngine, tDelayMean, tDelayRMS);
867 
868  // Time resolution of the interpolator
869  tDelay += CLHEP::RandGauss::shoot(&fRandEngine, 0, fTResInterpolator);
870 
871  // Propagation time
872  double tProp = CLHEP::RandGauss::shoot(fPropDelay, fPropDelayError) * r;
873 
874  double t = t0 + tProp + tDelay;
875 
876  // Get clock ticks
877  clock = clock.WithTime(t / 1e3); // WithTime takes microseconds
878 
879  if (fUltraVerbose) mf::LogInfo("CRT")
880  << "CRT TIMING: t0=" << t0
881  << ", tDelayMean=" << tDelayMean << ", tDelayRMS=" << tDelayRMS
882  << ", tDelay=" << tDelay << ", tDelay(interp)="
883  << tDelay << ", tProp=" << tProp << ", t=" << t
884  << ", ticks=" << clock.Ticks() << "\n";
885 
886  return (uint64_t)t;//clock.Ticks();
887  }//CRTDetSim::GetChannelTriggerTicks()
888 
889  //---------------------------------------------------------------
890  // function to clear member data at beginning of each art::event
891  void CRTDetSimAlg::ClearTaggers() {
892 
893  fTaggers.clear();
894  fHasFilledTaggers = false;
895 
896  fNsim_m = 0;
897  fNsim_d = 0;
898  fNsim_c = 0;
899  fNchandat_m = 0;
900  fNchandat_d = 0;
901  fNchandat_c = 0;
902  fNmissthr_c = 0;
903  fNmissthr_d = 0;
904  fNmissthr_m = 0;
905  fNmiss_strcoin_c = 0;
906  fNdual_m = 0;
907 
908  fRegions.clear();
909  fRegCounts.clear();
910  }
911 
912  //----------------------------------------------------------------
913  // function to make fill CRTData products a bit easer
914  CRTData CRTDetSimAlg::FillCRTData(uint8_t mac, uint32_t entry, uint64_t t0, uint64_t t1, uint16_t adc[64]){
915  CRTData dat;
916  dat.fMac5 = mac;
917  dat.fEntry = entry;
918  dat.fTs0 = t0;
919  dat.fTs1 = t1;
920  for(size_t chan=0; chan<64; chan++) {
921  dat.fAdc[chan] = adc[chan];
922  }
923 
924  return dat;
925  }
926 
927  //------------------------------------------------------------------
928  ChanData CRTDetSimAlg::FillChanData(int channel, uint16_t adc, uint64_t ts) {
929  ChanData dat;
930  dat.channel = channel;
931  dat.adc = adc;
932  dat.ts = ts;
933 
934  return dat;
935  }
936 
937  //---------------------------------------------------------------------------
938  void CRTDetSimAlg::FillAdcArr(const vector<ChanData>& data, uint16_t arr[64]) {
939  for(int chan=0; chan<64; chan++)
940  arr[chan] = 0;
941  for(auto const& dat : data) {
942  if(dat.channel<0 || dat.channel>63){
943  std::cout << "ChanData channel out of bounds!" << std::endl;
944  return;
945  }
946  arr[dat.channel] = dat.adc;
947  }
948  return;
949  }
950 
951  }//namespace crt
952 }//namespace icarus
953 
954 #endif
icarus::crt::CRTData::fMac5
uint8_t fMac5
FEB address following numbering convention common for both data and MC.
Definition: ICARUS/CRT/CRTData.hh:28
z
process_name opflash particleana ie ie ie z
Definition: run_opflash_electron.fcl:92
icarus::crt::CRTDetSimAlg::fTDelayRMSGausSigma
double fTDelayRMSGausSigma
Time delay RMS fit: Gaussian width.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:95
icarus::crt::CRTDetSimAlg::fTDelayNorm
double fTDelayNorm
Time delay fit: Gaussian normalization.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:89
second
second
Definition: detectorproperties_icarus.fcl:8
icarus::crt::CRTDetSimAlg::fLayerCoincidenceWindowM
double fLayerCoincidenceWindowM
Time window for two layer coincidence between MINOS modules [ns].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:116
icarus::crt::CRTDetSimAlg::GetTransAtten
pair< double, double > GetTransAtten(const double pos)
Definition: CRTDetSimAlg.cc:779
icarus::crt::CRTDetSimAlg::fLayerCoincidenceWindowC
double fLayerCoincidenceWindowC
Time window for two layer coincidence in a CERN module [ns].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:115
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
icarus::crt::CRTDetSimAlg::fApplyCoincidenceM
bool fApplyCoincidenceM
Whether or not to apply coincidence between hits in adjacent layers.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:113
icarus::crt::CRTDetSimAlg::fDeadTime
double fDeadTime
Dead Time inherent in the front-end electronics.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:119
icarus::crt::CRTDetSimAlg::FillAdcArr
void FillAdcArr(const vector< ChanData > &data, uint16_t arr[32])
Definition: CRTDetSimAlg.cc:938
Tagger::data
vector< pair< ChanData, AuxDetIDE > > data
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:65
icarus::crt::CRTDetSimAlg::fQThresholdC
uint16_t fQThresholdC
ADC charge threshold for CERN system [ADC].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:104
icarus::crt::CRTDetSimAlg::fStripCoincidenceWindow
double fStripCoincidenceWindow
Time window for two-fiber coincidence [ns].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:110
geo::AuxDetGeo::SensitiveVolume
AuxDetSensitiveGeo const & SensitiveVolume(size_t sv) const
Definition: AuxDetGeo.h:171
make_calib_df.dat
tuple dat
Definition: make_calib_df.py:20
icarus::crt::CRTDetSimAlg::fQThresholdD
uint16_t fQThresholdD
ADC charge threshold for DC system [ADC].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:106
icarus::crt::CRTDetSimAlg::fQPed
double fQPed
ADC offset for the single-peak peak mean [ADC].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:100
p
pdgs p
Definition: selectors.fcl:22
geo::AuxDetGeo::WorldToLocal
void WorldToLocal(const double *world, double *auxdet) const
Transform point from world frame to local auxiliary detector frame.
Definition: AuxDetGeo.h:136
icarus::crt::CRTDetSimAlg::fQThresholdM
uint16_t fQThresholdM
ADC charge threshold for MINOS system [ADC].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:105
icarus::crt::CRTData::fEntry
uint32_t fEntry
Hit index for given FEB in an event (starts from 0 for each event).
Definition: ICARUS/CRT/CRTData.hh:29
icarus::crt::CRTDetSimAlg::fQ0
double fQ0
Average energy deposited for mips in 1cm for charge scaling [GeV].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:99
detinfo::ElecClock::Ticks
constexpr int Ticks() const noexcept
Current clock tick (that is, the number of tick Time() falls in).
Definition: ElecClock.h:205
icarus::crt::CRTDetSimAlg::fPropDelayError
double fPropDelayError
Delay in pulse arrival time, uncertainty [ns/m].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:109
icarus::crt::CRTDetSimAlg::FillChanData
ChanData FillChanData(int channel, uint16_t adc, uint64_t ts)
Definition: CRTDetSimAlg.cc:928
icarus::crt::CRTDetSimAlg::fTDelayOffset
double fTDelayOffset
Time delay fit: Gaussian baseline offset.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:92
icarus::crt::CRTDetSimAlg::fTResInterpolator
double fTResInterpolator
Interpolator time resolution [ns].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:107
icarus::crt::CRTDetSimAlg::fTDelayRMSExpScale
double fTDelayRMSExpScale
Time delay RMS fit: Exponential scale.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:98
detinfo::ElecClock::WithTime
constexpr ElecClock WithTime(double const time) const noexcept
Definition: ElecClock.h:108
path
BEGIN_PROLOG triggeremu_data_config_icarus settings PMTADCthresholds sequence::icarus_stage0_multiTPC_TPC physics sequence::icarus_stage0_EastHits_TPC physics sequence::icarus_stage0_WestHits_TPC physics producers purityana0 caloskimCalorimetryCryoE physics caloskimCalorimetryCryoW physics path
Definition: magic_raw_to_triggeremu_icarus_Run1_data.fcl:108
icarus::crt::CRTDetSimAlg::fTDelaySigma
double fTDelaySigma
Time delay fit: Gaussian width.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:91
make_calib_df.t0
int t0
Definition: make_calib_df.py:42
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
icarus::crt::CRTDetSimAlg::fBiasTime
double fBiasTime
Hard cut off for follow-up hits after primary trigger to bias ADC level.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:120
y
process_name opflash particleana ie ie y
Definition: run_opflash_electron.fcl:91
icarus::crt::CRTDetSimAlg::fApplyCoincidenceC
bool fApplyCoincidenceC
Whether or not to apply coincidence between hits in adjacent layers.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:112
icarus::crt::CRTDetSimAlg::fUseBirks
bool fUseBirks
Whether or not to apply Birks&#39; quenching to light output.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:121
icarus::crt::CRTDetSimAlg::CreateData
vector< pair< CRTData, vector< AuxDetIDE > > > CreateData()
Definition: CRTDetSimAlg.cc:77
icarus::crt::CRTDetSimAlg::CRTDetSimAlg
CRTDetSimAlg(fhicl::ParameterSet const &p, CLHEP::HepRandomEngine &fRandEngine)
Definition: CRTDetSimAlg.cc:16
icarus::crt::CRTDetSimAlg::fQRMS
double fQRMS
ADC single-pe spectrum width [ADC].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:102
geo::AuxDetSensitiveGeo::TotalVolume
const TGeoVolume * TotalVolume() const
Definition: AuxDetSensitiveGeo.h:97
icarus::ns::util::end
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
icarus::crt::CRTDetSimAlg::fGlobalT0Offset
double fGlobalT0Offset
Time delay fit: Gaussian normalization.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:88
icarus::crt::CRTData
Single hit (self trigger) of a CRT board.
Definition: ICARUS/CRT/CRTData.hh:26
icarus::crt::CRTDetSimAlg::FillCRTData
CRTData FillCRTData(uint8_t mac, uint32_t entry, uint64_t t0, uint64_t ts1, uint16_t adc[32])
Definition: CRTDetSimAlg.cc:914
icarus::crt::CRTDetSimAlg::fPropDelay
double fPropDelay
Delay in pulse arrival time [ns/m].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:108
detinfo::DetectorClocksData::TriggerClock
ElecClock const & TriggerClock() const noexcept
Borrow a const Trigger clock with time set to Trigger time [us].
Definition: DetectorClocksData.h:414
icarus::crt::CRTDetSimAlg::fApplyStripCoinC
bool fApplyStripCoinC
Whether or not to apply coincence between fibers in a strip (c modules only)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:111
icarus::crt::CRTDetSimAlg::fTDelayRMSGausShift
double fTDelayRMSGausShift
Time delay RMS fit: Gaussian x shift.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:94
icarus::crt::CRTDetSimAlg::fTDelayShift
double fTDelayShift
Time delay fit: Gaussian x shift.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:90
icarus::crt::CRTDetSimAlg::fTDelayRMSGausNorm
double fTDelayRMSGausNorm
Time delay RMS fit: Gaussian normalization.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:93
icarus::crt::TimeOrderCRTData
bool TimeOrderCRTData(std::pair< ChanData, AuxDetIDE > crtdat1, std::pair< ChanData, AuxDetIDE > crtdat2)
Definition: CRTDetSimAlg.cc:9
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
util::absDiff
constexpr auto absDiff(A const &a, B const &b)
Returns the absolute value of the difference between two values.
Definition: NumericUtils.h:43
icarus::crt::CRTDetSimAlg::fTDelayRMSExpShift
double fTDelayRMSExpShift
Time delay RMS fit: Exponential x shift.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:97
icarus::crt::CRTDetSimAlg::fTDelayRMSExpNorm
double fTDelayRMSExpNorm
Time delay RMS fit: Exponential normalization.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:96
icarus::crt::CRTData::fAdc
uint16_t fAdc[64]
ADC readout for each channel. CAEN (Bern) CRT FEBs use only indices 0-31.
Definition: ICARUS/CRT/CRTData.hh:32
icarus::crt::CRTDetSimAlg::FillTaggers
void FillTaggers(detinfo::DetectorClocksData const &clockData, const uint32_t adid, const uint32_t adsid, const vector< AuxDetIDE > &ides)
Definition: CRTDetSimAlg.cc:470
geo::AuxDetGeo::TotalVolume
const TGeoVolume * TotalVolume() const
Definition: AuxDetGeo.h:106
detinfo::DetectorClocksData
Contains all timing reference information for the detector.
Definition: DetectorClocksData.h:283
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:563
icarus::crt::CRTDetSimAlg::reconfigure
void reconfigure(fhicl::ParameterSet const &p)
Definition: CRTDetSimAlg.cc:30
sim::AuxDetIDE
MC truth information to make RawDigits and do back tracking.
Definition: AuxDetSimChannel.h:27
geo::AuxDetSensitiveGeo::WorldToLocal
void WorldToLocal(const double *world, double *auxdet) const
Transform point from world frame to local auxiliary detector frame.
Definition: AuxDetSensitiveGeo.h:126
icarus::crt::CRTDetSimAlg::GetLongAtten
double GetLongAtten(const double dist)
Definition: CRTDetSimAlg.cc:830
icarus::crt::CRTDetSimAlg::fUseEdep
bool fUseEdep
Use the true G4 energy deposited, assume mip if false.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:118
detinfo::ElecClock
Class representing the time measured by an electronics clock.
Definition: ElecClock.h:91
crt
process_name crt
Definition: run_genericCRT.fcl:82
icarus::crt::CRTDetSimAlg::fLayerCoincidenceWindowD
double fLayerCoincidenceWindowD
Time window for two layer coincidence in a DC module [ns].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:117
r
esac echo uname r
Definition: condor_lar_pubs.sh:698
icarus::crt::CRTDetSimAlg::GetChannelTriggerTicks
uint64_t GetChannelTriggerTicks(detinfo::ElecClock &clock, double t0, float npeMean, float r)
Definition: CRTDetSimAlg.cc:851
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
icarus::crt::CRTData::fTs0
uint64_t fTs0
Absolute hit timestamp [ns].
Definition: ICARUS/CRT/CRTData.hh:30
icarus::crt::CRTDetSimAlg::fApplyCoincidenceD
bool fApplyCoincidenceD
Whether or not to apply coincidence between hits in adjacent layers.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTDetSimAlg.h:114
geo::origin
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
icarus::crt::CRTData::fTs1
uint64_t fTs1
Trigger time, not well defined as of Apr 14, 2021.
Definition: ICARUS/CRT/CRTData.hh:31
Generated by   doxygen 1.8.5