ReadoutSim_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief TODO
3 /// \author Ryan Patterson - rbpatter@caltech.edu
4 ////////////////////////////////////////////////////////////////////////
5 
6 #include <algorithm>
7 #include "sys/stat.h"
8 
9 // Framework includes
18 
19 #include "cetlib/search_path.h"
20 #include "fhiclcpp/ParameterSet.h"
22 
23 #include "CLHEP/Random/RandFlat.h"
27 
28 // NOvA includes
30 #include "Geometry/Geometry.h"
32 #include "MCCheater/BackTracker.h"
33 #include "NovaDAQConventions/DAQConventions.h"
34 #include "RawData/RawTrigger.h"
36 #include "ReadoutSim/NoiseMaker.h"
39 #include "Simulation/Simulation.h"
41 #include "Utilities/func/MathUtil.h"
43 #include "ReadoutSim/PulseShaper.h"
45 
46 // ROOT includes
47 #include "TFile.h"
48 #include "TH2.h"
49 #include "TH3.h"
50 #include "TMath.h"
51 #include "TSystem.h"
52 #include "TGraph.h"
53 #include "TString.h"
54 
55 #include "TH1D.h"
56 
57 #include <set>
58 #include <utility>
59 #include <iomanip>
60 #include <iostream>
61 #include <vector>
62 
63 namespace rsim
64 {
66  {
67  public:
68  explicit ReadoutSim(const fhicl::ParameterSet& pset);
69  virtual ~ReadoutSim();
70 
71  void produce(art::Event& evt);
72  void beginRun(art::Run& run);
73 
74  protected:
75  // internal methods
76  double GetSmearedPE(int nPE);
77  double PE2ADC(double smearedPE, double gain);
78 
79  //double PerfectASICCurve(double t);
80  bool CreateRawDigits(const std::vector<sim::PhotonSignal>& hitlist,
81  std::unique_ptr< std::vector<rawdata::RawDigit> > &rawdigitlist,
82  rawdata::RawTrigger rawtrig);
83 
84  // Draw a trace and save it as a TGraph using the TFileService
85  void DrawASICOutput(double *ASIC_Output, int offset, int stride,
86  int plane, int cell, std::string name="trace_");
87 
88  // configuration settings
89  std::string fPhotonModuleLabel; ///< photons made in the module with this label
90  std::string fCopyNoTruthHitsLabel; ///< hits with no truth (i.e., overlay) from this label will be copied into the digit vector after simulation is done
91 
93 
95 
96  double fScaleEmptyCellNoiseRate; ///< scale factor for empty cell noise
98 
99  std::string fEmptyCellNoiseFile_FD; ///< File holding FD noise PDF
100  std::string fEmptyCellNoiseFile_ND; ///< File holding ND noise PDF
101  std::string fEmptyCellNoiseFile_NDOS; ///< File holding NDOS noise PDF
102  std::string fEmptyCellNoiseFile_TB; ///< File holding TB noise PDF
103 
104  // internal utilities
106  rsim::NoiseMaker* fNoiseMaker; ///< NoiseMaker object for making correlated noise
107  rsim::ExcessNoiseMaker* fExcessNoiseMaker; ///< ExcessNoiseMaker object for modeling
108  rsim::IPulseShaper* fPulseShaper; ///< PulseShaper object for generating traces
109 
110  TH3D* fEmptyCellNoisePDF; ///< for fast generation of empty-cell noise hits
112 
113  std::vector< TH1* > fGainPDF;
114  std::vector< std::vector< float > > fGainMap;
123  std::vector< std::vector< double > > fBaselineMap;
124 
125  int LoadEmptyCellNoisePDF(const char* fname);
126  void PopulateRawTrig(rawdata::RawTrigger &rawtrig, art::Event const &evt);
127  TH1D *fADC;
128  TH1D *fMaxFebADC;
131  TH1D *fNoiseADC;
132 
133  bool fSag;
134  double fSagFrac;
135 
136  /// Vector of plane, cell pairs for which to write traces to the TFile
137  /// for debugging purposes.
138  std::vector<std::pair<int, int> > fTracePlaneCells;
139  };
140 
141  //----------------------------------------------------------------------
143  fPhotonModuleLabel(pset.get<std::string>("PhotonModuleLabel")),
144  fCopyNoTruthHitsLabel(pset.get<std::string>("CopyNoTruthHitsFromLabel")),
145  fPSetNDOS(pset.get<fhicl::ParameterSet>("ndos")),
146  fPSetND (pset.get<fhicl::ParameterSet>("nd")),
147  fPSetFD (pset.get<fhicl::ParameterSet>("fd")),
148  fPSetTB (pset.get<fhicl::ParameterSet>("tb")),
149  fParams(0),
150  fFPGAAlgo(0),
151  fNoiseMaker(0),
153  fPulseShaper(0),
156  fGainPDF(), fUseRandomGains(true),
158  fFEBFlashThreshold(10000.0),
159  fADC(0), fPlaneCellPE(0), fPlaneCellADC(0)
160  {
161 
162  // get the random number seed, use a random default if not specified
163  // in the configuration file.
164  unsigned int seed = pset.get< unsigned int >("Seed", sim::GetRandomNumberSeed());
165 
166  createEngine( seed );
167 
168  produces< std::vector<rawdata::RawDigit> >();
169  produces< std::vector<rawdata::RawTrigger> >();
170 
171  std::vector<int> tracePlanes = pset.get<std::vector<int> >("TracePlanes");
172  std::vector<int> traceCells = pset.get<std::vector<int> >("TraceCells");
173  if(traceCells.size() != tracePlanes.size())
174  throw cet::exception("InvalidTracePlaneCell") <<
175  "FCL configurable TracePlanes and TraceCells must have " <<
176  "same number of elements \n";
177 
178  for(size_t i = 0; i < traceCells.size(); ++i){
179  int plane = tracePlanes[i];
180  int cell = traceCells[i];
181  fTracePlaneCells.push_back(std::make_pair(plane, cell));
182  }
183 
184  }
185 
186  //----------------------------------------------------------------------
188  {
189  if (fParams) delete fParams;
190  if (fFPGAAlgo) delete fFPGAAlgo;
191  if (fNoiseMaker) delete fNoiseMaker;
194  if (fPulseShaper) delete fPulseShaper;
195  }
196 
197  //----------------------------------------------------------------------
198 
200  {
201  // Figure out which parameter set we're using
202  fhicl::ParameterSet pset;
204  switch(geom->DetId()){
205  case novadaq::cnv::kNDOS:
206  pset = fPSetNDOS;
207  break;
209  pset = fPSetND;
210  break;
212  pset = fPSetFD;
213  break;
215  pset = fPSetTB;
216  break;
217  default:
218  assert(0 && "Unknown detector");
219  }
220 
221  fNanosliceVersion = pset.get<int>("FPGA_DCS_NanosliceVersion");
222  fUseRandomGains = pset.get<bool>("UseRandomGains");
223  fVaryNoiseByThreshold = pset.get<bool>("VaryNoiseByThreshold");
224  fVaryBaseline = pset.get<bool>("VaryBaseline");
225  fUseNewExcessNoise = pset.get<bool>("UseNewExcessNoise");
226  fFEBFlashThreshold = pset.get<double>("FEBFlashThreshold");
227  fLookbackSamples = pset.get<int>("FPGA_DCS_LookbackSamples");
228  fUseLegacyShaper = pset.get<bool>("UseLegacyShaper");
229 
230  // Noise histograms
231  fScaleEmptyCellNoiseRate = pset.get<double>("ScaleEmptyCellNoiseRate");
232  fGeneratePhysicsCellNoise = pset.get<bool>("GeneratePhysicsCellNoise");
233 
234  // Create helpers
235  if(!fParams) fParams = new CommonParameters(pset);
237  if(!fNoiseMaker) fNoiseMaker = new rsim::NoiseMaker(pset);
239  if(!fPulseShaper)
240  {
241  if (fUseLegacyShaper) fPulseShaper = new rsim::LegacyPulseShaper(pset);
242  else fPulseShaper = new rsim::PulseShaper(pset);
243  }
244 
245  int runNumber = run.run();
246  if(runNumber >= 1000000 || runNumber < 10000) {
247  int forceGain = pset.get<int>("ForceGain");
248  if (forceGain == 0) {
249  throw cet::exception("NoGainSetting")
250  << "Ideal MC run number detected, but ForceGain == 0. Set ForceGain to the value of gain you want to simulate.\n"
251  << __FILE__ << ":" << __LINE__ << "\n";
252  }
253  }
254 
255  // constructor decides if initialized value is a path or an environment variable
256  cet::search_path sp("FW_SEARCH_PATH");
257 
258  std::string noiseFile;
259  sp.find_file(fParams->fEmptyCellNoiseFile, noiseFile);
260 
261  std::string gainFile;
262  sp.find_file(fParams->fGainFile, gainFile);
263 
264  struct stat sb;
265 
266  if(noiseFile.empty() || stat(noiseFile.c_str(), &sb) != 0){
267  // failed to resolve the file name
268  throw cet::exception("NoCellNoiseFile")
269  << "noise file " << noiseFile << " not found!\n"
270  << __FILE__ << ":" << __LINE__ << "\n";
271  }
272 
273  if(gainFile.empty() || stat(gainFile.c_str(), &sb) != 0){
274  // failed to resolve the file name
275  throw cet::exception("NoGainFile")
276  << "gain file " << gainFile << " not found!\n"
277  << __FILE__ << ":" << __LINE__ << "\n";
278  }
279 
280  if(fEmptyCellNoisePDF == 0){
281  if( LoadEmptyCellNoisePDF(noiseFile.c_str()) ){
282  mf::LogError("ReadoutSimBeginJob") << "failed to load cell noise pdf";
283  }
284  }
285 
286  if (fGainPDF.size() == 0){
287  if (gSystem->AccessPathName(gainFile.c_str())){
288  throw cet::exception("FileNotOpen")
289  << "Unable to open " << gainFile << "\n"
290  << __FILE__ << ":" << __LINE__ << "\n";
291  }
292 
293  TFile f(gainFile.c_str());
294  for (int ipix = 0; ipix < 32; ++ipix)
295  {
296  TString gainName = "hGains_PIX";
297  gainName += ipix;
298  TH1D *htemp = (TH1D*)f.Get(gainName);
299  if (!htemp){
300  throw cet::exception("FatalRootError") << "Unable to get gain histogram from " << gainFile << "\n"
301  << __FILE__ << ":" << __LINE__ << "\n";
302  }
303  gainName = "GainPDF_PIX";
304  gainName += ipix;
305  fGainPDF.push_back( (TH1D*)htemp->Clone(gainName) );
306  fGainPDF.back()->SetDirectory(0);
307  }
308  f.Close();
309 
310  //first time through, fill with nominal gains
311  for (unsigned int plane = 0; plane < geom->NPlanes(); plane++){
312  std::vector< float > plane_gains;
313  for (unsigned int cell = 0; cell < geom->Plane(plane)->Ncells(); cell++){
314  plane_gains.push_back(fParams->fGain);
315  }
316  fGainMap.push_back(plane_gains);
317  }
318  }
319 
320  if (fVaryBaseline) {
321  //Random baseline
322  fBaselineMap.clear();
323  double baselineMean = pset.get<double>("BaselineMean");
324  double baselineSigma = pset.get<double>("BaselineSigma");
326  CLHEP::HepRandomEngine &engine = rng->getEngine();
327  CLHEP::RandGaussQ gauss(engine);
328  for (unsigned int plane = 0; plane < geom->NPlanes(); plane++){
329  std::vector< double > plane_baselines;
330  for (unsigned int cell = 0; cell < geom->Plane(plane)->Ncells(); cell++){
331  double baseline = gauss.fire(baselineMean, baselineSigma);
332  plane_baselines.push_back(baseline);
333  }
334  fBaselineMap.push_back(plane_baselines);
335  }
336  }
337 
339 
340  fSag = pset.get<bool>("Sag");
341  fSagFrac = pset.get<double>("SagFrac");
342 
343  // Histograms
344  if(fPlaneCellPE) return; // Don't make histograms twice
345 
347 
348  fADC = tfs->make<TH1D>("ADC", ";ADC", 4095, 0, 4095);
349  fMaxFebADC = tfs->make<TH1D>("MaxFebADC", ";MaxFebADC", 2000, 0, 20000);
350  fNoiseADC = tfs->make<TH1D>("NoiseADC",";ADC3-ADC0", 150, 0, 300);
351  fPlaneCellPE = tfs->make<TH2D>("planeCellPE", ";Plane;Cell",
352  geom->NPlanes(), 0, geom->NPlanes(),
353  geom->Plane(0)->Ncells(), 0, geom->Plane(0)->Ncells());
354  fPlaneCellADC = tfs->make<TH2D>("planeCellADC", ";Plane;Cell",
355  geom->NPlanes(), 0., geom->NPlanes(),
356  geom->Plane(0)->Ncells(), 0, geom->Plane(0)->Ncells());
357 
358  }
359 
360  //----------------------------------------------------------------------
362  {
363  std::unique_ptr< std::vector<rawdata::RawDigit> > rawdigitlist( new std::vector<rawdata::RawDigit>);
364  std::unique_ptr< std::vector<rawdata::RawTrigger> > rawtriglist( new std::vector<rawdata::RawTrigger>);
365 
367 
368  //needs to be here for CMap to be initialized to get pixel-by-pixel gains
369  //random gains will change event-by-event, but that shouldn't be an issue
370  if (fUseRandomGains)
371  {
372  // get the RandomNumberGenerator service
375 
376  CLHEP::HepRandomEngine &engine = rng->getEngine();
377  RandHisto randHisto(engine);
378  //double gainScale(0.0), dummy1(0.0), dummy2(0.0);
379  for (unsigned int plane = 0; plane < geom->NPlanes(); plane++) {
380  for (unsigned int cell = 0; cell < geom->Plane(plane)->Ncells(); cell++) {
381  const daqchannelmap::lchan lc = cmap->Map()->encodeLChan(geom->DetId(), plane, cell);
382  const daqchannelmap::dchan dc = cmap->Map()->encodeDChan(lc);
383  const int pix = cmap->Map()->getPixel(dc);
384 
385  //randHisto.GetRandom( fGainPDF[pix], gainScale, dummy1, dummy2 );
386  //for now, just use the mean of the distibution for each pixel position
387  double gainScale = fGainPDF[pix]->GetMean(1);
388  fGainMap[plane][cell] = gainScale*fParams->fGain;
389  }
390  }
391  }
392 
393  // Load hit list (PhotonSignal objects)
395  evt.getByLabel(fPhotonModuleLabel, hitlistHandle);
396 
397  std::vector<sim::PhotonSignal> hitlist;
398  for(unsigned int i = 0; i < hitlistHandle->size(); ++i){
399  if((*hitlistHandle)[i].NPhoton() > 0)
400  hitlist.push_back((*hitlistHandle)[i]);
401  }
402 
403  // Need to ensure that hits on the same channel are adjacent in the vector.
404  // The efficiency optimizations in CreateRawDigits assume this to be true
405  // so it's required for correctness.
406  sim::SortByPlaneAndCell(hitlist);
407 
408  // Make a "trigger" happen first. (We'll need the trigger
409  // time when setting hits' TDCs.)
410  // Write the raw trigger into the event.
411  rawdata::RawTrigger rawtrig;
412  PopulateRawTrig(rawtrig, evt);
413  rawtriglist->push_back(rawtrig);
414  evt.put(std::move(rawtriglist));
415 
416  // Create raw digit list from hit list (i.e., do all the work)
417  int failed = CreateRawDigits(hitlist, rawdigitlist, rawtrig);
418  if (failed){
419  mf::LogWarning("FailedDigitCreation") << "failed to create RawDigits, \n"
420  << "putting empty RawDigit list into event";
421  evt.put(std::move(rawdigitlist));
422  return;
423  }
424 
426  for(unsigned int i = 0; i < rawdigitlist->size(); ++i){
427  TString debugOut = cmap->GetPlane(&rawdigitlist->at(i));
428  debugOut += " ";
429  debugOut += cmap->GetCell(&rawdigitlist->at(i));
430  debugOut += " ";
431  for (int iADC = 0; iADC < (int)rawdigitlist->at(i).NADC(); ++ iADC) {
432  debugOut += rawdigitlist->at(i).ADC(iADC);
433  debugOut += " ";
434  }
435  LOG_DEBUG("ReadoutSim") << debugOut;
436 
437  if (fNanosliceVersion == 3) {
438  fADC->Fill(rawdigitlist->at(i).ADC(3) - rawdigitlist->at(i).ADC(0));
439  fPlaneCellADC->Fill(cmap->GetPlane(&rawdigitlist->at(i)), cmap->GetCell(&rawdigitlist->at(i)), rawdigitlist->at(i).ADC(3) - rawdigitlist->at(i).ADC(0) );
440  }
441  }
442 
443  // Ensure all the digits are flagged as MC
444  for(unsigned int i = 0; i < rawdigitlist->size(); ++i)
445  rawdigitlist->at(i).SetMC();
446 
447  // If we want to copy unsimulated hits from a preceding DAQ run, do that now before storing them.
448  // (Primary use case: cosmic overlay hits in light level variations where we're resimulating the
449  // photons in the MC hits but don't want to drop the data cosmic hits.)
450  if (!fCopyNoTruthHitsLabel.empty())
451  {
453  evt.getByLabel(fCopyNoTruthHitsLabel, digitcol);
454 
455  if (digitcol->empty())
456  mf::LogWarning("ReadoutSim") << "Found no digits in fCopyNoTruthHitsLabel label '" << fCopyNoTruthHitsLabel << "' to copy!";
457 
458  for (const auto & srcDigit : *digitcol)
459  {
460  // only want digits NOT from simulation
461  if (!srcDigit.IsMC())
462  rawdigitlist->push_back(srcDigit);
463  }
464  }
465 
466  evt.put(std::move(rawdigitlist));
467 
468  try{
470  // This is the first point in MC production that we have all the products
471  // that BackTracker wants to build maps of. Ensure they're built
472  // correctly for all downstream BackTracker users.
473  bt->Rebuild(evt);
474  }
475  catch(...){
476  // There's no BackTracker service in this job. That's fine, it just
477  // means we don't have to do anything here.
478  }
479 
480  return;
481  }
482 
483 
484  //---------------------------------------------
485  //-------------- CreateRawDigits --------------
486  //---------------------------------------------
488  (const std::vector<sim::PhotonSignal>& hitlist,
489  std::unique_ptr< std::vector<rawdata::RawDigit> > &rawdigitlist,
490  rawdata::RawTrigger rawtrig)
491  {
492  // This will return (via the argument) a vector of RawDigits
493  // that holds the simulated readout for the hits in hitlist.
494  // There will be one RawDigit object per active channel in the
495  // returned vector.
496  //
497  // The steps:
498  // (0) Init (e.g., build template ASIC output curve)
499  // (1) For each cell, lay down scaled/jittered ASIC output curves
500  // on an array of clockticks (all clockticks for simplicity)
501  // (2) Pick out the ones the MUXer will actually pass through
502  // (depends on detector and channel)
503  // (3) Pass these digitizations to the appropriate digital signal
504  // processor. For readibility, also pass a time array even
505  // though it can always be back-computed from the array indices.
506  //
507  // The DSP steps are discussed in the DSP code.
508  //
509 
510  // get the RandomNumberGenerator service
512  CLHEP::HepRandomEngine &engine = rng->getEngine();
514  CLHEP::RandFlat flat(engine);
515  CLHEP::RandGaussQ gauss(engine);
516  CLHEP::RandExponential expo(engine);
517 
518  RandHisto randHisto(engine);
519 
522 
523  // Things we'll need:
524  //const int nHits = hitlist.size();
525 
526  // Not all of the clockticks result in a digitization, thanks to the MUXing.
527  // Which ones get seen by the FPGA depends on the channel and the detector.
528  // Thus, we must pass to the FPGA algorithm two things: (1) the first clocktick
529  // to consider, and (2) the number of clockticks that elapse between relevant
530  // digitizations. (The algorithm may also have other parameters that can be
531  // passed in to create algo variants without duplicating code.)
532 
533  // Loop over hits.
534  // In order to handle all the hits on each channel together, simply
535  // loop twice. The inner loop is to look for more hits on the
536  // same channel. Because the list had been sorted, this code
537  // runs with O(N) complexity.
538 
539  // If APD sag is enabled, create a sagged baseline trace for every FEB that
540  // has physics in it. Indexing here is first by plane and then by "FEB
541  // number", which for these purposes is just cell/32.
542  //key is (plane, cell)
543  std::map< std::pair<int, int>, std::list< std::pair<double, double> > > ADCPulseMap;
544  //key is (plane, feb)
545  std::map< std::pair<int, int>, std::list< std::pair<double, double> > > SaggedPulseMap;
546 
547  std::vector<bool> channelHasHits[geom->NPlanes()];
548  for (unsigned int plane=0;plane<geom->NPlanes();plane++) {
549  for (unsigned int cell=0;cell<geom->Plane(plane)->Ncells();cell++) {
550  channelHasHits[plane].push_back(false);
551  }
552  }
553 
554  for(const sim::PhotonSignal& phot: hitlist){
555  const int plane = phot.Plane();
556  const int cell = phot.Cell();
557 
558  std::pair<int, int> hitCoord(plane, cell);
559  std::pair<int, int> febCoord(plane, cell/32);
560 
561  int nPE = phot.NPhoton();
562  double smearedPE = GetSmearedPE(nPE);
563  double ADC = PE2ADC(smearedPE, fGainMap[plane][cell]);
564  double time = phot.TimeMean();
565 
566  channelHasHits[plane][cell] = true;
567 
568  if (fSag) {
569  //sagged pulses originating from the same channel as the unsagged version will cancel to give the original ADC after clustering
570  SaggedPulseMap[febCoord].push_back( std::make_pair(time, -fSagFrac*ADC) );
571  ADCPulseMap[hitCoord].push_back( std::make_pair(time, (1.0 + fSagFrac)*ADC) );
572  }
573  else ADCPulseMap[hitCoord].push_back( std::make_pair(time, ADC) );
574  }
575 
576  for (auto&& it: ADCPulseMap) it.second.sort();
577  for (auto&& it: SaggedPulseMap) it.second.sort();
578 
579  //cluster the sagged hits within 15 ns to determine if enough charge was collected to trigger an FEB flash
580  std::vector<bool> febHasFlashHit[geom->NPlanes()];
581  for (unsigned int plane=0;plane<geom->NPlanes();plane++) {
582  for (unsigned int feb=0;feb<geom->Plane(plane)->Ncells()/32;feb++) {
583  std::list< std::pair<double,double> > SaggedPulses = SaggedPulseMap[ std::make_pair(plane, feb) ];
584  double meanTime(0.0);
585  double ADC(0.0);
586  double maxADC(0.0);
587  for (auto it = SaggedPulses.begin(); it != SaggedPulses.end(); ++it) {
588  double currTime = it->first;
589  double currADC = fabs(it->second);
590  if (ADC == 0) {
591  meanTime = currTime;
592  ADC = currADC;
593  }
594  else if ( std::fabs(meanTime - currTime) < 15.0 ) {
595  meanTime = meanTime*ADC + currTime*currADC;
596  ADC += currADC;
597  meanTime /= ADC;
598  }
599  else {
600  if (ADC > maxADC) maxADC = ADC;
601  meanTime = currTime;
602  ADC = currADC;
603  }
604  }
605  if (ADC > maxADC) maxADC = ADC;
606  maxADC /= fSagFrac;
607  fMaxFebADC->Fill(maxADC);
608  if (maxADC > fFEBFlashThreshold) febHasFlashHit[plane].push_back(true);
609  else febHasFlashHit[plane].push_back(false);
610  }
611  }
612 
613  // For keeping track of which hits were treated already (by
614  // having had a previous hit on the same channel)...
615  //bool alreadyDone[nHits];
616  //for (int iH=0; iH<nHits; iH++) alreadyDone[iH] = false;
617 
618  // Keeping it simple here. Channels without physics activity need
619  // noise manually applied. That's most of them. Here's how I'll
620  // keep track.
621 
622  // Each FEB, in each plane, has a readout phase set at random. Precalculate
623  // them all here.
624  //
625  // One theoretical problem with this method is that each new spill gets a
626  // new phasing pattern, but I'm hard pressed to think why you'd be doing a
627  // study that could notice this.
628  std::vector<std::vector<int>> phases;
629  phases.resize(geom->NPlanes());
630  for(unsigned int i = 0; i < geom->NPlanes(); ++i){
631  for(unsigned int j = 0; j < geom->Plane(i)->Ncells()/32; ++j){
632  phases[i].push_back(flat.fireInt(fParams->fNumClockticksPerDigitization));
633  }
634  }
635 
636  //will need this for making noise on empty cells
637  if (fEmptyCellNoisePDF == 0) return 1; // return if error
638 
639  //loop over all planes and cells
640  for (unsigned int plane=0;plane<geom->NPlanes();plane++) {
641  for (unsigned int cell=0;cell<geom->Plane(plane)->Ncells();cell++) {
642  const daqchannelmap::lchan lchan = cmap->Map()->encodeLChan(geom->DetId(), plane, cell);
643  const daqchannelmap::dchan dchan = cmap->Map()->encodeDChan(lchan);
644  const int pix = cmap->Map()->getPixel(dchan);
645  // Within the FEB the phase of the digitization proceeds round-robin
646  // around the pixels: 0,8,16,24 together then 1,9,17,25, etc...
647  const int phase = phases[plane][cell/32];
648  int firstDigitizationClocktick = ((pix+phase)%fParams->fNumClockticksPerDigitization);
649 
650  //Set plane and cell so the proper threshold is used
652  //Get threshold for noise in this channel
653  double thresh = fFPGAAlgo->ThresholdADC();
654  //Get the baseline
655  double baseline = fParams->fASICBaselineInADCCounts;
656  if (fVaryBaseline) baseline = fBaselineMap[plane][cell];
657  //Get the noise scaling
658  double noiseFactor = 1.0;
660 
661  //normally only simulate a trace if a channel has a real hit on it
662  //if any channel on the feb has a hit larger than the flasher threshold, simulate a trace for all channels on that feb
663  if ( channelHasHits[plane][cell] || febHasFlashHit[plane][cell/32] ) {
664  std::list< std::pair<double,double> > ADCPulses = ADCPulseMap[ std::make_pair(plane, cell) ];
665  std::list< std::pair<double,double> > SaggedPulses = SaggedPulseMap[ std::make_pair(plane, cell/32) ];
666 
667  fPulseShaper->SetPhase(firstDigitizationClocktick);
668  fPulseShaper->SetBaseline(baseline);
669  //construct the trace
670  // This holds the real-valued "voltages" coming out of the ASIC
671  // (post-shaping). These will be digitized (i.e., integer-ized)
672  // after all the hits on the channel have contributed.
673  std::vector<double> ASIC_Output = fPulseShaper->CreateTrace( ADCPulses, SaggedPulses );
674  // Both "current" and "voltage" noise are being handled by the NoiseMaker
675  // object. See doc-4401 for background. This eliminates the need for
676  // explicit sprinkling of dark noise.
677  //
678  // An annoyance: the noise formalism expects that I will be needing noise
679  // only every digitization cycle, not every clocktick. The correlations
680  // will be wrong if I draw a noise value every clocktick. So, I'll only
681  // put noise on the to-be-used clockticks, which is kind of weird and
682  // will make the DumpASICOutput routine less useful. I can maybe clean
683  // this up later.
684  //
685  // In the same loop, pass the ASIC output through the ADC. This amounts to:
686  // (a) lifting the output by a baseline
687  // (b) flooring the real-valued numbers to make integers
688  // (c) truncating the range to [0, (1<<12)-1]
689  //
690  for (int iClocktick=firstDigitizationClocktick;
691  iClocktick<fParams->fNumClockticksInSpill;
692  iClocktick += fParams->fNumClockticksPerDigitization) {
693  if (fGeneratePhysicsCellNoise) ASIC_Output[iClocktick] += noiseFactor*fNoiseMaker->GetNoise(&gauss)/fParams->fADCMaxPE * ( (1<<12) - 1 );
694  ASIC_Output[iClocktick] = int(ASIC_Output[iClocktick]);
695  if ( ASIC_Output[iClocktick] < 0 ) ASIC_Output[iClocktick] = 0 ;
696  else if ( ASIC_Output[iClocktick] > (1<<12)-1 ) ASIC_Output[iClocktick] = (1<<12)-1;
697  } //end iClocktick
698 
699  // Draw traces if requested in fcl
700  if(!fTracePlaneCells.empty()) {
701  for(auto tracePlaneCell:fTracePlaneCells) {
702  if(plane != (unsigned int)tracePlaneCell.first) continue;
703  // Draw everything in the same FEB, important in the sag case
704  if(cell/32 != (unsigned int)tracePlaneCell.second/32) continue;
705 
706  DrawASICOutput(&ASIC_Output[0], firstDigitizationClocktick,
708  plane, cell, "trace_");
709 
710  } // for tracePlaneCell
711  }
712 
713  // Process the now-complete signal through the FPGA.
714  //
715  // The list of FPGA algorithms can grow. Add new ones by following the
716  // examples already present.
717  //
718  // Different detectors may have different signal processing
719  // algorithms. At the moment, though, all the detectors are controlled by
720  // a single string parameter
721 
722  // Call the requested algorithm to turn the ASIC output into a list of
723  // RawDigits to be added to the outgoing list
724  std::vector<rawdata::RawDigit> local_digitlist =
725  fFPGAAlgo->ASICToDigits(&ASIC_Output[0], firstDigitizationClocktick);
726  LOG_DEBUG("ReadoutSim") << " == size of digilist == "
727  << local_digitlist.size();
728 
729  // Store any rawdigits (if there are actually triggered hits on the
730  // channel)
731  for (unsigned int iD=0; iD<local_digitlist.size(); iD++) {
732  rawdata::RawDigit digit = local_digitlist[iD];
733 
734  if ( ! digit.NADC() ) {
735  throw cet::exception("NoADC")
736  << "ReadoutSim::CreateRawDigits() -- FPGA algorithm returned "
737  << "a digit with no ADC values. This should never happen. Fatal!\n"
738  << __FILE__ << ":" << __LINE__ << "\n";
739  return 1;
740  }
741  // Finalize the digit
742  digit.SetChannel(lchan);
743  digit.SetDaqChannel(dchan);
744 
745  // Store the digit!
746  LOG_DEBUG("ReadoutSim") << " -- digit -- " << plane << " " << cell
747  << " " << digit.TDC() << " " << digit.ADC(0);
748 
749  rawdigitlist->push_back(digit);
750  } // end iD
751 
752  } //end for channels with hits
753  else {
754  if( fScaleEmptyCellNoiseRate == 0) continue;
755  double rate = fScaleEmptyCellNoiseRate*fEmptyCellNoisePDFIntegral; // prob. per cell per ns
756  double mean = rate*fParams->fClocktick*fParams->fNumClockticksPerDigitization; //prob. per cell per digitization
757  double mean_arrival = 0;
758  if (mean != 0)
759  mean_arrival = 1.0/mean;
760 
761  //start at the offset first digitization clocktick
762  for (Long64_t iTick = firstDigitizationClocktick;
765  //std::floor( expo.fire(mean_arrival) ) is the number of digitization periods until a noise hit
766  //It's rounded down since any hit within a period should be assigned to the first clocktick of the period
767  iTick += std::floor( expo.fire(mean_arrival) )*fParams->fNumClockticksPerDigitization;
768  if (iTick < fParams->fNumClockticksInSpill) {
769  // random adc values
770  double adc0, adc1, adc2;
771  randHisto.GetRandom( fEmptyCellNoisePDF, adc0, adc1, adc2 );
772  // make the noise RawDigit
773  rawdata::RawDigit digit;
775 
776  //need to store peak time, not baseline time
777  digit.SetTDC((iTick+fLookbackSamples*fParams->fNumClockticksPerDigitization)<<2); // shift by two bits to move 16 MHz-->64 MHz
778  //digit.SetTDC(iTick<<2); // shift by two bits to move 16 MHz-->64 MHz
779  // GetRandom will have given us a value between binLowEdge and
780  // binHighEdge. All that bin's contents represent binLowEdge, so
781  // rounding down is correct.
782  if (fNanosliceVersion == 0) {
783  if (fParams->fUseNewEmptyCellNoise) digit.SetADC(0,uint16_t(adc2));
784  else digit.SetADC(0,uint16_t(adc0));
785  }
786  else {
787  digit.SetADC(0,uint16_t(baseline));
788  digit.SetADC(1,uint16_t(adc0 + baseline));
789  digit.SetADC(2,uint16_t(adc1 + baseline));
790  digit.SetADC(3,uint16_t(adc2 + baseline));
791  }
792  digit.SetChannel(lchan);
793  digit.SetDaqChannel(dchan);
794 
795  // Store the noise digit!
796  // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
797  if (digit.ADC() >= thresh) {
798  fNoiseADC->Fill(digit.ADC());
799  rawdigitlist->push_back(digit);
800  }
801  }
802  }
803  } //end for channels without hits
804  } //end for cell
805  } //end for plane
806  return 0;
807  }
808 
809  //---------------------------------------------
810  //-------------- PE2ADC --------------
811  //---------------------------------------------
812  double ReadoutSim::PE2ADC(double smearedPE, double gain)
813  {
814  //Converts photoelectrons to ADC
815 
816  return smearedPE / fParams->fADCMaxPE * ( (1<<12) - 1 ) * (gain/fParams->fGain);
817  }
818 
819  //------------------------------------------------
820  //-------------- GetSmearedPE------ --------------
821  //------------------------------------------------
822  double ReadoutSim::GetSmearedPE(int nPE)
823  {
824  // Smear PEs to account for excess noise such that the mean PE remains the same
825  // and the spread is sqrt(fAPDExcessNoiseFactor*nPE)
826 
827  // For smaller numbers of PE, use the ExcessNoiseMaker
828  // For larger numbers, approximate with the log-normal distribution
829 
830  // get the RandomNumberGenerator service
832  CLHEP::HepRandomEngine &engine = rng->getEngine();
833  CLHEP::RandFlat flat(engine);
834  CLHEP::RandGaussQ gauss(engine);
835 
836  if (fUseNewExcessNoise && nPE < 250) return fExcessNoiseMaker->DrawSmearedPE(nPE, &flat);
837  else {
838  if ( nPE == 0 ) {
839  throw cet::exception("NPEEq0")
840  << "ReadoutSim::GetASICScaleSpread() -- nPE==0. This should never happen.\n"
841  << __FILE__ << ":" << __LINE__ << "\n";
842  return 1;
843  }
844 
845  // This is the variance we need, to be combined with the preexisting poisson variance
846  const double var = (fParams->fAPDExcessNoiseFactor-1)/nPE;
847  // These paramters define a log-normal distribution with mean 1 and
848  // variance given by 'var'
849  const double mu = -log(1+var)/2;
850  const double sigma = sqrt(log(1+var));
851 
852  // Generate a log-normal number with parameters mu and sigma
853  return nPE*exp(mu+sigma*gauss.fire(0, 1));
854  }
855  }
856 
857  //---------------------------------------------
858  //----------- LoadEmptyCellNoisePDF -----------
859  //---------------------------------------------
861 
862  // Load PDF for generating noise on physics-less cell.
863  // I'd love for these files to stay as root files. Please
864  // don't change that unless absolutely necessary.
865 
866  if (gSystem->AccessPathName(fname)) {
867  cet::exception("FileNotOpen") << "ReadoutSim: Unable to open " << fname
868  << " Cannot get noise PDFs!";
869  return 1;
870  }
871 
872  TFile f(fname);
873 
874  TH3D *htemp = (TH3D*)f.Get("noise/h3_adc0adc1adc2");
875  if (!htemp) {
876  std::cerr << "ReadoutSim: Histogram h3_adc0adc1adc2 not found in " << fname << std::endl;
877  std::cerr << " Cannot get fine timing! " << std::endl;
878  return 1;
879  }
880  fEmptyCellNoisePDF = (TH3D*)htemp->Clone();
881  fEmptyCellNoisePDF->SetDirectory(0);
882 
883  f.Close();
884 
885  // Apply threshold to noise histogram so that it is consistent with the
886  // parameters we're currently simulating.
887  //const double thresh = fFPGAAlgo->ThresholdADC();
888  //for(int i = 0; i < fEmptyCellNoisePDF->GetNbinsX()+2; ++i){
889  //for(int j = 0; j < fEmptyCellNoisePDF->GetNbinsY()+2; ++j){
890  //for(int k = 0; k < fEmptyCellNoisePDF->GetNbinsZ()+2; ++k){
891  // if(fEmptyCellNoisePDF->GetXaxis()->GetBinCenter(i) < thresh)
892  // fEmptyCellNoisePDF->SetBinContent(i, j, k, 0);
893  //}
894  //}
895  //}
896 
898 
899  return 0;
900  }
901 
902 
904  // We may need to flesh out more of the trigger block in the
905  // future (although not all of the items make sense in MC).
906 
907  // Say this is MC:
908  rawtrig.fTriggerMask_MCBit = 1;
909 
910  // Units of this in data are 500ns
912 
913  // Set the trigger time such that it increments by 0.1 seconds
914  // every event (to mimic calibration triggers, though not beam
915  // triggers). Also, add one "day" per subrun and one "year" per
916  // run, just to keep the times in different files different.
917  // I'm not making these FHICL parameters because we will want
918  // to do this properly (e.g., have each MC file representing a
919  // specific timestamp drawn to match real data).
920  //
921  // The TDC ticks run at 64 MHz
922 
923  // using unsigned long long rather than uint64_t to ensure consistent behavior on OSX and LINUX
924  static const unsigned long long ticksPerYear = (unsigned long long)(64000000) * 31556926;
925  static const unsigned long long ticksPerDay = (unsigned long long)(64000000) * 86400;
926  static const unsigned long long ticksPerEvent = (unsigned long long)(64000000) / 10; // 10 Hz trigger
927 
929  evt.run() * ticksPerYear +
930  evt.subRun() * ticksPerDay +
931  evt.event() * ticksPerEvent;
932 
933  }
934 
935  //---------------------------------------------
936  //----------- DrawASICOutput -----------
937  //---------------------------------------------
938  void ReadoutSim::DrawASICOutput(double *ASIC_Output, int offset, int stride,
939  int plane, int cell, std::string name){
940 
941 
942  // Make this function sort of noisy so that it doesn't get accidentally
943  // get called in production or when folks don't expect it is being called.
944  std::cout << "Saving a trace for plane: " << plane
945  << " cell: " << cell << std::endl;
946 
947  std::vector<double> traceTimes;
948  std::vector<double> trace;
949 
950  for(int iClocktick = offset;
951  iClocktick < fParams->fNumClockticksInSpill;
952  iClocktick += stride)
953  {
954 
955  traceTimes.push_back(iClocktick * fParams->fClocktick);
956  trace.push_back(ASIC_Output[iClocktick]);
957  }
958 
959 
961  TGraph* traceGraph = tfs->make<TGraph>(traceTimes.size(),
962  &traceTimes[0], &trace[0]);
963  std::stringstream traceName;
964  traceName<<name<<"_plane" <<plane<<"_feb"<<cell/32<< "_cell"<<cell;
965  traceGraph->SetName(traceName.str().c_str());
966  traceGraph->SetTitle(";Time [ns]; ADC");
967  traceGraph->Write();
968  delete traceGraph;
969  }
970 
971 
972 
974 
975 }//namespace
double fAPDExcessNoiseFactor
APD&#39;s "excess noise factor" (see .cxx for more)
void SortByPlaneAndCell(std::vector< sim::PhotonSignal > &c)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
const XML_Char * name
Definition: expat.h:151
int fNumClockticksPerDigitization
How many ticks between each ADC.
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
double DrawSmearedPE(int nPE, CLHEP::RandFlat *flat)
virtual std::vector< double > CreateTrace(std::list< std::pair< double, double > > &ADCPulses, std::list< std::pair< double, double > > &SaggedPulses)=0
virtual std::vector< rawdata::RawDigit > ASICToDigits(double *ASIC_Output, int offset)=0
set< int >::iterator it
int32_t TDC() const
The time of the last baseline sample.
Definition: RawDigit.h:94
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::vector< TH1 * > fGainPDF
void beginRun(art::Run &run)
double GetNominalThreshold()
Definition: NoiseMaker.h:37
void SetTDC(int32_t iTDC)
Definition: RawDigit.cxx:116
ReadoutSim(const fhicl::ParameterSet &pset)
rsim::IPulseShaper * fPulseShaper
PulseShaper object for generating traces.
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
TH3D * fEmptyCellNoisePDF
for fast generation of empty-cell noise hits
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::string fEmptyCellNoiseFile_NDOS
File holding NDOS noise PDF.
void SetBaseline(double Baseline)
Definition: IPulseShaper.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
base_engine_t & createEngine(seed_t seed)
void SetChannel(uint32_t iChan)
Definition: RawDigit.h:99
OStream cerr
Definition: OStream.cxx:7
const daqchannelmap::DAQChannelMap * Map() const
Definition: CMap.h:57
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
IFPGAAlgorithm * CreateFPGAAlgorithm(const fhicl::ParameterSet &pset)
Create and configure the correct algorithm based on pset.
void PopulateRawTrig(rawdata::RawTrigger &rawtrig, art::Event const &evt)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
std::string fGainFile
File containing the distribution of gains on the detector.
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
RunNumber_t run() const
Definition: Run.h:47
std::string find_file(std::string const &filename) const
lchan encodeLChan(int detId, plane_t plane, cell_t cell) const
std::string fEmptyCellNoiseFile_FD
File holding FD noise PDF.
rsim::ExcessNoiseMaker * fExcessNoiseMaker
ExcessNoiseMaker object for modeling.
std::string fPhotonModuleLabel
photons made in the module with this label
std::string fEmptyCellNoiseFile_TB
File holding TB noise PDF.
Definition: Run.h:31
void SetVersion(uint8_t v)
Definition: RawDigit.h:101
rsim::IFPGAAlgorithm * fFPGAAlgo
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double fScaleEmptyCellNoiseRate
scale factor for empty cell noise
Common configuration params for SimpleReadout, FPGAAlgorithms, NoiseMaker.
int LoadEmptyCellNoisePDF(const char *fname)
void SetADC(uint32_t i, int16_t iADC)
Definition: RawDigit.cxx:108
base_engine_t & getEngine() const
double GetNoise(CLHEP::RandGaussQ *gauss)
Definition: NoiseMaker.h:32
void SetPlaneAndCell(int plane, int cell)
Far Detector at Ash River, MN.
if(dump)
fhicl::ParameterSet fPSetND
unsigned int seed
Definition: runWimpSim.h:102
std::string fCopyNoTruthHitsLabel
hits with no truth (i.e., overlay) from this label will be copied into the digit vector after simulat...
bool fUseNewEmptyCellNoise
For backwards compatibility with the old noise model until the 2nd analysis is over.
Prototype Near Detector on the surface at FNAL.
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
int evt
Interface implemented by all FPGA algorithms.
Near Detector in the NuMI cavern.
std::vector< std::vector< float > > fGainMap
void SetPhase(int Phase)
Definition: IPulseShaper.h:27
unsigned short GetPlane(const rawdata::RawDigit *dig)
Definition: CMap.cxx:285
double fGain
APD gain (electrons per photoelectron)
void SetDaqChannel(uint32_t iChan)
Definition: RawDigit.h:100
bool CreateRawDigits(const std::vector< sim::PhotonSignal > &hitlist, std::unique_ptr< std::vector< rawdata::RawDigit > > &rawdigitlist, rawdata::RawTrigger rawtrig)
void produce(art::Event &evt)
const double j
Definition: BetheBloch.cxx:29
fhicl::ParameterSet fPSetTB
double ThresholdADC() const
EventNumber_t event() const
Definition: Event.h:67
std::vector< std::vector< double > > fBaselineMap
fhicl::ParameterSet fPSetNDOS
static constexpr Double_t gauss
Definition: Munits.h:360
rsim::NoiseMaker * fNoiseMaker
NoiseMaker object for making correlated noise.
double sigma(TH1F *hist, double percentile)
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
void Rebuild(const art::Event &evt)
double fADCMaxPE
maximum signal size that the ADC can handle without saturating (in PE)
std::vector< std::pair< int, int > > fTracePlaneCells
unsigned int GetRandomNumberSeed()
Definition: Simulation.cxx:13
pixel_t getPixel(dchan daqchan) const
Decode the pixel id from a dchan.
fhicl::ParameterSet fPSetFD
T * make(ARGS...args) const
std::string fEmptyCellNoiseFile_ND
File holding ND noise PDF.
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
uint32_t NADC() const
Definition: RawDigit.h:81
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
unsigned long long fTriggerTimingMarker_TimeStart
Definition: RawTrigger.h:38
def lc(target="")
Definition: g4zmq.py:63
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Can be used as either a member holding configurations, or a mix-in.
void geom(int which=0)
Definition: geom.C:163
T trace(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Definition: trace.hpp:19
int fNumClockticksInSpill
number of clockticks in a spill
double GetSmearedPE(int nPE)
assert(nhit_max >=nhit_nbins)
unsigned short GetCell(const rawdata::RawDigit *dig)
Definition: CMap.cxx:327
void GetRandom(TH1 *histo, double &x, double &y, double &z)
Definition: RandHisto.cxx:8
uint32_t fTriggerRange_TriggerLength
Definition: RawTrigger.h:40
dchan encodeDChan(int detID, diblock_t diblock, dcm_id_t dcm, feb_t feb, pixel_t pixel) const
CommonParameters * fParams
unsigned int NPlanes() const
double PE2ADC(double smearedPE, double gain)
uint32_t dchan
< DAQ Channel Map Package
RunNumber_t run() const
Definition: Event.h:77
std::string fEmptyCellNoiseFile
File containing the empty cell noise templates.
Encapsulate the geometry of one entire detector (near, far, ndos)
double fASICBaselineInADCCounts
ASIC baseline output level, in ADC counts.
long fireInt(long n)
double fClocktick
digitization clock period, ns
void DrawASICOutput(double *ASIC_Output, int offset, int stride, int plane, int cell, std::string name="trace_")
enum BeamMode string