ImprovedTransport_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 /// \brief This slightly less simple photon transport uses a template for photon collection in
3 /// time and in distance along the fiber.
4 ///
5 /// \author aurisano@fnal.gov
6 /// \date
7 ////////////////////////////////////////////////////////////////////////
8 //Framework Includes
14 #include "fhiclcpp/ParameterSet.h"
18 #include "cetlib/search_path.h"
22 
23 //NOvA includes
24 #include "Geometry/Geometry.h"
26 #include "Simulation/FLSHit.h"
27 #include "Simulation/FLSHitList.h"
29 #include "Simulation/Simulation.h"
30 #include "Utilities/func/MathUtil.h"
34 
35 //stdlib Includes
36 #include <string>
37 #include <list>
38 #include <vector>
39 #include <cmath>
40 
41 //ROOT Includes
42 #include "TH1D.h"
43 #include "TH2D.h"
44 #include "TProfile.h"
45 #include "TFile.h"
46 #include "TMath.h"
47 #include "TRandom3.h"
48 
49 //System Includes
50 #include <sys/stat.h>
51 
52 ///photon transport simulation
53 namespace photrans
54 {
55  /// a struct to hold collection rate info
56  struct RateInfo
57  {
58  double rate;
59  double z;
60  double zWidth;
61  TH1D timePDF;
62  };
63 
64  /// a class for transporting photons in a roughly realistic way
66  {
67  public:
68  explicit ImprovedTransport(fhicl::ParameterSet const &pset);
69  virtual ~ImprovedTransport();
70 
71  void beginRun(art::Run& run);
72  void produce(art::Event& evt);
73 
74  protected:
75  double FiberTransmission(double x);
76  double ScintTime(double dEdx);
77  bool LoadHistos();
78  void StepAlongHit(sim::FLSHit& hit);
79  /// \param meanPE The photons entering the fiber
80  /// \param stepTime When they entered
81  /// \param dist How far to the readout
82  /// \param lambda Accumulates total mean expected at readout into here
83  void CreatePhotoElectrons(double meanPE, double stepTime, double dist, double dEdx, bool isScint, RateInfo& info, double& lambda);
84  void ClusterPhotoElectrons(int planeId, int cellId, int trackId);
85 
86  /// Get correction for the varying efficiency across the cross
87  /// section of of the cell caused by fiber position.
88  double GetPosCorr( const sim::FLSHit& hit );
89 
90  private:
91  //fcl parameters
92  std::string fFilterModuleLabel; ///< module that filtered the fls hits
93  fhicl::ParameterSet fPSetNDOS, fPSetND, fPSetFD, fPSetTB;
94 
95  double fPhotonsPerMeV; ///< blue photon production scale
96  double fXViewFactor; ///< scale factor for the x-view
97  double fYViewFactor; ///< scale factor for the y-view
98  double fLightLevelScale; ///< shift overall light yield by multiplicative factor
99  double fCerenkovEff; ///< Absorption and re-emission efficiency for Cerenkov photons
100  double fQuantumEff; ///< APD quantum efficiency
101  double fFiberIndexOfRefraction; ///< n for fiber core
102  double fEmissionTau; ///< exponential time constant for green photon emission
103  std::vector<double> fAttenPars; ///< from Leon's fit
104  double fTimeClustering; ///< groups photons closer than this many ns apart
105  bool fWriteEmptySignals; ///< Write an empty PhotonSignal if there's 0 PE
106  double fStep; ///< step size for walking along the track (cm)
107  bool fApplyBirks; ///< Use EdepBirks from FLSHits instead of calculating here
108  std::string fCollectionRateFile; ///< Relative path to the collection rate histo file
109  bool fApplyFiberSmearing; ///< Turn on fiber smearing
110  std::string fSmearingHistoFile; ///< Relative path to the fiber smearing histo file
111  std::string fBrightnessMode; ///< correct the fiber brightness by cell, by module, or not at all
112  std::string fBrightnessFile; ///< Relative path to the brightness correction file
113  bool fApplyTweak; ///< Turn on the data driven attenuation correction for each view
114  std::string fTweakFile; ///< Relative path to the data driven attenuation correction for each view
115  bool fSimulateDrift; ///< whether to vary light level as a function of time
116  int fDriftReference; ///< drift reference time (ie. where the slope is fixed to 1)
117  double fDriftGradient; ///< gradient for light level detector aging systematic [fraction/year]
118 
119  bool fUseScintTime; ///< Turn on sampling from PDF of scintillator decay time, see arXiv:1704.02291 for more details
120  double fCherenkovTau; ///< Time constant for directly excited PC to PPO scintillation
121  double fScintBetaDEDX; ///< Approximate dE/dx of betas used to determine PDF
122  std::vector<double> fScintBetaProb; ///< Cumulative probability for each exponential component for betas
123  std::vector<double> fScintBetaTau; ///< Time constants of each exponential component for betas
124  double fScintAlphaDEDX; ///< Approximate dE/dx of alphas used to determine PDF
125  std::vector<double> fScintAlphaProb; ///< Cumulative probability for each exponential component for alphas
126  std::vector<double> fScintAlphaTau; ///< Time constants of each exponential component for alphas
127 
128  /// Relative path to the cell efficiency maps. If left empty, do not
129  /// apply this correction.
131 
132  int msglevel;
133 
134  //non-fcl members
135  double fFiberSpeedOfLight; ///< Speed of light in the fiber
136  TH2D fCollectionRateHisto; ///< Collection rate for given dT vs. dZ
137  TH1D fFiberSmearingHisto; ///< Smearing histo for photons bouncing around in fiber
138  TH1D fTweak_X; ///< A data/mc correction for the X-View
139  TH1D fTweak_Y; ///< A data/mc correction for the Y-View
140  TH1D fTweak_XMC; ///< A data/mc correction for the X-View in the muon catcher
141  TH1D fTweak_YMC; ///< A data/mc correction for the Y-View in the muon catcher
142 
143  std::vector< std::vector< double > > fBrightnessMap; ///< The brightness correction by cell
144  //Add by Pengfei -- begin
149  //Add by Pengfei -- end
150 
151  // Efficiency maps, read from fPosCorrFile and applied in GetPosCorr().
152  //
153  // First the horizontals, middle 14 cells of each extrusion, which are
154  // the bigger cells, with the fibers randomly on the bottom, and with
155  // the fiber together on the bottom.
156  std::vector<TH2D *> fBigHorizMap;
157  std::vector<TH2D *> fBigHorizClumpMap;
158 
159  // Outer two cells of horizontals, the small ones.
160  std::vector<TH2D *> fSmallHorizMap;
161  std::vector<TH2D *> fSmallHorizClumpMap;
162 
163  // Same for the top cell in modules with air bubbles.
164  std::vector<TH2D *> fSmallHorizBubbleMap;
165  std::vector<TH2D *> fSmallHorizBubbleClumpMap;
166 
167  // Now the verticals, with fibers randomly distribted, and with fibers
168  // in the corners towards the bottom of cells, because of slack.
169  std::vector<TH2D *> fBigVertMap;
170  std::vector<TH2D *> fBigVertCornerMap;
171 
172  // Same for the small verticals.
173  std::vector<TH2D *> fSmallVertMap;
174  std::vector<TH2D *> fSmallVertCornerMap;
175 
176  // Index for the map being used now. Changes for each hit.
177  // This is for saving in the output so we know what map we used.
178  // We certainly aren't going to use more than 65535 maps, and we
179  // want to save this number for every hit, so use a 16 bit integer.
180  uint16_t fPosCorrMapIndex = -1;
181 
182  // Compute the global index of the fiber efficiency map we used
183  // and store it in fPosCorrMapIndex.
184  void SetPosCorrMapIndex(const std::vector<TH2D *> * const maps,
185  const unsigned int i);
186 
187  /// Reads in the efficiency maps from fPosCorrFile.
188  void LoadPosCorrHists();
189 
190  // Used to choose efficiency cell maps
191  TRandom3 fPosCorrHasher;
192  uint64_t fPosCorrEvtId;
193 
194  art::ServiceHandle<geo::Geometry> fGeo; ///< Handle to the geometry service
195  std::list< double > fPETimes; ///< List of photo-electron times (each representing one observed PE)
196  std::list< RateInfo > fRateInfo; ///< List of RateInfo object
197 
198  std::unique_ptr<std::vector<sim::PhotonSignal> > fSignals; ///< Produced vector of PhotonSignals
205 
206  // Histo for testing
207  TH1D* fScintTime;
208  };
209 
210  //............................................................
212  : fFilterModuleLabel (pset.get< std::string>("FilterModuleLabel"))
213  , fPSetNDOS (pset.get<fhicl::ParameterSet>("ndos"))
214  , fPSetND (pset.get<fhicl::ParameterSet>("nd"))
215  , fPSetFD (pset.get<fhicl::ParameterSet>("fd"))
216  , fPSetTB (pset.get<fhicl::ParameterSet>("tb"))
217  , fSignals (new std::vector<sim::PhotonSignal>)
218  , fFirstTime(true)
219  {
220  // get the random number seed, use a random default if not specified
221  // in the configuration file.
222  unsigned int seed = pset.get< unsigned int >("Seed", sim::GetRandomNumberSeed());
223 
224  createEngine( seed );
225 
226  // get the RandomNumberGenerator service
228  CLHEP::HepRandomEngine &engine = rng->getEngine();
229  fPoisson = new CLHEP::RandPoissonQ(engine);
230  fFlat = new CLHEP::RandFlat(engine);
231  fRandHisto = new RandHisto(engine);
232  fGauss = new CLHEP::RandGaussQ(engine);
233  fExp = new CLHEP::RandExponential(engine);
234 
235  produces< std::vector<sim::PhotonSignal> >();
236  //Add by Pengfei -- begin
237  produces<BrightnessLevel, art::InRun> ();
238  //Add by Pengfei -- end
239  }
240 
241  //............................................................
243  {
244  delete fPoisson;
245  delete fFlat;
246  delete fRandHisto;
247  delete fGauss;
248  }
249 
251  {
252  if (fFirstTime) {
254  fScintTime = tfs->make<TH1D>("hScintTime", ";time;photons", 1000, 0, 5000);
255 
256  // Figure out which parameter set we're using
257  fhicl::ParameterSet pset;
259  switch(geom->DetId()){
260  case novadaq::cnv::kNDOS:
261  pset = fPSetNDOS;
262  break;
264  pset = fPSetND;
265  break;
267  pset = fPSetFD;
268  break;
270  pset = fPSetTB;
271  break;
272  default:
273  assert(0 && "Unknown detector");
274  }
275 
276  fPhotonsPerMeV = pset.get< double >("PhotonsPerMeV");
277  fXViewFactor = pset.get< double >("XViewFactor");
278  fYViewFactor = pset.get< double >("YViewFactor");
279  fLightLevelScale = pset.get< double >("LightLevelScale");
280  fCerenkovEff = pset.get< double >("CerenkovEff");
281  fQuantumEff = pset.get< double >("QuantumEff");
282  fFiberIndexOfRefraction = pset.get< double >("FiberIndexOfRefraction");
283  fEmissionTau = pset.get< double >("EmissionTau");
284  fAttenPars = pset.get< std::vector<double> >("AttenPars");
285  fTimeClustering = pset.get< double >("TimeClustering");
286  fWriteEmptySignals = pset.get< bool >("WriteEmptySignals");
287  fStep = pset.get< double >("Step");
288  fApplyBirks = pset.get< bool >("ApplyBirks");
289  fCollectionRateFile = pset.get< std::string >("CollectionRateFile");
290  fSmearingHistoFile = pset.get< std::string >("SmearingHistoFile");
291  fBrightnessFile = pset.get< std::string >("BrightnessFile");
292  fBrightnessMode = pset.get< std::string >("BrightnessMode");
293  fUseScintTime = pset.get< bool >("UseScintTime");
294  fCherenkovTau = pset.get< double >("CherenkovTau");
295  fScintBetaDEDX = pset.get< double >("ScintBetaDEDX");
296  fScintBetaProb = pset.get< std::vector<double> >("ScintBetaProb");
297  fScintBetaTau = pset.get< std::vector<double> >("ScintBetaTau");
298  fScintAlphaDEDX = pset.get< double >("ScintAlphaDEDX");
299  fScintAlphaProb = pset.get< std::vector<double> >("ScintAlphaProb");
300  fScintAlphaTau = pset.get< std::vector<double> >("ScintAlphaTau");
301  //Add by Pengfei -- begin
302  fBrightnessMapName = pset.get< std::string >("BrightnessMapName");
303  fBrightnessValueName = pset.get< std::string >("BrightnessValueName");
304  fBrightnessValueMapName = pset.get< std::string >("BrightnessValueMapName");
305  //Add by Pengfei -- end
306  fApplyTweak = pset.get< bool >("ApplyTweak");
307  fTweakFile = pset.get< std::string >("TweakFile");
308  fSimulateDrift = pset.get< bool >("SimulateDrift");
309  fDriftReference = pset.get< int >("DriftReference");
310  fDriftGradient = pset.get< double >("DriftGradient");
311  fPosCorrFile = pset.get< std::string >("PosCorrFile");
312  fApplyFiberSmearing = pset.get< bool >("ApplyFiberSmearing");
313  msglevel = pset.get< int >("MessageLevel");
314 
315  // constructor decides if initialized value is a path or an environment variable
316  cet::search_path sp("FW_SEARCH_PATH");
317 
318  const double c = TMath::Ccgs()*1e-9; // c in cm/ns
320 
321  sp.find_file(pset.get< std::string >("CollectionRateFile"), fCollectionRateFile);
322  // Check the collection rate file
323  struct stat sb;
324  if ( fCollectionRateFile.empty() || stat(fCollectionRateFile.c_str(), &sb)!=0 ){
325  // failed to resolve the file name
326  throw cet::exception("NoCollectionRateFile") << "collection rate file " << fCollectionRateFile << " not found!\n" << __FILE__ << ":" << __LINE__ << "\n";
327  }
328 
329  sp.find_file(pset.get< std::string >("SmearingHistoFile"), fSmearingHistoFile);
330  if ( fSmearingHistoFile.empty() || stat(fSmearingHistoFile.c_str(), &sb)!=0 ){
331  // failed to resolve the file name
332  throw cet::exception("NoSmearingHistoFile") << "smearing histo file " << fSmearingHistoFile << " not found!\n" << __FILE__ << ":" << __LINE__ << "\n";
333  }
334 
335  if (fBrightnessMode != "None") {
336  fBrightnessFile = pset.get< std::string >("BrightnessFile");
337  sp.find_file(pset.get< std::string >("BrightnessFile"), fBrightnessFile);
338  if ( fBrightnessFile.empty() || stat(fBrightnessFile.c_str(), &sb)!=0 ){
339  // failed to resolve the file name
340  throw cet::exception("NoBrightnessFile") << "fiber brightness histo file " << fBrightnessFile << " not found!\n" << __FILE__ << ":" << __LINE__ << "\n";
341  }
342  }
343 
344  if (fApplyTweak) {
345  sp.find_file(pset.get< std::string >("TweakFile"), fTweakFile);
346  if ( fTweakFile.empty() || stat(fTweakFile.c_str(), &sb)!=0 ) {
347  // failed to resolve the file name
348  throw cet::exception("NoTweakFile") << "tweak file " << fTweakFile << " not found!\n" << __FILE__ << ":" << __LINE__ << "\n";
349  }
350  }
351 
352  sp.find_file(pset.get< std::string >("PosCorrFile"), fPosCorrFile);
353  if(fPosCorrFile.empty()){
354  if(msglevel >= 10)
355  std::cout << "Fiber position efficiency correction disabled\n";
356  }
357  else{
358  if(stat(fPosCorrFile.c_str(), &sb) != 0)
359  throw cet::exception("NoPosCorrFile") << "Efficiency map file "
360  << fPosCorrFile << " not found!\n" << __FILE__ << ":" << __LINE__
361  << "\n";
362  }
363 
364  if ( !LoadHistos() )
365  {
366  throw cet::exception("NoCollectionRateHisto")
367  << "can't load histos\n"
368  << __FILE__ << ":" << __LINE__ << "\n";
369  return;
370  }
371 
372  //Add by Pengfei -- begin
373  if (fBrightnessMode == "ByCell"){
374  std::unique_ptr<BrightnessLevel> pfb(new BrightnessLevel(fBrightnessFile, fBrightnessMapName, fBrightnessValueName,fBrightnessValueMapName));
375  run.put(std::move(pfb));
376  }
377  else if (fBrightnessMode == "ByBin"){
378  std::unique_ptr<BrightnessLevel> pfb(new BrightnessLevel(fBrightnessFile, fBrightnessMapName, fBrightnessValueName));
379  run.put(std::move(pfb));
380  }
381  //Add by Pengfei -- end
382 
383  //don't reinitialize after the first time through
384  fFirstTime = false;
385  }
386  }
387 
388 
389  //............................................................
391  {
392  // get list of hits
394  evt.getByLabel(fFilterModuleLabel, flslistHandle);
395 
396  //clear output vector
397  fSignals->clear();
398 
399  // Reasonably unique event number (unique if run, subrun < 65536)
400  // for use with cell efficiency maps.
401  fPosCorrEvtId = ((uint64_t)evt.run() << 48)
402  +((uint64_t)evt.subRun() << 32)
403  +((uint64_t)evt.event() );
404 
405  for (unsigned int ilist = 0; ilist < flslistHandle->size(); ++ilist) {
406  const std::vector<sim::FLSHit> hits = flslistHandle->at(ilist).fHits;
407  for (unsigned int ihit = 0; ihit < hits.size(); ++ihit) {
408  sim::FLSHit hit = hits[ihit];
409  //step along the hit and fill fPETimes (clustering done at end of StepAlongHit)
410  StepAlongHit( hit );
411  }
412  }
413 
414  // make a new collection into which we will put the sim::PhotonSignals. That new collection
415  // will then be put into the event. If we tried to put fSignals in the event, the std::move
416  // call changes ownership of the vector and the next time we try to access fSignals would
417  // cause a seg fault
418  std::unique_ptr<std::vector<sim::PhotonSignal>> pscol(new std::vector<sim::PhotonSignal>);
419  pscol->swap(*fSignals);
420 
421  //put vector<sim::PhotonSignal> into the event
422  evt.put(std::move( pscol ));
423 
424  return;
425  }
426 
427  // Helper function for GetPosCorr(). Given an efficiency map and two
428  // points, return the mean efficiency along the line defined by the points.
429  static double line_efficiency(TH2D * themap,
430  const double Xa, const double Ya,
431  const double Xb, const double Yb)
432  {
433  double sum = 0, dem = 0;
434 
435  // This is a dumb method which just takes little steps across and finds
436  // the average, but the performance is fine, so we'll go with it.
437  const double stepsize = 0.01; // cm
438 
439  const double dy = Yb - Ya;
440  const double dx = Xb - Xa;
441  const double seglen = hypot(dx, dy);
442  const int nstep = std::max(1, int(seglen/stepsize));
443  const double thisstepsize = seglen/nstep;
444  const double thisstepstart = thisstepsize/2;
445 
446  for(int step = 0; step < nstep; step++){
447  const double s = thisstepstart + step*thisstepsize;
448  const double sfrac = seglen == 0? 0: s/seglen;
449 
450  // Flip the x coordinate here because of an artifact in how I generated
451  // the efficiency maps. It only matters for horizontals where we want
452  // the fibers to be at the bottom (and air to be at the top). Local +x
453  // is global -y, where fibers should fall, but I generated the maps the
454  // other way around.
455  const double x = -(Xa + sfrac*dx);
456 
457  const double y = Ya + sfrac*dy;
458 
459  double eff = themap->GetBinContent(themap->FindBin(x, y));
460 
461  // There are some zero bins in the corners since the bins are square and
462  // were evaluated at the bin centers. Use any non-zero adjacent bin in
463  // this case.
464  int binx = 0, biny = 0;
465  if(eff <= 0){
466  int binz = 0;
467  themap->GetBinXYZ(themap->FindBin(x, y), binx, biny, binz);
468  eff = themap->GetBinContent(binx+1, biny);
469  }
470  if(eff <= 0) eff = themap->GetBinContent(binx-1, biny);
471  if(eff <= 0) eff = themap->GetBinContent(binx, biny+1);
472  if(eff <= 0) eff = themap->GetBinContent(binx, biny-1);
473 
474  if(eff <= 0) std::cerr << "GetPosCorr should not happen: " << eff
475  << " efficiency with " << themap->GetName()
476  << " point " << -x << " " << y
477  << " from path " << Xa << " " << Ya << " " << Xb << " " << Yb
478  << " length " << seglen << "\n";
479 
480  // Handle the case of zero step size, even though in all cases I've seen,
481  // these particles also have zero energy, so it doesn't matter.
482  sum += thisstepsize == 0? eff: thisstepsize*eff;
483  dem += thisstepsize;
484  }
485 
486  // If the length is zero, return the single efficiency stored in "sum". If
487  // somehow denominator is zero, prevent a NaN, but that shouldn't happen.
488  return seglen == 0? sum: dem == 0? 0: sum/dem;
489  }
490 
491  void load_hset(TFile * f, const char * name, const char * humanname,
492  std::vector<TH2D *> & hold, const bool verbose)
493  {
494  for(int n = 1; ; n++){
495  TH2D * h = dynamic_cast<TH2D *>(f->Get(Form("%s_%d", name, n)));
496  if(h == NULL) break;
497  h->SetDirectory(NULL);
498  hold.push_back(h);
499  }
500 
501  if(verbose) std::cout << "ImprovedTransport: using "
502  << hold.size() << " " << humanname << "\n";
503 
504  if(hold.empty())
505  throw cet::exception("NoPosCorrHists")
506  << "ImprovedTransport: No " << humanname << "\n";
507  }
508 
510  {
511  TFile * f = new TFile(fPosCorrFile.c_str());
512  if(f == NULL || f->IsZombie())
513  throw cet::exception("NoPosCorrHistFile")
514  << "ImprovedTransport: Unable to open "
515  << fPosCorrFile << "\n" << __FILE__ << ":" << __LINE__ << "\n";
516 
517  const bool verbose = msglevel >= 10;
518 
519  load_hset(f, "big_horiz",
520  "big horizontal cell maps",
521  fBigHorizMap, verbose);
522  load_hset(f, "big_horiz_clump",
523  "big horizontal cell maps with clumped fibers",
524  fBigHorizClumpMap, verbose);
525  load_hset(f, "small_horiz",
526  "small horizontal cell maps",
527  fSmallHorizMap, verbose);
528  load_hset(f, "small_horiz_clump",
529  "small horizontal cell maps with clumped fibers",
530  fSmallHorizClumpMap, verbose);
531  load_hset(f, "small_horiz_bubble",
532  "small horizontal cell maps with bubbles",
533  fSmallHorizBubbleMap, verbose);
534  load_hset(f, "small_horiz_bubble_clump",
535  "small horizontal cell maps with bubbles and clumped fibers",
536  fSmallHorizBubbleClumpMap, verbose);
537  load_hset(f, "big_vert",
538  "big vertical cell maps",
539  fBigVertMap, verbose);
540  load_hset(f, "big_vert_corner",
541  "big vertical cell maps with fibers in corners",
542  fBigVertCornerMap, verbose);
543  load_hset(f, "small_vert",
544  "small vertical cell maps",
545  fSmallVertMap, verbose);
546  load_hset(f, "small_vert_corner",
547  "small vertical cell maps with fibers in corners",
548  fSmallVertCornerMap, verbose);
549 
550  delete f;
551  }
552 
553  // Assigns a global index to each map, in the order that the maps are kept
554  // in the class definition. Of course, in order to trace back what map was
555  // used, you also need to know what ROOT file you loaded.
557  const std::vector<TH2D *> * const maps, const unsigned int i)
558  {
559  fPosCorrMapIndex = 0;
560  if(maps == &fBigHorizMap){ fPosCorrMapIndex += i; return; }
561  fPosCorrMapIndex += fBigHorizMap.size();
562  if(maps == &fBigHorizClumpMap){ fPosCorrMapIndex += i; return; }
564  if(maps == &fSmallHorizMap){ fPosCorrMapIndex += i; return; }
566  if(maps == &fSmallHorizClumpMap){ fPosCorrMapIndex += i; return; }
568  if(maps == &fSmallHorizBubbleMap){ fPosCorrMapIndex += i; return; }
570  if(maps == &fSmallHorizBubbleClumpMap){ fPosCorrMapIndex += i; return; }
572  if(maps == &fBigVertMap){ fPosCorrMapIndex += i; return; }
573  fPosCorrMapIndex += fBigVertMap.size();
574  if(maps == &fBigVertCornerMap){ fPosCorrMapIndex += i; return; }
576  if(maps == &fSmallVertMap){ fPosCorrMapIndex += i; return; }
578  if(maps == &fSmallVertCornerMap){ fPosCorrMapIndex += i; return; }
580  }
581 
583  {
584  // If no file is specified for the efficiency maps -> disable
585  if(fPosCorrFile.empty()) return 1;
586 
587  // Y is the large dimension of the cell
588  const double Xa = hit.GetEntryX();
589  const double Ya = hit.GetEntryY();
590  const double Xb = hit.GetExitX();
591  const double Yb = hit.GetExitY();
592  const int planeId = hit.GetPlaneID();
593  const int cellId = hit.GetCellID();
594  const geo::CellUniqueId uniqueID = hit.GetCellUniqueId();
595 
596  const bool horizontal = planeId%2 == 0;
597 
598  const bool largecell = cellId%16 != 0 && cellId%16 != 15;
599 
600  const bool topcell = horizontal && cellId%32 == 31;
601 
603 
604  // This does not test whether we used a geometry that includes air
605  // bubbles, so you need to only use this code if you are. Preferably,
606  // the size of the bubbles should match between the geometry and
607  // the efficiency maps too, which also is not enforced here.
608  const bool hasbubble = topcell &&
609  (geom->DetId() == novadaq::cnv::kFARDET ||
610  geom->DetId() == novadaq::cnv::kTESTBEAM);
611 
612  if(!fPosCorrFile.empty() && fBigHorizMap.empty()) LoadPosCorrHists();
613 
614  // 64 bits of seed are settable. Try to make a decently unique
615  // combination of the event number and cell number with this space. This
616  // gets us the same map within each art::Event given the same cell, but a
617  // different map if either the event or the cell is different. So
618  // multiple energy depositions in a cell in one event are treated sensibly
619  // (although not perfectly: this takes the efficiency to be constant along
620  // the length of the cell)
621  fPosCorrHasher.SetSeed(fPosCorrEvtId ^ uniqueID);
622 
623  // TODO: make probability into a fcl parameter
624  const double clumpfrac = geom->DetId() == novadaq::cnv::kFARDET? 0.16
625  :geom->DetId() == novadaq::cnv::kNEARDET? 0.08
626  : 0.07 /* TESTBEAM (and any others...) */;
627  const bool clump = horizontal && fPosCorrHasher.Rndm() < clumpfrac;
628 
629  bool sagged = false;
630  if(!horizontal){
631  const double Z = (hit.GetEntryZ() + hit.GetExitZ())/2;
632  int dum1 = 0, dum2 = 0;
633  const double frombottom =
634  Z + geom->IdToCell(uniqueID, &dum1, &dum2)->HalfL();
635 
636  // Characteristic distance from bottom where fiber starts to
637  // be free-floating instead sagged against the walls as per
638  // doc-39723. We don't know this number precisely at all, but it
639  // should probably be slightly larger for small cells, so we do
640  // that by the ratio of perimeters of the cells.
641  // TODO: make into a fcl parameter.
642  const double chardist = 100.0 * (1 + 0.012*(!largecell));
643 
644  sagged = fPosCorrHasher.Rndm() < exp(-frombottom/chardist);
645  }
646 
647  // If the vertical fiber is sagged, guess that it will be in opposite
648  // corners half the time, and against the walls elsewhere otherwise.
649  const bool saggedcorner = sagged? fPosCorrHasher.Rndm() < 0.5: false;
650 
651  std::vector<TH2D *> * maps =
652  &(horizontal?
653  (hasbubble?
655  : // no bubble
656  largecell?
658  : // no bubble small cell
660  : // vertical
661  (largecell?
662  (sagged?
663  // For sagged fibers not in the corner, use a *horizontal* map
664  // since those are also have the fibers against the walls. Of
665  // course, it is always the same side wall, so this isn't perfect.
666  (saggedcorner? fBigVertCornerMap: fBigHorizMap /* sic */)
667  : fBigVertMap)
668  : // small cell
669  (sagged?
670  (saggedcorner? fSmallVertCornerMap: fSmallHorizMap /* sic */)
671  : fSmallVertMap)));
672 
673  const unsigned int mapi = fPosCorrHasher.Integer(maps->size());
674  TH2D * themap = (*maps)[mapi];
675  SetPosCorrMapIndex(maps, mapi);
676 
677  // Want this correction to be 1 on average, which this approximates.
678  // It's not important to get this correct to three (or more) digits.
679  // This scale doesn't take into account thresholds, for instance.
680  // It will be necessary to calibrate at a higher level in any case.
681  const double scale = 1/0.134;
682 
683  const double eff = line_efficiency(themap, Xa, Ya, Xb, Yb);
684  #if 0
685  printf("GetPosCorr gave %f for path %f %f %f %f len %f when %s %s %s %s\n",
686  eff, Xa, Ya, Xb, Yb, hypot(Yb-Ya, Xb-Xa),
687  horizontal?"horizontal":"vertical",
688  sagged?"sagged":"notsagged",
689  hasbubble?"bubble":"nobubble",
690  largecell?"large":"small");
691  #endif
692  return eff * scale;
693  }
694 
695  //............................................................
697  {
698  //clear photo-electron times for this hit
699  fPETimes.clear();
700 
701  const double MeVPerGeV = 1000.;
702  //const double green_scale = MeVPerGeV * fPhotonsPerMeV * fQuantumEff;
703 
704  //for cleaner-looking code below...
705  double Xa = hit.GetEntryX();
706  double Ya = hit.GetEntryY();
707  double Za = hit.GetEntryZ();
708  double Ta = hit.GetEntryT();
709 
710  double Xb = hit.GetExitX();
711  double Yb = hit.GetExitY();
712  double Zb = hit.GetExitZ();
713  double Tb = hit.GetExitT();
714 
715  int planeId = hit.GetPlaneID();
716  int cellId = hit.GetCellID();
717  int trackId = hit.GetTrackID();
718 
719  const double positionCorr = GetPosCorr(hit);
720 
721  double brightnessCorr = fBrightnessMap[planeId][cellId];
722 
723  double seglength = util::pythag(Xb - Xa, Yb - Ya, Zb - Za);
724 
725  //cell half-width
726  double DZ = fGeo->Plane((unsigned int)(planeId))->Cell ((unsigned int)(cellId ))->HalfL();
727  //pigtail for cell
728  double pigtail = fGeo->GetPigtail(fGeo->Plane(planeId)->Cell(cellId)->Id());
729 
730  double factor;
731  if (fGeo->Plane(planeId)->View() == geo::kX) factor = fXViewFactor;
732  else factor = fYViewFactor;
733 
734  factor *= fLightLevelScale; // Nominally this is one, but can be modified for LL shifts
735 
736  // If we're simulating drift, randomly throw the light level
737  // not ready to go just yet -- j.
738  if (fSimulateDrift) {
740  int time = rh->TStart() + (0.5*rh->Duration()); // set to midpoint of run;
741  double yDiff = 3.17098e-8 * (time - fDriftReference); // Amount of drift in years
742  //std::cout << "Drifting " << yDiff << " years with respect to nominal." << std::endl;
743  double drift = 1 - (fDriftGradient*yDiff);
744  //std::cout << "Calculated drift of " << drift
745  // << " from nominal, factor corrected from " << factor << " to ";
746  factor *= drift;
747  //std::cout << factor << std::endl;
748  }
749 
750  bool inMC = false;
751  if (fGeo->DetId() == novadaq::cnv::kNEARDET) {
752  if (planeId >= int(fGeo->FirstPlaneInMuonCatcher())) inMC = true;
753  }
754 
755  //find out how many steps to take
756  int nsteps = ceil( seglength/fStep );
757 
758  //Get approximate dEdx for scintillation time scaling
759  //If segment length is zero, dEdx won't make sense.
760  //Default to zero - the beta PDF is almost always going to be more relevant than the alpha one, I think
761  double dEdx = (seglength > 0) ? MeVPerGeV*hit.GetEdep()/seglength : 0.0;
762 
763  // The total number of PE that would have been produced are accumulated in
764  // here, whether or not fluctuations allowed a PhotonSignal to be created
765  // for them to live in.
766  double lambda = 0;
767 
768  for (int istep = 0; istep < nsteps; ++istep) {
769  //at midpoint of the step
770  double stepZ = Za + ((Zb-Za)/nsteps)*(istep+0.5);
771  double stepT = Ta + ((Tb-Ta)/nsteps)*(istep+0.5);
772  double stepE = (fApplyBirks) ? hit.GetEdepBirks()/nsteps : hit.GetEdep()/nsteps;
773  double stepCerenkov = hit.GetNCerenkov()/nsteps;
774 
775  double tweakCorr(1.0);
776  if (fApplyTweak) {
777  //Tweak is currently only available for non-muon catcher cells
778  if (fGeo->DetId() == novadaq::cnv::kFARDET) {
779  if (fGeo->Plane(planeId)->View() == geo::kX) tweakCorr = fTweak_X.Interpolate(DZ - stepZ);
780  else tweakCorr = fTweak_Y.Interpolate(DZ - stepZ);
781  }
782  if (fGeo->DetId() == novadaq::cnv::kNEARDET) {
783  if (inMC) {
784  if (fGeo->Plane(planeId)->View() == geo::kX) tweakCorr = fTweak_XMC.Interpolate(DZ - stepZ);
785  else tweakCorr = fTweak_YMC.Interpolate(DZ - stepZ);
786  }
787  else {
788  if (fGeo->Plane(planeId)->View() == geo::kX) tweakCorr = fTweak_X.Interpolate(DZ - stepZ);
789  else tweakCorr = fTweak_Y.Interpolate(DZ - stepZ);
790  }
791  }
792  }
793 
794  //green_scale*stepE*collectionRate = N Photo-electrons (before attenuation)
795  std::list<RateInfo>::iterator rate_it;
796  for (rate_it = fRateInfo.begin(); rate_it != fRateInfo.end(); ++rate_it) {
797  double collectionRate = rate_it->rate;
798  double currZ = stepZ + rate_it->z;
799  double fraction(1.0);
800 
801  //is this currZ inside the cell?
802  if ( std::fabs(currZ) - 0.5*rate_it->zWidth > DZ ) continue;
803  else if ( std::fabs(currZ) + 0.5*rate_it->zWidth < DZ ) fraction = 1.0;
804  else fraction = (DZ - ( std::fabs(currZ) - 0.5*rate_it->zWidth ))/rate_it->zWidth;
805  collectionRate *= fraction;
806 
807  double scintillationPhotons = MeVPerGeV*fPhotonsPerMeV*stepE;
808  double cerenkovPhotons = fCerenkovEff*stepCerenkov;
809  double muGreenScint = fQuantumEff*collectionRate*brightnessCorr*tweakCorr*
810  positionCorr*factor*scintillationPhotons;
811  double muGreenCkv = fQuantumEff*collectionRate*brightnessCorr*tweakCorr*
812  positionCorr*factor*cerenkovPhotons;
813  //short and long length to the electronics
814  double dist0 = DZ - currZ + pigtail;
815  double dist1 = 3*DZ + currZ + pigtail;
816 
817  CreatePhotoElectrons( 0.5*muGreenScint, stepT, dist0, dEdx, true, *rate_it, lambda );
818  CreatePhotoElectrons( 0.5*muGreenScint, stepT, dist1, dEdx, true, *rate_it, lambda );
819  CreatePhotoElectrons( 0.5*muGreenCkv, stepT, dist0, dEdx, false, *rate_it, lambda );
820  CreatePhotoElectrons( 0.5*muGreenCkv, stepT, dist1, dEdx, false, *rate_it, lambda );
821  }
822  }
823 
824  // Keep track of how many we have so far so we know which ones get newly
825  // added by ClusterPhotoElectrons().
826  unsigned int start = fSignals->size();
827 
828  ClusterPhotoElectrons(planeId, cellId, trackId);
829 
830  // Divide the poission lambda evenly among all PhotonSignals. When someone
831  // adds them up downstream they'll get the right overall answer for this
832  // FLSHit.
833  const int N = fSignals->size()-start;
834  if(N > 0){
835  for(unsigned int i = start; i < fSignals->size(); ++i){
836  (*fSignals)[i].SetPoissonLambda(lambda/N);
837  }
838  }
839  else{
840  if(fWriteEmptySignals){
841  // Create a PhotonSignal with zero PE, just to have somewhere to attach
842  // the PoissonLambda to.
843  sim::PhotonSignal photon(hit.GetEntryT(), 0);
844  photon.SetPlaneCell(planeId, cellId);
845  photon.SetTrackId(trackId);
846  photon.SetPoissonLambda(lambda);
847  fSignals->push_back(photon);
848  }
849  }
850 
851  return;
852  }
853 
854  //............................................................
855  void ImprovedTransport::CreatePhotoElectrons(double meanPE, double stepTime, double dist, double dEdx, bool isScint, RateInfo& info, double& lambda)
856  {
857  double expect = meanPE*FiberTransmission(dist);
858  lambda += expect;
859 
860  int nPE = fPoisson->fire( expect );
861  double rand0(0), rand1(0), rand2(0);
862  for (int iPE=0; iPE < nPE; ++iPE) {
863  double time = stepTime;
864  time += dist/fFiberSpeedOfLight;
865 
866  fRandHisto->GetRandom( &info.timePDF, rand0, rand1, rand2 );
867  time += rand0;
868 
869  time += fExp->fire( fEmissionTau );
870 
871  double scintTime(0);
872  if (fUseScintTime)
873  {
874  if (isScint) scintTime = ScintTime(dEdx);
875  else scintTime = fExp->fire( fCherenkovTau );
876  }
877  time += scintTime;
878 
879  double fiberspread = 0;
880  if (fApplyFiberSmearing) {
881  double weight(-1.0), uweight(0.0);
882  while (weight < uweight) {
883  // Draw a Dt/z value then multiply by "dist" to get time (in ns)
884  fRandHisto->GetRandom( &fFiberSmearingHisto, rand0, rand1, rand2 );
885  fiberspread = dist*rand0;
886 
887  // Use fiber attenuation to deweight longer times (<==>longer distances)
888  // This is an approximation, but should be pretty close.
889  double Delta_dist = fiberspread*fFiberSpeedOfLight;
890  weight = FiberTransmission(dist+Delta_dist)/FiberTransmission(dist);
891  uweight = fFlat->fire(0,1);
892  }
893  }
894  time += fiberspread;
895  fPETimes.push_back( time );
896  fScintTime->Fill(time - stepTime);
897  }
898  return;
899  }
900 
901 
902  //............................................................
903  void ImprovedTransport::ClusterPhotoElectrons(int planeId, int cellId, int trackId)
904  {
905  if (fPETimes.empty()) return;
906  //cluster photons into PhotonSignals fTimeClustering long
907 
908  sim::PhotonSignal photon;
909  photon.SetPlaneCell(planeId, cellId);
910  photon.SetTrackId(trackId);
911 
912  // TODO: Would like to record which map is used, but having trouble
913  // getting it to work.
914  // photon.SetPosCorrI(fPosCorrFile.empty()?-1:fPosCorrMapIndex);
915 
916  double meanTime(0);
917  int nPE(0);
918 
919  fPETimes.sort();
920  std::list< double >::iterator it;
921  for (it = fPETimes.begin(); it != fPETimes.end(); ++it) {
922  double currTime = *it;
923  if (nPE == 0) {
924  meanTime = currTime;
925  nPE = 1;
926  }
927  else if ( fabs(meanTime - currTime) < fTimeClustering ) {
928  //add to cluster
929  meanTime = meanTime*nPE + currTime;
930  ++nPE;
931  meanTime /= nPE;
932  }
933  else {
934  photon.SetTimeMean(meanTime);
935  photon.SetNPhoton(nPE);
936  fSignals->push_back(photon);
937 
938  //Set variables with current iteration
939  nPE = 1;
940  meanTime = currTime;
941  }
942  }
943  //there is one PhotonSignal that hasn't been saved yet
944  photon.SetTimeMean(meanTime);
945  photon.SetNPhoton(nPE);
946  fSignals->push_back(photon);
947 
948  return;
949  }
950 
951  //............................................................
953  {
954  // just read a file from disk for now. should eventually be in the DB.
955  // use FileInPath to get the histogram file
956  struct stat sb;
957  if ( stat(fCollectionRateFile.c_str(), &sb) !=0 ) {
958  throw cet::exception("NoCollectionRateFile")
959  << "ImprovedTransport: Unable to open "
960  << fCollectionRateFile << "\n"
961  << __FILE__ << ":" << __LINE__ << "\n";
962  return false;
963  }
964 
965  if ( stat(fSmearingHistoFile.c_str(), &sb) !=0 ) {
966  throw cet::exception("NoSmearingHistoFile") << "ImprovedTransport: Unable to open " << fSmearingHistoFile << "\n" << __FILE__ << ":" << __LINE__ << "\n";
967  return false;
968  }
969 
970  TFile f1(fCollectionRateFile.c_str());
971  TH2D *htemp = (TH2D*)f1.Get("dT_dZ_CollectionRate")->Clone();
972  if (!htemp) {
973  std::cerr << "ImprovedTransport: Histogram dT_dZ_CollectionRate not found in " << fCollectionRateFile << std::endl;
974  return false;
975  }
976 
977  fCollectionRateHisto = *htemp;
978  fCollectionRateHisto.SetDirectory(0);
979  f1.Close();
980 
981  TFile f2(fSmearingHistoFile.c_str());
982  TH1D *htemp2 = (TH1D*)f2.Get("Dt_per_z")->Clone();
983  if (!htemp2) {
984  std::cerr << "ImprovedTransport: Histogram Dt_per_z_distribution not found in " << fSmearingHistoFile << std::endl;
985  return false;
986  }
987 
988  fFiberSmearingHisto = *htemp2;
989  fFiberSmearingHisto.SetDirectory(0);
990  f2.Close();
991 
992  for (unsigned int plane = 0; plane < fGeo->NPlanes(); plane++){
993  std::vector< double > plane_brightness;
994  for (unsigned int cell = 0; cell < fGeo->Plane(plane)->Ncells(); cell++){
995  plane_brightness.push_back(1.0);
996  }
997  fBrightnessMap.push_back(plane_brightness);
998  }
999 
1000  if (fBrightnessMode != "None") {
1001  TFile f3(fBrightnessFile.c_str());
1002  TH1D* hBrightness1D(0);
1003  TH2D* hBrightness2D(0);
1004  if (fBrightnessMode == "ByCell") {
1005  hBrightness2D = (TH2D*)f3.Get("BrightnessByCell");
1006  if (!hBrightness2D) {
1007  std::cerr << "ImprovedTransport: Histogram BrightnessByCell not found in " << fBrightnessFile << std::endl;
1008  return false;
1009  }
1010  }
1011  else if (fBrightnessMode == "ByModule") {
1012  hBrightness2D = (TH2D*)f3.Get("BrightnessByModule");
1013  if (!hBrightness2D) {
1014  std::cerr << "ImprovedTransport: Histogram BrightnessByModule not found in " << fBrightnessFile << std::endl;
1015  return false;
1016  }
1017  }
1018  else if (fBrightnessMode == "ByPlane") {
1019  hBrightness1D = (TH1D*)f3.Get("BrightnessByPlane");
1020  if (!hBrightness1D) {
1021  std::cerr << "ImprovedTransport: Histogram BrightnessByPlane not found in " << fBrightnessFile << std::endl;
1022  return false;
1023  }
1024  }
1025  //Add by Pengfei -- begin
1026  else if (fBrightnessMode == "ByBin") {
1028  }
1029  //Add by Pengfei -- end
1030  else {
1031  std::cerr << "ImprovedTransport: " << fBrightnessMode << " is not a recognized brightness correction mode. Valid choices are None, ByCell, ByModule, or ByPlane " << std::endl;
1032  return false;
1033  }
1034 
1035  for (unsigned int plane = 0; plane < fGeo->NPlanes(); plane++){
1036  for (unsigned int cell = 0; cell < fGeo->Plane(plane)->Ncells(); cell++){
1037  double brightness(1.0);
1038  if (fBrightnessMode == "ByCell") brightness = hBrightness2D->GetBinContent(plane+1, cell+1);
1039  else if (fBrightnessMode == "ByModule") brightness = hBrightness2D->GetBinContent(plane+1, cell/32 + 1);
1040  else if (fBrightnessMode == "ByPlane") brightness = hBrightness1D->GetBinContent(plane+1);
1041  //Add by Pengfei -- begin
1042  else if (fBrightnessMode == "ByBin") {
1044  }
1045  //Add by Pengfei -- end
1047  }
1048  }
1049  f3.Close();
1050  }
1051 
1052  if (fApplyTweak) {
1053  TFile f4(fTweakFile.c_str());
1054  TH1D* tmpTweak(0);
1055  tmpTweak = (TH1D*)f4.Get("Tweak_X")->Clone();
1056  if (!tmpTweak) {
1057  std::cerr << "Histogram Tweak_X not found in " << fTweakFile << std::endl;
1058  return false;
1059  }
1060  fTweak_X = *tmpTweak;
1061  fTweak_X.SetDirectory(0);
1062 
1063  tmpTweak = 0;
1064  tmpTweak = (TH1D*)f4.Get("Tweak_Y")->Clone();
1065  if (!tmpTweak) {
1066  std::cerr << "Histogram Tweak_Y not found in " << fTweakFile <<std::endl;
1067  return false;
1068  }
1069  fTweak_Y = *tmpTweak;
1070  fTweak_Y.SetDirectory(0);
1071 
1072  f4.Close();
1073  }
1074 
1075  fRateInfo.clear();
1076  std::stringstream hName;
1077  int nbinsY = fCollectionRateHisto.GetNbinsY();
1078  double ylo = fCollectionRateHisto.GetYaxis()->GetXmin();
1079  double yhi = fCollectionRateHisto.GetYaxis()->GetXmax();
1080  for (int iZ = 1; iZ <= fCollectionRateHisto.GetNbinsX(); ++iZ) {
1081  RateInfo info_pos;
1082  info_pos.rate = 0.0;
1083  info_pos.z = fCollectionRateHisto.GetXaxis()->GetBinCenter(iZ);
1084  info_pos.zWidth = fCollectionRateHisto.GetXaxis()->GetBinWidth(iZ);
1085  hName.str("timePDF_Pos_");
1086  hName << iZ;
1087  info_pos.timePDF.SetName(hName.str().c_str());
1088  info_pos.timePDF.SetBins(nbinsY, ylo, yhi);
1089  info_pos.timePDF.SetDirectory(0);
1090 
1091  RateInfo info_neg;
1092  info_neg.rate = 0.0;
1093  info_neg.z = -fCollectionRateHisto.GetXaxis()->GetBinCenter(iZ);
1094  info_neg.zWidth = fCollectionRateHisto.GetXaxis()->GetBinWidth(iZ);
1095  hName.str("timePDF_Neg_");
1096  hName << iZ;
1097  info_neg.timePDF.SetName(hName.str().c_str());
1098  info_neg.timePDF.SetBins(nbinsY, ylo, yhi);
1099  info_neg.timePDF.SetDirectory(0);
1100 
1101  for (int iT = 1; iT <= fCollectionRateHisto.GetNbinsY(); ++iT) {
1102  double collectionRate = 0.5*fCollectionRateHisto.GetBinContent(iZ,iT);
1103  info_pos.rate += collectionRate;
1104  info_pos.timePDF.SetBinContent(iT, collectionRate);
1105 
1106  info_neg.rate += collectionRate;
1107  info_neg.timePDF.SetBinContent(iT, collectionRate);
1108  }
1109  fRateInfo.push_back(info_pos);
1110  fRateInfo.push_back(info_neg);
1111  }
1112  return true;
1113  }
1114 
1115  //............................................................
1117  {
1118  // currently, all fibers are the same
1119 
1120  // fAttenPars holds a series of (c,L) pairs that form the
1121  // function:
1122  // transmission = c0*exp(-x/L0) + c1*exp(-x/L1) + ...
1123  //
1124  // Thus, we should have an even number of parameters
1125 
1126  double answer = 0;
1127 
1128  int nPar = fAttenPars.size();
1129  if (!nPar) {
1130  mf::LogWarning("AttenPars") << "No AttenPars entries. Fiber transmission set to 1!";
1131  answer = 1;
1132  }
1133  else if (nPar%2) {
1134  mf::LogWarning("AttenPars") << "Odd number of AttenPars entries. Fiber transmission set to 1!";
1135  answer = 1;
1136  }
1137  else {
1138  answer = 0;
1139  for (Int_t iP=0;iP<nPar;iP+=2) answer += fAttenPars[iP]*exp(-x/fAttenPars[iP+1]);
1140  }
1141  return answer;
1142  }
1143 
1144  //............................................................
1145  double ImprovedTransport::ScintTime(double dEdx)
1146  {
1147  // sample from a sum of multiple exponentials to get the scintilator decay time for each photon
1148  // if dEdx is less than that used for measuring the beta distribution, use the beta distribution
1149  // if dEdx is greater than that used for measuring the alpha distribution, use the alpha distribution
1150  // if dEdx is between the two, smoothly transition between the distributions
1151 
1152  int nProbBeta = fScintBetaProb.size();
1153  int nProbAlpha = fScintAlphaProb.size();
1154  int nTauBeta = fScintBetaTau.size();
1155  int nTauAlpha = fScintAlphaTau.size();
1156 
1157  double time = 0;
1158  if (!nProbBeta) {
1159  mf::LogWarning("ScintTime") << "No ScintBetaProb entries. No scintillation delays applied!";
1160  time = 0;
1161  }
1162  else if (!nProbAlpha) {
1163  mf::LogWarning("ScintTime") << "No ScintAlphaProb entries. No scintillation delays applied!";
1164  time = 0;
1165  }
1166  else if (nProbBeta != nTauBeta) {
1167  mf::LogWarning("ScintTime") << "Number of entries in ScintBetaProb and ScintBetaTau not the same. No scintillation delays applied!";
1168  time = 0;
1169  }
1170  else if (nProbAlpha != nTauAlpha) {
1171  mf::LogWarning("ScintTime") << "Number of entries in ScintBetaAlpha and ScintBetaAlpha not the same. No scintillation delays applied!";
1172  time = 0;
1173  }
1174  else if (fScintBetaDEDX >= fScintAlphaDEDX) {
1175  mf::LogWarning("ScintTime") << "Beta dEdx must be strictly smaller than alpha dEdx. No scintillation delays applied!";
1176  time = 0;
1177  }
1178  else {
1179  //actually sample the distribution
1180  double probBeta = (dEdx - fScintBetaDEDX)/(fScintBetaDEDX - fScintAlphaDEDX) + 1;
1181  double u = fFlat->fire(0,1);
1182  bool useBeta = (u < probBeta);
1183  u = fFlat->fire(0,1);
1184  if (useBeta) {
1185  for (int iC = 0; iC < nProbBeta; ++iC) {
1186  if (u <= fScintBetaProb[iC]) {
1187  time = fExp->fire(fScintBetaTau[iC]);
1188  break;
1189  }
1190  }
1191  }
1192  else {
1193  for (int iC = 0; iC < nProbAlpha; ++iC) {
1194  if (u <= fScintAlphaProb[iC]) {
1195  time = fExp->fire(fScintAlphaTau[iC]);
1196  break;
1197  }
1198  }
1199  }
1200  }
1201  return time;
1202  }
1203 
1205 } //namespace
double fCerenkovEff
Absorption and re-emission efficiency for Cerenkov photons.
geo::CellUniqueId GetCellUniqueId() const
Unique Cell ID.
Definition: FLSHit.h:41
bool fApplyFiberSmearing
Turn on fiber smearing.
T max(const caf::Proxy< T > &a, T b)
const XML_Char XML_Encoding * info
Definition: expat.h:530
Float_t f4
float GetEntryX() const
Entry point of the particle (position, time and energy)
Definition: FLSHit.h:48
ImprovedTransport(fhicl::ParameterSet const &pset)
const XML_Char * name
Definition: expat.h:151
bool fApplyTweak
Turn on the data driven attenuation correction for each view.
double HalfL() const
Definition: CellGeo.cxx:198
void ClusterPhotoElectrons(int planeId, int cellId, int trackId)
SubRunNumber_t subRun() const
Definition: Event.h:72
TH2D fCollectionRateHisto
Collection rate for given dT vs. dZ.
void SetTrackId(int t)
Definition: PhotonSignal.h:25
double fStep
step size for walking along the track (cm)
std::string fCollectionRateFile
Relative path to the collection rate histo file.
std::vector< TH2D * > fSmallHorizBubbleClumpMap
TH2 * rh
Definition: drawXsec.C:5
set< int >::iterator it
int GetPlaneID() const
Plane ID.
Definition: FLSHit.h:37
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::enable_if_t< std::is_arithmetic< T >::value, T > hypot(T x, T y)
Definition: hypot.h:60
double fYViewFactor
scale factor for the y-view
TH1D fTweak_XMC
A data/mc correction for the X-View in the muon catcher.
TH1D fTweak_Y
A data/mc correction for the Y-View.
const Var weight
double fXViewFactor
scale factor for the x-view
double GetPosCorr(const sim::FLSHit &hit)
int GetCellID() const
Cell ID.
Definition: FLSHit.h:39
float GetExitT() const
Definition: FLSHit.h:57
const int nsteps
double fTimeClustering
groups photons closer than this many ns apart
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
double fLightLevelScale
shift overall light yield by multiplicative factor
void LoadPosCorrHists()
Reads in the efficiency maps from fPosCorrFile.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
a struct to hold collection rate info
bool fWriteEmptySignals
Write an empty PhotonSignal if there&#39;s 0 PE.
double fFiberSpeedOfLight
Speed of light in the fiber.
base_engine_t & createEngine(seed_t seed)
float GetNCerenkov() const
Get total N Cerenkov photons.
Definition: FLSHit.h:35
OStream cerr
Definition: OStream.cxx:7
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
art::ServiceHandle< geo::Geometry > fGeo
Handle to the geometry service.
A single unit of energy deposition in the liquid scintillator.
Definition: FLSHit.h:19
Float_t f2
const PlaneGeo * Plane(unsigned int i) const
This slightly less simple photon transport uses a template for photon collection in time and in dista...
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
float GetExitX() const
Exit point of the particle (position, time and energy)
Definition: FLSHit.h:54
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
double fFiberIndexOfRefraction
n for fiber core
a class for transporting photons in a roughly realistic way
TH1D fTweak_YMC
A data/mc correction for the Y-View in the muon catcher.
double fScintBetaDEDX
Approximate dE/dx of betas used to determine PDF.
Definition: Run.h:31
double fCherenkovTau
Time constant for directly excited PC to PPO scintillation.
void SetPlaneCell(int plane, int cell)
Definition: PhotonSignal.h:24
std::unique_ptr< std::vector< sim::PhotonSignal > > fSignals
Produced vector of PhotonSignals.
double dist
Definition: runWimpSim.h:113
const CellUniqueId & Id() const
Definition: CellGeo.h:55
std::list< RateInfo > fRateInfo
List of RateInfo object.
Float_t Z
Definition: plot.C:38
std::vector< double > fScintAlphaTau
Time constants of each exponential component for alphas.
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
const XML_Char * s
Definition: expat.h:262
Float_t f3
base_engine_t & getEngine() const
Double_t scale
Definition: plot.C:25
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
std::vector< double > fAttenPars
from Leon&#39;s fit
double getBrightnessLevel(unsigned int plane, unsigned int cell) const
Far Detector at Ash River, MN.
std::vector< double > fScintBetaProb
Cumulative probability for each exponential component for betas.
void hits()
Definition: readHits.C:15
double dy[NP][NC]
double dx[NP][NC]
double fDriftGradient
gradient for light level detector aging systematic [fraction/year]
unsigned int seed
Definition: runWimpSim.h:102
Prototype Near Detector on the surface at FNAL.
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::string fBrightnessMode
correct the fiber brightness by cell, by module, or not at all
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
printf("%d Experimental points found\n", nlines)
void SetPosCorrMapIndex(const std::vector< TH2D * > *const maps, const unsigned int i)
Near Detector in the NuMI cavern.
void SetNPhoton(int npe)
Definition: PhotonSignal.h:23
int Duration()
in units of seconds
Definition: RunHistory.h:389
std::string fTweakFile
Relative path to the data driven attenuation correction for each view.
bool fSimulateDrift
whether to vary light level as a function of time
Float_t f1
EventNumber_t event() const
Definition: Event.h:67
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
Definition: FillTruth.h:15
std::list< double > fPETimes
List of photo-electron times (each representing one observed PE)
float GetEdepBirks() const
Get total Energy with Birks suppression deposited into the cell for the whole FLSHit.
Definition: FLSHit.h:33
Definition: run.py:1
double fEmissionTau
exponential time constant for green photon emission
OStream cout
Definition: OStream.cxx:6
double fScintAlphaDEDX
Approximate dE/dx of alphas used to determine PDF.
double GetPigtail(const CellUniqueId &id) const
Return length of fiber in cm from end of cell to apd.
bool fUseScintTime
Turn on sampling from PDF of scintillator decay time, see arXiv:1704.02291 for more details...
float GetEntryZ() const
Definition: FLSHit.h:50
std::vector< std::vector< double > > fBrightnessMap
The brightness correction by cell.
#define expect(expr, value)
Definition: lz4.cxx:117
bool fApplyBirks
Use EdepBirks from FLSHits instead of calculating here.
float GetEntryT() const
Definition: FLSHit.h:51
void CreatePhotoElectrons(double meanPE, double stepTime, double dist, double dEdx, bool isScint, RateInfo &info, double &lambda)
unsigned int GetRandomNumberSeed()
Definition: Simulation.cxx:13
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
std::vector< TH2D * > fSmallVertCornerMap
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
TH1D fFiberSmearingHisto
Smearing histo for photons bouncing around in fiber.
int GetTrackID() const
Track ID.
Definition: FLSHit.h:45
std::vector< double > fScintAlphaProb
Cumulative probability for each exponential component for alphas.
static double line_efficiency(TH2D *themap, const double Xa, const double Ya, const double Xb, const double Yb)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
Definition: event.h:1
std::vector< double > fScintBetaTau
Time constants of each exponential component for betas.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
std::string fSmearingHistoFile
Relative path to the fiber smearing histo file.
assert(nhit_max >=nhit_nbins)
double fPhotonsPerMeV
blue photon production scale
void GetRandom(TH1 *histo, double &x, double &y, double &z)
Definition: RandHisto.cxx:8
std::vector< TH2D * > fSmallHorizClumpMap
std::string fBrightnessFile
Relative path to the brightness correction file.
float GetEdep() const
Get total Energy deposited into the cell for the whole FLSHit.
Definition: FLSHit.h:31
unsigned int NPlanes() const
float GetExitY() const
Definition: FLSHit.h:55
Double_t sum
Definition: plot.C:31
double fQuantumEff
APD quantum efficiency.
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
float GetExitZ() const
Definition: FLSHit.h:56
void load_hset(TFile *f, const char *name, const char *humanname, std::vector< TH2D * > &hold, const bool verbose)
const CellGeo * IdToCell(const CellUniqueId &id, int *iplane, int *icell) const
void SetTimeMean(double t)
Definition: PhotonSignal.h:22
std::string fFilterModuleLabel
module that filtered the fls hits
int fDriftReference
drift reference time (ie. where the slope is fixed to 1)
Encapsulate the geometry of one entire detector (near, far, ndos)
fvar< T > ceil(const fvar< T > &x)
Definition: ceil.hpp:11
TH1D fTweak_X
A data/mc correction for the X-View.
std::vector< TH2D * > fSmallHorizBubbleMap
float GetEntryY() const
Definition: FLSHit.h:49
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.