ImprovedTransport_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////
/// \brief   This slightly less simple photon transport uses a template for photon collection in 
///          time and in distance along the fiber.
///
/// \author  aurisano@fnal.gov
/// \date
////////////////////////////////////////////////////////////////////////
//Framework Includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "cetlib/search_path.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/FLSHit.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/PhotonSignal.h"
#include "Simulation/Simulation.h"
#include "Utilities/func/MathUtil.h"
#include "Utilities/func/RandHisto.h"
#include "PhotonTransport/BrightnessLevel.h"
#include "RunHistory/service/RunHistoryService.h"

//stdlib Includes
#include <string>
#include <list>
#include <vector>
#include <cmath>

//ROOT Includes
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile.h"
#include "TFile.h"
#include "TMath.h"
#include "TRandom3.h"

//System Includes
#include <sys/stat.h>

///photon transport simulation
namespace photrans
{
  /// a struct to hold collection rate info
  struct RateInfo
  {
    double rate;
    double z;
    double zWidth;
    TH1D   timePDF;
  };

  /// a class for transporting photons in a roughly realistic way
  class ImprovedTransport : public art::EDProducer
  {
  public:
    explicit ImprovedTransport(fhicl::ParameterSet const &pset);
    virtual ~ImprovedTransport();
    
    void beginRun(art::Run& run);
    void produce(art::Event& evt);
    
  protected:
    double FiberTransmission(double x);
    double ScintTime(double dEdx);
    bool   LoadHistos();
    void   StepAlongHit(sim::FLSHit& hit);
    /// \param meanPE   The photons entering the fiber
    /// \param stepTime When they entered
    /// \param dist     How far to the readout
    /// \param lambda   Accumulates total mean expected at readout into here
    void   CreatePhotoElectrons(double meanPE, double stepTime, double dist, double dEdx, bool isScint, RateInfo& info, double& lambda);
    void   ClusterPhotoElectrons(int planeId, int cellId, int trackId);

    /// Get correction for the varying efficiency across the cross
    /// section of of the cell caused by fiber position.
    double GetPosCorr( const sim::FLSHit& hit );

  private:
    //fcl parameters
    std::string         fFilterModuleLabel;      ///< module that filtered the fls hits
    fhicl::ParameterSet fPSetNDOS, fPSetND, fPSetFD, fPSetTB;

    double              fPhotonsPerMeV;          ///< blue photon production scale
    double              fXViewFactor;            ///< scale factor for the x-view
    double              fYViewFactor;            ///< scale factor for the y-view
    double              fLightLevelScale;        ///< shift overall light yield by multiplicative factor
    double              fCerenkovEff;            ///< Absorption and re-emission efficiency for Cerenkov photons
    double              fQuantumEff;             ///< APD quantum efficiency
    double              fFiberIndexOfRefraction; ///< n for fiber core
    double              fEmissionTau;            ///< exponential time constant for green photon emission 
    std::vector<double> fAttenPars;              ///< from Leon's fit
    double              fTimeClustering;         ///< groups photons closer than this many ns apart
    bool                fWriteEmptySignals;      ///< Write an empty PhotonSignal if there's 0 PE
    double              fStep;                   ///< step size for walking along the track (cm)
    bool                fApplyBirks;             ///< Use EdepBirks from FLSHits instead of calculating here
    std::string         fCollectionRateFile;     ///< Relative path to the collection rate histo file
    bool                fApplyFiberSmearing;     ///< Turn on fiber smearing
    std::string         fSmearingHistoFile;      ///< Relative path to the fiber smearing histo file
    std::string         fBrightnessMode;         ///< correct the fiber brightness by cell, by module, or not at all
    std::string         fBrightnessFile;         ///< Relative path to the brightness correction file
    bool                fApplyTweak;             ///< Turn on the data driven attenuation correction for each view
    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
    int                 fDriftReference;         ///< drift reference time (ie. where the slope is fixed to 1)
    double              fDriftGradient;          ///< gradient for light level detector aging systematic [fraction/year]

    bool                fUseScintTime;           ///< Turn on sampling from PDF of scintillator decay time, see arXiv:1704.02291 for more details
    double              fCherenkovTau;           ///< Time constant for directly excited PC to PPO scintillation
    double              fScintBetaDEDX;          ///< Approximate dE/dx of betas used to determine PDF
    std::vector<double> fScintBetaProb;          ///< Cumulative probability for each exponential component for betas
    std::vector<double> fScintBetaTau;           ///< Time constants of each exponential component for betas
    double              fScintAlphaDEDX;         ///< Approximate dE/dx of alphas used to determine PDF
    std::vector<double> fScintAlphaProb;         ///< Cumulative probability for each exponential component for alphas
    std::vector<double> fScintAlphaTau;          ///< Time constants of each exponential component for alphas

    /// Relative path to the cell efficiency maps.  If left empty, do not
    /// apply this correction.
    std::string         fPosCorrFile;

    int                 msglevel;
    
    //non-fcl members
    double                                fFiberSpeedOfLight;    ///< Speed of light in the fiber
    TH2D                                  fCollectionRateHisto;  ///< Collection rate for given dT vs. dZ
    TH1D                                  fFiberSmearingHisto;   ///< Smearing histo for photons bouncing around in fiber
    TH1D                                  fTweak_X;              ///< A data/mc correction for the X-View
    TH1D                                  fTweak_Y;              ///< A data/mc correction for the Y-View 
    TH1D                                  fTweak_XMC;            ///< A data/mc correction for the X-View in the muon catcher
    TH1D                                  fTweak_YMC;            ///< A data/mc correction for the Y-View in the muon catcher

    std::vector< std::vector< double > >  fBrightnessMap;        ///< The brightness correction by cell
    //Add by Pengfei -- begin
    std::string fBrightnessMapName;
    std::string fBrightnessValueName;
    std::string fBrightnessValueMapName; //Added by RJN
    BrightnessLevel*                      fBrightnessLevel;
    //Add by Pengfei -- end

    // Efficiency maps, read from fPosCorrFile and applied in GetPosCorr().
    //
    // First the horizontals, middle 14 cells of each extrusion, which are
    // the bigger cells, with the fibers randomly on the bottom, and with
    // the fiber together on the bottom.
    std::vector<TH2D *> fBigHorizMap;
    std::vector<TH2D *> fBigHorizClumpMap;

    // Outer two cells of horizontals, the small ones.
    std::vector<TH2D *> fSmallHorizMap;
    std::vector<TH2D *> fSmallHorizClumpMap;

