OpDetDigitizerDUNE_module.cc
Go to the documentation of this file.
1 //=========================================================
2 // OpDetDigitizerDUNE_module.cc
3 // This module produces digitized output
4 // (creating OpDetWaveform)
5 // from photon detectors taking SimPhotonsLite as input.
6 //
7 // Gleb Sinev, Duke, 2015
8 // Based on OpMCDigi_module.cc
9 //=========================================================
10 
11 #ifndef OpDetDigitizerDUNE_h
12 #define OpDetDigitizerDUNE_h
13 
14 // Framework includes
15 
23 #include "cetlib_except/exception.h"
24 #include "fhiclcpp/ParameterSet.h"
26 #include "art_root_io/TFileService.h"
27 #include "art_root_io/TFileDirectory.h"
28 
29 
30 
31 // ART extensions
32 #include "nurandom/RandomUtils/NuRandomService.h"
33 
34 // LArSoft includes
35 
48 
49 // CLHEP includes
50 
51 #include "CLHEP/Random/RandGauss.h"
52 #include "CLHEP/Random/RandExponential.h"
53 #include "CLHEP/Random/RandFlat.h"
54 
55 // C++ includes
56 
57 #include <vector>
58 #include <map>
59 #include <cmath>
60 #include <memory>
61 
62 // ROOT includes
63 
64 #include "TTree.h"
65 
66 
67 namespace opdet {
68 
69  class FocusList
70  {
71  public:
72  FocusList(int nSamples, int padding)
73  : fNSamples(nSamples), fPadding(padding) {}
74 
75  void AddRange(int from, int to)
76  {
77  from -= fPadding;
78  to += fPadding;
79 
80  if(from < 0) from = 0;
81  if(to >= fNSamples) to = fNSamples-1;
82 
83  for(unsigned int i = 0; i < ranges.size(); ++i){
84  std::pair<int, int>& r = ranges[i];
85  // Completely nested, discard
86  if(from >= r.first && to <= r.second) return;
87  // Extend end
88  if(from >= r.first && from <= r.second){
89  r.second = to;
90  return;
91  }
92  // Extend front
93  if(to >= r.first && to <= r.second){
94  r.first = from;
95  return;
96  }
97  }
98  // Discontiguous, add
99  ranges.emplace_back(from, to);
100  }
101 
102  std::vector<std::pair<int, int>> ranges;
103 
104  protected:
106  int fPadding;
107  };
108 
110 
111  public:
112 
114  // Should the destructor be empty?
115  // virtual ~OpDetDigitizerDUNE();
116 
117  void produce(art::Event&);
118 
119  private:
120 
121  // The parameters read from the FHiCL file
122  std::vector<std::string> vInputModules;
123  std::set<std::string> fInputModules; // Input tag for OpDet collection
124  double fSampleFreq; // Sampling frequency in MHz
125  double fTimeBegin; // Beginning of waveform in us
126  double fTimeEnd; // End of waveform in us
127  double fVoltageToADC; // Conversion factor mV to ADC counts
128  double fLineNoiseRMS; // Pedestal RMS in ADC counts
129  double fDarkNoiseRate; // In Hz
130  double fCrossTalk; // Probability of SiPM producing 2 PE signal
131  // in response to 1 photon
132  short fPedestal; // In ADC counts
133  bool fDefaultSimWindow; // Set the start time to -1 drift window and
134  // the end time to the end time
135  // of the TPC readout
136  bool fFullWaveformOutput; // Output full waveforms -- produces large
137  // output. Mostly for debug purposes
138  size_t fReadoutWindow; // In ticks
139  size_t fPreTrigger; // In ticks
140 
141  int fPadding; // In ticks
142 
143  bool fDigiTree_SSP_LED; // To create a analysis Tree for SSP LED
144  bool fUseSDPs; // = pset.get< bool >("UseSDPs", true);
145 
146  double fQEOverride;
148 
149  //-----------------------------------------------------
150  // Trigger analysis variables
151  std::vector<double> t_photon; // vitor
152  std::vector<int> op_photon;
153 
154  TTree *arvore2;
155  //-----------------------------------------------------
156 
157  // Threshold algorithm
158  std::unique_ptr< pmtana::AlgoSSPLeadingEdge > fThreshAlg;
159 
160  // Random number engines
161  std::unique_ptr< CLHEP::RandGauss > fRandGauss;
162  std::unique_ptr< CLHEP::RandExponential > fRandExponential;
163  std::unique_ptr< CLHEP::RandFlat > fRandFlat;
164 
165  // Function that adds n pulses to a waveform
166  void AddPulse(size_t timeBin, int scale,
167  std::vector< double >& waveform,
168  FocusList& fl) const;
169 
170  // Functional response to one photoelectron (time in ns)
171  double Pulse1PE(double time) const;
172 
173  // Single photoelectron pulse parameters
174  double fPulseLength; // 1PE pulse length in us
175  double fPeakTime; // Time when the pulse reaches its maximum in us
176  double fMaxAmplitude; // Maximum amplitude of the pulse in mV
177  double fFrontTime; // Constant in the exponential function
178  // of the leading edge in us
179  double fBackTime; // Constant in the exponential function
180  // of the tail in us
181 
182  // Make sure the FHiCL parameters make sense
183  void CheckFHiCLParameters() const;
184 
185  std::vector< double > fSinglePEWaveform;
186  void CreateSinglePEWaveform();
187 
188  // Produce waveform on one of the optical detectors
189  void CreatePDWaveform(art::Ptr<sim::OpDetBacktrackerRecord> const& btr_p,
190  geo::Geometry const& geometry,
191  std::vector< std::vector< double > >& pdWaveforms,
192  std::vector<FocusList>& fls,
193  sim::OpDetDivRec& DivRec,
194  bool Reflected);
195 
196  // Vary the pedestal
197  void AddLineNoise(std::vector< std::vector< double > >&,
198  const std::vector<FocusList>& fls) const;
199 
200  void AddDarkNoise(std::vector< std::vector< double > >&,
201  std::vector<FocusList>& fls) const;
202 
203  unsigned short CrossTalk() const;
204 
205  // Create a vector of shorts from a vector of doubles
206  // rounding it properly
207  std::vector< short > VectorOfDoublesToVectorOfShorts(std::vector< double > const&) const;
208 
209  // Make several shorter waveforms out of a long one using a hit finder,
210  // recording also when they start in the long waveform
211  std::map< size_t, std::vector< short > > SplitWaveform(std::vector< short > const&,
212  const FocusList&);
213 
214  double GetDriftWindow(detinfo::DetectorPropertiesData const& detProp) const;
215 
216  // Convert time to ticks or the other way around
217  // without any checks
218  double TickToTime(size_t tick) const;
219  size_t TimeToTick(double time) const;
220 
221 
222  bool Detected(int opDet, int &readoutChannel, bool Reflected);
223  };
224 
225 }
226 
227 #endif
228 
229 namespace opdet {
230 
232 
233 }
234 
235 namespace opdet {
236 
237  //---------------------------------------------------------------------------
238  // Constructor
240  : EDProducer{pset}
241  , vInputModules(pset.get< std::vector<std::string> >("InputModules"))
242  , fInputModules(vInputModules.begin(), vInputModules.end())
243  {
244 
245 
246  // Read the fcl-file
247 // auto tempvec =
248 // =
249  fVoltageToADC = pset.get< double >("VoltageToADC" );
250  fLineNoiseRMS = pset.get< double >("LineNoiseRMS" );
251  fDarkNoiseRate = pset.get< double >("DarkNoiseRate" );
252  fCrossTalk = pset.get< double >("CrossTalk" );
253  fPedestal = pset.get< short >("Pedestal" );
254  fDefaultSimWindow = pset.get< bool >("DefaultSimWindow" );
255  fFullWaveformOutput = pset.get< bool >("FullWaveformOutput");
256  fReadoutWindow = pset.get< size_t >("ReadoutWindow" );
257  fPreTrigger = pset.get< size_t >("PreTrigger" );
258 
259  fPadding = pset.get< int >("Padding" );
260 
261  fPulseLength = pset.get< double >("PulseLength" );
262  fPeakTime = pset.get< double >("PeakTime" );
263  fMaxAmplitude = pset.get< double >("MaxAmplitude" );
264  fFrontTime = pset.get< double >("FrontTime" );
265  fBackTime = pset.get< double >("BackTime" );
266  fDigiTree_SSP_LED = pset.get< bool >("SSP_LED_DigiTree" );
267  fUseSDPs = pset.get< bool >("UseSDPs", true );
268 
269  if (!fUseSDPs) {
270  throw art::Exception(art::errors::UnimplementedFeature) << "SimPhotonsLite is now deprecated in favor SDPs. If you do not have SDPs because your input file is old, use an older version of dunetpc to run this digitizer";
271  }
272 
273  // Replace QE with effective area. Get area from OpDet geometry
274  double tempQE = pset.get< double >("QEOverride", 0 );
275  double tempRefQE = pset.get< double >("QERefOverride", 0 );
276 
277  fQEOverride = 0;
278  fRefQEOverride = 0;
279 
280  if (tempQE > 0) {
281  mf::LogInfo("OpDetDigitizerDUNE") << "Overriding QE. All the functions of OpDetResponseService being ignored.";
282 
283  // Correct out the prescaling applied during simulation
284  auto const *LarProp = lar::providerFrom<detinfo::LArPropertiesService>();
285  fQEOverride = tempQE / LarProp->ScintPreScale();
286 
287  if (fQEOverride > 1.0001 ) {
288  mf::LogError("OpDetDigitizerDUNE") << "Quantum efficiency set as OpDetDigitizerDUNE.QEOverride, " << tempQE
289  << " is too large. It is larger than the prescaling applied during simulation, "
290  << LarProp->ScintPreScale()
291  << " (fQEOverride = "
292  << fQEOverride
293  << "). Final QE must be equal to or smaller than the QE applied at simulation time.";
294  std::abort();
295  }
296  }
297  if (tempRefQE > 0) {
298  mf::LogInfo("OpDetDigitizerDUNE") << "Overriding QE. All the functions of OpDetResponseService being ignored.";
299 
300  // Correct out the prescaling applied during simulation
301  auto const *LarProp = lar::providerFrom<detinfo::LArPropertiesService>();
302  fRefQEOverride = tempRefQE / LarProp->ScintPreScale();
303 
304  if (fRefQEOverride > 1.0001 ) {
305  mf::LogError("OpDetDigitizerDUNE") << "Quantum efficiency set as OpDetDigitizerDUNE.QERefOverride, " << tempRefQE
306  << " is too large. It is larger than the prescaling applied during simulation, "
307  << LarProp->ScintPreScale()
308  << " (fRefQEOverride = "
309  << fRefQEOverride
310  << "). Final QE must be equal to or smaller than the QE applied at simulation time.";
311  std::abort();
312  }
313  }
314 
315  fThreshAlg = std::make_unique< pmtana::AlgoSSPLeadingEdge >
316  (pset.get< fhicl::ParameterSet >("algo_threshold"));
317 
318  if(fDigiTree_SSP_LED){
320  arvore2 = tfs->make<TTree>("PhotonData", "Photon_analysis");
321  arvore2->Branch("photon_opCh",&op_photon);
322  arvore2->Branch("photon_pulse",&t_photon);
323  }
324 
325  // This module produces (infrastructure piece)
326  produces< std::vector< raw::OpDetWaveform > >();
327  produces<std::vector<sim::OpDetDivRec> > ();
328 
329 
330  // Obtaining parameters from the DetectorClocksService
331  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
332  fSampleFreq = clockData.OpticalClock().Frequency();
333 
334  if (fDefaultSimWindow)
335  {
336  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob(clockData);
337 
338  // Assume the readout starts at -1 drift window
339  fTimeBegin = -1*GetDriftWindow(detProp);
340 
341  // Take the TPC readout window size and convert
342  // to us with the electronics clock frequency
343  fTimeEnd = detProp.ReadOutWindowSize() / clockData.TPCClock().Frequency();
344  }
345  else
346  {
347  fTimeBegin = pset.get< double >("TimeBegin");
348  fTimeEnd = pset.get< double >("TimeEnd" );
349  }
350 
352 
353  // Initializing random number engines
354  unsigned int seed = pset.get< unsigned int >("Seed", sim::GetRandomNumberSeed());
355  auto& engine = createEngine(seed);
356  fRandGauss = std::make_unique< CLHEP::RandGauss >(engine);
357  fRandExponential = std::make_unique< CLHEP::RandExponential >(engine);
358  fRandFlat = std::make_unique< CLHEP::RandFlat >(engine);
359 
360  // Creating a single photoelectron waveform
361  // Hardcoded, probably need to read them from the FHiCL file
362  //fPulseLength = 4.0;
363  //fPeakTime = 0.260;
364  //fMaxAmplitude = 0.12;
365  //fFrontTime = 0.009;
366  //fBackTime = 0.476;
368 
369  }
370 
371  //---------------------------------------------------------------------------
373  {
374 
376 
377  // A pointer that will store produced OpDetWaveforms
378  // std::unique_ptr< std::vector< raw::OpDetWaveform > >
379  // std::make_unique< std::vector< raw::OpDetWaveform > >());
380  auto wave_forms_p = std::make_unique< std::vector< raw::OpDetWaveform > >();
381  auto bt_DivRec_p = std::make_unique< std::vector< sim::OpDetDivRec > >();
382 
383 
384  // Total number of ticks in the full waveform
385  // Including one pretrigger window before the waveform
386  // and one readout window - pretrigger window after the waveform
387  unsigned int nSamples = (fTimeEnd - fTimeBegin)*fSampleFreq
388  + fReadoutWindow;
389 
390  // Geometry service
392 
393  auto const btr_handles = evt.getMany<std::vector<sim::OpDetBacktrackerRecord>>();
394 
395  if (btr_handles.size() == 0)
396  throw art::Exception(art::errors::ProductNotFound)<<"No OpDetBacktrackerRecords retrieved.";
397 
398  for (auto btr_handle: btr_handles) {
399  // Do some checking before we proceed
400  if (!btr_handle.isValid()) continue;
401  if (!fInputModules.count(btr_handle.provenance()->moduleLabel())) continue;
402 
403  bool Reflected = (btr_handle.provenance()->productInstanceName() == "Reflected");
404 
405 
406  std::vector<art::Ptr<sim::OpDetBacktrackerRecord>> btr_vec;
407  art::fill_ptr_vector(btr_vec, btr_handle);
408 
409 
410  // For every optical detector:
411  for (auto const& btr : btr_vec)
412  {
413  int opDet = btr->OpDetNum();
414  //unsigned int opDet = btr->OpDetNum();
415 
416  // Get number of channels in this optical detector
417  unsigned int nChannelsPerOpDet = geometry->NOpHardwareChannels(opDet);
418 
419  std::vector<FocusList> fls(nChannelsPerOpDet, FocusList(nSamples, fPadding));
420  //std::vector<sim::OpDetDivRec> DivRec;
421  sim::OpDetDivRec DivRec(opDet);
422  //DivRec.chans.resize(nChannelsPerOpDet);
423 
424  // This vector stores waveforms created for each optical channel
425  std::vector< std::vector< double > > pdWaveforms(nChannelsPerOpDet,
426  std::vector< double >(nSamples, static_cast< double >(fPedestal)));
427 
428  CreatePDWaveform(btr, *geometry, pdWaveforms, fls, DivRec, Reflected);
429  //DivRec comes out with all of the ticks filled correctly, with each channel filled in it's map.
430  //Break here to investigate div recs as they are made and compare them to btrs
431 
432 
433  // Generate dark noise //I will not at this time include dark noise in my split backtracking records.
434  if (fDarkNoiseRate > 0.0) AddDarkNoise(pdWaveforms, fls);
435 
436  // Uncomment to undo the effect of FocusLists. Replaces the accumulated
437  // lists with ones asserting we need to look at the whole trace.
438  // for(FocusList& fl: fls){
439  // fl.ranges.clear();
440  // fl.ranges.emplace_back(0, nSamples-1);
441  // }
442 
443  // Vary the pedestal
444  if (fLineNoiseRMS > 0.0) AddLineNoise(pdWaveforms, fls);
445 
446  // Loop over all the created waveforms, split them into shorter
447  // waveforms and use them to initialize OpDetWaveforms
448  for (unsigned int hardwareChannel = 0;
449  hardwareChannel < nChannelsPerOpDet; ++hardwareChannel)
450  {
451  for(const std::pair<int, int>& p: fls[hardwareChannel].ranges){
452  // It's a shame we copy here. We could actually avoid by making the
453  // functions below take a begin()/end() pair.
454  const std::vector<double> sub(pdWaveforms[hardwareChannel].begin()+p.first,
455  pdWaveforms[hardwareChannel].begin()+p.second+1);
456 
457  std::vector< short > waveformOfShorts = VectorOfDoublesToVectorOfShorts(sub);
458 
459  std::map< size_t, std::vector < short > > mapTickWaveform =
461  SplitWaveform(waveformOfShorts, fls[hardwareChannel]) :
462  std::map< size_t, std::vector< short > >{ std::make_pair(0,
463  waveformOfShorts) };
464 
465  unsigned int opChannel = geometry->OpChannel(opDet, hardwareChannel);
466 
467  for (auto const& pairTickWaveform : mapTickWaveform)
468  {
469  double timeStamp =
470  static_cast< double >(TickToTime(pairTickWaveform.first+p.first));
471 
472  raw::OpDetWaveform adcVec(timeStamp, opChannel,
473  pairTickWaveform.second.size());
474 
475  for (short const& value : pairTickWaveform.second){
476  adcVec.emplace_back(value);
477  }
478 
479  // wave_forms_p->emplace_back(std::move(adcVec));
480  wave_forms_p->push_back(std::move(adcVec));
481  //art::Ptr< raw::OpDetWaveform > wave_form_p = make_wave_form_ptr(wave_forms_p->size()-1);
482  //bt_wave_assns->addSingle(btr, wave_form_p, DivRec/*DIVREC*/);
483  }
484  }
485  }
486  bt_DivRec_p->push_back(std::move(DivRec));
487  }
488  }
489 
490  if(fDigiTree_SSP_LED){
491  arvore2->Fill();
492  t_photon.clear();
493  op_photon.clear();
494  }
495 
496 
497  // Push the OpDetWaveforms into the event
498  evt.put(std::move(wave_forms_p));
499  if(fUseSDPs){ evt.put(std::move(bt_DivRec_p));}
500 
501 
502  }
503 
504  //---------------------------------------------------------------------------
505  void OpDetDigitizerDUNE::AddPulse(size_t timeBin,
506  int scale, std::vector< double >& waveform,
507  FocusList& fl) const
508  {
509 
510  // How many bins will be changed
511  size_t pulseLength = fSinglePEWaveform.size();
512  if ((timeBin + fSinglePEWaveform.size()) > waveform.size())
513  pulseLength = (waveform.size() - timeBin);
514 
515  fl.AddRange(timeBin, timeBin+pulseLength-1);
516 
517  // Adding a pulse to the waveform
518  for (size_t tick = 0; tick != pulseLength; ++tick)
519  waveform[timeBin + tick] += scale*fSinglePEWaveform[tick];
520 
521  }
522 
523  //---------------------------------------------------------------------------
524  double OpDetDigitizerDUNE::Pulse1PE(double time) const
525  {
526 
527  if (time < fPeakTime) return
528  (fVoltageToADC*fMaxAmplitude*std::exp((time - fPeakTime)/fFrontTime));
529  else return
530  (fVoltageToADC*fMaxAmplitude*std::exp(-(time - fPeakTime)/fBackTime));
531 
532  }
533 
534  //---------------------------------------------------------------------------
536  {
537 
538  size_t length =
539  static_cast< size_t > (std::round(fPulseLength*fSampleFreq));
540  fSinglePEWaveform.resize(length);
541  for (size_t tick = 0; tick != length; ++tick)
543  Pulse1PE(static_cast< double >(tick)/fSampleFreq);
544 
545  }
546 
547  //---------------------------------------------------------------------------
550  geo::Geometry const& geometry,
551  std::vector< std::vector< double > >& pdWaveforms,
552  std::vector<FocusList>& fls,
553  sim::OpDetDivRec& DivRec,
554  bool Reflected)
555  {
556 
557  int const opDet = btr_p->OpDetNum();
558  //unsigned int const opDet = btr_p->OpDetNum();
559  // This is int because otherwise detectedLite doesn't work
560  int readoutChannel;
561  // For a group of photons arriving at the same time this is a map
562  // of < arrival time (in ns), number of photons >
563  // std::map< int, int > const& photonsMap = litePhotons.DetectedPhotons;
564  //sim::timePDclockSDPs_t time_sdps_vector = btr_p->timePDclockSDPsMap();
565  //int time_sdps_vector = btr_p->timePDclockSDPsMap();
566  auto time_sdps_vector = btr_p->timePDclockSDPsMap();
567  /*
568  for(auto& divchan : DivRec.chans){
569  if(divchan.tick_photons_frac.size()<time_sdps_vector.size())
570  divchan.tick_photons_frac.resize(time_sdps_vector.size(), 0.0);
571  }*/
572 
573  // For every pair of (arrival time, number of photons) in the map:
574  //for (auto const& pulse : photonsMap)
575  //for (auto const& time_sdps : time_sdps_vector)
576  for (size_t i = 0; i<time_sdps_vector.size(); ++i) //I is now the time bin.
577  {
578  auto time_sdps = time_sdps_vector[i];
579  // Converting ns to us
580  // double photonTime = static_cast< double >(time_sdps.first)/1000.0;
581  //Really need to do something better here. This conversion is fragile, and makes the populating of the DivRecs confusing.
582  double photonTime = time_sdps.first/1000.0;//This should be done with the timing service
583  //for (int i = 0; i < pulse.second; ++i)
584  //for (size_t i = 0; i < time_sdps.second.size(); ++i)
585  for (auto const& sdp : time_sdps.second)
586  {
587  int tid = sdp.trackID;
588  for(int j=0; j<sdp.numPhotons;++j)
589  {
590  if ((photonTime >= fTimeBegin) && (photonTime < fTimeEnd))
591  {
592  // Sample a random subset according to QE
593  if ( Detected(opDet, readoutChannel, Reflected) )
594  {
595  unsigned int hardwareChannel =
596  geometry.HardwareChannelFromOpChannel(readoutChannel);
597  // Convert the time of the pulse to ticks
598  size_t timeBin = TimeToTick(photonTime);
599  // Add 1 pulse to the waveform
600  AddPulse(timeBin, CrossTalk(), pdWaveforms.at(hardwareChannel), fls[hardwareChannel]);
601 
602  unsigned int opChannel = geometry.OpChannel(opDet, hardwareChannel);
603  //Set/find tick. Set/find Channel
604  sim::OpDet_Time_Chans::stored_time_t tmp_time=time_sdps.first;
605  DivRec.AddPhoton(opChannel, tid, tmp_time);
606  if(fDigiTree_SSP_LED){
607  op_photon.emplace_back(opChannel);
608  t_photon.emplace_back(photonTime); //vitor: devo usar o time ou o tick?
609  }
610  }//else{
611  /* unsigned int hardwareChannel =
612  geometry.HardwareChannelFromOpChannel(readoutChannel);
613  unsigned int opChannel = geometry.OpChannel(opDet, hardwareChannel);
614  DivRec.tick_chans[i].DivChans.IfNotInit(opChannel);
615  }*/
616  }
617  // else
618  // remove this as it fills up logfiles for cosmic-ray runs
619  //mf::LogInfo("OpDetDigitizerDUNE")
620  // << "Throwing away an out-of-time photon at " << photonTime << '\n';
621  }
622  }//end of this sdp.
623  /*double total=0.0; //This section is no longer needed. I store the un-normalized values, and normalize them at retrieval.
624  for(auto it=DivRec.tick_chans[i].DivChans.bit(); it!= DivRec.tick_chans[i].DivChans.lit(); ++it)
625  { //With a bit of work, we could avoid this loop by remembering this number from above.
626  total+=it->second.tick_photons_frac;
627  }
628  for(auto chan : DivRec.tick_chans[i].DivChans.chans())
629  {
630  DivRec.tick_chans[i].DivChans.NormPhotonFrac(chan, total);
631  }
632  total=0.0; // This isn't actually needed.*/
633  }//End this time tick.
634 
635  }
636 
637  //---------------------------------------------------------------------------
639  AddLineNoise(std::vector< std::vector< double > >& waveforms,
640  const std::vector<FocusList>& fls) const
641  {
642  int i = 0;
643  for(auto& waveform : waveforms){
644  for(unsigned int j = 0; j < fls[i].ranges.size(); ++j){
645  const std::pair<int, int>& p = fls[i].ranges[j];
646  for(int k = p.first; k <= p.second; ++k){
647  waveform[k] += fRandGauss->fire(0, fLineNoiseRMS);
648  }
649  }
650 
651  ++i;
652  }
653  }
654 
655  //---------------------------------------------------------------------------
657  AddDarkNoise(std::vector< std::vector< double > >& waveforms,
658  std::vector<FocusList>& fls) const
659  {
660  int i = 0;
661  for (auto& waveform : waveforms)
662  {
663  // Multiply by 10^6 since fDarkNoiseRate is in Hz
664  double darkNoiseTime = static_cast< double >(fRandExponential->
665  fire(1.0/fDarkNoiseRate)*1000000.0) + fTimeBegin;
666  while (darkNoiseTime < fTimeEnd)
667  {
668  size_t timeBin = TimeToTick(darkNoiseTime);
669  AddPulse(timeBin, CrossTalk(), waveform, fls[i]);
670  // Find next time to simulate a single PE pulse
671  darkNoiseTime += static_cast< double >
672  (fRandExponential->fire(1.0/fDarkNoiseRate)*1000000.0);
673  }
674 
675  ++i;
676  }
677  }
678 
679  //---------------------------------------------------------------------------
680  unsigned short OpDetDigitizerDUNE::CrossTalk() const
681  {
682  // Sometimes this should produce 3 or more PEs (not implemented)
683  if (fCrossTalk <= 0.0) return 1;
684  else if (fRandFlat->fire(1.0) > fCrossTalk) return 1;
685  else return 2;
686  }
687 
688  //---------------------------------------------------------------------------
690  (std::vector< double > const& vectorOfDoubles) const
691  {
692  // Don't bother to round properly, it's faster this way
693  return std::vector<short>(vectorOfDoubles.begin(), vectorOfDoubles.end());
694 
695  /*
696  std::vector< short > vectorOfShorts;
697  vectorOfShorts.reserve(vectorOfDoubles.size());
698 
699  for (short const& value : vectorOfDoubles)
700  vectorOfShorts.emplace_back(static_cast< short >(std::round(value)));
701 
702  return vectorOfShorts;
703  */
704  }
705 
706  //---------------------------------------------------------------------------
707  std::map< size_t, std::vector< short > > OpDetDigitizerDUNE::
708  SplitWaveform(std::vector< short > const& waveform,
709  const FocusList& fls)
710  {
711 
712  std::map< size_t, std::vector< short > > mapTickWaveform;
713 
714  ::pmtana::PedestalMean_t ped_mean (waveform.size(),0);
715  ::pmtana::PedestalSigma_t ped_sigma(waveform.size(),0);
716 
717  fThreshAlg->Reconstruct(waveform,ped_mean,ped_sigma);
718 
719  std::vector< pmtana::pulse_param > pulses;
720  for (size_t pulseCounter = 0; pulseCounter < fThreshAlg->GetNPulse();
721  ++pulseCounter)
722  pulses.emplace_back(fThreshAlg->GetPulse(pulseCounter));
723 
724  // We have to refine this algorithm later
725  for (pmtana::pulse_param const& pulse : pulses)
726  {
727  if (pulse.t_end <= pulse.t_start)
728  // Can I call it a logic error?
730  << "Pulse ends before or at the same time it starts!\n";
731 
733  waveform.begin() + static_cast< size_t >(pulse.t_start);
735  waveform.begin() + static_cast< size_t >(pulse.t_end );
736  mapTickWaveform.emplace(static_cast< size_t >(pulse.t_start),
737  std::vector< short >(window_start, window_end));
738  // Don't forget to check that the time output by the (new) algortihm
739  // is < fTimeEnd and > fTimeBegin!
740 
741  }
742 
743  return mapTickWaveform;
744 
745  }
746 
747  //---------------------------------------------------------------------------
749  {
750 
751  double driftWindow;
752 
753  double maxDrift = 0.0;
754  for (geo::TPCGeo const& tpc :
755  art::ServiceHandle< geo::Geometry >()->IterateTPCs())
756  if (maxDrift < tpc.DriftDistance()) maxDrift = tpc.DriftDistance();
757 
758  driftWindow = maxDrift/detProp.DriftVelocity();
759 
760  return driftWindow;
761 
762  }
763 
764  //---------------------------------------------------------------------------
766  {
767 
768  if (tick > fPreTrigger)
769  return (static_cast< double >(tick - fPreTrigger)/fSampleFreq
770  + fTimeBegin);
771  else
772  return (static_cast< double >(fPreTrigger - tick)/fSampleFreq*(-1.0)
773  + fTimeBegin);
774 
775  }
776 
777  //---------------------------------------------------------------------------
779  {
780 
781  return static_cast< size_t >(std::round((time - fTimeBegin)*fSampleFreq
782  + fPreTrigger));
783 
784  }
785 
786  //---------------------------------------------------------------------------
788  {
789 
790  // Are all these logic errors?
791 
792  if (fLineNoiseRMS < 0.0)
794  << "fLineNoiseRMS: " << fLineNoiseRMS << '\n'
795  << "Line noise RMS should be non-negative!\n";
796 
797  if (fDarkNoiseRate < 0.0)
799  << "fDarkNoiseRate: " << fDarkNoiseRate << '\n'
800  << "Dark noise rate should be non-negative!\n";
801 
804  << "PreTrigger: " << fPreTrigger << " and "
805  << "ReadoutWindow: " << fReadoutWindow << '\n'
806  << "Pretrigger window has to be shorter than readout window!\n";
807 
808  if (fTimeBegin >= fTimeEnd)
810  << "TimeBegin: " << fTimeBegin << " and "
811  << "TimeEnd: " << fTimeEnd << '\n'
812  << "TimeBegin should be less than TimeEnd!\n";
813 
814  }
815 
816  //---------------------------------------------------------------------------
817  bool OpDetDigitizerDUNE::Detected(int OpDet, int &readoutChannel, bool Reflected) {
818 
819  if (Reflected && !fRefQEOverride) {
820  cet::exception("OpDetDigitizerDUNE") << "Cannot handle reflected light if QERefOverride isn't set";
821  }
822  if ( (!Reflected && fQEOverride > 0) || Reflected ) {
823  // Find the Optical Detector using the geometry service
825 
826  // Here OpDet must be opdet since we are introducing
827  // channel mapping here.
828  float NOpHardwareChannels = geom->NOpHardwareChannels(OpDet);
829  int hardwareChannel = (int) ( CLHEP::RandFlat::shoot(1.0) * NOpHardwareChannels );
830  readoutChannel = geom->OpChannel(OpDet, hardwareChannel);
831 
832  // Check QE
833  if (!Reflected) {
834  if ( CLHEP::RandFlat::shoot(1.0) > fQEOverride ) return false;
835  }
836  else {
837  if ( CLHEP::RandFlat::shoot(1.0) > fRefQEOverride ) return false;
838  }
839  return true;
840  }
841  else {
843  // Service for determining optical detector responses
844  return odResponse->detectedLite(OpDet, readoutChannel);
845  }
846  }
847 
848 
849 }
base_engine_t & createEngine(seed_t seed)
Store parameters for running LArG4.
bool Detected(int opDet, int &readoutChannel, bool Reflected)
void AddLineNoise(std::vector< std::vector< double > > &, const std::vector< FocusList > &fls) const
std::vector< double > PedestalSigma_t
FocusList(int nSamples, int padding)
void AddRange(int from, int to)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
struct vector vector
unsigned int NOpHardwareChannels(int opDet) const
std::map< size_t, std::vector< short > > SplitWaveform(std::vector< short > const &, const FocusList &)
void AddPulse(size_t timeBin, int scale, std::vector< double > &waveform, FocusList &fl) const
size_t TimeToTick(double time) const
intermediate_table::const_iterator const_iterator
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
double GetDriftWindow(detinfo::DetectorPropertiesData const &detProp) const
std::vector< std::pair< int, int > > ranges
OpDetDigitizerDUNE(fhicl::ParameterSet const &)
std::vector< double > fSinglePEWaveform
Simulation objects for optical detectors.
int OpDetNum() const
Returns the readout Optical Detector this object describes.
std::unique_ptr< CLHEP::RandExponential > fRandExponential
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
double Pulse1PE(double time) const
unsigned int HardwareChannelFromOpChannel(int opChannel) const
Convert unique channel to hardware channel.
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const
Definition: DataViewImpl.h:479
def move(depos, offset)
Definition: depos.py:107
std::vector< std::string > vInputModules
std::unique_ptr< pmtana::AlgoSSPLeadingEdge > fThreshAlg
void AddDarkNoise(std::vector< std::vector< double > > &, std::vector< FocusList > &fls) const
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:75
p
Definition: test.py:223
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
void CreatePDWaveform(art::Ptr< sim::OpDetBacktrackerRecord > const &btr_p, geo::Geometry const &geometry, std::vector< std::vector< double > > &pdWaveforms, std::vector< FocusList > &fls, sim::OpDetDivRec &DivRec, bool Reflected)
std::unique_ptr< CLHEP::RandGauss > fRandGauss
std::unique_ptr< CLHEP::RandFlat > fRandFlat
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
The geometry of one entire detector, as served by art.
Definition: Geometry.h:196
virtual bool detectedLite(int OpChannel, int &newOpChannel) const
unsigned int OpChannel(int detNum, int hardwareChannel) const
Convert detector number and hardware channel to unique channel.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
unsigned int GetRandomNumberSeed()
Definition: sim.h:32
double TickToTime(size_t tick) const
Index NOpHardwareChannels(Index opDet)
std::vector< short > VectorOfDoublesToVectorOfShorts(std::vector< double > const &) const
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:72
Tools and modules for checking out the basics of the Monte Carlo.
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
std::vector< double > PedestalMean_t
timePDclockSDPs_t const & timePDclockSDPsMap() const
Returns all the deposited energy information as stored.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void AddPhoton(int opchan, int tid, OpDet_Time_Chans::stored_time_t pdTime)
Definition: OpDetDivRec.cxx:25
std::set< std::string > fInputModules
std::string sub(const std::string &a, const std::string &b)
Definition: TruthText.cxx:100