OpDetDigitizerDUNEDP_module.cc
Go to the documentation of this file.
1 //=========================================================
2 // OpDetDigitizerDUNEDP_module.cc
3 // This module produces digitized output
4 // (creating OpDetWaveform)
5 // from photon detectors taking SimPhotonsLite as input.
6 //
7 // J. Soto, CIEMAT, 2018
8 // Based on OpMCDigi_module.cc
9 //=========================================================
10 
11 #ifndef OpDetDigitizerDUNEDP_h
12 #define OpDetDigitizerDUNEDP_h 1
13 
14 // Framework includes
15 
22 #include "fhiclcpp/ParameterSet.h"
23 #include "art_root_io/TFileService.h"
24 #include "art_root_io/TFileDirectory.h"
25 #include "CLHEP/Random/RandFlat.h"
26 
27 
28 // ART extensions
29 #include "nurandom/RandomUtils/NuRandomService.h"
30 
31 // LArSoft includes
32 
43 
44 // CLHEP includes
45 
46 #include "CLHEP/Random/RandGauss.h"
47 #include "CLHEP/Random/RandExponential.h"
48 #include "CLHEP/Random/RandFlat.h"
49 
50 // C++ includes
51 
52 #include <vector>
53 #include <map>
54 #include <cmath>
55 #include <memory>
56 
57 // ROOT includes
58 
59 #include "TTree.h"
60 
61 
62 namespace opdet {
63 
64  class FocusList
65  {
66  public:
67  FocusList(int nSamples, int padding)
68  : fNSamples(nSamples), fPadding(padding) {}
69 
70  //Warning, this FocusList feature is actually not used to split the waveform,
71  //the use a standard hit finding algorithm instead.
72  void AddRange(int from, int to)
73  {
74  from -= fPadding;
75  to += fPadding;
76  if(from < 0) from = 0;
77  if(to >= fNSamples) to = fNSamples-1;
78 
79  for(unsigned int i = 0; i < ranges.size(); ++i){
80  std::pair<int, int>& r = ranges[i];
81  // Completely nested, discard
82  if(from >= r.first && to <= r.second) return;
83  // Extend end
84  if(from >= r.first && from <= r.second){
85  r.second = to;
86  return;
87  }
88  // Extend front
89  if(to >= r.first && to <= r.second){
90  r.first = from;
91  return;
92  }
93  }
94  // Discontiguous, add
95  ranges.emplace_back(from, to);
96  }
97 
98  std::vector<std::pair<int, int>> ranges;
99 
100  protected:
101  int fNSamples;
102  int fPadding;
103  };
104 
106 
107  public:
108 
110  // Should the destructor be empty?
111 // virtual ~OpDetDigitizerDUNEDP();
112 
113  void produce(art::Event&);
114 
115  private:
116 
117  // The parameters read from the FHiCL file
118  std::vector < std::string > fInputModule; // Input tag for OpDet collection
119  double fSampleFreq; // Sampling frequency in MHz
120  double fTimeBegin; // Beginning of waveform in us
121  double fTimeEnd; // End of waveform in us
122  double fVoltageToADC; // Conversion factor mV to ADC counts
123  double fLineNoiseRMS; // Pedestal RMS in ADC counts
124  double fDarkNoiseRate; // In Hz
125  double fCrossTalk; // Probability of SiPM producing 2 PE signal
126  // in response to 1 photon
127  double fPedestal; // In ADC counts
128  double fGain; // Gain of the PMT
129  bool fDefaultSimWindow; // Set the start time to -1 drift window and
130  // the end time to the end time
131  // of the TPC readout
132  bool fFullWaveformOutput; // Output full waveforms -- produces large
133  // output. Mostly for debug purposes
134  size_t fReadoutWindow; // In ticks
135  size_t fPreTrigger; // In ticks
136 
137  int fPadding; // In ticks
138 
139  bool fDigiTree_SSP_LED; // To create a analysis Tree for SSP LED
140 
141 //-----------------------------------------------------
142  // Trigger analysis variables
143  std::vector<double> t_photon; // vitor
144  std::vector<int> op_photon;
145 
146  TTree *arvore2;
147 //-----------------------------------------------------
148 
149  // Threshold algorithm
150  std::unique_ptr< pmtana::AlgoSSPLeadingEdge > fThreshAlg;
151 
152  // Random number engines
153  std::unique_ptr< CLHEP::RandGauss > fRandGauss;
154  std::unique_ptr< CLHEP::RandExponential > fRandExponential;
155  std::unique_ptr< CLHEP::RandFlat > fRandFlat;
156 
157  // Function that adds n pulses to a waveform
158  void AddPulse(size_t timeBin, int scale,
159  std::vector< double >& waveform,
160  FocusList& fl, double Gain) const;
161  double fQE;
162  bool fCustomQE;
163  std::vector<double> fQEVector;
165  std::vector<double> fGainVector;
166  double getQE(int OpDet);
167  double getGain(int OpDet);
168  double fGainsError;
169  double fSPEWidth;
170  // Make sure the FHiCL parameters make sense
171  void CheckFHiCLParameters() const;
172 
173  std::vector< double > fSinglePEWaveform;
174  void CreateSinglePEWaveform();
175  bool fNegativeSignal; // negative signal if true (as real pmts)
176 
177  // Produce waveform on one of the optical detectors
178  void CreatePDWaveform(sim::SimPhotonsLite const&,
180  geo::Geometry const&,
181  std::vector< std::vector< double > >&,
182  std::vector<FocusList>&);
183 
184  // Produce waveform on one of the optical detectors
185  void CreatePDWaveform(sim::SimPhotons const&,
187  geo::Geometry const&,
188  std::vector< std::vector< double > >&,
189  std::vector<FocusList>&);
190 
191  // Vary the pedestal
192  void AddLineNoise(std::vector< std::vector< double > >&,
193  const std::vector<FocusList>& fls) const;
194 
195  void AddDarkNoise(std::vector< std::vector< double > >&,
196  std::vector<FocusList>& fls, int opDet);
197 
198  unsigned short CrossTalk() const;
199 
200  // Create a vector of shorts from a vector of doubles
201  // rounding it properly
202  std::vector< short > VectorOfDoublesToVectorOfShorts
203  (std::vector< double > const&) const;
204 
205  // Make several shorter waveforms out of a long one using a hit finder,
206  // recording also when they start in the long waveform
207  std::map< size_t, std::vector< short > >
208  SplitWaveform(std::vector< short > const&,
209  const FocusList&);
210 
211  double GetDriftWindow(detinfo::DetectorPropertiesData const& detProp) const;
212 
213  // Convert time to ticks or the other way around
214  // without any checks
215  double TickToTime(size_t tick) const;
216  size_t TimeToTick(double time) const;
217 
218  int PMTSaturationFunction(int);
219  double SumOfElements(std::vector<double>);
222  TTree *fRunInfo;
223 
224  int Run;
225  int SubRun;
226  int Event;
227  unsigned int nSamples;
228  double SampleSize;
229  std::vector<short*> adc_value;
230  unsigned int nOpDet;
231  };
232 
233 }
234 
235 #endif
236 
237 namespace opdet {
238 
240 
241 }
242 
243 namespace opdet {
244 
245  //---------------------------------------------------------------------------
246  // Constructor
248  : EDProducer{pset}
249  {
250 
251  // This module produces (infrastructure piece)
252  produces< std::vector< raw::OpDetWaveform > >();
253 
254  // Read the fcl-file
255  fInputModule = pset.get<std::vector<std::string>>("InputModule",{"largeant"});
256  fVoltageToADC = pset.get< double >("VoltageToADC" );
257  fLineNoiseRMS = pset.get< double >("LineNoiseRMS" );
258  fCrossTalk = pset.get< double >("CrossTalk" );
259  fPedestal = pset.get< double >("Pedestal" );
260  fDefaultSimWindow = pset.get< bool >("DefaultSimWindow" );
261  fFullWaveformOutput = pset.get< bool >("FullWaveformOutput");
262  fReadoutWindow = pset.get< size_t >("ReadoutWindow" );
263  fPreTrigger = pset.get< size_t >("PreTrigger" );
264  fPadding = pset.get< int >("Padding" );
265  fGain = pset.get< double >("Gain" );
266  fDigiTree_SSP_LED = pset.get< bool >("SSP_LED_DigiTree" );
267  fNegativeSignal = pset.get< bool >("NegativeSignal" );
268  fCustomQE = pset.get< bool >("CustomQEperOpDet",false );
269  fCustomGain = pset.get< bool >("CustomGainperOpDet",false );
270  fGainsError = pset.get< double >("GainsError",0.0 ); //Factor of gain smearing in units of the gain value
271  fSPEWidth = pset.get< double >("SPEWidth",0.0 ); //Smear the SPE amplitude by the SPE width, in units of the gain value.
272  fExportWaveformTree = pset.get< bool >("ExportWaveformTree", false);
273 
274  if(fCustomQE) fQEVector = pset.get<std::vector<double>>("QEVector");
275  if(fCustomGain) fGainVector = pset.get<std::vector<double>>("GainVector");
276 
277  fThreshAlg = std::make_unique< pmtana::AlgoSSPLeadingEdge >
278  (pset.get< fhicl::ParameterSet >("algo_threshold"));
279 
280  if(fDigiTree_SSP_LED){
282  arvore2 = tfs->make<TTree>("PhotonData", "Photon_analysis");
283  arvore2->Branch("photon_opCh",&op_photon);
284  arvore2->Branch("photon_pulse",&t_photon);
285  }
286 //
288 
289  // Obtaining parameters from the DetectorClocksService
290  auto const clockData = art::ServiceHandle< detinfo::DetectorClocksService const>()->DataForJob();
291  fSampleFreq = clockData.OpticalClock().Frequency();
292  fDarkNoiseRate = odp->DarkRate();
293  fQE = odp->QE();
294 
295 // fSampleFreq = odp->SampleFreq(); //Sample Freq in MHz
296 
297  if (fDefaultSimWindow)
298  {
299  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob(clockData);
300  // Assume the readout starts at -1 drift window
301  fTimeBegin = -1*GetDriftWindow(detProp);
302 
303  // Take the TPC readout window size and convert
304  // to us with the electronics clock frequency
305  fTimeEnd = detProp.ReadOutWindowSize() / clockData.TPCClock().Frequency();
306 
307  fPreTrigger = 0;//Since we are using negative times as the pretrigger, PreTrigger is zero.
308  //fPreTrigger=GetDriftWindow()*clockData.OpticalClock().Frequency();
309  fReadoutWindow = (fTimeEnd- fTimeBegin)*clockData.OpticalClock().Frequency();
310  }
311  else
312  {
313  fTimeBegin = odp->TimeBegin();
314  fTimeEnd = odp->TimeEnd();
315 
316  }
317 
319 
320  // Initializing random number engines
321  unsigned int seed = pset.get< unsigned int >("Seed", sim::GetRandomNumberSeed());
322  auto& engine = createEngine(seed);
323  fRandGauss = std::make_unique< CLHEP::RandGauss >(engine);
324  fRandExponential = std::make_unique< CLHEP::RandExponential >(engine);
325  fRandFlat = std::make_unique< CLHEP::RandFlat >(engine);
326 
327 
329  nOpDet=geometry->NOpDets();
330 
331  if(fGainsError!=0)
332  {
333  std::cout << "Smearing gains vector by " << fGainsError << " percent. " << std::endl;
334  std::vector<double> newGains;
335  for (unsigned int pm=0; pm<nOpDet; pm++) newGains.push_back(fRandGauss->fire(getGain(pm), fGainsError*getGain(pm)));
336  std::cout << "Old gains vector: "; for (unsigned int pm=0; pm<nOpDet; pm++) std::cout << getGain(pm) << " "; std::cout << std::endl;
337  fCustomGain=true; fGainVector=newGains;
338  std::cout << "New gains vector: "; for (unsigned int pm=0; pm<nOpDet; pm++) std::cout << getGain(pm) << " "; std::cout << std::endl;
339  }
340 
341  fSinglePEWaveform = odp->SinglePEWaveform();
342  if(fSampleFreq!=250)
343  {
344  int rebin=(int)(1.0*250/fSampleFreq);
345  std::vector <double> aux;
346 
347  std::cout << "\tbefore: SinglePEWaveform: " << fSinglePEWaveform.size() <<"bins, rebin: "<< rebin<<std::endl;
348  if (rebin<(int)fSinglePEWaveform.size())
349  {
350  for(unsigned int i=0; i<fSinglePEWaveform.size();i++)
351  {
352  if(i%rebin==0) aux.push_back(fSinglePEWaveform[i]);
353  else aux[aux.size()-1]+=fSinglePEWaveform[i];
354  }
355  }
356  else
357  {
358  aux.resize(1);aux[0]=1.0;
359  }
360  fSinglePEWaveform=aux;
361  }
362  double normalization =0;
363  for(unsigned int i=0; i<fSinglePEWaveform.size();i++) normalization+=fSinglePEWaveform[i];
364  std::cout << "Generating waveforms of " << fTimeEnd - fTimeBegin << "us = "<< (fTimeEnd - fTimeBegin)*fSampleFreq <<" Samples"<< std::endl;
365  std::cout << "\tTimeBegin" << fTimeBegin <<" "<< std::endl;
366  std::cout << "\tfTimeEnd" << fTimeEnd <<" "<< std::endl;
367  std::cout << "\tSampleFreq" << fSampleFreq <<" "<< std::endl;
368  std::cout << "\tReadoutWindow" << fReadoutWindow <<" "<< std::endl;
369  std::cout << "\tfPreTrigger" << fPreTrigger <<" "<< std::endl;
370  std::cout << "\fSinglePEWaveform: " << fSinglePEWaveform.size() <<"bins, norm: "<< normalization<<std::endl;
371 
373  SampleSize = 1000.0/fSampleFreq;
374 
376 
377  if(fExportWaveformTree)
378  {
379  adc_value.resize(nOpDet);
380  for (unsigned int i=0; i<nOpDet; i++)adc_value[i] = (short*)malloc(sizeof(short)*nSamples);
381 
382  double *vaux; vaux = (double*)malloc(sizeof(double)*nOpDet);
383  fRunInfo = tfs->make<TTree>("RunInfo","MonteCarlo Run Info");
384  fRunInfo->Branch("SampleSize" , &SampleSize , "SampleSize/D" );
385  fRunInfo->Branch("nOpDet" , &nOpDet , Form("nOpDet/I") );
386  fRunInfo->Branch("Gains" , vaux , Form("Gains[nOpDet]/D") );
387  for (unsigned int i=0; i<nOpDet; i++) vaux[i]=getGain(i);
388  for (unsigned int i=0; i<nOpDet; i++) std::cout << " Gain [ " << i << " ] = " << vaux[i] << " ADCxticks/PE" << std::endl;
389  fRunInfo->Fill();
390 
391 
392  fWaveformTree = tfs->make<TTree>("WaveformTree","Waveforms Tree");
393  fWaveformTree->Branch("Run" , &Run , "Run/I" );
394  fWaveformTree->Branch("SubRun" , &SubRun , "SubRun/I" );
395  fWaveformTree->Branch("Event" , &Event , "Event/I" );
396  fWaveformTree->Branch("NSamples" , &nSamples , "NSamples/I" );
397  for (unsigned int i=0; i<nOpDet; i++) fWaveformTree->Branch(Form("adc_channel_%i",i), adc_value[i] , Form("adc_value_%i[NSamples]/S",i) );
398 
399  }
400 
401  }
402 
403  //---------------------------------------------------------------------------
405  {
407  {
408  Run = evt.run();
409  SubRun = evt.subRun();
410  Event = evt.event();
411  }
413 
414  // A pointer that will store produced OpDetWaveforms
415  std::unique_ptr< std::vector< raw::OpDetWaveform > >
416  pulseVecPtr(std::make_unique< std::vector< raw::OpDetWaveform > >());
417 
419  bool fUseLitePhotons = lgp->UseLitePhotons();
420 
421  if (!fUseLitePhotons)
423  << "Sorry, but for now only Lite Photon digitization is implemented!"
424  << '\n';
425 
426  // Geometry service
428 
429  // Service for determining optical detector responses
431 
432 // std::vector<art::Handle< std::vector< sim::SimPhotonsLite > > > litePhotonHandle;
433 
434  int modulecounter=0;
435 
436 
437 // unsigned int nChannelsPerOpDet=1;
438 
439  unsigned int nChannelsPerOpDet = geometry->NOpHardwareChannels(nOpDet);
440 
441 
442  // This vector stores waveforms created for each optical channel
443 
444  for(unsigned int opDet=0; opDet<nOpDet; opDet++)
445  {
446 
447  std::vector< std::vector< double > > pdWaveforms(nChannelsPerOpDet, std::vector< double >(nSamples, static_cast< double >(fPedestal)));
448  std::vector<std::vector<FocusList>> fls(nOpDet,std::vector<FocusList>(nChannelsPerOpDet, FocusList(nSamples, fPadding)));
449 
450  for(auto mod : fInputModule){
451  // For every module collection
452  //std::cout << "ANALYZING " << mod<<std::endl;
453  modulecounter++;
454 
456  if(fUseLitePhotons) litePhotonHandle = evt.getHandle< std::vector< sim::SimPhotonsLite > >(mod);
457 
459  if(!fUseLitePhotons) PhotonHandle = evt.getHandle< std::vector< sim::SimPhotons > >(mod);
460 
461  if(fUseLitePhotons)
462  {
463 
464  // For every optical detector:
465  for (auto const& litePhotons : (*litePhotonHandle))
466  {
467  // OpChannel in SimPhotonsLite is actually the photon detector number
468  //unsigned int opDet = litePhotons.OpChannel;
469  if(opDet == (unsigned)litePhotons.OpChannel)
470  {
471  // Get number of channels in this optical detector
472  //unsigned int nChannelsPerOpDet = geometry->NOpHardwareChannels(opDet);
473  // it is one by default in the dual phase geometry, but let's keep it to have compatible functions with single phase.
474 
475  CreatePDWaveform(litePhotons, *odResponse, *geometry, pdWaveforms, fls[opDet]);
476  if((unsigned)modulecounter<fInputModule.size()) continue;//==fInputModule.size()
477 
478  if(fExportWaveformTree) for (unsigned int hardwareChannel = 0;
479  hardwareChannel < nChannelsPerOpDet; ++hardwareChannel) fls[opDet][hardwareChannel].AddRange(0,SampleSize);
480 
481 
482  // Generate dark noise
483  if (fDarkNoiseRate > 0.0) AddDarkNoise(pdWaveforms, fls[opDet], opDet);
484 
485  // Uncomment to undo the effect of FocusLists. Replaces the accumulated
486  // lists with ones asserting we need to look at the whole trace.
487  // for(FocusList& fl: fls){
488  // fl.ranges.clear();
489  // fl.ranges.emplace_back(0, nSamples-1);
490  // }
491 
492  // Vary the pedestal
493 
494  if (fLineNoiseRMS > 0.0) AddLineNoise(pdWaveforms, fls[opDet]);
495 
496  // Loop over all the created waveforms, split them into shorter
497  // waveforms and use them to initialize OpDetWaveforms
498 
499  for (unsigned int hardwareChannel = 0;
500  hardwareChannel < nChannelsPerOpDet; ++hardwareChannel)
501  {
502  for(const std::pair<int, int>& p: fls[opDet][hardwareChannel].ranges){
503  // It's a shame we copy here. We could actually avoid by making the
504  // functions below take a begin()/end() pair.
505  const std::vector<double> sub(pdWaveforms[hardwareChannel].begin()+p.first,
506  pdWaveforms[hardwareChannel].begin()+p.second+1);
507 
508  std::vector< short > waveformOfShorts =
510  //std::cout << "waveformOfShorts " << waveformOfShorts.size()<< std::endl;
511 
512  std::map< size_t, std::vector < short > > mapTickWaveform =
514  SplitWaveform(waveformOfShorts, fls[opDet][hardwareChannel]) :
515  std::map< size_t, std::vector< short > >{ std::make_pair(0,
516  waveformOfShorts) };
517 
518  //std::cout << "mapTickWaveform " << mapTickWaveform.size()<< std::endl;
519  unsigned int opChannel = geometry->OpChannel(opDet, hardwareChannel);
520  for (auto const& pairTickWaveform : mapTickWaveform)
521  {
522  double timeStamp =
523  static_cast< double >(TickToTime(pairTickWaveform.first+p.first));
524  //std::cout << "\tp " << p<< std::endl;
525  //std::cout << "\ttimeStamp " << timeStamp<< ", pairTickWaveform.first " << pairTickWaveform.first <<", p.first " << p.first<< std::endl;
526  // std::cout << "\tpairTickWaveform.second.size() " << pairTickWaveform.second.size()<< std::endl;
527 
528  raw::OpDetWaveform adcVec(timeStamp, opChannel,
529  pairTickWaveform.second.size());
530  int counter=0;
531  for (short const& value : pairTickWaveform.second)
532  {
533  //std::cout <<"\t\tvalue " << value<< std::endl;
534  adcVec.emplace_back(value); counter++;
535  }
536  //std::cout <<"\t\tcounter " << counter<< std::endl;
537 
538  pulseVecPtr->emplace_back(std::move(adcVec));
539 
540  }//endloop per
541  }//endloop per Focus List <fls>
542  }//endloop per Hardware channel
543  }//endif pmt
544  }//endloop per OpDet
545  }//endif litephotons
546  else //if SimPhotons
547  {
548  // For every optical detector:
549  for (auto const& Photons : (*PhotonHandle))
550  {
551  if(opDet == (unsigned)Photons.OpChannel())
552  {
553  CreatePDWaveform(Photons, *odResponse, *geometry, pdWaveforms, fls[opDet]);
554  if((unsigned)modulecounter<fInputModule.size()) continue;//==fInputModule.size()
555 
556  if (fDarkNoiseRate > 0.0) AddDarkNoise(pdWaveforms, fls[opDet], opDet);
557 
558  if (fLineNoiseRMS > 0.0) AddLineNoise(pdWaveforms, fls[opDet]);
559 
560  // Loop over all the created waveforms, split them into shorter
561  // waveforms and use them to initialize OpDetWaveforms
562 
563  for (unsigned int hardwareChannel = 0;
564  hardwareChannel < nChannelsPerOpDet; ++hardwareChannel)
565  {
566  for(const std::pair<int, int>& p: fls[opDet][hardwareChannel].ranges){
567  // It's a shame we copy here. We could actually avoid by making the
568  // functions below take a begin()/end() pair.
569  const std::vector<double> sub(pdWaveforms[hardwareChannel].begin()+p.first,
570  pdWaveforms[hardwareChannel].begin()+p.second+1);
571 
572  std::vector< short > waveformOfShorts =
574 
575  std::map< size_t, std::vector < short > > mapTickWaveform =
577  SplitWaveform(waveformOfShorts, fls[opDet][hardwareChannel]) :
578  std::map< size_t, std::vector< short > >{ std::make_pair(0,
579  waveformOfShorts) };
580 
581  //std::cout << "mapTickWaveform " << mapTickWaveform.size()<< std::endl;
582  unsigned int opChannel = geometry->OpChannel(opDet, hardwareChannel);
583  for (auto const& pairTickWaveform : mapTickWaveform)
584  {
585  double timeStamp =
586  static_cast< double >(TickToTime(pairTickWaveform.first+p.first));
587  //std::cout << "\tp " << p<< std::endl;
588  //std::cout << "\ttimeStamp " << timeStamp<< ", pairTickWaveform.first " << pairTickWaveform.first <<", p.first " << p.first<< std::endl;
589  // std::cout << "\tpairTickWaveform.second.size() " << pairTickWaveform.second.size()<< std::endl;
590 
591  raw::OpDetWaveform adcVec(timeStamp, opChannel,
592  pairTickWaveform.second.size());
593  int counter=0;
594  for (short const& value : pairTickWaveform.second)
595  {
596  //std::cout <<"\t\tvalue " << value<< std::endl;
597  adcVec.emplace_back(value); counter++;
598  }
599  //std::cout <<"\t\tcounter " << counter<< std::endl;
600 
601  pulseVecPtr->emplace_back(std::move(adcVec));
602  }//endloop per
603  }//endloop per Focus List <fls>
604  }//endloop per Hardware channel
605  }//endif pmt
606  }//endloop per OpDet
607  } //endif litePhotons
608  if(fUseLitePhotons) litePhotonHandle.clear();
609  if(!fUseLitePhotons) PhotonHandle.clear();
610  }//endloop per Input Module (S1 and S2 light)
611 
612  }//endloop per pmt
614  {
615  arvore2->Fill();
616  t_photon.clear();
617  op_photon.clear();
618  }
619 
621  {
622  for(unsigned int j=0; j < pulseVecPtr->size(); j++)
623  {
624  unsigned int opDet= pulseVecPtr->at(j).ChannelNumber();
625 // std::cout <<" one - pmt - " << j << "\t" << opDet << "\t"<< pulseVecPtr->at(j).Waveform().size() << " " << pulseVecPtr->at(j).Waveform()[0] << std::endl;
626  for (unsigned int i = 0; i< pulseVecPtr->at(j).Waveform().size() ; i++) adc_value[opDet][i] = pulseVecPtr->at(j).Waveform()[i];
627  }
628  fWaveformTree->Fill();
629 
630  }
631 
632  // Push the OpDetWaveforms into the event
633  evt.put(std::move(pulseVecPtr));
634 
635 
636 
637 
638  }
639 
640  //---------------------------------------------------------------------------
641  void OpDetDigitizerDUNEDP::AddPulse(size_t timeBin,
642  int scale, std::vector< double >& waveform,
643  FocusList& fl, double Gain) const
644  {
645  if(timeBin>waveform.size()) return; //photon out of time
646 
647  size_t pulseLength = fSinglePEWaveform.size()+1;
648  if ((timeBin + fSinglePEWaveform.size()) > waveform.size())
649  pulseLength = (waveform.size() - timeBin);
650 
651  fl.AddRange(timeBin, timeBin+pulseLength-1);
652 
653 
654  double interbinallignment = fRandFlat->fire(1.0);//shifting the deposited charge among consecutive bins.
655  double customGain = fRandGauss->fire(Gain, fSPEWidth*Gain); //smear the SPE amplitude by the SPEwidth
656  // Adding a pulse to the waveform
657  for (size_t tick = 0; tick != pulseLength; ++tick)
658  {
659  double bincharge;
660  if (tick==0) bincharge=interbinallignment*scale*customGain*fSinglePEWaveform[tick];
661  else bincharge = interbinallignment*scale*customGain*fSinglePEWaveform[tick] + (1-interbinallignment)*scale*customGain*fSinglePEWaveform[tick-1];
662 
663  if(!fNegativeSignal) waveform[timeBin + tick] += bincharge;
664  else waveform[timeBin + tick] -= (double)bincharge;
665 
666  }
667  }
668 
669  //---------------------------------------------------------------------------
671  (sim::SimPhotonsLite const& litePhotons,
672  opdet::OpDetResponseInterface const& odResponse,
673  geo::Geometry const& geometry,
674  std::vector< std::vector< double > >& pdWaveforms,
675  std::vector<FocusList>& fls)
676  {
677 
678  unsigned int const opDet = litePhotons.OpChannel;
679  // This is int because otherwise detectedLite doesn't work
680  int readoutChannel;
681  // For a group of photons arriving at the same time this is a map
682  // of < arrival time (in ns), number of photons >
683  std::map< int, int > const& photonsMap = litePhotons.DetectedPhotons;
684 
685  int counter=0;
686  // For every pair of (arrival time, number of photons) in the map:
687 
688  for (auto const& pulse : photonsMap)
689  {
690  // Converting ns to us
691  double photonTime = static_cast< double >(pulse.first)/1000.0;
692 
693  int NumberOfPEs = PMTSaturationFunction(pulse.second);
694 // std::cout << "Adding " <<NumberOfPEs << "PEs @ " << photonTime << " in ch" << opDet << std::endl;
695  for (int i = 0; i < NumberOfPEs; ++i)
696  {
697  if ((photonTime >= fTimeBegin) && (photonTime < fTimeEnd))
698  {
699  // Sample a random subset according to QE
700  if (CLHEP::RandFlat::shoot(1.0) <getQE(opDet))
701  {
702  odResponse.detectedLite(opDet, readoutChannel);
703  unsigned int hardwareChannel = geometry.HardwareChannelFromOpChannel(readoutChannel);
704  // Convert the time of the pulse to ticks
705  size_t timeBin = TimeToTick(photonTime);
706  // Add 1 pulse to the waveform
707  AddPulse(timeBin, CrossTalk(), pdWaveforms.at(hardwareChannel), fls[hardwareChannel], getGain(opDet));
708  counter++;
709  unsigned int opChannel = geometry.OpChannel(opDet, hardwareChannel);
710  if(fDigiTree_SSP_LED){
711  op_photon.emplace_back(opChannel);
712  t_photon.emplace_back(photonTime);
713  }
714  }
715  //else std::cout << "photon not detected "<< std::endl;
716  }
717  }
718  }
719 // std::cout << "Created waveform for channel " << opDet << " with " << counter << " photons." << std::endl;
720 
721  }
722 
723  //---------------------------------------------------------------------------
725  (sim::SimPhotons const& Photons,
726  opdet::OpDetResponseInterface const& odResponse,
727  geo::Geometry const& geometry,
728  std::vector< std::vector< double > >& pdWaveforms,
729  std::vector<FocusList>& fls)
730  {
731 
732  unsigned int const opDet = Photons.OpChannel();
733  int readoutChannel;
734  // For a group of photons arriving at the same time this is a map
735  // of < arrival time (in ns), number of photons >
736  std::map< int, int > photonsMap;
737  for(unsigned int i=0; i<Photons.size();i++)photonsMap[Photons[i].Time]=0;
738  for(unsigned int i=0; i<Photons.size();i++)photonsMap[Photons[i].Time]++;
739  int counter=0;
740  // For every pair of (arrival time, number of photons) in the map:
741 
742  for (auto const& pulse : photonsMap)
743  {
744  // Converting ns to us
745  double photonTime = static_cast< double >(pulse.first)/1000.0;
746 
747  int NumberOfPEs = PMTSaturationFunction(pulse.second);
748 // std::cout << "Adding " <<NumberOfPEs << "PEs @ " << photonTime << " in ch" << opDet << std::endl;
749  for (int i = 0; i < NumberOfPEs; ++i)
750  {
751  if ((photonTime >= fTimeBegin) && (photonTime < fTimeEnd))
752  {
753  // Sample a random subset according to QE
754  if (CLHEP::RandFlat::shoot(1.0) <getQE(opDet))
755  {
756  odResponse.detectedLite(opDet, readoutChannel);
757  unsigned int hardwareChannel = geometry.HardwareChannelFromOpChannel(readoutChannel);
758  // Convert the time of the pulse to ticks
759  size_t timeBin = TimeToTick(photonTime);
760  // Add 1 pulse to the waveform
761  AddPulse(timeBin, CrossTalk(), pdWaveforms.at(hardwareChannel), fls[hardwareChannel], getGain(opDet));
762  counter++;
763  unsigned int opChannel = geometry.OpChannel(opDet, hardwareChannel);
764  if(fDigiTree_SSP_LED){
765  op_photon.emplace_back(opChannel);
766  t_photon.emplace_back(photonTime);
767  }
768  }
769  //else std::cout << "photon not detected "<< std::endl;
770  }
771  }
772  }
773 // std::cout << "Created waveform for channel " << opDet << " with " << counter << " photons." << std::endl;
774 
775  }
776 
777 
778  //---------------------------------------------------------------------------
780  AddLineNoise(std::vector< std::vector< double > >& waveforms,
781  const std::vector<FocusList>& fls) const
782  {
783  int i = 0;
784  for(auto& waveform : waveforms){
785  for(unsigned int j = 0; j < fls[i].ranges.size(); ++j){
786  const std::pair<int, int>& p = fls[i].ranges[j];
787  for(int k = p.first; k <= p.second; ++k){
788  waveform[k] += fRandGauss->fire(0, fLineNoiseRMS);
789  }
790  }
791 
792  ++i;
793  }
794  }
795  //---------------------------------------------------------------------------
797  AddDarkNoise(std::vector< std::vector< double > >& waveforms,
798  std::vector<FocusList>& fls, int optDet)
799  {
800  int i = 0;
801  for (auto& waveform : waveforms)
802  {
803  // Multiply by 10^6 since fDarkNoiseRate is in Hz
804  double darkNoiseTime = static_cast< double >(fRandExponential->
805  fire(1.0/fDarkNoiseRate)*1000000.0) + fTimeBegin;
806  while (darkNoiseTime < fTimeEnd)
807  {
808  size_t timeBin = TimeToTick(darkNoiseTime);
809  AddPulse(timeBin, CrossTalk(), waveform, fls[i], getGain(optDet) );
810  // Find next time to simulate a single PE pulse
811  darkNoiseTime += static_cast< double >
812  (fRandExponential->fire(1.0/fDarkNoiseRate)*1000000.0);
813  }
814 
815  ++i;
816  }
817  }
818 
819  //---------------------------------------------------------------------------
820  unsigned short OpDetDigitizerDUNEDP::CrossTalk() const
821  {
822  // Sometimes this should produce 3 or more PEs (not implemented)
823  if (fCrossTalk <= 0.0) return 1;
824  else if (fRandFlat->fire(1.0) > fCrossTalk) return 1;
825  else return 2;
826  }
827 
828  //---------------------------------------------------------------------------
830  (std::vector< double > const& vectorOfDoubles) const
831  {
832  // Don't bother to round properly, it's faster this way
833  return std::vector<short>(vectorOfDoubles.begin(), vectorOfDoubles.end());
834 
835  /*
836  std::vector< short > vectorOfShorts;
837  vectorOfShorts.reserve(vectorOfDoubles.size());
838  int i=0;
839  for (short const& value : vectorOfDoubles)
840  {
841  vectorOfShorts.emplace_back(static_cast< short >(std::round(value)));
842  }
843  return vectorOfShorts;
844  // */
845  }
846 
847  //---------------------------------------------------------------------------
848  std::map< size_t, std::vector< short > > OpDetDigitizerDUNEDP::
849  SplitWaveform(std::vector< short > const& waveform,
850  const FocusList& fls)
851  {
852 
853  std::map< size_t, std::vector< short > > mapTickWaveform;
854 
855  ::pmtana::PedestalMean_t ped_mean (waveform.size(),0);
856  ::pmtana::PedestalSigma_t ped_sigma(waveform.size(),0);
857 
858 
859  fThreshAlg->Reconstruct(waveform,ped_mean,ped_sigma);
860 
861  std::vector< pmtana::pulse_param > pulses;
862  for (size_t pulseCounter = 0; pulseCounter < fThreshAlg->GetNPulse();
863  ++pulseCounter)
864  pulses.emplace_back(fThreshAlg->GetPulse(pulseCounter));
865 
866  // We have to refine this algorithm later
867  for (pmtana::pulse_param const& pulse : pulses)
868  {
869  if (pulse.t_end <= pulse.t_start)
870  // Can I call it a logic error?
872  << "Pulse ends before or at the same time it starts!\n";
873 
875  waveform.begin() + static_cast< size_t >(pulse.t_start);
877  waveform.begin() + static_cast< size_t >(pulse.t_end );
878  mapTickWaveform.emplace(static_cast< size_t >(pulse.t_start),
879  std::vector< short >(window_start, window_end));
880  // Don't forget to check that the time output by the (new) algortihm
881  // is < fTimeEnd and > fTimeBegin!
882 
883  }
884 
885  return mapTickWaveform;
886 
887  }
888 
889  //---------------------------------------------------------------------------
891  {
892 
893  double driftWindow;
894 
895  double maxDrift = 0.0;
896  for (geo::TPCGeo const& tpc :
897  art::ServiceHandle< geo::Geometry >()->IterateTPCs())
898  if (maxDrift < tpc.DriftDistance()) maxDrift = tpc.DriftDistance();
899 
900  driftWindow = maxDrift / detProp.DriftVelocity();
901 
902  return driftWindow;
903 
904  }
905 
906  //---------------------------------------------------------------------------
908  {
909 
910 //std::cout << "tick "<< tick<< " " << (tick > fPreTrigger) << " " << (static_cast< double >(tick - fPreTrigger)/fSampleFreq
911  // + fTimeBegin) << " " << (static_cast< double >(fPreTrigger - tick)/fSampleFreq*(-1.0)
912 // + fTimeBegin)<< std::endl;
913  if (tick > fPreTrigger)
914  return (static_cast< double >(tick - fPreTrigger)/fSampleFreq
915  + fTimeBegin);
916  else
917  return (static_cast< double >(fPreTrigger - tick)/fSampleFreq*(-1.0)
918  + fTimeBegin);
919 
920  }
921 
922  //---------------------------------------------------------------------------
924  {
925 
926  return static_cast< size_t >(std::round((time - fTimeBegin)*fSampleFreq
927  + fPreTrigger));
928 
929  }
930 
931  //---------------------------------------------------------------------------
933  {
934 
935  // Are all these logic errors?
936 
937  if (fLineNoiseRMS < 0.0)
939  << "fLineNoiseRMS: " << fLineNoiseRMS << '\n'
940  << "Line noise RMS should be non-negative!\n";
941 
942  if (fDarkNoiseRate < 0.0)
944  << "fDarkNoiseRate: " << fDarkNoiseRate << '\n'
945  << "Dark noise rate should be non-negative!\n";
946 
949  << "PreTrigger: " << fPreTrigger << " and "
950  << "ReadoutWindow: " << fReadoutWindow << '\n'
951  << "Pretrigger window has to be shorter than readout window!\n";
952 
953  if (fTimeBegin >= fTimeEnd)
955  << "TimeBegin: " << fTimeBegin << " and "
956  << "TimeEnd: " << fTimeEnd << '\n'
957  << "TimeBegin should be less than TimeEnd!\n";
958 
959  }
960 
961  int OpDetDigitizerDUNEDP::PMTSaturationFunction(int photons){ return photons;}
962  double OpDetDigitizerDUNEDP::SumOfElements( std::vector<double> aa)
963  {
964  double sum=0;
965  for (auto& n : aa) sum+=n;
966  return sum;
967  }
968  double OpDetDigitizerDUNEDP::getQE(int OpDet)
969  {
970  if(!fCustomQE) return fQE;
971  else return fQEVector[OpDet];
972  }
973 
975  {
976  if(!fCustomGain) return fGain;
977  else return fGainVector[OpDet];
978  }
979 
980 }
base_engine_t & createEngine(seed_t seed)
Store parameters for running LArG4.
void CreatePDWaveform(sim::SimPhotonsLite const &, opdet::OpDetResponseInterface const &, geo::Geometry const &, std::vector< std::vector< double > > &, std::vector< FocusList > &)
EventNumber_t event() const
Definition: DataViewImpl.cc:85
void AddDarkNoise(std::vector< std::vector< double > > &, std::vector< FocusList > &fls, int opDet)
std::vector< double > PedestalSigma_t
int OpChannel() const
Returns the optical channel number this object is associated to.
Definition: SimPhotons.h:254
void AddPulse(size_t timeBin, int scale, std::vector< double > &waveform, FocusList &fl, double Gain) const
std::unique_ptr< CLHEP::RandFlat > fRandFlat
FocusList(int nSamples, int padding)
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
void AddRange(int from, int to)
double GetDriftWindow(detinfo::DetectorPropertiesData const &detProp) const
Geometry information for a single TPC.
Definition: TPCGeo.h:38
std::unique_ptr< CLHEP::RandGauss > fRandGauss
struct vector vector
unsigned int NOpHardwareChannels(int opDet) const
intermediate_table::const_iterator const_iterator
size_t TimeToTick(double time) const
std::map< int, int > DetectedPhotons
Number of photons detected at each given time: time tick -> photons.
Definition: SimPhotons.h:117
std::vector< std::pair< int, int > > ranges
double SumOfElements(std::vector< double >)
Simulation objects for optical detectors.
std::unique_ptr< pmtana::AlgoSSPLeadingEdge > fThreshAlg
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
unsigned int HardwareChannelFromOpChannel(int opChannel) const
Convert unique channel to hardware channel.
std::void_t< T > n
def move(depos, offset)
Definition: depos.py:107
double TickToTime(size_t tick) const
int OpChannel
Optical detector channel associated to this data.
Definition: SimPhotons.h:114
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:75
p
Definition: test.py:223
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
RunNumber_t run() const
Definition: DataViewImpl.cc:71
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
std::unique_ptr< CLHEP::RandExponential > fRandExponential
virtual bool detectedLite(int OpChannel, int &newOpChannel) const
std::vector< short > VectorOfDoublesToVectorOfShorts(std::vector< double > const &) const
std::vector< std::string > fInputModule
void AddLineNoise(std::vector< std::vector< double > > &, const std::vector< FocusList > &fls) 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
Collection of photons which recorded on one channel.
Definition: SimPhotons.h:136
Compact representation of photons on a channel.
Definition: SimPhotons.h:103
void clear()
Definition: Handle.h:236
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:72
OpDetDigitizerDUNEDP(fhicl::ParameterSet const &)
Tools and modules for checking out the basics of the Monte Carlo.
TCEvent evt
Definition: DataStructs.cxx:7
def rebin(a, args)
Definition: arrays.py:2
std::vector< double > PedestalMean_t
std::map< size_t, std::vector< short > > SplitWaveform(std::vector< short > const &, const FocusList &)
bool UseLitePhotons() const
QTextStream & endl(QTextStream &s)
std::string sub(const std::string &a, const std::string &b)
Definition: TruthText.cxx:100