    // Same for the top cell in modules with air bubbles.
    std::vector<TH2D *> fSmallHorizBubbleMap;
    std::vector<TH2D *> fSmallHorizBubbleClumpMap;

    // Now the verticals, with fibers randomly distribted, and with fibers
    // in the corners towards the bottom of cells, because of slack.
    std::vector<TH2D *> fBigVertMap;
    std::vector<TH2D *> fBigVertCornerMap;

    // Same for the small verticals.
    std::vector<TH2D *> fSmallVertMap;
    std::vector<TH2D *> fSmallVertCornerMap;

    // Index for the map being used now.  Changes for each hit.
    // This is for saving in the output so we know what map we used.
    // We certainly aren't going to use more than 65535 maps, and we
    // want to save this number for every hit, so use a 16 bit integer.
    uint16_t fPosCorrMapIndex = -1;

    // Compute the global index of the fiber efficiency map we used
    // and store it in fPosCorrMapIndex.
    void SetPosCorrMapIndex(const std::vector<TH2D *> * const maps,
                            const unsigned int i);

    /// Reads in the efficiency maps from fPosCorrFile.
    void LoadPosCorrHists();

    // Used to choose efficiency cell maps
    TRandom3 fPosCorrHasher;
    uint64_t fPosCorrEvtId;

    art::ServiceHandle<geo::Geometry>     fGeo;                  ///< Handle to the geometry service
    std::list< double > fPETimes;                                ///< List of photo-electron times (each representing one observed PE)
    std::list< RateInfo >                 fRateInfo;             ///< List of RateInfo object
   
    std::unique_ptr<std::vector<sim::PhotonSignal> >   fSignals; ///< Produced vector of PhotonSignals
    CLHEP::RandPoissonQ*    fPoisson;
    CLHEP::RandFlat*        fFlat;
    RandHisto*              fRandHisto;
    CLHEP::RandGaussQ*      fGauss;
    CLHEP::RandExponential* fExp;
    bool                    fFirstTime;

    // Histo for testing
    TH1D* fScintTime;
  };

  //............................................................
  ImprovedTransport::ImprovedTransport(fhicl::ParameterSet const &pset)
    : fFilterModuleLabel     (pset.get< std::string>("FilterModuleLabel"))
    , fPSetNDOS              (pset.get<fhicl::ParameterSet>("ndos"))
    , fPSetND                (pset.get<fhicl::ParameterSet>("nd"))
    , fPSetFD                 (pset.get<fhicl::ParameterSet>("fd"))
    , fPSetTB                 (pset.get<fhicl::ParameterSet>("tb"))
    , fSignals               (new std::vector<sim::PhotonSignal>)
    , fFirstTime(true)
  {
    // 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 );

    // get the RandomNumberGenerator service
    art::ServiceHandle<art::RandomNumberGenerator> rng;
    CLHEP::HepRandomEngine &engine = rng->getEngine();
    fPoisson = new CLHEP::RandPoissonQ(engine);
    fFlat = new CLHEP::RandFlat(engine);
    fRandHisto = new RandHisto(engine);
    fGauss = new CLHEP::RandGaussQ(engine);
    fExp   = new CLHEP::RandExponential(engine);
    
    produces< std::vector<sim::PhotonSignal> >();
    //Add by Pengfei -- begin
    produces<BrightnessLevel, art::InRun> ();
    //Add by Pengfei -- end
  }
  
  //............................................................
  ImprovedTransport::~ImprovedTransport()
  {
    delete fPoisson;
    delete fFlat;
    delete fRandHisto;
    delete fGauss;
  }

