SimWireDUNE10kt_module.cc
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////
2 //
3 // SimWireDUNE10kt class designed to simulate signal on a wire in the TPC
4 //
5 // katori@fnal.gov
6 //
7 // jti3@fnal.gov
8 // - Revised to use sim::RawDigit instead of rawdata::RawDigit, and to
9 // - save the electron clusters associated with each digit.
10 //
11 /////////////////////////////////////////////////////////////////////////
12 
13 #include <vector>
14 #include <string>
15 #include <algorithm>
16 #include <sstream>
17 #include <fstream>
18 #include <bitset>
19 
20 extern "C" {
21 #include <sys/types.h>
22 #include <sys/stat.h>
23 }
24 
30 #include "art_root_io/TFileService.h"
31 #include "art_root_io/TFileDirectory.h"
32 #include "fhiclcpp/ParameterSet.h"
34 
35 // art extensions
36 #include "nurandom/RandomUtils/NuRandomService.h"
37 
39 #include "lardataobj/RawData/raw.h"
43 
48 
49 #include "TMath.h"
50 #include "TComplex.h"
51 #include "TString.h"
52 #include "TH2.h"
53 #include "TH1D.h"
54 #include "TFile.h"
55 
56 #include "CLHEP/Random/RandFlat.h"
57 #include "CLHEP/Random/RandGaussQ.h"
58 
59 ///Detector simulation of raw signals on wires
60 namespace detsim {
61 
62  // Base class for creation of raw signals on wires.
64 
65  public:
66 
67  explicit SimWireDUNE10kt(fhicl::ParameterSet const& pset);
68 
69  // read/write access to event
70  void produce (art::Event& evt);
71  void beginJob();
72  void endJob();
73  void reconfigure(fhicl::ParameterSet const& p);
74 
75  private:
76 
77  void GenNoise(std::vector<float>& array);
78 
79  std::string fDriftEModuleLabel;///< module making the ionization electrons
80  raw::Compress_t fCompression; ///< compression type to use
81  unsigned int fNoiseOn; ///< noise turned on or off for debugging; default is on
82  unsigned int fNoiseModel; ///< noise model
83  float fNoiseFact; ///< noise scale factor
84  float fNoiseWidth; ///< exponential noise width (kHz)
85  float fLowCutoff; ///< low frequency filter cutoff (kHz)
86  float fNoiseFactZ; ///< noise scale factor for Z (collection) plane
87  float fNoiseWidthZ; ///< exponential noise width (kHz) for Z (collection) plane
88  float fLowCutoffZ; ///< low frequency filter cutoff (kHz) for Z (collection) plane
89  float fNoiseFactU; ///< noise scale factor for U plane
90  float fNoiseWidthU; ///< exponential noise width (kHz) for U plane
91  float fLowCutoffU; ///< low frequency filter cutoff (kHz) for U plane
92  float fNoiseFactV; ///< noise scale factor for V plane
93  float fNoiseWidthV; ///< exponential noise width (kHz) for V plane
94  float fLowCutoffV; ///< low frequency filter cutoff (kHz) for V plane
95  unsigned int fZeroThreshold; ///< Zero suppression threshold
96  int fNearestNeighbor; ///< Maximum distance between hits above threshold before they are separated into different blocks
97  int fNTicks; ///< number of ticks of the clock
98  double fSampleRate; ///< sampling rate in ns
99  unsigned int fNSamplesReadout; ///< number of ADC readout samples in 1 readout frame
100  unsigned int fNTimeSamples; ///< number of ADC readout samples in all readout frames (per event)
101  unsigned int fNoiseArrayPoints; ///< number of points in randomly generated noise array
102 
103  std::vector<double> fChargeWork;
104  //std::vector< std::vector<float> > fNoise;///< noise on each channel for each time
105  std::vector< std::vector<float> > fNoiseZ;///< noise on each channel for each time for Z (collection) plane
106  std::vector< std::vector<float> > fNoiseU;///< noise on each channel for each time for U plane
107  std::vector< std::vector<float> > fNoiseV;///< noise on each channel for each time for V plane
108 
109  TH1D* fNoiseDist; ///< distribution of noise counts
110 
111  //define max ADC value - if one wishes this can
112  //be made a fcl parameter but not likely to ever change
113  const float adcsaturation = 4095;
114 
115  bool fSaveEmptyChannel; // switch for saving channels with all zero entries
116  float fCollectionPed; ///< ADC value of baseline for collection plane
117  float fInductionPed; ///< ADC value of baseline for induction plane
118  CLHEP::HepRandomEngine& fEngine;
119  }; // class SimWireDUNE10kt
120 
121  //-------------------------------------------------
123  : EDProducer{pset}
124  // create a default random engine; obtain the random seed from NuRandomService,
125  // unless overridden in configuration with key "Seed"
127  {
128 
129  this->reconfigure(pset);
130 
131  produces< std::vector<raw::RawDigit> >();
133  produces< std::vector<raw::RawDigit> >("preSpill");
134  produces< std::vector<raw::RawDigit> >("postSpill");
135  }
136 
138  TString compression(pset.get< std::string >("CompressionType"));
139  if(compression.Contains("Huffman",TString::kIgnoreCase)) fCompression = raw::kHuffman;
140  if(compression.Contains("ZeroSuppression",TString::kIgnoreCase)) fCompression = raw::kZeroSuppression;
141  }
142 
143  //-------------------------------------------------
145  {
146  fDriftEModuleLabel= p.get< std::string >("DriftEModuleLabel");
147 
148 
149  fNoiseFactZ = p.get< double >("NoiseFactZ");
150  fNoiseWidthZ = p.get< double >("NoiseWidthZ");
151  fLowCutoffZ = p.get< double >("LowCutoffZ");
152  fNoiseFactU = p.get< double >("NoiseFactU");
153  fNoiseWidthU = p.get< double >("NoiseWidthU");
154  fLowCutoffU = p.get< double >("LowCutoffU");
155  fNoiseFactV = p.get< double >("NoiseFactV");
156  fNoiseWidthV = p.get< double >("NoiseWidthV");
157  fLowCutoffV = p.get< double >("LowCutoffV");
158  fZeroThreshold = p.get< unsigned int >("ZeroThreshold");
159  fNearestNeighbor = p.get< int >("NearestNeighbor");
160  fNoiseArrayPoints = p.get< unsigned int >("NoiseArrayPoints");
161  fNoiseOn = p.get< unsigned int >("NoiseOn");
162  fNoiseModel = p.get< unsigned int >("NoiseModel");
163  fCollectionPed = p.get< float >("CollectionPed");
164  fInductionPed = p.get< float >("InductionPed");
165  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
166  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob(clockData);
167  fSampleRate = sampling_rate(clockData);
168  fNSamplesReadout = detProp.ReadOutWindowSize();
169  fNTimeSamples = detProp.NumberTimeSamples();
170  fSaveEmptyChannel = p.get< bool >("SaveEmptyChannel");
171  return;
172  }
173 
174  //-------------------------------------------------
176  {
177 
178  // get access to the TFile service
180 
181  fNoiseDist = tfs->make<TH1D>("Noise", ";Noise (ADC);", 1000, -10., 10.);
182 
184  fNTicks = fFFT->FFTSize();
185 
186  if ( fNTicks%2 != 0 )
187  MF_LOG_DEBUG("SimWireDUNE10kt") << "Warning: FFTSize not a power of 2. "
188  << "May cause issues in (de)convolution.\n";
189 
190  if ( (int)fNSamplesReadout > fNTicks )
191  mf::LogError("SimWireDUNE10kt") << "Cannot have number of readout samples "
192  << "greater than FFTSize!";
193 
194  fChargeWork.resize(fNTicks, 0.);
196 
197  //Generate noise if selected to be on
198  if(fNoiseOn && fNoiseModel==1){
199 
200  //fNoise.resize(geo->Nchannels());
201  fNoiseZ.resize(fNoiseArrayPoints);
202  fNoiseU.resize(fNoiseArrayPoints);
203  fNoiseV.resize(fNoiseArrayPoints);
204 
205  // GenNoise() will further resize each channel's
206  // fNoise vector to fNoiseArrayPoints long.
207 
208  for(unsigned int p = 0; p < fNoiseArrayPoints; ++p){
212 
213  GenNoise(fNoiseZ[p]);
214  for(int i = 0; i < fNTicks; ++i)
215  fNoiseDist->Fill(fNoiseZ[p][i]);
216 
220 
221  GenNoise(fNoiseU[p]);
222  for(int i = 0; i < fNTicks; ++i)
223  fNoiseDist->Fill(fNoiseU[p][i]);
224 
225 
229 
230 
231  GenNoise(fNoiseV[p]);
232  for(int i = 0; i < fNTicks; ++i)
233  fNoiseDist->Fill(fNoiseV[p][i]);
234 
235  }// end loop over wires
236  }
237  return;
238 
239  }
240 
241  //-------------------------------------------------
243  {
244  }
245 
246  //-------------------------------------------------
248  {
249  // get the geometry to be able to figure out signal types and chan -> plane mappings
251  unsigned int signalSize = fNTicks;
252 
253  // vectors for working
254  std::vector<short> adcvec(signalSize, 0);
255  std::vector<short> adcvecPreSpill(signalSize, 0);
256  std::vector<short> adcvecPostSpill(signalSize, 0);
257  std::vector<const sim::SimChannel*> chanHandle;
258  evt.getView(fDriftEModuleLabel,chanHandle);
259 
260  //Get fIndShape and fColShape from SignalShapingService, on the fly
262 
263  // make a vector of const sim::SimChannel* that has same number
264  // of entries as the number of channels in the detector
265  // and set the entries for the channels that have signal on them
266  // using the chanHandle
267  std::vector<const sim::SimChannel*> channels(geo->Nchannels());
268  for(size_t c = 0; c < chanHandle.size(); ++c){
269  channels[chanHandle[c]->Channel()] = chanHandle[c];
270  }
271 
272  // whether or not to do the prespill and postspill digitization:
273  const bool prepost = (fNSamplesReadout != fNTimeSamples);
274 
275  // make an unique_ptr of sim::SimDigits that allows ownership of the produced
276  // digits to be transferred to the art::Event after the put statement below
277  std::unique_ptr< std::vector<raw::RawDigit> > digcol(new std::vector<raw::RawDigit>);
278  std::unique_ptr< std::vector<raw::RawDigit> > digcolPreSpill(prepost? new std::vector<raw::RawDigit>: nullptr);
279  std::unique_ptr< std::vector<raw::RawDigit> > digcolPostSpill(prepost? new std::vector<raw::RawDigit>: nullptr);
280 
281  unsigned int chan = 0;
282  fChargeWork.clear();
283  fChargeWork.resize(fNTicks, 0.);
284 
285  std::vector<double> fChargeWorkPreSpill, fChargeWorkPostSpill;
286 
288 
289  // Add all channels
290  CLHEP::RandFlat flat(fEngine);
291 
293 
294  digcol->reserve(geo->Nchannels());
295  if (prepost) {
296  digcolPreSpill->reserve(geo->Nchannels());
297  digcolPostSpill->reserve(geo->Nchannels());
298  }
299 
300  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
301  for(chan = 0; chan < geo->Nchannels(); chan++) {
302 
303  fChargeWork.clear();
304  // fChargeWork.resize(fNTicks, 0.);
305  fChargeWork.resize(fNTimeSamples, 0.);
306 
307 
308  if (prepost) {
309  fChargeWorkPreSpill.clear();
310  fChargeWorkPostSpill.clear();
311  fChargeWorkPreSpill.resize(fNTicks,0);
312  fChargeWorkPostSpill.resize(fNTicks,0);
313  }
314 
315  // get the sim::SimChannel for this channel
316  const sim::SimChannel* sc = channels[chan];
317 
318  const geo::View_t view = geo->View(chan);
319 
320  if( sc ){
321 
322  // loop over the tdcs and grab the number of electrons for each
323  for(size_t t = 0; t < fChargeWork.size(); ++t)
324  fChargeWork[t] = sc->Charge(t);
325 
326 
327  // Convolve charge with appropriate response function
328  if(prepost) {
329  fChargeWorkPreSpill.assign(fChargeWork.begin(), fChargeWork.begin()+fNSamplesReadout);
330  fChargeWorkPostSpill.assign(fChargeWork.begin()+2*fNSamplesReadout, fChargeWork.end());
331 
332  fChargeWork.erase( fChargeWork.begin()+2*fNSamplesReadout, fChargeWork.end() );
333  fChargeWork.erase( fChargeWork.begin(), fChargeWork.begin()+fNSamplesReadout );
334 
335  fChargeWorkPreSpill.resize(fNTicks,0);
336  fChargeWorkPostSpill.resize(fNTicks,0);
337  sss->Convolute(clockData, chan, fChargeWorkPreSpill);
338  sss->Convolute(clockData, chan, fChargeWorkPostSpill);
339  }
340  fChargeWork.resize(fNTicks,0);
341  sss->Convolute(clockData, chan,fChargeWork);
342 
343 
344  }
345 
346  float ped_mean = fCollectionPed;
347  geo::SigType_t sigtype = geo->SignalType(chan);
348  if (sigtype == geo::kInduction){
349  ped_mean = fInductionPed;
350  }
351  else if (sigtype == geo::kCollection){
352  ped_mean = fCollectionPed;
353  }
354  // noise was already generated for each wire in the event
355  // raw digit vec is already in channel order
356  // pick a new "noise channel" for every channel - this makes sure
357  // the noise has the right coherent characteristics to be on one channel
358 
359 
360  int noisechan = nearbyint(flat.fire()*(1.*(fNoiseArrayPoints-1)+0.1));
361  int noisechanpre = nearbyint(flat.fire()*(1.*(fNoiseArrayPoints-1)+0.1));
362  int noisechanpost = nearbyint(flat.fire()*(1.*(fNoiseArrayPoints-1)+0.1));
363 
364  // optimize for speed -- access vectors as arrays
365 
366  double *fChargeWork_a=0;
367  double *fChargeWorkPreSpill_a=0;
368  double *fChargeWorkPostSpill_a=0;
369  short *adcvec_a=0;
370  short *adcvecPreSpill_a=0;
371  short *adcvecPostSpill_a=0;
372  float *noise_a_U=0;
373  float *noise_a_V=0;
374  float *noise_a_Z=0;
375  float *noise_a_Upre=0;
376  float *noise_a_Vpre=0;
377  float *noise_a_Zpre=0;
378  float *noise_a_Upost=0;
379  float *noise_a_Vpost=0;
380  float *noise_a_Zpost=0;
381 
382  if (signalSize>0) {
383  fChargeWork_a = fChargeWork.data();
384  adcvec_a = adcvec.data();
385  if (prepost) {
386  fChargeWorkPreSpill_a = fChargeWorkPreSpill.data();
387  fChargeWorkPostSpill_a = fChargeWorkPostSpill.data();
388  adcvecPreSpill_a = adcvecPreSpill.data();
389  adcvecPostSpill_a = adcvecPostSpill.data();
390  }
391  if (fNoiseOn && fNoiseModel==1) {
392  noise_a_U=(fNoiseU[noisechan]).data();
393  noise_a_V=(fNoiseV[noisechan]).data();
394  noise_a_Z=(fNoiseZ[noisechan]).data();
395  if (prepost) {
396  noise_a_Upre=(fNoiseU[noisechanpre]).data();
397  noise_a_Vpre=(fNoiseV[noisechanpre]).data();
398  noise_a_Zpre=(fNoiseZ[noisechanpre]).data();
399  noise_a_Upost=(fNoiseU[noisechanpost]).data();
400  noise_a_Vpost=(fNoiseV[noisechanpost]).data();
401  noise_a_Zpost=(fNoiseZ[noisechanpost]).data();
402  }
403  }
404  }
405 
406  float tmpfv=0; // this is here so we do our own rounding from floats to short ints (saves CPU time)
407  float tnoise=0;
408  float tnoisepre=0;
409  float tnoisepost=0;
410 
411  if (view != geo::kU && view != geo::kV && view != geo::kZ) {
412  mf::LogError("SimWireDUNE10kt") << "ERROR: CHANNEL NUMBER " << chan << " OUTSIDE OF PLANE";
413  }
414 
415  //std::cout << "Xin " << fNoiseOn << " " << fNoiseModel << std::endl;
416 
417  if(fNoiseOn && fNoiseModel==1) {
418  for(unsigned int i = 0; i < signalSize; ++i){
419  if(view==geo::kU) { tnoise = noise_a_U[i]; }
420  else if (view==geo::kV) { tnoise = noise_a_V[i]; }
421  else { tnoise = noise_a_Z[i]; }
422  tmpfv = tnoise + fChargeWork_a[i];
423  //allow for ADC saturation
424  if ( tmpfv > adcsaturation - ped_mean)
425  tmpfv = adcsaturation- ped_mean;
426  //don't allow for "negative" saturation
427  if ( tmpfv < 0 - ped_mean)
428  tmpfv = 0- ped_mean;
429  adcvec_a[i] = (tmpfv >=0) ? (short) (tmpfv+0.5) : (short) (tmpfv-0.5);
430  }
431  if (prepost) {
432  for(unsigned int i = 0; i < signalSize; ++i){
433  if(view==geo::kU) { tnoisepre = noise_a_Upre[i]; tnoisepost = noise_a_Upost[i]; }
434  else if(view==geo::kV) { tnoisepre = noise_a_Vpre[i]; tnoisepost = noise_a_Vpost[i]; }
435  else { tnoisepre = noise_a_Zpre[i]; tnoisepost = noise_a_Zpost[i]; }
436 
437  tmpfv = tnoisepre + fChargeWorkPreSpill_a[i];
438  //allow for ADC saturation
439  if ( tmpfv > adcsaturation - ped_mean){
440  tmpfv = adcsaturation- ped_mean;
441  }
442  //don't allow for "negative" saturation
443  if ( tmpfv < 0- ped_mean ){
444  tmpfv = 0- ped_mean;
445  }
446  adcvecPreSpill_a[i] = (tmpfv >=0) ? (short) (tmpfv+0.5) : (short) (tmpfv-0.5);
447 
448  tmpfv = tnoisepost + fChargeWorkPostSpill_a[i];
449  //allow for ADC saturation
450  if ( tmpfv > adcsaturation - ped_mean){
451  tmpfv = adcsaturation- ped_mean;
452  }
453  //don't allow for "negative" saturation
454  if ( tmpfv < 0 - ped_mean){
455  tmpfv = 0- ped_mean;
456  }
457  adcvecPostSpill_a[i] = (tmpfv >=0) ? (short) (tmpfv+0.5) : (short) (tmpfv-0.5);
458  }
459  }
460  }else if (fNoiseOn && fNoiseModel==2){
461  float fASICGain = sss->GetASICGain(chan);
462  double fShapingTime = sss->GetShapingTime(chan);
463  std::map< double, int > fShapingTimeOrder;
464  fShapingTimeOrder = { {0.5, 0}, {1.0, 1}, {2.0, 2}, {3.0, 3} };
465  DoubleVec fNoiseFactVec;
466  auto tempNoiseVec = sss->GetNoiseFactVec();
467  if ( fShapingTimeOrder.find( fShapingTime ) != fShapingTimeOrder.end() ){
468  size_t i = 0;
469  fNoiseFactVec.resize(2);
470  for (auto& item : tempNoiseVec) {
471  fNoiseFactVec[i] = item.at( fShapingTimeOrder.find( fShapingTime )->second );
472  fNoiseFactVec[i] *= fASICGain/4.7;
473  ++i;
474  }
475  }
476  else {//Throw exception...
477  throw cet::exception("SimWireMicroBooNE")
478  << "\033[93m"
479  << "Shaping Time received from signalservices_microboone.fcl is not one of allowed values"
480  << std::endl
481  << "Allowed values: 0.5, 1.0, 2.0, 3.0 usec"
482  << "\033[00m"
483  << std::endl;
484  }
485  // std::cout << "Xin " << fASICGain << " " << fShapingTime << " " << fNoiseFactVec[0] << " " << fNoiseFactVec[1] << std::endl;
486 
487  CLHEP::RandGaussQ rGauss_Ind(fEngine, 0.0, fNoiseFactVec[0]);
488  CLHEP::RandGaussQ rGauss_Col(fEngine, 0.0, fNoiseFactVec[1]);
489 
490 
491  for(unsigned int i = 0; i < signalSize; ++i){
492  if(view==geo::kU) { tnoise = rGauss_Ind.fire(); }
493  else if (view==geo::kV) { tnoise = rGauss_Ind.fire(); }
494  else { tnoise = rGauss_Col.fire(); }
495  tmpfv = tnoise + fChargeWork_a[i];
496  //allow for ADC saturation
497  if ( tmpfv > adcsaturation - ped_mean)
498  tmpfv = adcsaturation- ped_mean;
499  //don't allow for "negative" saturation
500  if ( tmpfv < 0 - ped_mean)
501  tmpfv = 0- ped_mean;
502  adcvec_a[i] = (tmpfv >=0) ? (short) (tmpfv+0.5) : (short) (tmpfv-0.5);
503  }
504  if (prepost) {
505  for(unsigned int i = 0; i < signalSize; ++i){
506  if(view==geo::kU) { tnoisepre = rGauss_Ind.fire(); tnoisepost = rGauss_Ind.fire(); }
507  else if(view==geo::kV) { tnoisepre = rGauss_Ind.fire(); tnoisepost = rGauss_Ind.fire(); }
508  else { tnoisepre = rGauss_Col.fire(); tnoisepost = rGauss_Col.fire(); }
509 
510  tmpfv = tnoisepre + fChargeWorkPreSpill_a[i];
511  //allow for ADC saturation
512  if ( tmpfv > adcsaturation - ped_mean){
513  tmpfv = adcsaturation- ped_mean;
514  }
515  //don't allow for "negative" saturation
516  if ( tmpfv < 0 - ped_mean){
517  tmpfv = 0- ped_mean;
518  }
519  adcvecPreSpill_a[i] = (tmpfv >=0) ? (short) (tmpfv+0.5) : (short) (tmpfv-0.5);
520 
521  tmpfv = tnoisepost + fChargeWorkPostSpill_a[i];
522  //allow for ADC saturation
523  if ( tmpfv > adcsaturation - ped_mean){
524  tmpfv = adcsaturation- ped_mean;
525  }
526  //don't allow for "negative" saturation
527  if ( tmpfv < 0 - ped_mean){
528  tmpfv = 0- ped_mean;
529  }
530  adcvecPostSpill_a[i] = (tmpfv >=0) ? (short) (tmpfv+0.5) : (short) (tmpfv-0.5);
531  }
532  }
533  }
534  else { // no noise, so just round the values to nearest short ints and store them
535  for(unsigned int i = 0; i < signalSize; ++i){
536  tmpfv = fChargeWork_a[i];
537  //allow for ADC saturation
538  if ( tmpfv > adcsaturation- ped_mean )
539  tmpfv = adcsaturation- ped_mean;
540  //don't allow for "negative" saturation
541  if ( tmpfv < 0 - ped_mean)
542  tmpfv = 0- ped_mean;
543  adcvec_a[i] = (tmpfv >=0) ? (short) (tmpfv+0.5) : (short) (tmpfv-0.5);
544  }
545  if (prepost) {
546  for(unsigned int i = 0; i < signalSize; ++i){
547  tmpfv = fChargeWorkPreSpill_a[i];
548  //allow for ADC saturation
549  if ( tmpfv > adcsaturation - ped_mean)
550  tmpfv = adcsaturation- ped_mean;
551  //don't allow for "negative" saturation
552  if ( tmpfv < 0 - ped_mean) {
553  tmpfv = 0- ped_mean; }
554  adcvecPreSpill_a[i] = (tmpfv >=0) ? (short) (tmpfv+0.5) : (short) (tmpfv-0.5);
555  tmpfv = fChargeWorkPostSpill_a[i];
556  //allow for ADC saturation
557  if ( tmpfv > adcsaturation - ped_mean)
558  tmpfv = adcsaturation- ped_mean;
559  //don't allow for "negative" saturation
560  if ( tmpfv < 0 - ped_mean) {
561  tmpfv = 0- ped_mean; }
562  adcvecPostSpill_a[i] = (tmpfv >=0) ? (short) (tmpfv+0.5) : (short) (tmpfv-0.5);
563  }
564  }
565  }
566 
567  // resize the adcvec to be the correct number of time samples,
568  // just drop the extra samples
569 
570  // compress the adc vector using the desired compression scheme,
571  // if raw::kNone is selected nothing happens to adcvec
572  // This shrinks adcvec, if fCompression is not kNone.
573 
574  //std::cout << "Channel view: " << view << std::endl;
575 
576  adcvec.resize(fNSamplesReadout);
578  raw::RawDigit rd(chan, fNSamplesReadout, adcvec, fCompression);
579  //rd.SetPedestal(ped_mean);
580  adcvec.resize(signalSize); // Then, resize adcvec back to full length. Do not initialize to zero (slow)
581 
582  if(fSaveEmptyChannel || adcvec[1]>0)
583  digcol->push_back(rd); // add this digit to the collection
584 
585  // add the digit to the collection (in-place constructon)
586  // digcol->emplace_back(chan, fNSamplesReadout, adcvec, fCompression);
587  // adcvec.resize(signalSize); // Then, resize adcvec back to full length. Do not initialize to zero (slow)
588 
589  // do all this for the prespill and postspill samples if need be
590  if (prepost) {
591  adcvecPreSpill.resize(fNSamplesReadout);
592  adcvecPostSpill.resize(fNSamplesReadout);
595  raw::RawDigit rdPreSpill(chan, fNSamplesReadout, adcvecPreSpill, fCompression);
596  raw::RawDigit rdPostSpill(chan, fNSamplesReadout, adcvecPostSpill, fCompression);
597  // rdPreSpill.SetPedestal(ped_mean);
598  // rdPostSpill.SetPedestal(ped_mean);
599  adcvecPreSpill.resize(signalSize);
600  adcvecPostSpill.resize(signalSize);
601 
602  if(fSaveEmptyChannel || adcvecPreSpill[1]>0)
603  digcolPreSpill->push_back(rdPreSpill);
604  if(fSaveEmptyChannel || adcvecPostSpill[1]>0)
605  digcolPostSpill->push_back(rdPostSpill);
606  // raw::Compress(adcvecPreSpill, fCompression, fZeroThreshold, fNearestNeighbor);
607  // raw::Compress(adcvecPostSpill, fCompression, fZeroThreshold, fNearestNeighbor);
608  // digcolPreSpill->emplace_back(chan, fNSamplesReadout, adcvecPreSpill, fCompression);
609  // digcolPostSpill->emplace_back(chan, fNSamplesReadout, adcvecPostSpill, fCompression);
610  // adcvecPreSpill.resize(signalSize);
611  // adcvecPostSpill.resize(signalSize);
612  }
613 
614  }// end loop over channels
615 
616  evt.put(std::move(digcol));
617  if(prepost) {
618  evt.put(std::move(digcolPreSpill), "preSpill");
619  evt.put(std::move(digcolPostSpill), "postSpill");
620  }
621 
622  return;
623  }
624 
625  //-------------------------------------------------
626  void SimWireDUNE10kt::GenNoise(std::vector<float>& noise)
627  {
628  CLHEP::RandFlat flat(fEngine);
629 
630  noise.clear();
631  noise.resize(fNTicks,0.0);
632  // noise in frequency space
633  std::vector<TComplex> noiseFrequency(fNTicks/2+1, 0.);
634 
635  double pval = 0.;
636  double lofilter = 0.;
637  double phase = 0.;
638  double rnd[2] = {0.};
639 
640  // width of frequencyBin in kHz
641  double binWidth = 1.0/(fNTicks*fSampleRate*1.0e-6);
642  for(int i=0; i< fNTicks/2+1; ++i){
643  // exponential noise spectrum
644  pval = fNoiseFact*exp(-(double)i*binWidth/fNoiseWidth);
645  // low frequency cutoff
646  lofilter = 1.0/(1.0+exp(-(i-fLowCutoff/binWidth)/0.5));
647  // randomize 10%
648  flat.fireArray(2,rnd,0,1);
649  pval *= lofilter*(0.9+0.2*rnd[0]);
650  // random pahse angle
651  phase = rnd[1]*2.*TMath::Pi();
652 
653  TComplex tc(pval*cos(phase),pval*sin(phase));
654  noiseFrequency[i] += tc;
655  }
656 
657  // std::cout << "filled noise freq" << std::endl;
658 
659  // inverse FFT MCSignal
661  fFFT->DoInvFFT(noiseFrequency, noise);
662 
663  noiseFrequency.clear();
664 
665  // multiply each noise value by fNTicks as the InvFFT
666  // divides each bin by fNTicks assuming that a forward FFT
667  // has already been done.
668  for(unsigned int i = 0; i < noise.size(); ++i) noise[i] *= 1.*fNTicks;
669  }
670 }
671 
unsigned int fNoiseArrayPoints
number of points in randomly generated noise array
float fNoiseWidthV
exponential noise width (kHz) for V plane
intermediate_table::iterator iterator
base_engine_t & createEngine(seed_t seed)
Huffman Encoding.
Definition: RawTypes.h:10
std::string fDriftEModuleLabel
module making the ionization electrons
double GetShapingTime(Channel channel) const override
enum raw::_compress Compress_t
float fLowCutoffU
low frequency filter cutoff (kHz) for U plane
Collection of charge vs time digitized from a single readout channel.
Definition: RawDigit.h:69
Energy deposited on a readout channel by simulated tracks.
Definition: SimChannel.h:140
float fNoiseFactV
noise scale factor for V plane
unsigned int fNoiseModel
noise model
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
std::string string
Definition: nybbler.cc:12
Planes which measure V.
Definition: geo_types.h:130
std::vector< DoubleVec > GetNoiseFactVec() const override
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
float fLowCutoff
low frequency filter cutoff (kHz)
TH1D * fNoiseDist
distribution of noise counts
Detector simulation of raw signals on wires.
float fInductionPed
ADC value of baseline for induction plane.
Planes which measure Z direction.
Definition: geo_types.h:132
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
std::vector< std::vector< float > > fNoiseV
noise on each channel for each time for V plane
unsigned int fNoiseOn
noise turned on or off for debugging; default is on
void reconfigure(fhicl::ParameterSet const &p)
Service to provide microboone-specific signal shaping for simulation (convolution) and reconstruction...
no compression
Definition: RawTypes.h:9
art framework interface to geometry description
void DoInvFFT(std::vector< TComplex > &input, std::vector< T > &output)
Definition: LArFFT.h:120
unsigned int Nchannels() const
Returns the number of TPC readout channels in the detector.
std::vector< std::vector< float > > fNoiseZ
noise on each channel for each time for Z (collection) plane
float fCollectionPed
ADC value of baseline for collection plane.
Zero Suppression algorithm.
Definition: RawTypes.h:11
unsigned int fNTimeSamples
number of ADC readout samples in all readout frames (per event)
Planes which measure U.
Definition: geo_types.h:129
int FFTSize() const
Definition: LArFFT.h:69
int fNearestNeighbor
Maximum distance between hits above threshold before they are separated into different blocks...
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
double fSampleRate
sampling rate in ns
float fNoiseFactZ
noise scale factor for Z (collection) plane
float fNoiseWidth
exponential noise width (kHz)
Signal from induction planes.
Definition: geo_types.h:145
double Charge(TDC_t tdc) const
Returns the total number of ionization electrons on this channel in the specified TDC...
Definition: SimChannel.cxx:134
enum geo::_plane_sigtype SigType_t
def move(depos, offset)
Definition: depos.py:107
T get(std::string const &key) const
Definition: ParameterSet.h:271
p
Definition: test.py:223
SimWireDUNE10kt(fhicl::ParameterSet const &pset)
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
pval
Definition: tracks.py:168
float fLowCutoffZ
low frequency filter cutoff (kHz) for Z (collection) plane
auto array(Array const &a)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:228
float fNoiseWidthU
exponential noise width (kHz) for U plane
void Convolute(detinfo::DetectorClocksData const &clockData, Channel channel, std::vector< T > &func) const
std::vector< double > fChargeWork
void GenNoise(std::vector< float > &array)
float fNoiseWidthZ
exponential noise width (kHz) for Z (collection) plane
#define MF_LOG_DEBUG(id)
float fNoiseFact
noise scale factor
std::vector< std::vector< float > > fNoiseU
noise on each channel for each time for U plane
void Compress(std::vector< short > &adc, raw::Compress_t compress)
Compresses a raw data buffer.
Definition: raw.cxx:19
CLHEP::HepRandomEngine & fEngine
raw::Compress_t fCompression
compression type to use
std::size_t getView(std::string const &moduleLabel, std::string const &productInstanceName, std::string const &processName, std::vector< ELEMENT const * > &result) const
Definition: DataViewImpl.h:500
Tools and modules for checking out the basics of the Monte Carlo.
TCEvent evt
Definition: DataStructs.cxx:7
double GetASICGain(Channel channel) const override
unsigned int fNSamplesReadout
number of ADC readout samples in 1 readout frame
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
LArSoft geometry interface.
Definition: ChannelGeo.h:16
unsigned int fZeroThreshold
Zero suppression threshold.
int fNTicks
number of ticks of the clock
float fNoiseFactU
noise scale factor for U plane
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
QTextStream & endl(QTextStream &s)
void produce(art::Event &evt)
Signal from collection planes.
Definition: geo_types.h:146
std::vector< double > DoubleVec
float fLowCutoffV
low frequency filter cutoff (kHz) for V plane