PhotonTransport_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   This slightly less simple photon transport doesn't propagate
///          individual photons (thus "simple"), but the model is gaining
///          some realism.
///
/// \author  rbpatter@caltech.edu
////////////////////////////////////////////////////////////////////////

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Persistency/Common/Ptr.h"

#include "cetlib/search_path.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandExponential.h"
#include "CLHEP/Random/RandPoissonQ.h"
#include "CLHEP/Random/RandGaussQ.h"

// NOvA includes
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/PhotonSignal.h"
#include "Simulation/Simulation.h"
#include "Utilities/func/MathUtil.h"
#include "Utilities/func/RandHisto.h"

// ROOT includes
#include "TFile.h"
#include "TH1.h"
#include "TMath.h"

#include <cmath>
#include <sys/stat.h>

namespace photrans{
  /// a class for transporting photons in a roughly realistic way
  class PhotonTransport: public art::EDProducer
  {
   public:
     explicit PhotonTransport(fhicl::ParameterSet const &pset);
     virtual ~PhotonTransport();

     void produce(art::Event& evt);

   protected:
     double FiberTransmission(double x);
     int LoadSmearingHisto();
     /// The length (in cm) of the fibre from the end of the cell to the APD

   private:

     double              fPhotonsPerMeV;          ///< blue photon production scale
     double              fQuantumEff;             ///< APD quantum efficiency
     double              fCollectionEff;          ///< fiber collection efficiency
     double              fFiberIndexOfRefraction; ///< n for fiber core
     double              fEmissionTau;            ///< exponential time constant for green photon emission
     double              fTimeSpread;             ///< spread in time due to everything
     std::vector<double> fAttenPars;              ///< from Leon's fit
     double              fTimeClustering;         ///< groups photons closer than this many ns apart
     double              fStep;                   ///< step size for walking along the track (cm)
     bool                fApplyFiberSmearing;     ///< apply fiber smearing
     bool                fNoAttenNoPoissonMode;   ///< run special no-attenuation, no-Poisson-fluct mode?
     bool                fApplyBirks;             ///< Use EdepBirks from FLSHits instead of calculating here

     TH1D                fFiberSmearingHisto;
     int                 msglevel;

     std::string         fFilterModuleLabel;      ///< module that filtered the fls hits
     std::string         fSmearingHistoFile;      ///< relative path to the smearing histo file
   };