  void ImprovedTransport::beginRun(art::Run& run /*run*/)
  {
    if (fFirstTime) {
      art::ServiceHandle<art::TFileService> tfs;
      fScintTime = tfs->make<TH1D>("hScintTime", ";time;photons", 1000, 0, 5000);

      // Figure out which parameter set we're using
      fhicl::ParameterSet pset;
      art::ServiceHandle<geo::Geometry> geom;
      switch(geom->DetId()){
      case novadaq::cnv::kNDOS:
        pset = fPSetNDOS;
        break;
      case novadaq::cnv::kNEARDET:
        pset = fPSetND;
        break;
      case novadaq::cnv::kFARDET:
        pset = fPSetFD;
        break;
      case novadaq::cnv::kTESTBEAM:
        pset = fPSetTB;
        break;
      default:
        assert(0 && "Unknown detector");
      }

      fPhotonsPerMeV          = pset.get< double              >("PhotonsPerMeV");
      fXViewFactor            = pset.get< double              >("XViewFactor");
      fYViewFactor            = pset.get< double              >("YViewFactor");
      fLightLevelScale        = pset.get< double              >("LightLevelScale");
      fCerenkovEff            = pset.get< double              >("CerenkovEff");
      fQuantumEff             = pset.get< double              >("QuantumEff");
      fFiberIndexOfRefraction = pset.get< double              >("FiberIndexOfRefraction");
      fEmissionTau            = pset.get< double              >("EmissionTau");
      fAttenPars              = pset.get< std::vector<double> >("AttenPars");
      fTimeClustering         = pset.get< double              >("TimeClustering"); 
      fWriteEmptySignals      = pset.get< bool                >("WriteEmptySignals");
      fStep                   = pset.get< double              >("Step");
      fApplyBirks             = pset.get< bool                >("ApplyBirks"); 
      fCollectionRateFile     = pset.get< std::string         >("CollectionRateFile");
      fSmearingHistoFile      = pset.get< std::string         >("SmearingHistoFile");
      fBrightnessFile         = pset.get< std::string         >("BrightnessFile");
      fBrightnessMode         = pset.get< std::string         >("BrightnessMode");
      fUseScintTime           = pset.get< bool                >("UseScintTime");
      fCherenkovTau           = pset.get< double              >("CherenkovTau");
      fScintBetaDEDX          = pset.get< double              >("ScintBetaDEDX");
      fScintBetaProb          = pset.get< std::vector<double> >("ScintBetaProb");
      fScintBetaTau           = pset.get< std::vector<double> >("ScintBetaTau");
      fScintAlphaDEDX         = pset.get< double              >("ScintAlphaDEDX");
      fScintAlphaProb         = pset.get< std::vector<double> >("ScintAlphaProb");
      fScintAlphaTau          = pset.get< std::vector<double> >("ScintAlphaTau");
      //Add by Pengfei -- begin
      fBrightnessMapName      = pset.get< std::string         >("BrightnessMapName");
      fBrightnessValueName    = pset.get< std::string         >("BrightnessValueName");
      fBrightnessValueMapName = pset.get< std::string         >("BrightnessValueMapName");
      //Add by Pengfei -- end
      fApplyTweak             = pset.get< bool                >("ApplyTweak");
      fTweakFile              = pset.get< std::string         >("TweakFile");
      fSimulateDrift          = pset.get< bool                >("SimulateDrift");
      fDriftReference         = pset.get< int                 >("DriftReference");
      fDriftGradient          = pset.get< double              >("DriftGradient");
      fPosCorrFile            = pset.get< std::string         >("PosCorrFile");
      fApplyFiberSmearing     = pset.get< bool                >("ApplyFiberSmearing");    
      msglevel                = pset.get< int                 >("MessageLevel");

      // constructor decides if initialized value is a path or an environment variable
      cet::search_path sp("FW_SEARCH_PATH");
      
      const double c = TMath::Ccgs()*1e-9; // c in cm/ns
      fFiberSpeedOfLight = c/fFiberIndexOfRefraction;
      
      sp.find_file(pset.get< std::string >("CollectionRateFile"), fCollectionRateFile);
      // Check the collection rate file
      struct stat sb;
      if ( fCollectionRateFile.empty() || stat(fCollectionRateFile.c_str(), &sb)!=0 ){
        // failed to resolve the file name
        throw cet::exception("NoCollectionRateFile") << "collection rate file " << fCollectionRateFile << " not found!\n" << __FILE__ << ":" << __LINE__ << "\n";
      }
      
      sp.find_file(pset.get< std::string >("SmearingHistoFile"), fSmearingHistoFile);
      if ( fSmearingHistoFile.empty() || stat(fSmearingHistoFile.c_str(), &sb)!=0 ){
        // failed to resolve the file name
        throw cet::exception("NoSmearingHistoFile") << "smearing histo file " << fSmearingHistoFile << " not found!\n" << __FILE__ << ":" << __LINE__ << "\n";
      }

      if (fBrightnessMode != "None") {
        fBrightnessFile  = pset.get< std::string >("BrightnessFile");
        sp.find_file(pset.get< std::string >("BrightnessFile"), fBrightnessFile);
        if ( fBrightnessFile.empty() || stat(fBrightnessFile.c_str(), &sb)!=0 ){
          // failed to resolve the file name      
          throw cet::exception("NoBrightnessFile") << "fiber brightness histo file " << fBrightnessFile << " not found!\n" << __FILE__ << ":" << __LINE__ << "\n";
        }
      }
      
      if (fApplyTweak) {
        sp.find_file(pset.get< std::string >("TweakFile"), fTweakFile);
        if ( fTweakFile.empty() || stat(fTweakFile.c_str(), &sb)!=0 ) {
          // failed to resolve the file name                                                                                                     
          throw cet::exception("NoTweakFile") << "tweak file " << fTweakFile << " not found!\n" << __FILE__ << ":" << __LINE__ << "\n";
        }
      }
      
      sp.find_file(pset.get< std::string >("PosCorrFile"), fPosCorrFile);
      if(fPosCorrFile.empty()){
        if(msglevel >= 10)
          std::cout << "Fiber position efficiency correction disabled\n";
      }
      else{
        if(stat(fPosCorrFile.c_str(), &sb) != 0)
          throw cet::exception("NoPosCorrFile") << "Efficiency map file "
            << fPosCorrFile << " not found!\n" << __FILE__ << ":" << __LINE__
            << "\n";
      }

      if ( !LoadHistos() ) 
        {
          throw cet::exception("NoCollectionRateHisto") 
          << "can't load histos\n"
          << __FILE__ << ":" << __LINE__ << "\n";
        return;
        }

        //Add by Pengfei -- begin
      if (fBrightnessMode == "ByCell"){
        std::unique_ptr<BrightnessLevel> pfb(new BrightnessLevel(fBrightnessFile, fBrightnessMapName, fBrightnessValueName,fBrightnessValueMapName));
        run.put(std::move(pfb));
      }
      else if (fBrightnessMode == "ByBin"){
        std::unique_ptr<BrightnessLevel> pfb(new BrightnessLevel(fBrightnessFile, fBrightnessMapName, fBrightnessValueName));
        run.put(std::move(pfb));
      }
        //Add by Pengfei -- end

      //don't reinitialize after the first time through
      fFirstTime = false;
    }
  }


  //............................................................
  void ImprovedTransport::produce(art::Event& evt) 
  {
    // get list of hits
    art::Handle< std::vector<sim::FLSHitList> > flslistHandle;
    evt.getByLabel(fFilterModuleLabel, flslistHandle);
    
    //clear output vector
    fSignals->clear();
    
    // Reasonably unique event number (unique if run, subrun < 65536)
    // for use with cell efficiency maps.
    fPosCorrEvtId = ((uint64_t)evt.run()    << 48)
                   +((uint64_t)evt.subRun() << 32)
                   +((uint64_t)evt.event()       );

    for (unsigned int ilist = 0; ilist < flslistHandle->size(); ++ilist) {
      const std::vector<sim::FLSHit> hits = flslistHandle->at(ilist).fHits;
      for (unsigned int ihit = 0; ihit < hits.size(); ++ihit) {
        sim::FLSHit hit = hits[ihit];
        //step along the hit and fill fPETimes (clustering done at end of StepAlongHit)
        StepAlongHit( hit );
      }
    }
    
    // make a new collection into which we will put the sim::PhotonSignals.  That new collection
    // will then be put into the event.  If we tried to put fSignals in the event, the std::move
    // call changes ownership of the vector and the next time we try to access fSignals would
    // cause a seg fault
    std::unique_ptr<std::vector<sim::PhotonSignal>> pscol(new std::vector<sim::PhotonSignal>);
    pscol->swap(*fSignals);
    
    //put vector<sim::PhotonSignal> into the event
    evt.put(std::move( pscol ));
    
    return;
  }

