PhotonTransport_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief This slightly less simple photon transport doesn't propagate
3 /// individual photons (thus "simple"), but the model is gaining
4 /// some realism.
5 ///
6 /// \author rbpatter@caltech.edu
7 ////////////////////////////////////////////////////////////////////////
8 
9 // Framework includes
16 
17 #include "cetlib/search_path.h"
18 #include "fhiclcpp/ParameterSet.h"
20 
21 #include "CLHEP/Random/RandFlat.h"
25 
26 // NOvA includes
27 #include "Geometry/Geometry.h"
29 #include "Simulation/FLSHitList.h"
31 #include "Simulation/Simulation.h"
32 #include "Utilities/func/MathUtil.h"
34 
35 // ROOT includes
36 #include "TFile.h"
37 #include "TH1.h"
38 #include "TMath.h"
39 
40 #include <cmath>
41 #include <sys/stat.h>
42 
43 namespace photrans{
44  /// a class for transporting photons in a roughly realistic way
46  {
47  public:
48  explicit PhotonTransport(fhicl::ParameterSet const &pset);
49  virtual ~PhotonTransport();
50 
51  void produce(art::Event& evt);
52 
53  protected:
54  double FiberTransmission(double x);
55  int LoadSmearingHisto();
56  /// The length (in cm) of the fibre from the end of the cell to the APD
57 
58  private:
59 
60  double fPhotonsPerMeV; ///< blue photon production scale
61  double fQuantumEff; ///< APD quantum efficiency
62  double fCollectionEff; ///< fiber collection efficiency
63  double fFiberIndexOfRefraction; ///< n for fiber core
64  double fEmissionTau; ///< exponential time constant for green photon emission
65  double fTimeSpread; ///< spread in time due to everything
66  std::vector<double> fAttenPars; ///< from Leon's fit
67  double fTimeClustering; ///< groups photons closer than this many ns apart
68  double fStep; ///< step size for walking along the track (cm)
69  bool fApplyFiberSmearing; ///< apply fiber smearing
70  bool fNoAttenNoPoissonMode; ///< run special no-attenuation, no-Poisson-fluct mode?
71  bool fApplyBirks; ///< Use EdepBirks from FLSHits instead of calculating here
72 
74  int msglevel;
75 
76  std::string fFilterModuleLabel; ///< module that filtered the fls hits
77  std::string fSmearingHistoFile; ///< relative path to the smearing histo file
78  };
79 
80  //............................................................
82  : fPhotonsPerMeV (pset.get< double >("PhotonsPerMeV"))
83  , fQuantumEff (pset.get< double >("QuantumEff"))
84  , fCollectionEff (pset.get< double >("CollectionEff"))
85  , fFiberIndexOfRefraction(pset.get< double >("FiberIndexOfRefraction"))
86  , fEmissionTau (pset.get< double >("EmissionTau"))
87  , fTimeSpread (pset.get< double >("TimeSpread"))
88  , fAttenPars (pset.get< std::vector<double> >("AttenPars"))
89  , fTimeClustering (pset.get< double >("TimeClustering"))
90  , fStep (pset.get< double >("Step"))
91  , fApplyFiberSmearing (pset.get< bool >("ApplyFiberSmearing"))
92  , fNoAttenNoPoissonMode (pset.get< bool >("NoAttenNoPoissonMode"))
93  , fApplyBirks (pset.get< bool >("ApplyBirks"))
94  , msglevel (pset.get< int >("MessageLevel"))
95  , fFilterModuleLabel (pset.get< std::string >("FilterModuleLabel"))
96  {
97  // get the random number seed, use a random default if not specified
98  // in the configuration file.
99  unsigned int seed = pset.get< unsigned int >("Seed", sim::GetRandomNumberSeed());
100 
101  createEngine( seed );
102 
103  produces< std::vector<sim::PhotonSignal> >();
104 
105  // constructor decides if initialized value is a path or an environment variable
106  cet::search_path sp("FW_SEARCH_PATH");
107 
108  sp.find_file(pset.get< std::string >("SmearingHistoFile"), fSmearingHistoFile);
109 
110  // Open the GDML file
111  struct stat sb;
112  if ( fSmearingHistoFile.empty() || stat(fSmearingHistoFile.c_str(), &sb)!=0 ) {
113  // failed to resolve the file name
114  throw cet::exception("NoSmearingHisto")
115  << "smearing histogram file " << fSmearingHistoFile << " not found!\n"
116  << __FILE__ << ":" << __LINE__ << "\n";
117  }
118 
119  }
120 
121  //............................................................
123  {
124  }
125 
126  //............................................................
128  {
129 
130  if (fApplyFiberSmearing&&!fFiberSmearingHisto.Integral()) {
131  // Load up the distribution of Delta(t) (per z-distance). Note:
132  // this distribution was generated using an offline fiber MC. The
133  // exact distribution depends on the choice of green photon starting
134  // positions within the fiber. The version stamped Oct 25 2008
135  // assumes that blue photons originate 2 cm away with a +/- 45deg
136  // angular spread in Z (arbitrary!) and with a 0.03 cm absorption
137  // length in the fiber.
138  if ( LoadSmearingHisto() ) {
139  throw cet::exception("NoSmearingHisto")
140  << "can't load smearing histo\n"
141  << __FILE__ << ":" << __LINE__ << "\n";
142  }
143  }
144 
146 
147  // get the RandomNumberGenerator service
149  CLHEP::HepRandomEngine &engine = rng->getEngine();
151  CLHEP::RandGaussQ gauss(engine);
152  CLHEP::RandExponential exponential(engine);
153  CLHEP::RandFlat flat(engine);
154  RandHisto histoDstrb(engine);
155 
156 
157  // get list of hits
159  evt.getByLabel(fFilterModuleLabel, flslistHandle);
160  // make a art::PtrVector of the FLSHitLists
162  for(unsigned int i = 0; i < flslistHandle->size(); ++i){
163  art::Ptr<sim::FLSHitList> flshl(flslistHandle, i);
164  hitlist.push_back(flshl);
165  }
166 
167  // make output vector
168  std::unique_ptr<std::vector<sim::PhotonSignal> > photlist(new std::vector<sim::PhotonSignal>);
169 
170  const double MeVPerGeV = 1000.;
171  const double c = TMath::Ccgs()*1e-9; // c in cm/ns
172  const double cn = c/fFiberIndexOfRefraction;
173 
174  const double green_scale = MeVPerGeV * fPhotonsPerMeV * fCollectionEff * fQuantumEff;
175 
176  int nHit = 0;
177 
178  /// Loop over multiple lists of hits in event
180  while ( iter_hitlist != hitlist.end() ) {
181 
182  /// loop over multiple hits in current hitlist
183  std::vector<sim::FLSHit> hits = (*iter_hitlist)->fHits;
184 
185  std::vector<sim::FLSHit>::iterator iter_hit = hits.begin();
186  while ( iter_hit != hits.end() ) {
187 
188  nHit++;
189 
190  // for cleaner-looking code below...
191  double Xa = (*iter_hit).GetEntryX();
192  double Ya = (*iter_hit).GetEntryY();
193  double Za = (*iter_hit).GetEntryZ();
194  double Xb = (*iter_hit).GetExitX();
195  double Yb = (*iter_hit).GetExitY();
196  double Zb = (*iter_hit).GetExitZ();
197  // double segtime = (*iter_hit).fT;
198 
199  double Ta = (*iter_hit).GetEntryT();
200  double Tb = (*iter_hit).GetExitT();
201 
202  int planeId = (*iter_hit).GetPlaneID();
203  int cellId = (*iter_hit).GetCellID();
204  int trackId = (*iter_hit).GetTrackID();
205 
206 
207  double seglength = util::pythag(Xb - Xa, Yb - Ya, Zb - Za);
208 
209  //Calculate the velocity:
210  double v = ((Tb-Ta>0) ? seglength/(Tb-Ta) : c);
211  if (v>c) v=c;
212 
213 
214  // cell half-width
215  double DZ = fGeo->Plane((unsigned int)(planeId))->Cell ((unsigned int)(cellId ))->HalfL();
216 
217  // walk along track segment
218  std::list<double> currtimes;
219  currtimes.clear();
220  double muTot = 0;
221  for (double currS=0; currS<=seglength; currS+=fStep) {
222 
223  // last step is truncated
224  double dS = fStep;
225  if ( dS > seglength-currS ) dS = seglength-currS;
226 
227  // current fraction along segment, at step midpoint
228  double u = (seglength) ? (currS + dS / 2) / seglength : 0;
229 
230  // current location (only Z component is needed currently)
231  double currZ = Za + (Zb - Za) * u;
232 
233  // green photon travel distances for this step
234  // "pigtail" is the additional distance outside the cell
235  // --> short path and long path through fiber:
236  double pigtail = fGeo->GetPigtail(fGeo->Plane(planeId)->Cell(cellId)->Id());
237  double dist0 = DZ - currZ + pigtail;
238  double dist1 = 3*DZ + currZ + pigtail;
239 
240  // Mean number of green photons for this step
241  // If segment has length zero, dump all the energy right here. The upstream
242  // code should be examined to figure out why seglength is sometimes zero.
243  double muGreenT = green_scale * (seglength ? dS / seglength : 1);
244  if (fApplyBirks) muGreenT *= iter_hit->GetEdepBirks();
245  else muGreenT *= iter_hit->GetEdep();
246 
247  double muGreen0 = 0.5 * muGreenT * FiberTransmission(dist0);
248  double muGreen1 = 0.5 * muGreenT * FiberTransmission(dist1);
249 
250  int nGreen0 = poisson.fire(muGreen0);
251  int nGreen1 = poisson.fire(muGreen1);
252 
253  muTot += muGreen0+muGreen1;
254 
255  // If running in this special mode, just pass on the mu's,
256  // only via the short route and all rounded to an integer
257  // (unfortunately). Only half the mu-ness is sent, and
258  // the transmission at dist=0 is applied (in case it is
259  // far from unity).
260  if (fNoAttenNoPoissonMode) {
261  nGreen0 = Int_t( 0.5 * muGreenT * FiberTransmission(0) + 0.5 );
262  nGreen1 = 0;
263  }
264 
265  // make photons
266  // just collect a list of times so that we can later cluster the output
267  for (Int_t iP=0; iP<nGreen0+nGreen1; iP++) {
268 
269  // time is:
270  // track start time + shift to current track position + Gaussian spread +
271  // exponential delay + green photon time-of-flight +
272  // fiber time-of-flight spread
273 
274  double dist = ( iP<nGreen0 ? dist0 : dist1 );
275  double fiberspread = 0;
276 
277  if (fApplyFiberSmearing) {
278  double weight = -1;
279  double uweight = 0;
280  while (weight<uweight) {
281  // Draw a Dt/z value then multiply by "dist" to get time (in ns)
282  double rand0(0.0), rand1(0.0), rand2(0.0);
283  histoDstrb.GetRandom( &fFiberSmearingHisto, rand0, rand1, rand2 );
284  fiberspread = dist*rand0;
285 
286  // Use fiber attenuation to deweight longer times (<==>longer distances)
287  // This is an approximation, but should be pretty close.
288  double Delta_dist = fiberspread * cn ;
289  weight = FiberTransmission(dist+Delta_dist)/FiberTransmission(dist);
290  uweight = flat.fire(0,1);
291  }
292  }
293 
294  double time = (Ta
295  + u * seglength / v
296  + gauss.fire(0,fTimeSpread)
297  + exponential.fire(fEmissionTau)
298  + dist / cn
299  + fiberspread);
300 
301  currtimes.push_back(time);
302 
303  } // iP loop end
304  } // z stepping end
305 
306  // sort times for this [track,plane,cell] hit and group photons together
307  // at whatever resolution desired to save time and (mostly) space
308  const unsigned int nTimes = currtimes.size();
309  if (nTimes) {
310 
311  currtimes.sort();
312  std::list<double>::iterator it = currtimes.begin();
313 
314  // here's our photon template
315  sim::PhotonSignal photon;
316  photon.SetPlaneCell(planeId, cellId);
317  photon.SetTrackId(trackId);
318 
319  // assemble clusters
320  double tmin = (*it);
321  double sum = tmin;
322  int count = 1;
323  it++;
324  while ( it != currtimes.end() ) {
325  if ( (*it) <= tmin + fTimeClustering){
326  // add to cluster
327  count++;
328  sum += (double)(*it);
329  }
330  else {
331  // make a signal with the average time of this cluster
332  photon.SetTimeMean(sum/count);
333  photon.SetNPhoton(count);
334  // Share truth equally. This is an approximation
335  photon.SetPoissonLambda((muTot*count)/nTimes);
336  photlist->push_back(photon);
337  // init next cluster
338  count = 1;
339  sum = tmin = (*it);
340  }
341  it++;
342  } // end while loop
343  // make final signal with leftovers
344  photon.SetTimeMean(sum/count);
345  photon.SetNPhoton(count);
346  // Share truth equally. This is an approximation
347  photon.SetPoissonLambda((muTot*count)/nTimes);
348  photlist->push_back(photon);
349  } // if (currtimes.size())
350 
351  iter_hit++;
352 
353  } // end of hit loop
354 
355  iter_hitlist++;
356 
357  } // end of hitlist loop
358 
359  mf::LogDebug("PhotonTransport") << "Read " << nHit << " FLS hits";
360  mf::LogDebug("PhotonTransport") << "Created " << photlist->size() << " photon signals";
361 
362  evt.put(std::move(photlist));
363 
364  return;
365 
366  }
367 
368  //............................................................
370  // currently, all fibers are the same
371 
372  // fAttenPars holds a series of (c,L) pairs that form the
373  // function:
374  // transmission = c0*exp(-x/L0) + c1*exp(-x/L1) + ...
375  //
376  // Thus, we should have an even number of parameters
377 
378  double answer = 0;
379 
380  int nPar = fAttenPars.size();
381  if (!nPar) {
382  std::cerr << "In PhotonTransport: No AttenPars entries. Fiber transmission set to 1!" << std::endl;
383  answer = 1;
384  }
385  else if (nPar%2) {
386  std::cerr << "In PhotonTransport: Odd number of AttenPars entries. Fiber transmission set to 1!" << std::endl;
387  answer = 1;
388  }
389  else {
390  answer = 0;
391  for (Int_t iP=0;iP<nPar;iP+=2) {
392  answer += fAttenPars[iP]*exp(-x/fAttenPars[iP+1]);
393  }
394  }
395 
396  return answer;
397  }
398 
399 
400  //............................................................
402  // just read a file from disk for now. should eventually be in the DB.
403  // use FileInPath to get the histogram file
404  struct stat sb;
405  if ( stat(fSmearingHistoFile.c_str(), &sb) !=0 ) {
406  throw cet::exception("NoSmearingHistoFile")
407  << "PhotonTransport: Unable to open "
409  << "\nCannot apply fiber smearing as requested!\n"
410  << __FILE__ << ":" << __LINE__ << "\n";
411  return 1;
412  }
413 
414  TFile f(fSmearingHistoFile.c_str());
415  TH1D *htemp = (TH1D*)f.Get("Dt_per_z")->Clone();
416  if (!htemp) {
417  std::cerr << "PhotonTransport: Histogram Dt_per_z not found in "
419  std::cerr << " Cannot apply fiber smearing as requested! " << std::endl;
420  return 1;
421  }
422  fFiberSmearingHisto = *htemp;
423  fFiberSmearingHisto.SetDirectory(0);
424  f.Close();
425  return 0;
426  }
427 
429 
430 } // namespace
void SetPoissonLambda(double l)
Definition: PhotonSignal.h:26
double HalfL() const
Definition: CellGeo.cxx:198
void SetTrackId(int t)
Definition: PhotonSignal.h:25
a class for transporting photons in a roughly realistic way
set< int >::iterator it
bool fApplyFiberSmearing
apply fiber smearing
const Var weight
iterator begin()
Definition: PtrVector.h:223
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
double fTimeClustering
groups photons closer than this many ns apart
base_engine_t & createEngine(seed_t seed)
OStream cerr
Definition: OStream.cxx:7
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
const PlaneGeo * Plane(unsigned int i) const
This slightly less simple photon transport uses a template for photon collection in time and in dista...
DEFINE_ART_MODULE(TestTMapFile)
void SetPlaneCell(int plane, int cell)
Definition: PhotonSignal.h:24
double dist
Definition: runWimpSim.h:113
const CellUniqueId & Id() const
Definition: CellGeo.h:55
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double fEmissionTau
exponential time constant for green photon emission
base_engine_t & getEngine() const
void hits()
Definition: readHits.C:15
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
unsigned int seed
Definition: runWimpSim.h:102
T get(std::string const &key) const
Definition: ParameterSet.h:231
double fFiberIndexOfRefraction
n for fiber core
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
iterator end()
Definition: PtrVector.h:237
void SetNPhoton(int npe)
Definition: PhotonSignal.h:23
static constexpr Double_t gauss
Definition: Munits.h:360
std::vector< double > fAttenPars
from Leon&#39;s fit
data_t::const_iterator const_iterator
Definition: PtrVector.h:61
double GetPigtail(const CellUniqueId &id) const
Return length of fiber in cm from end of cell to apd.
unsigned int GetRandomNumberSeed()
Definition: Simulation.cxx:13
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string fSmearingHistoFile
relative path to the smearing histo file
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
double fPhotonsPerMeV
The length (in cm) of the fibre from the end of the cell to the APD.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
double fTimeSpread
spread in time due to everything
void GetRandom(TH1 *histo, double &x, double &y, double &z)
Definition: RandHisto.cxx:8
PhotonTransport(fhicl::ParameterSet const &pset)
std::string fFilterModuleLabel
module that filtered the fls hits
bool fApplyBirks
Use EdepBirks from FLSHits instead of calculating here.
bool fNoAttenNoPoissonMode
run special no-attenuation, no-Poisson-fluct mode?
double fQuantumEff
APD quantum efficiency.
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
double fStep
step size for walking along the track (cm)
void SetTimeMean(double t)
Definition: PhotonSignal.h:22
Definition: fwd.h:28
Encapsulate the geometry of one entire detector (near, far, ndos)
double fCollectionEff
fiber collection efficiency