   //............................................................
  PhotonTransport::PhotonTransport(fhicl::ParameterSet const &pset)
    : fPhotonsPerMeV         (pset.get< double              >("PhotonsPerMeV")) 
    , fQuantumEff            (pset.get< double              >("QuantumEff")) 
    , fCollectionEff         (pset.get< double              >("CollectionEff")) 
    , fFiberIndexOfRefraction(pset.get< double              >("FiberIndexOfRefraction"))
    , fEmissionTau           (pset.get< double              >("EmissionTau")) 
    , fTimeSpread            (pset.get< double              >("TimeSpread")) 
    , fAttenPars             (pset.get< std::vector<double> >("AttenPars"))
    , fTimeClustering        (pset.get< double              >("TimeClustering")) 
    , fStep                  (pset.get< double              >("Step")) 
    , fApplyFiberSmearing    (pset.get< bool                >("ApplyFiberSmearing")) 
    , fNoAttenNoPoissonMode  (pset.get< bool                >("NoAttenNoPoissonMode")) 
    , fApplyBirks            (pset.get< bool                >("ApplyBirks")) 
    , msglevel               (pset.get< int                 >("MessageLevel"))
    , fFilterModuleLabel     (pset.get< std::string         >("FilterModuleLabel"))
  {
    // get the random number seed, use a random default if not specified    
    // in the configuration file.  
    unsigned int seed = pset.get< unsigned int >("Seed", sim::GetRandomNumberSeed());
  
    createEngine( seed );

    produces< std::vector<sim::PhotonSignal> >();

    // constructor decides if initialized value is a path or an environment variable
    cet::search_path sp("FW_SEARCH_PATH");

    sp.find_file(pset.get< std::string >("SmearingHistoFile"), fSmearingHistoFile);

    // Open the GDML file
    struct stat sb;
    if ( fSmearingHistoFile.empty() || stat(fSmearingHistoFile.c_str(), &sb)!=0 ) {
      // failed to resolve the file name
      throw cet::exception("NoSmearingHisto") 
        << "smearing histogram file " << fSmearingHistoFile << " not found!\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }

  }

  //............................................................
  PhotonTransport::~PhotonTransport()
  {
  }

  //............................................................
  void PhotonTransport::produce(art::Event& evt) 
  {

    if (fApplyFiberSmearing&&!fFiberSmearingHisto.Integral()) {
      // Load up the distribution of Delta(t) (per z-distance).  Note:
      // this distribution was generated using an offline fiber MC.  The
      // exact distribution depends on the choice of green photon starting
      // positions within the fiber.  The version stamped Oct 25 2008
      // assumes that blue photons originate 2 cm away with a +/- 45deg
      // angular spread in Z (arbitrary!) and with a 0.03 cm absorption
      // length in the fiber.
      if ( LoadSmearingHisto() ) {
        throw cet::exception("NoSmearingHisto") 
          << "can't load smearing histo\n"
          << __FILE__ << ":" << __LINE__ << "\n";
      }
    }

    art::ServiceHandle<geo::Geometry> fGeo;
 
    // get the RandomNumberGenerator service
    art::ServiceHandle<art::RandomNumberGenerator> rng;
    CLHEP::HepRandomEngine &engine = rng->getEngine();
    CLHEP::RandPoissonQ     poisson(engine);
    CLHEP::RandGaussQ       gauss(engine);
    CLHEP::RandExponential  exponential(engine);
    CLHEP::RandFlat         flat(engine);
    RandHisto               histoDstrb(engine);


    // get list of hits
    art::Handle< std::vector<sim::FLSHitList> > flslistHandle;
    evt.getByLabel(fFilterModuleLabel, flslistHandle);
    // make a art::PtrVector of the FLSHitLists
    art::PtrVector<sim::FLSHitList> hitlist;
    for(unsigned int i = 0; i < flslistHandle->size(); ++i){
      art::Ptr<sim::FLSHitList> flshl(flslistHandle, i);
      hitlist.push_back(flshl);
    }

    // make output vector
    std::unique_ptr<std::vector<sim::PhotonSignal> > photlist(new std::vector<sim::PhotonSignal>);

    const double MeVPerGeV = 1000.;
    const double c = TMath::Ccgs()*1e-9; // c in cm/ns
    const double cn = c/fFiberIndexOfRefraction;

    const double green_scale = MeVPerGeV * fPhotonsPerMeV * fCollectionEff * fQuantumEff;

    int nHit = 0;

    /// Loop over multiple lists of hits in event
    art::PtrVector<sim::FLSHitList>::const_iterator iter_hitlist = hitlist.begin();
    while ( iter_hitlist != hitlist.end() ) {
    
      /// loop over multiple hits in current hitlist
      std::vector<sim::FLSHit> hits = (*iter_hitlist)->fHits;

      std::vector<sim::FLSHit>::iterator iter_hit = hits.begin();
      while ( iter_hit != hits.end() ) {

        nHit++;

        // for cleaner-looking code below...
        double Xa       = (*iter_hit).GetEntryX();
        double Ya       = (*iter_hit).GetEntryY();
        double Za       = (*iter_hit).GetEntryZ();
        double Xb       = (*iter_hit).GetExitX();
        double Yb       = (*iter_hit).GetExitY();
        double Zb       = (*iter_hit).GetExitZ();
        //      double segtime  = (*iter_hit).fT;

        double Ta       = (*iter_hit).GetEntryT();
        double Tb       = (*iter_hit).GetExitT();

        int   planeId  = (*iter_hit).GetPlaneID();
        int   cellId   = (*iter_hit).GetCellID();
        int   trackId  = (*iter_hit).GetTrackID();


        double seglength = util::pythag(Xb - Xa, Yb - Ya, Zb - Za);
        
        //Calculate the velocity:
        double v = ((Tb-Ta>0) ? seglength/(Tb-Ta) : c);
        if (v>c) v=c;


        // cell half-width
        double DZ = fGeo->Plane((unsigned int)(planeId))->Cell ((unsigned int)(cellId ))->HalfL();

        // walk along track segment
        std::list<double> currtimes;
        currtimes.clear();
        double muTot = 0;
        for (double currS=0; currS<=seglength; currS+=fStep) {

          // last step is truncated
          double dS = fStep;
          if ( dS > seglength-currS ) dS = seglength-currS;

          // current fraction along segment, at step midpoint
          double u = (seglength) ? (currS + dS / 2) / seglength : 0;

          // current location (only Z component is needed currently)
          double currZ = Za + (Zb - Za) * u;

          // green photon travel distances for this step
          // "pigtail" is the additional distance outside the cell
          // --> short path and long path through fiber:
          double pigtail = fGeo->GetPigtail(fGeo->Plane(planeId)->Cell(cellId)->Id());
          double dist0 =   DZ - currZ + pigtail;
          double dist1 = 3*DZ + currZ + pigtail;
  
          // Mean number of green photons for this step
          // If segment has length zero, dump all the energy right here.  The upstream
          // code should be examined to figure out why seglength is sometimes zero.
          double muGreenT = green_scale * (seglength ? dS / seglength : 1);
          if (fApplyBirks) muGreenT *= iter_hit->GetEdepBirks();
          else muGreenT *= iter_hit->GetEdep();

          double muGreen0 = 0.5 * muGreenT * FiberTransmission(dist0);
          double muGreen1 = 0.5 * muGreenT * FiberTransmission(dist1);
  
          int nGreen0 = poisson.fire(muGreen0);
          int nGreen1 = poisson.fire(muGreen1);

          muTot += muGreen0+muGreen1;
  
          // If running in this special mode, just pass on the mu's,
          // only via the short route and all rounded to an integer
          // (unfortunately).  Only half the mu-ness is sent, and
          // the transmission at dist=0 is applied (in case it is
          // far from unity).
          if (fNoAttenNoPoissonMode) {
            nGreen0 = Int_t( 0.5 * muGreenT * FiberTransmission(0) + 0.5 );
            nGreen1 = 0;
          }

          // make photons
          // just collect a list of times so that we can later cluster the output
          for (Int_t iP=0; iP<nGreen0+nGreen1; iP++) {
  
            // time is:
            //   track start time + shift to current track position + Gaussian spread +
            //   exponential delay + green photon time-of-flight +
            //   fiber time-of-flight spread

            double dist = ( iP<nGreen0 ? dist0 : dist1 );
            double fiberspread = 0;

            if (fApplyFiberSmearing) {
              double weight = -1;
              double uweight = 0;
              while (weight<uweight) {
                // Draw a Dt/z value then multiply by "dist" to get time (in ns)
                double rand0(0.0), rand1(0.0), rand2(0.0);
                histoDstrb.GetRandom( &fFiberSmearingHisto, rand0, rand1, rand2 );
                fiberspread = dist*rand0;
              
                // Use fiber attenuation to deweight longer times (<==>longer distances)
                // This is an approximation, but should be pretty close.
                double Delta_dist = fiberspread * cn ;
                weight = FiberTransmission(dist+Delta_dist)/FiberTransmission(dist);
                uweight = flat.fire(0,1);
              }
            }

            double time = (Ta 
                           + u * seglength / v 
                           + gauss.fire(0,fTimeSpread)
                           + exponential.fire(fEmissionTau) 
                           + dist / cn 
                           + fiberspread);
          
            currtimes.push_back(time);
          
          } // iP loop end
        } // z stepping end
      
        // sort times for this [track,plane,cell] hit and group photons together
        // at whatever resolution desired to save time and (mostly) space
        const unsigned int nTimes = currtimes.size();
        if (nTimes) {

          currtimes.sort();
          std::list<double>::iterator it = currtimes.begin();

          // here's our photon template
          sim::PhotonSignal photon;
          photon.SetPlaneCell(planeId, cellId);
          photon.SetTrackId(trackId);

          // assemble clusters
          double tmin = (*it);
          double sum = tmin;
          int count = 1;
          it++;
          while ( it != currtimes.end() ) {
            if ( (*it) <= tmin + fTimeClustering){
              // add to cluster
              count++;
              sum += (double)(*it);
            }
            else {
              // make a signal with the average time of this cluster
              photon.SetTimeMean(sum/count);
              photon.SetNPhoton(count);
              // Share truth equally. This is an approximation
              photon.SetPoissonLambda((muTot*count)/nTimes);
              photlist->push_back(photon);
              // init next cluster
              count = 1;
              sum = tmin = (*it);
            }
            it++;
          } // end while loop
          // make final signal with leftovers
          photon.SetTimeMean(sum/count);
          photon.SetNPhoton(count);
          // Share truth equally. This is an approximation
          photon.SetPoissonLambda((muTot*count)/nTimes);
          photlist->push_back(photon);
        } // if (currtimes.size())

        iter_hit++;

      } // end of hit loop

      iter_hitlist++;

    } // end of hitlist loop

    mf::LogDebug("PhotonTransport") << "Read " << nHit << " FLS hits";
    mf::LogDebug("PhotonTransport") << "Created " << photlist->size() << " photon signals";
  
    evt.put(std::move(photlist));

    return;

  }

  //............................................................
  double PhotonTransport::FiberTransmission(double x) {
    // currently, all fibers are the same

    // fAttenPars holds a series of (c,L) pairs that form the
    // function:
    //    transmission = c0*exp(-x/L0) + c1*exp(-x/L1) + ...
    //
    // Thus, we should have an even number of parameters

    double answer = 0;

    int nPar = fAttenPars.size();
    if (!nPar) {
      std::cerr << "In PhotonTransport: No AttenPars entries.  Fiber transmission set to 1!" << std::endl;
      answer = 1;
    }
    else if (nPar%2) {
      std::cerr << "In PhotonTransport: Odd number of AttenPars entries.  Fiber transmission set to 1!" << std::endl;
      answer = 1;
    }
    else {
      answer = 0;    
      for (Int_t iP=0;iP<nPar;iP+=2) {
        answer += fAttenPars[iP]*exp(-x/fAttenPars[iP+1]);
      }
    }

    return answer;
  }


  //............................................................
  int PhotonTransport::LoadSmearingHisto() {
    // just read a file from disk for now.  should eventually be in the DB.
    // use FileInPath to get the histogram file
    struct stat sb;
    if ( stat(fSmearingHistoFile.c_str(), &sb) !=0 ) {
      throw cet::exception("NoSmearingHistoFile")
        << "PhotonTransport: Unable to open " 
        << fSmearingHistoFile
        << "\nCannot apply fiber smearing as requested!\n"
        << __FILE__ << ":" << __LINE__ << "\n";
      return 1;
    }

    TFile f(fSmearingHistoFile.c_str());
    TH1D *htemp = (TH1D*)f.Get("Dt_per_z")->Clone();
    if (!htemp) {
      std::cerr << "PhotonTransport: Histogram Dt_per_z not found in " 
                << fSmearingHistoFile << std::endl;
      std::cerr << "   Cannot apply fiber smearing as requested!  " << std::endl;
      return 1;
    }
    fFiberSmearingHisto = *htemp;
    fFiberSmearingHisto.SetDirectory(0);
    f.Close();
    return 0;
  }

  DEFINE_ART_MODULE(PhotonTransport)

} // namespace