  // Helper function for GetPosCorr().  Given an efficiency map and two
  // points, return the mean efficiency along the line defined by the points.
  static double line_efficiency(TH2D * themap,
                                const double Xa, const double Ya,
                                const double Xb, const double Yb)
  {
    double sum = 0, dem = 0;

    // This is a dumb method which just takes little steps across and finds
    // the average, but the performance is fine, so we'll go with it.
    const double stepsize = 0.01; // cm

    const double dy = Yb - Ya;
    const double dx = Xb - Xa;
    const double seglen = hypot(dx, dy);
    const int nstep = std::max(1, int(seglen/stepsize));
    const double thisstepsize = seglen/nstep;
    const double thisstepstart = thisstepsize/2;

    for(int step = 0; step < nstep; step++){
      const double s = thisstepstart + step*thisstepsize;
      const double sfrac = seglen == 0? 0: s/seglen;

      // Flip the x coordinate here because of an artifact in how I generated
      // the efficiency maps.  It only matters for horizontals where we want
      // the fibers to be at the bottom (and air to be at the top).  Local +x
      // is global -y, where fibers should fall, but I generated the maps the
      // other way around.
      const double x = -(Xa + sfrac*dx);

      const double y = Ya + sfrac*dy;

      double eff = themap->GetBinContent(themap->FindBin(x, y));

      // There are some zero bins in the corners since the bins are square and
      // were evaluated at the bin centers.  Use any non-zero adjacent bin in
      // this case.
      int binx = 0, biny = 0;
      if(eff <= 0){
        int binz = 0;
        themap->GetBinXYZ(themap->FindBin(x, y), binx, biny, binz);
        eff = themap->GetBinContent(binx+1, biny);
      }
      if(eff <= 0) eff = themap->GetBinContent(binx-1, biny);
      if(eff <= 0) eff = themap->GetBinContent(binx, biny+1);
      if(eff <= 0) eff = themap->GetBinContent(binx, biny-1);

      if(eff <= 0) std::cerr << "GetPosCorr should not happen: " << eff
          << " efficiency with " << themap->GetName()
          << " point " << -x << " " << y
          << " from path " << Xa << " " << Ya << " " << Xb << " " << Yb
          << " length " << seglen << "\n";

      // Handle the case of zero step size, even though in all cases I've seen,
      // these particles also have zero energy, so it doesn't matter.
      sum += thisstepsize == 0? eff: thisstepsize*eff;
      dem += thisstepsize;
    }

    // If the length is zero, return the single efficiency stored in "sum".  If
    // somehow denominator is zero, prevent a NaN, but that shouldn't happen.
    return seglen == 0? sum: dem == 0? 0: sum/dem;
  }

  void load_hset(TFile * f, const char * name, const char * humanname,
                 std::vector<TH2D *> & hold, const bool verbose)
  {
    for(int n = 1; ; n++){
      TH2D * h = dynamic_cast<TH2D *>(f->Get(Form("%s_%d", name, n)));
      if(h == NULL) break;
      h->SetDirectory(NULL);
      hold.push_back(h);
    }

    if(verbose) std::cout << "ImprovedTransport: using "
                          << hold.size() << " " << humanname << "\n";

    if(hold.empty())
      throw cet::exception("NoPosCorrHists")
        << "ImprovedTransport: No " << humanname << "\n";
  }

  void ImprovedTransport::LoadPosCorrHists()
  {
    TFile * f = new TFile(fPosCorrFile.c_str());
    if(f == NULL || f->IsZombie())
      throw cet::exception("NoPosCorrHistFile")
        << "ImprovedTransport: Unable to open " 
        << fPosCorrFile << "\n" << __FILE__ << ":" << __LINE__ << "\n";

    const bool verbose = msglevel >= 10;

    load_hset(f, "big_horiz",
                 "big horizontal cell maps",
                 fBigHorizMap, verbose);
    load_hset(f, "big_horiz_clump",
                 "big horizontal cell maps with clumped fibers",
                 fBigHorizClumpMap, verbose);
    load_hset(f, "small_horiz",
                 "small horizontal cell maps",
                 fSmallHorizMap, verbose);
    load_hset(f, "small_horiz_clump",
                 "small horizontal cell maps with clumped fibers",
                 fSmallHorizClumpMap, verbose);
    load_hset(f, "small_horiz_bubble",
                 "small horizontal cell maps with bubbles",
                 fSmallHorizBubbleMap, verbose);
    load_hset(f, "small_horiz_bubble_clump",
                 "small horizontal cell maps with bubbles and clumped fibers",
                 fSmallHorizBubbleClumpMap, verbose);
    load_hset(f, "big_vert",
                 "big vertical cell maps",
                 fBigVertMap, verbose);
    load_hset(f, "big_vert_corner",
                 "big vertical cell maps with fibers in corners",
                 fBigVertCornerMap, verbose);
    load_hset(f, "small_vert",
                 "small vertical cell maps",
                 fSmallVertMap, verbose);
    load_hset(f, "small_vert_corner",
                 "small vertical cell maps with fibers in corners",
                 fSmallVertCornerMap, verbose);

    delete f;
  }

  // Assigns a global index to each map, in the order that the maps are kept
  // in the class definition.  Of course, in order to trace back what map was
  // used, you also need to know what ROOT file you loaded.
  void ImprovedTransport::SetPosCorrMapIndex(
    const std::vector<TH2D *> * const maps, const unsigned int i)
  {
    fPosCorrMapIndex = 0;
    if(maps == &fBigHorizMap){ fPosCorrMapIndex += i; return; }
    fPosCorrMapIndex += fBigHorizMap.size();
    if(maps == &fBigHorizClumpMap){ fPosCorrMapIndex += i; return; }
    fPosCorrMapIndex += fBigHorizClumpMap.size();
    if(maps == &fSmallHorizMap){ fPosCorrMapIndex += i; return; }
    fPosCorrMapIndex += fSmallHorizMap.size();
    if(maps == &fSmallHorizClumpMap){ fPosCorrMapIndex += i; return; }
    fPosCorrMapIndex += fSmallHorizClumpMap.size();
    if(maps == &fSmallHorizBubbleMap){ fPosCorrMapIndex += i; return; }
    fPosCorrMapIndex += fSmallHorizBubbleMap.size();
    if(maps == &fSmallHorizBubbleClumpMap){ fPosCorrMapIndex += i; return; }
    fPosCorrMapIndex += fSmallHorizBubbleClumpMap.size();
    if(maps == &fBigVertMap){ fPosCorrMapIndex += i; return; }
    fPosCorrMapIndex += fBigVertMap.size();
    if(maps == &fBigVertCornerMap){ fPosCorrMapIndex += i; return; }
    fPosCorrMapIndex += fBigVertCornerMap.size();
    if(maps == &fSmallVertMap){ fPosCorrMapIndex += i; return; }
    fPosCorrMapIndex += fSmallVertMap.size();
    if(maps == &fSmallVertCornerMap){ fPosCorrMapIndex += i; return; }
    fPosCorrMapIndex += fSmallVertCornerMap.size();
  }

  double ImprovedTransport::GetPosCorr( const sim::FLSHit& hit )
  {
    // If no file is specified for the efficiency maps -> disable
    if(fPosCorrFile.empty()) return 1;

    // Y is the large dimension of the cell
    const double Xa         = hit.GetEntryX();
    const double Ya         = hit.GetEntryY();
    const double Xb         = hit.GetExitX();
    const double Yb         = hit.GetExitY();
    const int   planeId     = hit.GetPlaneID();
    const int   cellId      = hit.GetCellID();
    const geo::CellUniqueId uniqueID = hit.GetCellUniqueId();

    const bool horizontal = planeId%2 == 0;

    const bool largecell = cellId%16 != 0 && cellId%16 != 15;

    const bool topcell = horizontal && cellId%32 == 31;

    art::ServiceHandle<geo::Geometry> geom;

    // This does not test whether we used a geometry that includes air
    // bubbles, so you need to only use this code if you are.  Preferably,
    // the size of the bubbles should match between the geometry and
    // the efficiency maps too, which also is not enforced here.
    const bool hasbubble = topcell &&
      (geom->DetId() == novadaq::cnv::kFARDET ||
       geom->DetId() == novadaq::cnv::kTESTBEAM);

    if(!fPosCorrFile.empty() && fBigHorizMap.empty()) LoadPosCorrHists();

    // 64 bits of seed are settable.  Try to make a decently unique
    // combination of the event number and cell number with this space.  This
    // gets us the same map within each art::Event given the same cell, but a
    // different map if either the event or the cell is different.  So
    // multiple energy depositions in a cell in one event are treated sensibly
    // (although not perfectly: this takes the efficiency to be constant along
    // the length of the cell)
    fPosCorrHasher.SetSeed(fPosCorrEvtId ^ uniqueID);

    // TODO: make probability into a fcl parameter
    const double clumpfrac = geom->DetId() == novadaq::cnv::kFARDET? 0.16
                            :geom->DetId() == novadaq::cnv::kNEARDET? 0.08
                            : 0.07 /* TESTBEAM (and any others...) */;
    const bool clump = horizontal && fPosCorrHasher.Rndm() < clumpfrac;

    bool sagged = false;
    if(!horizontal){
      const double Z = (hit.GetEntryZ() + hit.GetExitZ())/2;
      int dum1 = 0, dum2 = 0;
      const double frombottom =
        Z + geom->IdToCell(uniqueID, &dum1, &dum2)->HalfL();

      // Characteristic distance from bottom where fiber starts to
      // be free-floating instead sagged against the walls as per
      // doc-39723. We don't know this number precisely at all, but it
      // should probably be slightly larger for small cells, so we do
      // that by the ratio of perimeters of the cells.
      // TODO: make into a fcl parameter.  
      const double chardist = 100.0 * (1 + 0.012*(!largecell));

      sagged = fPosCorrHasher.Rndm() < exp(-frombottom/chardist);
    }

    // If the vertical fiber is sagged, guess that it will be in opposite
    // corners half the time, and against the walls elsewhere otherwise.
    const bool saggedcorner = sagged? fPosCorrHasher.Rndm() < 0.5: false;

    std::vector<TH2D *> * maps =
      &(horizontal?
         (hasbubble?
            (clump?  fSmallHorizBubbleClumpMap: fSmallHorizBubbleMap)
         : // no bubble
          largecell?
            (clump?  fBigHorizClumpMap: fBigHorizMap)
         : // no bubble small cell
            (clump?  fSmallHorizClumpMap: fSmallHorizMap))
       : // vertical
         (largecell?
           (sagged?
             // For sagged fibers not in the corner, use a *horizontal* map
             // since those are also have the fibers against the walls.  Of
             // course, it is always the same side wall, so this isn't perfect.
             (saggedcorner? fBigVertCornerMap: fBigHorizMap /* sic */)
             : fBigVertMap)
           : // small cell
           (sagged?
             (saggedcorner? fSmallVertCornerMap: fSmallHorizMap /* sic */)
             : fSmallVertMap)));

    const unsigned int mapi = fPosCorrHasher.Integer(maps->size());
    TH2D * themap = (*maps)[mapi];
    SetPosCorrMapIndex(maps, mapi);

    // Want this correction to be 1 on average, which this approximates.
    // It's not important to get this correct to three (or more) digits.
    // This scale doesn't take into account thresholds, for instance.
    // It will be necessary to calibrate at a higher level in any case.
    const double scale = 1/0.134;

    const double eff = line_efficiency(themap, Xa, Ya, Xb, Yb);
    #if 0
    printf("GetPosCorr gave %f for path %f %f %f %f len %f when %s %s %s %s\n",
           eff, Xa, Ya, Xb, Yb, hypot(Yb-Ya, Xb-Xa),
           horizontal?"horizontal":"vertical",
           sagged?"sagged":"notsagged",
           hasbubble?"bubble":"nobubble",
           largecell?"large":"small");
    #endif
    return eff * scale;
  }
  
  //............................................................
  void ImprovedTransport::StepAlongHit( sim::FLSHit& hit )
  {
    //clear photo-electron times for this hit
    fPETimes.clear();
    
    const double MeVPerGeV = 1000.;
    //const double green_scale = MeVPerGeV * fPhotonsPerMeV * fQuantumEff;
    
    //for cleaner-looking code below...
    double Xa       = hit.GetEntryX();
    double Ya       = hit.GetEntryY();
    double Za       = hit.GetEntryZ();
    double Ta       = hit.GetEntryT();

    double Xb       = hit.GetExitX();
    double Yb       = hit.GetExitY();
    double Zb       = hit.GetExitZ();
    double Tb       = hit.GetExitT();
    
    int   planeId   = hit.GetPlaneID();
    int   cellId    = hit.GetCellID();
    int   trackId   = hit.GetTrackID();

    const double positionCorr = GetPosCorr(hit);

    double brightnessCorr = fBrightnessMap[planeId][cellId];
    
    double seglength = util::pythag(Xb - Xa, Yb - Ya, Zb - Za);

    //cell half-width
    double DZ = fGeo->Plane((unsigned int)(planeId))->Cell ((unsigned int)(cellId ))->HalfL();
    //pigtail for cell
    double pigtail = fGeo->GetPigtail(fGeo->Plane(planeId)->Cell(cellId)->Id());
    
    double factor;
    if (fGeo->Plane(planeId)->View() == geo::kX) factor = fXViewFactor;
    else factor = fYViewFactor;

    factor *= fLightLevelScale; // Nominally this is one, but can be modified for LL shifts

    // If we're simulating drift, randomly throw the light level
    // not ready to go just yet -- j.
    if (fSimulateDrift) {
      art::ServiceHandle<nova::dbi::RunHistoryService> rh;
      int time = rh->TStart() + (0.5*rh->Duration()); // set to midpoint of run;
      double yDiff = 3.17098e-8 * (time - fDriftReference); // Amount of drift in years
      //std::cout << "Drifting " << yDiff << " years with respect to nominal." << std::endl;
      double drift = 1 - (fDriftGradient*yDiff);
      //std::cout << "Calculated drift of " << drift
      //  << " from nominal, factor corrected from " << factor << " to ";
      factor *= drift;
      //std::cout << factor << std::endl;
    }

    bool inMC = false;
    if (fGeo->DetId() == novadaq::cnv::kNEARDET) {
      if (planeId >= int(fGeo->FirstPlaneInMuonCatcher())) inMC = true;
    }
    
    //find out how many steps to take
    int nsteps = ceil( seglength/fStep );

    //Get approximate dEdx for scintillation time scaling
    //If segment length is zero, dEdx won't make sense.  
    //Default to zero - the beta PDF is almost always going to be more relevant than the alpha one,  I think
    double dEdx = (seglength > 0) ? MeVPerGeV*hit.GetEdep()/seglength : 0.0;

    // The total number of PE that would have been produced are accumulated in
    // here, whether or not fluctuations allowed a PhotonSignal to be created
    // for them to live in.
    double lambda = 0;
    
    for (int istep = 0; istep < nsteps; ++istep) {
      //at midpoint of the step
      double stepZ = Za + ((Zb-Za)/nsteps)*(istep+0.5);
      double stepT = Ta + ((Tb-Ta)/nsteps)*(istep+0.5);
      double stepE = (fApplyBirks) ? hit.GetEdepBirks()/nsteps : hit.GetEdep()/nsteps;
      double stepCerenkov = hit.GetNCerenkov()/nsteps;

      double tweakCorr(1.0);
      if (fApplyTweak) {
        //Tweak is currently only available for non-muon catcher cells
        if (fGeo->DetId() == novadaq::cnv::kFARDET) {
          if (fGeo->Plane(planeId)->View() == geo::kX) tweakCorr = fTweak_X.Interpolate(DZ  - stepZ);
          else tweakCorr = fTweak_Y.Interpolate(DZ  - stepZ);
        }
        if (fGeo->DetId() == novadaq::cnv::kNEARDET) {
          if (inMC) {
            if (fGeo->Plane(planeId)->View() == geo::kX) tweakCorr = fTweak_XMC.Interpolate(DZ  - stepZ);
            else tweakCorr = fTweak_YMC.Interpolate(DZ  - stepZ);
          }
          else {
            if (fGeo->Plane(planeId)->View() == geo::kX) tweakCorr = fTweak_X.Interpolate(DZ  - stepZ);
            else tweakCorr = fTweak_Y.Interpolate(DZ  - stepZ);
          }
        }
      }

      //green_scale*stepE*collectionRate = N Photo-electrons (before attenuation)
      std::list<RateInfo>::iterator rate_it;
      for (rate_it = fRateInfo.begin(); rate_it != fRateInfo.end(); ++rate_it) {
        double collectionRate = rate_it->rate;
        double currZ = stepZ + rate_it->z;
        double fraction(1.0);

        //is this currZ inside the cell? 
        if ( std::fabs(currZ) - 0.5*rate_it->zWidth > DZ ) continue;
        else if ( std::fabs(currZ) + 0.5*rate_it->zWidth < DZ ) fraction = 1.0;
        else fraction = (DZ - ( std::fabs(currZ) - 0.5*rate_it->zWidth ))/rate_it->zWidth;
        collectionRate *= fraction;
        
        double scintillationPhotons = MeVPerGeV*fPhotonsPerMeV*stepE;
        double cerenkovPhotons      = fCerenkovEff*stepCerenkov;
        double muGreenScint = fQuantumEff*collectionRate*brightnessCorr*tweakCorr*
          positionCorr*factor*scintillationPhotons;
        double muGreenCkv = fQuantumEff*collectionRate*brightnessCorr*tweakCorr*
          positionCorr*factor*cerenkovPhotons;
        //short and long length to the electronics
        double dist0 = DZ  - currZ + pigtail;
        double dist1 = 3*DZ + currZ + pigtail;                      
        
        CreatePhotoElectrons( 0.5*muGreenScint, stepT, dist0, dEdx, true,  *rate_it, lambda );
        CreatePhotoElectrons( 0.5*muGreenScint, stepT, dist1, dEdx, true,  *rate_it, lambda );
        CreatePhotoElectrons( 0.5*muGreenCkv,   stepT, dist0, dEdx, false, *rate_it, lambda );
        CreatePhotoElectrons( 0.5*muGreenCkv,   stepT, dist1, dEdx, false, *rate_it, lambda );
      }
    }

    // Keep track of how many we have so far so we know which ones get newly
    // added by ClusterPhotoElectrons().
    unsigned int start = fSignals->size();

    ClusterPhotoElectrons(planeId, cellId, trackId);
    
    // Divide the poission lambda evenly among all PhotonSignals. When someone
    // adds them up downstream they'll get the right overall answer for this
    // FLSHit.
    const int N = fSignals->size()-start;
    if(N > 0){
      for(unsigned int i = start; i < fSignals->size(); ++i){
        (*fSignals)[i].SetPoissonLambda(lambda/N);
      }
    }
    else{
      if(fWriteEmptySignals){
        // Create a PhotonSignal with zero PE, just to have somewhere to attach
        // the PoissonLambda to.
        sim::PhotonSignal photon(hit.GetEntryT(), 0);
        photon.SetPlaneCell(planeId, cellId);
        photon.SetTrackId(trackId);
        photon.SetPoissonLambda(lambda);
        fSignals->push_back(photon);
      }
    }
    
    return;
  }
    
  //............................................................ 
  void ImprovedTransport::CreatePhotoElectrons(double meanPE, double stepTime, double dist, double dEdx, bool isScint, RateInfo& info, double& lambda)
  {
    double expect = meanPE*FiberTransmission(dist);
    lambda += expect;

    int nPE = fPoisson->fire( expect );
    double rand0(0), rand1(0), rand2(0);
    for (int iPE=0; iPE < nPE; ++iPE) {
      double time = stepTime;
      time += dist/fFiberSpeedOfLight;
      
      fRandHisto->GetRandom( &info.timePDF, rand0, rand1, rand2 );      
      time += rand0;
      
      time += fExp->fire( fEmissionTau );

      double scintTime(0);
      if (fUseScintTime)
        {
          if (isScint) scintTime = ScintTime(dEdx);
          else         scintTime = fExp->fire( fCherenkovTau );
        }
      time += scintTime;

      double fiberspread = 0;
      if (fApplyFiberSmearing) {
        double weight(-1.0), uweight(0.0);
        while (weight < uweight) {
          // Draw a Dt/z value then multiply by "dist" to get time (in ns)
          fRandHisto->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*fFiberSpeedOfLight;
          weight = FiberTransmission(dist+Delta_dist)/FiberTransmission(dist);
          uweight = fFlat->fire(0,1);
        }
      }
      time += fiberspread;
      fPETimes.push_back( time );
      fScintTime->Fill(time - stepTime);
    }
    return;
  }
  
  
  //............................................................ 
  void ImprovedTransport::ClusterPhotoElectrons(int planeId, int cellId, int trackId)
  {
    if (fPETimes.empty()) return;
    //cluster photons into PhotonSignals fTimeClustering long

    sim::PhotonSignal photon;
    photon.SetPlaneCell(planeId, cellId);
    photon.SetTrackId(trackId);

    // TODO: Would like to record which map is used, but having trouble
    // getting it to work.
    // photon.SetPosCorrI(fPosCorrFile.empty()?-1:fPosCorrMapIndex);

    double meanTime(0);
    int nPE(0);

    fPETimes.sort();
    std::list< double >::iterator it;
    for (it = fPETimes.begin(); it != fPETimes.end(); ++it) {
      double currTime = *it;
      if (nPE == 0) {
        meanTime = currTime;
        nPE = 1;
      }
      else if ( fabs(meanTime - currTime) < fTimeClustering ) {
        //add to cluster
        meanTime = meanTime*nPE + currTime;
        ++nPE;
        meanTime /= nPE;
      }
      else {
        photon.SetTimeMean(meanTime);
        photon.SetNPhoton(nPE);
        fSignals->push_back(photon);
        
        //Set variables with current iteration
        nPE = 1;
        meanTime = currTime;
      }
    }
    //there is one PhotonSignal that hasn't been saved yet
    photon.SetTimeMean(meanTime);
    photon.SetNPhoton(nPE);
    fSignals->push_back(photon);
    
    return;
  }
  
  //............................................................
  bool ImprovedTransport::LoadHistos() 
  {
    // 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(fCollectionRateFile.c_str(), &sb) !=0 ) {
      throw cet::exception("NoCollectionRateFile")
        << "ImprovedTransport: Unable to open " 
        << fCollectionRateFile << "\n"
        << __FILE__ << ":" << __LINE__ << "\n";
      return false;
    }
    
    if ( stat(fSmearingHistoFile.c_str(), &sb) !=0 ) {
      throw cet::exception("NoSmearingHistoFile") << "ImprovedTransport: Unable to open " << fSmearingHistoFile << "\n" << __FILE__ << ":" << __LINE__ << "\n";
      return false;
    }
    
    TFile f1(fCollectionRateFile.c_str());
    TH2D *htemp = (TH2D*)f1.Get("dT_dZ_CollectionRate")->Clone();
    if (!htemp) {
      std::cerr << "ImprovedTransport: Histogram dT_dZ_CollectionRate not found in " << fCollectionRateFile << std::endl;
      return false;
    }
    
    fCollectionRateHisto = *htemp;
    fCollectionRateHisto.SetDirectory(0);
    f1.Close();
    
    TFile f2(fSmearingHistoFile.c_str());
    TH1D *htemp2 = (TH1D*)f2.Get("Dt_per_z")->Clone();
    if (!htemp2) {
      std::cerr << "ImprovedTransport: Histogram Dt_per_z_distribution not found in " << fSmearingHistoFile << std::endl;
      return false;
    }
    
    fFiberSmearingHisto = *htemp2;
    fFiberSmearingHisto.SetDirectory(0);
    f2.Close();

    for (unsigned int plane = 0; plane < fGeo->NPlanes(); plane++){
      std::vector< double > plane_brightness;
      for (unsigned int cell = 0; cell < fGeo->Plane(plane)->Ncells(); cell++){
        plane_brightness.push_back(1.0);
      }
      fBrightnessMap.push_back(plane_brightness);
    }
    
    if (fBrightnessMode != "None") {
        TFile f3(fBrightnessFile.c_str());
        TH1D* hBrightness1D(0);
        TH2D* hBrightness2D(0);
        if (fBrightnessMode == "ByCell") {
            hBrightness2D = (TH2D*)f3.Get("BrightnessByCell");
            if (!hBrightness2D) {
                std::cerr << "ImprovedTransport: Histogram BrightnessByCell not found in " << fBrightnessFile << std::endl;
                return false;
            }
        }
        else if (fBrightnessMode == "ByModule") {
            hBrightness2D = (TH2D*)f3.Get("BrightnessByModule");
            if (!hBrightness2D) {
                std::cerr << "ImprovedTransport: Histogram BrightnessByModule not found in " << fBrightnessFile << std::endl;
                return false;
            }
        }
        else if (fBrightnessMode == "ByPlane") {
            hBrightness1D = (TH1D*)f3.Get("BrightnessByPlane");
            if (!hBrightness1D) {
                std::cerr << "ImprovedTransport: Histogram BrightnessByPlane not found in " << fBrightnessFile << std::endl;
                return false;
            }
        }
        //Add by Pengfei -- begin
        else if (fBrightnessMode == "ByBin") {
          fBrightnessLevel = new BrightnessLevel(fBrightnessFile, fBrightnessMapName, fBrightnessValueName,fBrightnessValueMapName);
        }
        //Add by Pengfei -- end
        else {
            std::cerr << "ImprovedTransport: " << fBrightnessMode << " is not a recognized brightness correction mode. Valid choices are None, ByCell, ByModule, or ByPlane " << std::endl;
            return false;
        }

        for (unsigned int plane = 0; plane < fGeo->NPlanes(); plane++){
            for (unsigned int cell = 0; cell < fGeo->Plane(plane)->Ncells(); cell++){
                double brightness(1.0);
                if (fBrightnessMode == "ByCell") brightness = hBrightness2D->GetBinContent(plane+1, cell+1);
                else if (fBrightnessMode == "ByModule") brightness = hBrightness2D->GetBinContent(plane+1, cell/32 + 1);
                else if (fBrightnessMode == "ByPlane")  brightness = hBrightness1D->GetBinContent(plane+1);
                //Add by Pengfei -- begin
                else if (fBrightnessMode == "ByBin") {
                    brightness = fBrightnessLevel->getBrightnessLevel(plane, cell);
                }
                //Add by Pengfei -- end
                fBrightnessMap[plane][cell] = brightness;
            }
        }
        f3.Close();
    }

    if (fApplyTweak) {
      TFile f4(fTweakFile.c_str());
      TH1D* tmpTweak(0);
      tmpTweak = (TH1D*)f4.Get("Tweak_X")->Clone();
      if (!tmpTweak) {
        std::cerr << "Histogram Tweak_X not found in " << fTweakFile << std::endl;
        return false;
      }
      fTweak_X = *tmpTweak;
      fTweak_X.SetDirectory(0);

      tmpTweak = 0;
      tmpTweak = (TH1D*)f4.Get("Tweak_Y")->Clone();
      if (!tmpTweak) {
        std::cerr << "Histogram Tweak_Y not found in " << fTweakFile <<std::endl;
        return false;
      }
      fTweak_Y = *tmpTweak;
      fTweak_Y.SetDirectory(0);

      f4.Close();
    }
    
    fRateInfo.clear();
    std::stringstream hName;
    int nbinsY = fCollectionRateHisto.GetNbinsY();
    double ylo = fCollectionRateHisto.GetYaxis()->GetXmin();
    double yhi = fCollectionRateHisto.GetYaxis()->GetXmax();
    for (int iZ = 1; iZ <= fCollectionRateHisto.GetNbinsX(); ++iZ) {
      RateInfo info_pos;
      info_pos.rate = 0.0;
      info_pos.z = fCollectionRateHisto.GetXaxis()->GetBinCenter(iZ);
      info_pos.zWidth = fCollectionRateHisto.GetXaxis()->GetBinWidth(iZ);
      hName.str("timePDF_Pos_");
      hName << iZ;
      info_pos.timePDF.SetName(hName.str().c_str());
      info_pos.timePDF.SetBins(nbinsY, ylo, yhi);
      info_pos.timePDF.SetDirectory(0);
      
      RateInfo info_neg;
      info_neg.rate = 0.0;
      info_neg.z = -fCollectionRateHisto.GetXaxis()->GetBinCenter(iZ);
      info_neg.zWidth = fCollectionRateHisto.GetXaxis()->GetBinWidth(iZ);
      hName.str("timePDF_Neg_");
      hName << iZ;
      info_neg.timePDF.SetName(hName.str().c_str());
      info_neg.timePDF.SetBins(nbinsY, ylo, yhi);
      info_neg.timePDF.SetDirectory(0);
      
      for (int iT = 1; iT <= fCollectionRateHisto.GetNbinsY(); ++iT) {
        double collectionRate = 0.5*fCollectionRateHisto.GetBinContent(iZ,iT);
        info_pos.rate += collectionRate;
        info_pos.timePDF.SetBinContent(iT, collectionRate);
        
        info_neg.rate += collectionRate;
        info_neg.timePDF.SetBinContent(iT, collectionRate);
      }
      fRateInfo.push_back(info_pos);
      fRateInfo.push_back(info_neg);
    }
    return true;
  }
  
  //............................................................
  double ImprovedTransport::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) {
      mf::LogWarning("AttenPars") << "No AttenPars entries.  Fiber transmission set to 1!";
      answer = 1;
    }
    else if (nPar%2) {
      mf::LogWarning("AttenPars") << "Odd number of AttenPars entries.  Fiber transmission set to 1!";
      answer = 1;
    }
    else {
      answer = 0;    
      for (Int_t iP=0;iP<nPar;iP+=2) answer += fAttenPars[iP]*exp(-x/fAttenPars[iP+1]);
    }
    return answer;
  }

  //............................................................
  double ImprovedTransport::ScintTime(double dEdx) 
  {
    // sample from a sum of multiple exponentials to get the scintilator decay time for each photon
    // if dEdx is less than that used for measuring the beta distribution, use the beta distribution
    // if dEdx is greater than that used for measuring the alpha distribution, use the alpha distribution
    // if dEdx is between the two, smoothly transition between the distributions
    
    int nProbBeta  = fScintBetaProb.size();
    int nProbAlpha = fScintAlphaProb.size();
    int nTauBeta   = fScintBetaTau.size();
    int nTauAlpha  = fScintAlphaTau.size();

    double time = 0;
    if (!nProbBeta) {
      mf::LogWarning("ScintTime") << "No ScintBetaProb entries.  No scintillation delays applied!";
      time = 0;
    }
    else if (!nProbAlpha) {
      mf::LogWarning("ScintTime") << "No ScintAlphaProb entries.  No scintillation delays applied!";
      time = 0;
    }
    else if (nProbBeta != nTauBeta) {
      mf::LogWarning("ScintTime") << "Number of entries in ScintBetaProb and ScintBetaTau not the same. No scintillation delays applied!";
      time = 0;
    }
    else if (nProbAlpha != nTauAlpha) {
      mf::LogWarning("ScintTime") << "Number of entries in ScintBetaAlpha and ScintBetaAlpha not the same. No scintillation delays applied!";
      time = 0;
    }
    else if (fScintBetaDEDX >= fScintAlphaDEDX) {
      mf::LogWarning("ScintTime") << "Beta dEdx must be strictly smaller than alpha dEdx. No scintillation delays applied!";
      time = 0;
    }
    else {
      //actually sample the distribution
       double probBeta = (dEdx - fScintBetaDEDX)/(fScintBetaDEDX - fScintAlphaDEDX) + 1;
       double u = fFlat->fire(0,1);
       bool useBeta = (u < probBeta);
       u = fFlat->fire(0,1);
       if (useBeta) {
         for (int iC = 0; iC < nProbBeta; ++iC) {
           if (u <= fScintBetaProb[iC]) {
             time = fExp->fire(fScintBetaTau[iC]);
             break;
           }
         }
       }
       else {
         for (int iC = 0; iC < nProbAlpha; ++iC) {
           if (u <= fScintAlphaProb[iC]) {
             time = fExp->fire(fScintAlphaTau[iC]);
             break;
           }
         }
       }
    }
    return time;      
  }

  DEFINE_ART_MODULE(ImprovedTransport)
} //namespace