OpFastScintillation.cxx

Go to the documentation of this file.

// Class adapted for LArSoft by Ben Jones, MIT 10/10/12
//
// This class is a physics process based on the standard Geant4 scintillation process.
//
// It has been stripped down and adapted to form the backbone of the LArG4 fast optical
// simulation.  Photons, instead of being produced and added to the geant4 particle stack,
// are logged and used to predict the visibility of this step to each PMT in the detector.
//
// The photonvisibilityservice looks up the visibility of the relevant xyz point, and if a
// photon is detected at a given PMT, one OnePhoton object is logged in the
// OpDetPhotonTable
//
// At the end of the event, the OpDetPhotonTable is read out by LArG4, and detected
// photons are stored in the event.
//
// This process can be used alongside the standard Cerenkov process, which does step
// geant4 opticalphotons.  Both the fast scintillation table and the geant4 sensitive
// detectors are read out by LArG4 to produce a combined SimPhoton collection.
//
// Added disclaimer : This code is gross.  Thats basically because it adheres to the
//  original, gross Geant4 implementation.

#include "Geant4/G4Alpha.hh"
#include "Geant4/G4DynamicParticle.hh"
#include "Geant4/G4Electron.hh"
#include "Geant4/G4ExceptionSeverity.hh"
#include "Geant4/G4Gamma.hh"
#include "Geant4/G4KaonMinus.hh"
#include "Geant4/G4KaonPlus.hh"
#include "Geant4/G4Material.hh"
#include "Geant4/G4MaterialPropertiesTable.hh"
#include "Geant4/G4MaterialPropertyVector.hh"
#include "Geant4/G4MaterialTable.hh"
#include "Geant4/G4MuonMinus.hh"
#include "Geant4/G4MuonPlus.hh"
#include "Geant4/G4ParticleChange.hh"
#include "Geant4/G4PhysicsVector.hh"
#include "Geant4/G4PionMinus.hh"
#include "Geant4/G4PionPlus.hh"
#include "Geant4/G4Poisson.hh"
#include "Geant4/G4Proton.hh"
#include "Geant4/G4Step.hh"
#include "Geant4/G4StepPoint.hh"
#include "Geant4/G4Track.hh"
#include "Geant4/G4VPhysicalVolume.hh"
#include "Geant4/G4ios.hh"
#include "Geant4/Randomize.hh"
#include "Geant4/globals.hh"

#include "larcore/CoreUtils/ServiceUtil.h"
#include "larcore/Geometry/Geometry.h"
#include "larcorealg/CoreUtils/counter.h"
#include "larcorealg/Geometry/OpDetGeo.h"
#include "larcorealg/Geometry/geo_vectors_utils.h"         // geo::vect::fillCoords()
#include "larcorealg/Geometry/geo_vectors_utils_TVector.h" // geo::vect::toTVector3()
#include "lardata/DetectorInfoServices/LArPropertiesService.h"
#include "lardataobj/Simulation/OpDetBacktrackerRecord.h"
#include "lardataobj/Simulation/SimPhotons.h"
#include "larsim/LegacyLArG4/IonizationAndScintillation.h"
#include "larsim/LegacyLArG4/OpDetPhotonTable.h"
#include "larsim/LegacyLArG4/OpFastScintillation.hh"
#include "larsim/LegacyLArG4/ParticleListAction.h"
#include "larsim/PhotonPropagation/PhotonVisibilityService.h"
#include "larsim/Simulation/LArG4Parameters.h"

// support libraries
#include "cetlib_except/exception.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "TGeoSphere.h"
#include "TMath.h"
#include "TRandom3.h"

#include "range/v3/view/enumerate.hpp"

#include <cassert>
#include <cmath>
#include <limits>

#include "boost/math/special_functions/ellint_1.hpp"
#include "boost/math/special_functions/ellint_3.hpp"

// Define a new policy *not* internally promoting RealType to double:
typedef boost::math::policies::policy<boost::math::policies::promote_double<false>>
  noLDoublePromote;

namespace {

  double finter_d(double* x, double* par)
  {
    double y1 = par[2] * TMath::Landau(x[0], par[0], par[1]);
    double y2 = TMath::Exp(par[3] + x[0] * par[4]);
    return TMath::Abs(y1 - y2);
  }

  double model_close(double* x, double* par)
  {
    // par0 = joining point
    // par1 = Landau MPV
    // par2 = Landau width
    // par3 = normalization
    // par4 = Expo cte
    // par5 = Expo tau
    // par6 = t_min
    double y1 = par[3] * TMath::Landau(x[0], par[1], par[2]);
    double y2 = TMath::Exp(par[4] + x[0] * par[5]);
    if (x[0] <= par[6] || x[0] > par[0]) y1 = 0.;
    if (x[0] < par[0]) y2 = 0.;
    return (y1 + y2);
  }

  double model_far(double* x, double* par)
  {
    // par1 = Landau MPV
    // par2 = Landau width
    // par3 = normalization
    // par0 = t_min
    double y = par[3] * TMath::Landau(x[0], par[1], par[2]);
    if (x[0] <= par[0]) y = 0.;
    return y;
  }

  double fast_acos(double x)
  {
    double negate = double(x < 0);
    x = std::abs(x);
    x -= double(x > 1.0) * (x - 1.0); // <- equivalent to min(1.0,x), but faster
    double ret = -0.0187293;
    ret = ret * x;
    ret = ret + 0.0742610;
    ret = ret * x;
    ret = ret - 0.2121144;
    ret = ret * x;
    ret = ret + 1.5707288;
    ret = ret * std::sqrt(1.0 - x);
    ret = ret - 2. * negate * ret;
    return negate * 3.14159265358979 + ret;
  }
}

namespace larg4 {

  OpFastScintillation::OpFastScintillation(const G4String& processName, G4ProcessType type)
    : G4VRestDiscreteProcess(processName, type)
    , fActiveVolumes{extractActiveVolumes(*(lar::providerFrom<geo::Geometry>()))}
    , bPropagate(!(art::ServiceHandle<sim::LArG4Parameters const>()->NoPhotonPropagation()))
    , fPVS(bPropagate ? art::ServiceHandle<phot::PhotonVisibilityService const>().get() : nullptr)
    , fUseNhitsModel(fPVS && fPVS->UseNhitsModel())
    // for now, limit to the active volume only if semi-analytic model is used
    , fOnlyActiveVolume(usesSemiAnalyticModel())
  {
    SetProcessSubType(25); // TODO: unhardcode
    fTrackSecondariesFirst = false;
    fFiniteRiseTime = false;
    YieldFactor = 1.0;
    ExcitationRatio = 1.0;

    const detinfo::LArProperties* larp = lar::providerFrom<detinfo::LArPropertiesService>();

    scintillationByParticleType = larp->ScintByParticleType();

    if (verboseLevel > 0) { G4cout << GetProcessName() << " is created " << G4endl; }

    BuildThePhysicsTable();
    emSaturation = NULL;

    if (bPropagate) {
      assert(fPVS);

      // Loading the position of each optical channel, neccessary for the
      // parametrizatiuons of Nhits and prop-time
      geo::GeometryCore const& geom = *(lar::providerFrom<geo::Geometry>());
      {
        auto log = mf::LogTrace("OpFastScintillation")
                   << "OpFastScintillation: active volume boundaries from " << fActiveVolumes.size()
                   << " volumes:";
        for (auto const& [iCryo, box] : fActiveVolumes | ranges::views::enumerate) {
          log << "\n - C:" << iCryo << ": " << box.Min() << " -- " << box.Max() << " cm";
        } // for
        log << "\n  (scintillation photons are propagated "
            << (fOnlyActiveVolume ? "only from active volumes" : "from anywhere") << ")";
      } // local scope

      if (usesSemiAnalyticModel() && (geom.Ncryostats() > 1U)) {
        if (fOnlyOneCryostat) {
          mf::LogWarning("OpFastScintillation")
            << std::string(80, '=') << "\nA detector with " << geom.Ncryostats()
            << " cryostats is configured"
            << " , and semi-analytic model is requested for scintillation photon propagation."
            << " THIS CONFIGURATION IS NOT SUPPORTED and it is open to bugs"
            << " (e.g. scintillation may be detected only in cryostat #0)."
            << "\nThis would be normally a fatal error, but it has been forcibly overridden."
            << "\n"
            << std::string(80, '=');
        }
        else {
          throw art::Exception(art::errors::Configuration)
            << "Photon propagation via semi-analytic model is not supported yet"
            << " on detectors with more than one cryostat.";
        }
      } // if

      geo::Point_t const Cathode_centre{geom.TPC().GetCathodeCenter().X(),
                                        fActiveVolumes[0].CenterY(),
                                        fActiveVolumes[0].CenterZ()};
      mf::LogTrace("OpFastScintillation") << "Cathode_centre: " << Cathode_centre << " cm";

      for (size_t const i : util::counter(fPVS->NOpChannels())) {
        geo::OpDetGeo const& opDet = geom.OpDetGeoFromOpDet(i);
        fOpDetCenter.push_back(opDet.GetCenter());

        if (dynamic_cast<TGeoSphere const*>(opDet.Shape()) != nullptr) { // sphere/dome
          fOpDetType.push_back(1);                                       // Dome PMTs
          fOpDetLength.push_back(-1);
          fOpDetHeight.push_back(-1);
        }
        else if (opDet.isBar()) {  // box
          fOpDetType.push_back(0); //Arapucas
          fOpDetLength.push_back(opDet.Length());
          fOpDetHeight.push_back(opDet.Height());
        }
        else {                     // disk
          fOpDetType.push_back(2); // Disk PMTs
          fOpDetLength.push_back(-1);
          fOpDetHeight.push_back(-1);
        }
      }

      if (fPVS->IncludePropTime()) {
        std::cout << "Using parameterisation of timings." << std::endl;
        //New VUV time parapetrization
        fPVS->LoadTimingsForVUVPar(fparameters,
                                   fstep_size,
                                   fmax_d,
                                   fmin_d,
                                   fvuv_vgroup_mean,
                                   fvuv_vgroup_max,
                                   finflexion_point_distance,
                                   fangle_bin_timing_vuv);

        // create vector of empty TF1s that will be replaces with the parameterisations
        // that are generated as they are required default TF1() constructor gives
        // function with 0 dimensions, can then check numDim to qucikly see if a
        // parameterisation has been generated
        const size_t num_params =
          (fmax_d - fmin_d) /
          fstep_size; // for d < fmin_d, no parameterisaton, a delta function is used instead
        size_t num_angles = std::round(90 / fangle_bin_timing_vuv);
        VUV_timing = std::vector(num_angles, std::vector(num_params, TF1()));

        // initialise vectors to contain range parameterisations sampled to in each case
        // when using TF1->GetRandom(xmin,xmax), must be in same range otherwise sampling
        // is regenerated, this is the slow part!
        VUV_max = std::vector(num_angles, std::vector(num_params, 0.0));
        VUV_min = std::vector(num_angles, std::vector(num_params, 0.0));

        // VIS time parameterisation
        if (fPVS->StoreReflected()) {
          // load parameters
          fPVS->LoadTimingsForVISPar(fdistances_refl,
                                     fradial_distances_refl,
                                     fcut_off_pars,
                                     ftau_pars,
                                     fvis_vmean,
                                     fangle_bin_timing_vis);
        }
      }
      if (usesSemiAnalyticModel()) {
        mf::LogVerbatim("OpFastScintillation")
          << "OpFastScintillation: using semi-analytic model for number of hits";

        // LAr absorption length in cm
        std::map<double, double> abs_length_spectrum =
          lar::providerFrom<detinfo::LArPropertiesService>()->AbsLengthSpectrum();
        std::vector<double> x_v, y_v;
        for (auto elem : abs_length_spectrum) {
          x_v.push_back(elem.first);
          y_v.push_back(elem.second);
        }
        fL_abs_vuv =
          std::round(interpolate(x_v, y_v, 9.7, false)); // 9.7 eV: peak of VUV emission spectrum

        // Load Gaisser-Hillas corrections for VUV semi-analytic hits
        std::cout << "Loading the GH corrections" << std::endl;
        fPVS->LoadVUVSemiAnalyticProperties(
          fIsFlatPDCorr, fIsDomePDCorr, fdelta_angulo_vuv, fradius);
        if (!fIsFlatPDCorr && !fIsDomePDCorr) {
          throw cet::exception("OpFastScintillation")
            << "Both isFlatPDCorr and isDomePDCorr parameters are false, at least one type of "
               "parameterisation is required for the semi-analytic light simulation."
            << "\n";
        }
        if (fIsFlatPDCorr) {
          fPVS->LoadGHFlat(fGHvuvpars_flat, fborder_corr_angulo_flat, fborder_corr_flat);
        }
        if (fIsDomePDCorr) {
          fPVS->LoadGHDome(fGHvuvpars_dome, fborder_corr_angulo_dome, fborder_corr_dome);
        }
        // cathode center coordinates required for corrections
        fcathode_centre = geom.TPC().GetCathodeCenter();
        fcathode_centre.SetY(fActiveVolumes[0].CenterY());
        fcathode_centre.SetZ(
          fActiveVolumes[0].CenterZ()); // to get full cathode dimension rather than just single tpc

        if (fPVS->StoreReflected()) {
          fStoreReflected = true;
          // Load corrections for VIS semi-anlytic hits
          std::cout << "Loading visible light corrections" << std::endl;
          fPVS->LoadVisSemiAnalyticProperties(fdelta_angulo_vis, fradius);
          if (fIsFlatPDCorr) {
            fPVS->LoadVisParsFlat(fvis_distances_x_flat, fvis_distances_r_flat, fvispars_flat);
          }
          if (fIsDomePDCorr) {
            fPVS->LoadVisParsDome(fvis_distances_x_dome, fvis_distances_r_dome, fvispars_dome);
          }

          // cathode dimensions
          fcathode_ydimension = fActiveVolumes[0].SizeY();
          fcathode_zdimension = fActiveVolumes[0].SizeZ();

          // set cathode plane struct for solid angle function
          fcathode_plane.h = fcathode_ydimension;
          fcathode_plane.w = fcathode_zdimension;
          fplane_depth = std::abs(fcathode_centre.X());
        }
        else
          fStoreReflected = false;
      }
    }
    tpbemission = lar::providerFrom<detinfo::LArPropertiesService>()->TpbEm();
    std::vector<double> parent;
    parent.reserve(tpbemission.size());
    for (auto iter = tpbemission.begin(); iter != tpbemission.end(); ++iter) {
      parent.push_back(iter->second);
    }
    fTPBEm = std::make_unique<CLHEP::RandGeneral>(parent.data(), parent.size());
  }

  OpFastScintillation::~OpFastScintillation()
  {
    if (theFastIntegralTable) theFastIntegralTable->clearAndDestroy();
    if (theSlowIntegralTable) theSlowIntegralTable->clearAndDestroy();
  }

  // AtRestDoIt
  // ----------
  //
  G4VParticleChange* OpFastScintillation::AtRestDoIt(const G4Track& aTrack, const G4Step& aStep)
  // This routine simply calls the equivalent PostStepDoIt since all the necessary
  // information resides in aStep.GetTotalEnergyDeposit()
  {
    return OpFastScintillation::PostStepDoIt(aTrack, aStep);
  }

  // PostStepDoIt
  // -------------
  //
  G4VParticleChange* OpFastScintillation::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
  // This routine is called for each tracking step of a charged particle in a
  // scintillator. A Poisson/Gauss-distributed number of photons is generated according to
  // the scintillation yield formula, distributed evenly along the track segment and
  // uniformly into 4pi.
  {
    aParticleChange.Initialize(aTrack);
    // Check that we are in a material with a properties table, if not just return
    const G4Material* aMaterial = aTrack.GetMaterial();
    G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
    if (!aMaterialPropertiesTable) return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);

    G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
    G4ThreeVector x0 = pPreStepPoint->GetPosition();
    G4ThreeVector p0 = aStep.GetDeltaPosition().unit();

    // We don't want to produce any trackable G4 secondaries
    aParticleChange.SetNumberOfSecondaries(0);

    // Retrieve the Scintillation Integral for this material new
    // G4PhysicsOrderedFreeVector allocated to hold CII's

    // Some explanation for later improvements to scint yield code:
    //
    // What does G4 do here?
    //  It produces light in 2 steps, fast (scnt=1) then slow (scnt=2)
    //
    // The ratio of slow photons to fast photons is related by the yieldratio parameter.
    //  G4's poisson fluctuating scheme is a bit different to ours - we should check that
    //  they are equivalent.
    //
    // G4 poisson fluctuates the number of initial photons then divides them with a
    // constant factor between fast + slow, whereas we poisson fluctuate separateyly the
    // fast and slow detection numbers.
    //

    // get the number of photons produced from the IonizationAndScintillation singleton
    larg4::IonizationAndScintillation::Instance()->Reset(&aStep);
    double MeanNumberOfPhotons =
      larg4::IonizationAndScintillation::Instance()->NumberScintillationPhotons();

    RecordPhotonsProduced(aStep, MeanNumberOfPhotons); //, stepEnergy);
    if (verboseLevel > 0) {
      G4cout << "\n Exiting from OpFastScintillation::DoIt -- NumberOfSecondaries = "
             << aParticleChange.GetNumberOfSecondaries() << G4endl;
    }
    return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
  }

  void OpFastScintillation::ProcessStep(const G4Step& step)
  {
    if (step.GetTotalEnergyDeposit() <= 0) return;

    OpDetPhotonTable::Instance()->AddEnergyDeposit(
      -1,
      -1,
      1.0,                                                 //scintillation yield
      (double)(step.GetTotalEnergyDeposit() / CLHEP::MeV), //energy in MeV
      (float)(step.GetPreStepPoint()->GetPosition().x() / CLHEP::cm),
      (float)(step.GetPreStepPoint()->GetPosition().y() / CLHEP::cm),
      (float)(step.GetPreStepPoint()->GetPosition().z() / CLHEP::cm),
      (float)(step.GetPostStepPoint()->GetPosition().x() / CLHEP::cm),
      (float)(step.GetPostStepPoint()->GetPosition().y() / CLHEP::cm),
      (float)(step.GetPostStepPoint()->GetPosition().z() / CLHEP::cm),
      (double)(step.GetPreStepPoint()->GetGlobalTime()),
      (double)(step.GetPostStepPoint()->GetGlobalTime()),
      ParticleListAction::GetCurrentTrackID(),
      step.GetTrack()->GetParticleDefinition()->GetPDGEncoding(),
      ParticleListAction::GetCurrentOrigTrackID(),
      step.GetPreStepPoint()->GetPhysicalVolume()->GetName());
  }

  bool OpFastScintillation::RecordPhotonsProduced(const G4Step& aStep, double MeanNumberOfPhotons)
  {
    // make sure that whatever happens afterwards, the energy deposition is stored
    art::ServiceHandle<sim::LArG4Parameters const> lgp;
    if (lgp->FillSimEnergyDeposits()) ProcessStep(aStep);

    // Get the pointer to the fast scintillation table
    OpDetPhotonTable* fst = OpDetPhotonTable::Instance();

    const G4Track* aTrack = aStep.GetTrack();

    G4StepPoint const* pPreStepPoint = aStep.GetPreStepPoint();

    const G4DynamicParticle* aParticle = aTrack->GetDynamicParticle();
    const G4Material* aMaterial = aTrack->GetMaterial();

    G4int materialIndex = aMaterial->GetIndex();
    G4int tracknumber = aStep.GetTrack()->GetTrackID();

    G4ThreeVector x0 = pPreStepPoint->GetPosition();
    G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
    G4double t0 = pPreStepPoint->GetGlobalTime();

    G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();

    G4MaterialPropertyVector* Fast_Intensity =
      aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
    G4MaterialPropertyVector* Slow_Intensity =
      aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");

    if (!Fast_Intensity && !Slow_Intensity) return 1;

    if (!bPropagate) return 0;

    // Get the visibility vector for this point
    assert(fPVS);

    G4int nscnt = 1;
    if (Fast_Intensity && Slow_Intensity) nscnt = 2;

    double Num = 0;
    double YieldRatio = 0;

    if (scintillationByParticleType) {
      // The scintillation response is a function of the energy deposited by particle
      // types.

      // Get the definition of the current particle
      G4ParticleDefinition* pDef = aParticle->GetDefinition();

      // Obtain the G4MaterialPropertyVectory containing the scintillation light yield as
      // a function of the deposited energy for the current particle type

      // Protons
      if (pDef == G4Proton::ProtonDefinition()) {
        YieldRatio = aMaterialPropertiesTable->GetConstProperty("PROTONYIELDRATIO");
      }
      // Muons
      else if (pDef == G4MuonPlus::MuonPlusDefinition() ||
               pDef == G4MuonMinus::MuonMinusDefinition()) {
        YieldRatio = aMaterialPropertiesTable->GetConstProperty("MUONYIELDRATIO");
      }
      // Pions
      else if (pDef == G4PionPlus::PionPlusDefinition() ||
               pDef == G4PionMinus::PionMinusDefinition()) {
        YieldRatio = aMaterialPropertiesTable->GetConstProperty("PIONYIELDRATIO");
      }
      // Kaons
      else if (pDef == G4KaonPlus::KaonPlusDefinition() ||
               pDef == G4KaonMinus::KaonMinusDefinition()) {
        YieldRatio = aMaterialPropertiesTable->GetConstProperty("KAONYIELDRATIO");
      }
      // Alphas
      else if (pDef == G4Alpha::AlphaDefinition()) {
        YieldRatio = aMaterialPropertiesTable->GetConstProperty("ALPHAYIELDRATIO");
      }
      // Electrons (must also account for shell-binding energy attributed to gamma from
      // standard PhotoElectricEffect)
      else if (pDef == G4Electron::ElectronDefinition() || pDef == G4Gamma::GammaDefinition()) {
        YieldRatio = aMaterialPropertiesTable->GetConstProperty("ELECTRONYIELDRATIO");
      }
      // Default for particles not enumerated/listed above
      else {
        YieldRatio = aMaterialPropertiesTable->GetConstProperty("ELECTRONYIELDRATIO");
      }
      // If the user has not specified yields for (p,d,t,a,carbon) then these unspecified
      // particles will default to the electron's scintillation yield
      if (YieldRatio == 0) {
        YieldRatio = aMaterialPropertiesTable->GetConstProperty("ELECTRONYIELDRATIO");
      }
    }

    geo::Point_t const ScintPoint = {x0[0] / CLHEP::cm, x0[1] / CLHEP::cm, x0[2] / CLHEP::cm};
    if (fOnlyActiveVolume && !isScintInActiveVolume(ScintPoint)) return 0;
    const phot::MappedCounts_t& Visibilities = fPVS->GetAllVisibilities(ScintPoint);

    phot::MappedCounts_t ReflVisibilities;

    // Store timing information in the object for use in propagationTime method
    if (fPVS->StoreReflected()) {
      ReflVisibilities = fPVS->GetAllVisibilities(ScintPoint, true);
      if (fPVS->StoreReflT0()) ReflT0s = fPVS->GetReflT0s(ScintPoint);
    }
    if (fPVS->IncludeParPropTime()) { ParPropTimeTF1 = fPVS->GetTimingTF1(ScintPoint); }

    for (G4int scnt = 1; scnt <= nscnt; scnt++) {
      G4double ScintillationTime = 0. * CLHEP::ns;
      G4double ScintillationRiseTime = 0. * CLHEP::ns;
      G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL;
      if (scnt == 1) {
        if (nscnt == 1) {
          if (Fast_Intensity) {
            ScintillationTime = aMaterialPropertiesTable->GetConstProperty("FASTTIMECONSTANT");
            if (fFiniteRiseTime) {
              ScintillationRiseTime =
                aMaterialPropertiesTable->GetConstProperty("FASTSCINTILLATIONRISETIME");
            }
            ScintillationIntegral =
              (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
          }
          if (Slow_Intensity) {
            ScintillationTime = aMaterialPropertiesTable->GetConstProperty("SLOWTIMECONSTANT");
            if (fFiniteRiseTime) {
              ScintillationRiseTime =
                aMaterialPropertiesTable->GetConstProperty("SLOWSCINTILLATIONRISETIME");
            }
            ScintillationIntegral =
              (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
          }
        } //endif nscnt=1
        else {
          if (YieldRatio == 0)
            YieldRatio = aMaterialPropertiesTable->GetConstProperty("YIELDRATIO");

          if (ExcitationRatio == 1.0) { Num = std::min(YieldRatio, 1.0) * MeanNumberOfPhotons; }
          else {
            Num = std::min(ExcitationRatio, 1.0) * MeanNumberOfPhotons;
          }
          ScintillationTime = aMaterialPropertiesTable->GetConstProperty("FASTTIMECONSTANT");
          if (fFiniteRiseTime) {
            ScintillationRiseTime =
              aMaterialPropertiesTable->GetConstProperty("FASTSCINTILLATIONRISETIME");
          }
          ScintillationIntegral =
            (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
        } //endif nscnt!=1
      }   //end scnt=1

      else {
        Num = MeanNumberOfPhotons - Num;
        ScintillationTime = aMaterialPropertiesTable->GetConstProperty("SLOWTIMECONSTANT");
        if (fFiniteRiseTime) {
          ScintillationRiseTime =
            aMaterialPropertiesTable->GetConstProperty("SLOWSCINTILLATIONRISETIME");
        }
        ScintillationIntegral =
          (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
      }

      if (!ScintillationIntegral) continue;

      if (!Visibilities && !usesSemiAnalyticModel()) continue;

      // detected photons from direct light
      std::map<size_t, int> DetectedNum;
      if (Visibilities && !usesSemiAnalyticModel()) {
        for (size_t const OpDet : util::counter(fPVS->NOpChannels())) {
          if (fOpaqueCathode && !isOpDetInSameTPC(ScintPoint, fOpDetCenter.at(OpDet))) continue;
          int const DetThis = std::round(G4Poisson(Visibilities[OpDet] * Num));
          if (DetThis > 0) DetectedNum[OpDet] = DetThis;
        }
      }
      else {
        detectedDirectHits(DetectedNum, Num, ScintPoint);
      }

      // detected photons from reflected light
      std::map<size_t, int> ReflDetectedNum;
      if (fPVS->StoreReflected()) {
        if (!usesSemiAnalyticModel()) {
          for (size_t const OpDet : util::counter(fPVS->NOpChannels())) {
            if (fOpaqueCathode && !isOpDetInSameTPC(ScintPoint, fOpDetCenter.at(OpDet))) continue;
            int const ReflDetThis = std::round(G4Poisson(ReflVisibilities[OpDet] * Num));
            if (ReflDetThis > 0) ReflDetectedNum[OpDet] = ReflDetThis;
          }
        }
        else {
          detectedReflecHits(ReflDetectedNum, Num, ScintPoint);
        }
      }

      std::vector<double> arrival_time_dist;
      // Now we run through each PMT figuring out num of detected photons
      for (size_t Reflected = 0; Reflected <= 1; ++Reflected) {
        // Only do the reflected loop if we have reflected visibilities
        if (Reflected && !fPVS->StoreReflected()) continue;

        std::map<size_t, int>::const_iterator itstart;
        std::map<size_t, int>::const_iterator itend;
        if (Reflected) {
          itstart = ReflDetectedNum.begin();
          itend = ReflDetectedNum.end();
        }
        else {
          itstart = DetectedNum.begin();
          itend = DetectedNum.end();
        }
        for (auto itdetphot = itstart; itdetphot != itend; ++itdetphot) {
          const size_t OpChannel = itdetphot->first;
          const int NPhotons = itdetphot->second;

          // Set up the OpDetBTR information
          sim::OpDetBacktrackerRecord tmpOpDetBTRecord(OpChannel);
          int thisG4TrackID = ParticleListAction::GetCurrentTrackID();
          double xyzPos[3];
          average_position(aStep, xyzPos);
          double Edeposited = 0;
          if (scintillationByParticleType) {
            // We use this when it is the only sensical information. It may be of limited
            // use to end users.
            Edeposited = aStep.GetTotalEnergyDeposit();
          }
          else if (emSaturation) {
            // If Birk Coefficient used, log VisibleEnergies.
            Edeposited =
              larg4::IonizationAndScintillation::Instance()->VisibleEnergyDeposit() / CLHEP::MeV;
          }
          else {
            // We use this when it is the only sensical information. It may be of limited
            // use to end users.
            Edeposited = aStep.GetTotalEnergyDeposit();
          }

          // Get the transport time distribution
          arrival_time_dist.resize(NPhotons);
          propagationTime(arrival_time_dist, x0, OpChannel, Reflected);

          // We need to split the energy up by the number of photons so that we never try
          // to write a 0 energy.
          Edeposited = Edeposited / double(NPhotons);

          // Loop through the photons
          for (G4int i = 0; i < NPhotons; ++i) {
            G4double Time = t0 + scint_time(aStep, ScintillationTime, ScintillationRiseTime) +
                            arrival_time_dist[i] * CLHEP::ns;

            // Always store the BTR
            tmpOpDetBTRecord.AddScintillationPhotons(thisG4TrackID, Time, 1, xyzPos, Edeposited);

            // Store as lite photon or as OnePhoton
            if (lgp->UseLitePhotons()) {
              fst->AddLitePhoton(OpChannel, static_cast<int>(Time), 1, Reflected);
            }
            else {
              // The sim photon in this case stores its production point and time
              TVector3 PhotonPosition(x0[0], x0[1], x0[2]);

              float PhotonEnergy = 0;
              if (Reflected)
                PhotonEnergy = reemission_energy() * CLHEP::eV;
              else
                PhotonEnergy = 9.7 * CLHEP::eV; // 9.7 eV peak of VUV emission spectrum

              // Make a photon object for the collection
              sim::OnePhoton PhotToAdd;
              PhotToAdd.InitialPosition = PhotonPosition;
              PhotToAdd.Energy = PhotonEnergy;
              PhotToAdd.Time = Time;
              PhotToAdd.SetInSD = false;
              PhotToAdd.MotherTrackID = tracknumber;

              fst->AddPhoton(OpChannel, std::move(PhotToAdd), Reflected);
            }
          }
          fst->AddOpDetBacktrackerRecord(tmpOpDetBTRecord, Reflected);
        }
      }
    }
    return 0;
  }

  // BuildThePhysicsTable for the scintillation process
  // --------------------------------------------------
  //
  void OpFastScintillation::BuildThePhysicsTable()
  {
    if (theFastIntegralTable && theSlowIntegralTable) return;

    const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
    G4int numOfMaterials = G4Material::GetNumberOfMaterials();

    // create new physics table
    if (!theFastIntegralTable)
      theFastIntegralTable = std::make_unique<G4PhysicsTable>(numOfMaterials);
    if (!theSlowIntegralTable)
      theSlowIntegralTable = std::make_unique<G4PhysicsTable>(numOfMaterials);

    // loop for materials
    for (G4int i = 0; i < numOfMaterials; i++) {
      G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector = new G4PhysicsOrderedFreeVector();
      G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector = new G4PhysicsOrderedFreeVector();

      // Retrieve vector of scintillation wavelength intensity for the material from the
      // material's optical properties table.
      G4Material* aMaterial = (*theMaterialTable)[i];

      G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();

      if (aMaterialPropertiesTable) {

        G4MaterialPropertyVector* theFastLightVector =
          aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");

        if (theFastLightVector) {
          // Retrieve the first intensity point in vector of (photon energy, intensity)
          // pairs
          G4double currentIN = (*theFastLightVector)[0];
          if (currentIN >= 0.0) {
            // Create first (photon energy, Scintillation Integral) pair
            G4double currentPM = theFastLightVector->Energy(0);
            G4double currentCII = 0.0;

            aPhysicsOrderedFreeVector->InsertValues(currentPM, currentCII);

            // Set previous values to current ones prior to loop
            G4double prevPM = currentPM;
            G4double prevCII = currentCII;
            G4double prevIN = currentIN;

            // loop over all (photon energy, intensity) pairs stored for this material
            for (size_t i = 1; i < theFastLightVector->GetVectorLength(); i++) {
              currentPM = theFastLightVector->Energy(i);
              currentIN = (*theFastLightVector)[i];

              currentCII = 0.5 * (prevIN + currentIN);

              currentCII = prevCII + (currentPM - prevPM) * currentCII;

              aPhysicsOrderedFreeVector->InsertValues(currentPM, currentCII);

              prevPM = currentPM;
              prevCII = currentCII;
              prevIN = currentIN;
            }
          }
        }

        G4MaterialPropertyVector* theSlowLightVector =
          aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
        if (theSlowLightVector) {
          // Retrieve the first intensity point in vector of (photon energy, intensity)
          // pairs
          G4double currentIN = (*theSlowLightVector)[0];
          if (currentIN >= 0.0) {
            // Create first (photon energy, Scintillation Integral) pair
            G4double currentPM = theSlowLightVector->Energy(0);
            G4double currentCII = 0.0;

            bPhysicsOrderedFreeVector->InsertValues(currentPM, currentCII);

            // Set previous values to current ones prior to loop
            G4double prevPM = currentPM;
            G4double prevCII = currentCII;
            G4double prevIN = currentIN;

            // loop over all (photon energy, intensity) pairs stored for this material
            for (size_t i = 1; i < theSlowLightVector->GetVectorLength(); i++) {
              currentPM = theSlowLightVector->Energy(i);
              currentIN = (*theSlowLightVector)[i];

              currentCII = 0.5 * (prevIN + currentIN);

              currentCII = prevCII + (currentPM - prevPM) * currentCII;

              bPhysicsOrderedFreeVector->InsertValues(currentPM, currentCII);

              prevPM = currentPM;
              prevCII = currentCII;
              prevIN = currentIN;
            }
          }
        }
      }
      // The scintillation integral(s) for a given material will be inserted in the
      // table(s) according to the position of the material in the material table.
      theFastIntegralTable->insertAt(i, aPhysicsOrderedFreeVector);
      theSlowIntegralTable->insertAt(i, bPhysicsOrderedFreeVector);
    }
  }

  // Called by the user to set the scintillation yield as a function of energy deposited
  // by particle type
  void OpFastScintillation::SetScintillationByParticleType(const G4bool scintType)
  {
    if (emSaturation) {
      G4Exception("OpFastScintillation::SetScintillationByParticleType",
                  "Scint02",
                  JustWarning,
                  "Redefinition: Birks Saturation is replaced by ScintillationByParticleType!");
      RemoveSaturation();
    }
    scintillationByParticleType = scintType;
  }

  G4double OpFastScintillation::GetMeanFreePath(const G4Track&,
                                                G4double,
                                                G4ForceCondition* condition)
  {
    *condition = StronglyForced;
    return DBL_MAX;
  }

  G4double OpFastScintillation::GetMeanLifeTime(const G4Track&, G4ForceCondition* condition)
  {
    *condition = Forced;
    return DBL_MAX;
  }

  G4double OpFastScintillation::scint_time(const G4Step& aStep,
                                           G4double ScintillationTime,
                                           G4double ScintillationRiseTime) const
  {
    G4StepPoint const* pPreStepPoint = aStep.GetPreStepPoint();
    G4StepPoint const* pPostStepPoint = aStep.GetPostStepPoint();
    G4double avgVelocity = (pPreStepPoint->GetVelocity() + pPostStepPoint->GetVelocity()) / 2.;
    G4double deltaTime = aStep.GetStepLength() / avgVelocity;
    if (ScintillationRiseTime == 0.0) {
      deltaTime = deltaTime - ScintillationTime * std::log(G4UniformRand());
    }
    else {
      deltaTime = deltaTime + sample_time(ScintillationRiseTime, ScintillationTime);
    }
    return deltaTime;
  }

  void OpFastScintillation::propagationTime(std::vector<double>& arrival_time_dist,
                                            G4ThreeVector x0,
                                            const size_t OpChannel,
                                            bool Reflected)
  {
    assert(fPVS);
    if (fPVS->IncludeParPropTime() && fPVS->IncludePropTime()) {
      throw cet::exception("OpFastScintillation")
        << "Cannot have both propagation time models simultaneously.";
    }
    else if (fPVS->IncludeParPropTime() &&
             !(ParPropTimeTF1 && (ParPropTimeTF1[OpChannel].GetNdim() == 1))) {
      // Warning: TF1::GetNdim()==1 will tell us if the TF1 is really defined or it is the
      // default one.  This will fix a segfault when using timing and interpolation.
      G4cout << "WARNING: Requested parameterized timing, but no function found. Not applying "
                "propagation time."
             << G4endl;
    }
    else if (fPVS->IncludeParPropTime()) {
      if (Reflected)
        throw cet::exception("OpFastScintillation")
          << "No parameterized propagation time for reflected light";
      for (size_t i = 0; i < arrival_time_dist.size(); ++i) {
        arrival_time_dist[i] = ParPropTimeTF1[OpChannel].GetRandom();
      }
    }
    else if (fPVS->IncludePropTime()) {
      // Get VUV photons arrival time distribution from the parametrization
      geo::Point_t const& opDetCenter = fOpDetCenter.at(OpChannel);
      if (!Reflected) {
        const G4ThreeVector OpDetPoint(
          opDetCenter.X() * CLHEP::cm, opDetCenter.Y() * CLHEP::cm, opDetCenter.Z() * CLHEP::cm);
        double distance_in_cm = (x0 - OpDetPoint).mag() / CLHEP::cm; // this must be in CENTIMETERS!
        double cosine = std::abs(x0[0] * CLHEP::cm - OpDetPoint[0] * CLHEP::cm) / distance_in_cm;
        double theta = fast_acos(cosine) * 180. / CLHEP::pi;
        int angle_bin = theta / fangle_bin_timing_vuv;
        getVUVTimes(arrival_time_dist, distance_in_cm, angle_bin); // in ns
      }
      else {
        TVector3 const ScintPoint(x0[0] / CLHEP::cm, x0[1] / CLHEP::cm, x0[2] / CLHEP::cm); // in cm
        getVISTimes(arrival_time_dist, ScintPoint, geo::vect::toTVector3(opDetCenter));     // in ns
      }
    }
  }

  G4double OpFastScintillation::sample_time(const G4double tau1, const G4double tau2) const
  {
    // tau1: rise time and tau2: decay time
    while (1) {
      // two random numbers
      G4double ran1 = G4UniformRand();
      G4double ran2 = G4UniformRand();
      //
      // exponential distribution as envelope function: very efficient
      //
      G4double d = (tau1 + tau2) / tau2;
      // make sure the envelope function is always larger than the bi-exponential
      G4double t = -1.0 * tau2 * std::log(1 - ran1);
      G4double g = d * single_exp(t, tau2);
      if (ran2 <= bi_exp(t, tau1, tau2) / g) return t;
    }
    return -1.0;
  }

  double OpFastScintillation::reemission_energy() const
  {
    return fTPBEm->fire() * ((*(--tpbemission.end())).first - (*tpbemission.begin()).first) +
           (*tpbemission.begin()).first;
  }

  void OpFastScintillation::average_position(G4Step const& aStep, double* xyzPos) const
  {
    CLHEP::Hep3Vector prePoint = (aStep.GetPreStepPoint())->GetPosition();
    CLHEP::Hep3Vector postPoint = (aStep.GetPostStepPoint())->GetPosition();
    xyzPos[0] = (((prePoint.getX() + postPoint.getX()) / 2.0) / CLHEP::cm);
    xyzPos[1] = (((prePoint.getY() + postPoint.getY()) / 2.0) / CLHEP::cm);
    xyzPos[2] = (((prePoint.getZ() + postPoint.getZ()) / 2.0) / CLHEP::cm);
  }

  // New Parametrization code parameterisation generation function
  void OpFastScintillation::generateParam(const size_t index, const size_t angle_bin)
  {
    // get distance
    double distance_in_cm = (index * fstep_size) + fmin_d;

    // time range
    const double signal_t_range = 5000.; // TODO: unhardcode

    // parameterisation TF1
    TF1 fVUVTiming;

    // For very short distances the time correction is just a shift
    double t_direct_mean = distance_in_cm / fvuv_vgroup_mean;
    double t_direct_min = distance_in_cm / fvuv_vgroup_max;

    // Defining the model function(s) describing the photon transportation timing vs
    // distance Getting the landau parameters from the time parametrization
    std::array<double, 3> pars_landau;
    interpolate3(pars_landau,
                 fparameters[0][0],
                 fparameters[2][angle_bin],
                 fparameters[3][angle_bin],
                 fparameters[1][angle_bin],
                 distance_in_cm,
                 true);
    // Deciding which time model to use (depends on the distance) defining useful times
    // for the VUV arrival time shapes
    if (distance_in_cm >= finflexion_point_distance) {
      double pars_far[4] = {t_direct_min, pars_landau[0], pars_landau[1], pars_landau[2]};
      // Set model: Landau
      fVUVTiming = TF1("fVUVTiming", model_far, 0, signal_t_range, 4);
      fVUVTiming.SetParameters(pars_far);
    }
    else {
      // Set model: Landau + Exponential
      fVUVTiming = TF1("fVUVTiming", model_close, 0, signal_t_range, 7);
      // Exponential parameters
      double pars_expo[2];
      // Getting the exponential parameters from the time parametrization
      pars_expo[1] =
        interpolate(fparameters[4][0], fparameters[5][angle_bin], distance_in_cm, true);
      pars_expo[0] =
        interpolate(fparameters[4][0], fparameters[6][angle_bin], distance_in_cm, true);
      pars_expo[0] *= pars_landau[2];
      pars_expo[0] = std::log(pars_expo[0]);
      // this is to find the intersection point between the two functions:
      TF1 fint = TF1("fint", finter_d, pars_landau[0], 4 * t_direct_mean, 5);
      double parsInt[5] = {
        pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
      fint.SetParameters(parsInt);
      double t_int = fint.GetMinimumX();
      double minVal = fint.Eval(t_int);
      // the functions must intersect - output warning if they don't
      if (minVal > 0.015) {
        std::cout << "WARNING: Parametrization of VUV light discontinuous for distance = "
                  << distance_in_cm << std::endl;
        std::cout << "WARNING: This shouldn't be happening " << std::endl;
      }
      double parsfinal[7] = {t_int,
                             pars_landau[0],
                             pars_landau[1],
                             pars_landau[2],
                             pars_expo[0],
                             pars_expo[1],
                             t_direct_min};
      fVUVTiming.SetParameters(parsfinal);
    }

    // set the number of points used to sample parameterisation for shorter distances,
    // peak is sharper so more sensitive sampling required
    int fsampling; // TODO: unhardcode
    if (distance_in_cm < 50) { fsampling = 10000; }
    else if (distance_in_cm < 100) {
      fsampling = 5000;
    }
    else {
      fsampling = 1000;
    }
    fVUVTiming.SetNpx(fsampling);

    // calculate max and min distance relevant to sample parameterisation max
    const size_t nq_max = 1;
    double xq_max[nq_max];
    double yq_max[nq_max];
    xq_max[0] = 0.995; // include 99.5%
    fVUVTiming.GetQuantiles(nq_max, yq_max, xq_max);
    double max = yq_max[0];
    // min
    double min = t_direct_min;

    // store TF1 and min/max, this allows identical TF1 to be used every time sampling the
    // first call of GetRandom generates the timing sampling and stores it in the TF1
    // object, this is the slow part all subsequent calls check if it has been generated
    // previously and are ~100+ times quicker
    VUV_timing[angle_bin][index] = fVUVTiming;
    VUV_max[angle_bin][index] = max;
    VUV_min[angle_bin][index] = min;
  }

  // VUV arrival times calculation function
  void OpFastScintillation::getVUVTimes(std::vector<double>& arrivalTimes,
                                        const double distance,
                                        const size_t angle_bin)
  {
    if (distance < fmin_d) {
      // times are fixed shift i.e. direct path only
      double t_prop_correction = distance / fvuv_vgroup_mean;
      for (size_t i = 0; i < arrivalTimes.size(); ++i) {
        arrivalTimes[i] = t_prop_correction;
      }
    }
    else { // distance >= fmin_d
      // determine nearest parameterisation in discretisation
      int index = std::round((distance - fmin_d) / fstep_size);
      // check whether required parameterisation has been generated, generating if not
      if (VUV_timing[angle_bin][index].GetNdim() == 0) { generateParam(index, angle_bin); }
      // randomly sample parameterisation for each photon
      for (size_t i = 0; i < arrivalTimes.size(); ++i) {
        arrivalTimes[i] = VUV_timing[angle_bin][index].GetRandom(VUV_min[angle_bin][index],
                                                                 VUV_max[angle_bin][index]);
      }
    }
  }

  // VIS arrival times calculation functions
  void OpFastScintillation::getVISTimes(std::vector<double>& arrivalTimes,
                                        const TVector3& ScintPoint,
                                        const TVector3& OpDetPoint)
  {
    // *************************************************************************************************
    //     Calculation of earliest arrival times and corresponding unsmeared distribution
    // *************************************************************************************************

    // set plane_depth for correct TPC:
    double plane_depth;
    if (ScintPoint[0] < 0) { plane_depth = -fplane_depth; }
    else {
      plane_depth = fplane_depth;
    }

    // calculate point of reflection for shortest path
    TVector3 bounce_point(plane_depth, ScintPoint[1], ScintPoint[2]);

    // calculate distance travelled by VUV light and by vis light
    double VUVdist = (bounce_point - ScintPoint).Mag();
    double Visdist = (OpDetPoint - bounce_point).Mag();

    // calculate times taken by VUV part of path
    int angle_bin_vuv = 0; // on-axis by definition
    getVUVTimes(arrivalTimes, VUVdist, angle_bin_vuv);

    // add visible direct path transport time
    for (size_t i = 0; i < arrivalTimes.size(); ++i) {
      arrivalTimes[i] += Visdist / fvis_vmean;
    }

    // *************************************************************************************************
    //      Smearing of arrival time distribution
    // *************************************************************************************************
    // calculate fastest time possible
    // vis part
    double vis_time = Visdist / fvis_vmean;
    // vuv part
    double vuv_time;
    if (VUVdist < fmin_d) { vuv_time = VUVdist / fvuv_vgroup_max; }
    else {
      // find index of required parameterisation
      const size_t index = std::round((VUVdist - fmin_d) / fstep_size);
      // find shortest time
      vuv_time = VUV_min[angle_bin_vuv][index];
    }
    // sum
    double fastest_time = vis_time + vuv_time;

    // calculate angle theta between bound_point and optical detector
    double cosine_theta = std::abs(OpDetPoint[0] - bounce_point[0]) / Visdist;
    double theta = fast_acos(cosine_theta) * 180. / CLHEP::pi;

    // determine smearing parameters using interpolation of generated points:
    // 1). tau = exponential smearing factor, varies with distance and angle
    // 2). cutoff = largest smeared time allowed, preventing excessively large times
    //     caused by exponential distance to cathode
    double distance_cathode_plane = std::abs(plane_depth - ScintPoint[0]);
    // angular bin
    size_t theta_bin = theta / fangle_bin_timing_vis;
    // radial distance from centre of TPC (y,z plane)
    double r = std::sqrt(std::pow(ScintPoint.Y() - fcathode_centre.Z(), 2) +
                         std::pow(ScintPoint.Y() - fcathode_centre.Z(), 2));

    // cut-off and tau
    // cut-off
    // interpolate in d_c for each r bin
    std::vector<double> interp_vals(fcut_off_pars[theta_bin].size(), 0.0);
    for (size_t i = 0; i < fcut_off_pars[theta_bin].size(); i++) {
      interp_vals[i] =
        interpolate(fdistances_refl, fcut_off_pars[theta_bin][i], distance_cathode_plane, true);
    }
    // interpolate in r
    double cutoff = interpolate(fradial_distances_refl, interp_vals, r, true);

    // tau
    // interpolate in x for each r bin
    std::vector<double> interp_vals_tau(ftau_pars[theta_bin].size(), 0.0);
    for (size_t i = 0; i < ftau_pars[theta_bin].size(); i++) {
      interp_vals_tau[i] =
        interpolate(fdistances_refl, ftau_pars[theta_bin][i], distance_cathode_plane, true);
    }
    // interpolate in r
    double tau = interpolate(fradial_distances_refl, interp_vals_tau, r, true);

    if (tau < 0) { tau = 0; } // failsafe if tau extrapolate goes wrong

    // apply smearing:
    for (size_t i = 0; i < arrivalTimes.size(); ++i) {
      double arrival_time_smeared;
      // if time is already greater than cutoff, do not apply smearing
      if (arrivalTimes[i] >= cutoff) { continue; }
      // otherwise smear
      else {
        unsigned int counter = 0;
        // loop until time generated is within cutoff limit; most are within single
        // attempt, very few take more than two
        do {
          // don't attempt smearings too many times
          if (counter >= 10) {                      // TODO: unhardcode
            arrival_time_smeared = arrivalTimes[i]; // don't smear
            break;
          }
          else {
            // generate random number in appropriate range
            double x = gRandom->Uniform(0.5, 1.0); // TODO: unhardcode
            // apply the exponential smearing
            arrival_time_smeared =
              arrivalTimes[i] + (arrivalTimes[i] - fastest_time) * (std::pow(x, -tau) - 1);
          }
          counter++;
        } while (arrival_time_smeared > cutoff);
      }
      arrivalTimes[i] = arrival_time_smeared;
    }
  }

  // ---------------------------------------------------------------------------
  bool OpFastScintillation::usesSemiAnalyticModel() const
  {
    return fUseNhitsModel;
  }

  // ---------------------------------------------------------------------------
  void OpFastScintillation::detectedDirectHits(std::map<size_t, int>& DetectedNum,
                                               const double Num,
                                               geo::Point_t const& ScintPoint) const
  {
    for (size_t const OpDet : util::counter(fPVS->NOpChannels())) {
      if (!isOpDetInSameTPC(ScintPoint, fOpDetCenter.at(OpDet))) continue;

      // set detector struct for solid angle function
      const OpFastScintillation::OpticalDetector op{fOpDetHeight.at(OpDet),
                                                    fOpDetLength.at(OpDet),
                                                    fOpDetCenter.at(OpDet),
                                                    fOpDetType.at(OpDet)};
      const int DetThis = VUVHits(Num, ScintPoint, op);
      if (DetThis > 0) { DetectedNum[OpDet] = DetThis; }
    }
  }

  void OpFastScintillation::detectedReflecHits(std::map<size_t, int>& ReflDetectedNum,
                                               const double Num,
                                               geo::Point_t const& ScintPoint) const
  {
    // 1). calculate total number of hits of VUV photons on reflective foils via solid
    //     angle + Gaisser-Hillas corrections:

    // set plane_depth for correct TPC:
    double plane_depth;
    if (ScintPoint.X() < 0.) { plane_depth = -fplane_depth; }
    else {
      plane_depth = fplane_depth;
    }

    // get scintpoint coords relative to centre of cathode plane
    std::array<double, 3> ScintPoint_relative = {std::abs(ScintPoint.X() - plane_depth),
                                                 std::abs(ScintPoint.Y() - fcathode_centre.Y()),
                                                 std::abs(ScintPoint.Z() - fcathode_centre.Z())};
    // calculate solid angle of cathode from the scintillation point
    double solid_angle_cathode = Rectangle_SolidAngle(fcathode_plane, ScintPoint_relative);

    // calculate distance and angle between ScintPoint and hotspot; vast majority of hits
    // in hotspot region directly infront of scintpoint, therefore consider attenuation
    // for this distance and on axis GH instead of for the centre coordinate
    double distance_cathode = std::abs(plane_depth - ScintPoint.X());
    // calculate hits on cathode plane via geometric acceptance
    double cathode_hits_geo = std::exp(-1. * distance_cathode / fL_abs_vuv) *
                              (solid_angle_cathode / (4. * CLHEP::pi)) * Num;
    // determine Gaisser-Hillas correction including border effects use flat correction
    double r = std::sqrt(std::pow(ScintPoint.Y() - fcathode_centre.Y(), 2) +
                         std::pow(ScintPoint.Z() - fcathode_centre.Z(), 2));
    double pars_ini[4] = {0, 0, 0, 0};
    double s1 = 0;
    double s2 = 0;
    double s3 = 0;
    if (fIsFlatPDCorr) {
      pars_ini[0] = fGHvuvpars_flat[0][0];
      pars_ini[1] = fGHvuvpars_flat[1][0];
      pars_ini[2] = fGHvuvpars_flat[2][0];
      pars_ini[3] = fGHvuvpars_flat[3][0];
      s1 = interpolate(fborder_corr_angulo_flat, fborder_corr_flat[0], 0, true);
      s2 = interpolate(fborder_corr_angulo_flat, fborder_corr_flat[1], 0, true);
      s3 = interpolate(fborder_corr_angulo_flat, fborder_corr_flat[2], 0, true);
    }
    else
      std::cout
        << "Error: flat optical detector VUV correction required for reflected semi-analytic hits."
        << std::endl;

    // add border correction
    pars_ini[0] = pars_ini[0] + s1 * r;
    pars_ini[1] = pars_ini[1] + s2 * r;
    pars_ini[2] = pars_ini[2] + s3 * r;
    pars_ini[3] = pars_ini[3];

    // calculate corrected number of hits
    double GH_correction = Gaisser_Hillas(distance_cathode, pars_ini);
    const double cathode_hits_rec = GH_correction * cathode_hits_geo;

    // detemine hits on each PD
    const std::array<double, 3> hotspot = {plane_depth, ScintPoint.Y(), ScintPoint.Z()};
    for (size_t const OpDet : util::counter(fPVS->NOpChannels())) {
      if (!isOpDetInSameTPC(ScintPoint, fOpDetCenter.at(OpDet))) continue;

      // set detector struct for solid angle function
      const OpFastScintillation::OpticalDetector op{fOpDetHeight.at(OpDet),
                                                    fOpDetLength.at(OpDet),
                                                    fOpDetCenter.at(OpDet),
                                                    fOpDetType.at(OpDet)};

      int const ReflDetThis = VISHits(ScintPoint, op, cathode_hits_rec, hotspot);
      if (ReflDetThis > 0) { ReflDetectedNum[OpDet] = ReflDetThis; }
    }
  }

  // VUV semi-analytic hits calculation
  int OpFastScintillation::VUVHits(const double Nphotons_created,
                                   geo::Point_t const& ScintPoint_v,
                                   OpticalDetector const& opDet) const
  {
    // the interface has been converted into geo::Point_t, the implementation not yet
    std::array<double, 3U> ScintPoint;
    std::array<double, 3U> OpDetPoint;
    geo::vect::fillCoords(ScintPoint, ScintPoint_v);
    geo::vect::fillCoords(OpDetPoint, opDet.OpDetPoint);

    // distance and angle between ScintPoint and OpDetPoint
    double distance = dist(&ScintPoint[0], &OpDetPoint[0], 3);
    double cosine = std::abs(ScintPoint[0] - OpDetPoint[0]) / distance;
    // double theta = std::acos(cosine) * 180. / CLHEP::pi;
    double theta = fast_acos(cosine) * 180. / CLHEP::pi;

    // calculate solid angle:
    double solid_angle = 0;
    // Arapucas/Bars (rectangle)
    if (opDet.type == 0) {
      // get scintillation point coordinates relative to arapuca window centre
      std::array<double, 3> ScintPoint_rel = {std::abs(ScintPoint[0] - OpDetPoint[0]),
                                              std::abs(ScintPoint[1] - OpDetPoint[1]),
                                              std::abs(ScintPoint[2] - OpDetPoint[2])};
      // calculate solid angle
      solid_angle = Rectangle_SolidAngle(Dims{opDet.h, opDet.w}, ScintPoint_rel);
    }
    // PMTs (dome)
    else if (opDet.type == 1) {
      solid_angle = Omega_Dome_Model(distance, theta);
    }
    // PMTs (disk approximation)
    else if (opDet.type == 2) {
      // offset in z-y plane
      double d = dist(&ScintPoint[1], &OpDetPoint[1], 2);
      // drift distance (in x)
      double h = std::abs(ScintPoint[0] - OpDetPoint[0]);
      // Solid angle of a disk
      solid_angle = Disk_SolidAngle(d, h, fradius);
    }
    else {
      std::cout << "Error: Invalid optical detector type. 0 = rectangular, 1 = dome, 2 = disk"
                << std::endl;
    }

    // calculate number of photons hits by geometric acceptance: accounting for solid
    // angle and LAr absorption length
    double hits_geo =
      std::exp(-1. * distance / fL_abs_vuv) * (solid_angle / (4 * CLHEP::pi)) * Nphotons_created;

    // apply Gaisser-Hillas correction for Rayleigh scattering distance and angular
    // dependence offset angle bin
    const size_t j = (theta / fdelta_angulo_vuv);

    // determine GH parameters, accounting for border effects radial distance from centre
    // of detector (Y-Z)

    // FIXME (KJK): Original implementation only accounts for the Y component
    double r = std::abs(ScintPoint[1] - fcathode_centre.Y());
    double pars_ini[4] = {0, 0, 0, 0};
    double s1 = 0;
    double s2 = 0;
    double s3 = 0;
    // flat PDs
    if ((opDet.type == 0 || opDet.type == 2) && fIsFlatPDCorr) {
      pars_ini[0] = fGHvuvpars_flat[0][j];
      pars_ini[1] = fGHvuvpars_flat[1][j];
      pars_ini[2] = fGHvuvpars_flat[2][j];
      pars_ini[3] = fGHvuvpars_flat[3][j];
      s1 = interpolate(fborder_corr_angulo_flat, fborder_corr_flat[0], theta, true);
      s2 = interpolate(fborder_corr_angulo_flat, fborder_corr_flat[1], theta, true);
      s3 = interpolate(fborder_corr_angulo_flat, fborder_corr_flat[2], theta, true);
    }
    // dome PDs
    else if (opDet.type == 1 && fIsDomePDCorr) {
      pars_ini[0] = fGHvuvpars_dome[0][j];
      pars_ini[1] = fGHvuvpars_dome[1][j];
      pars_ini[2] = fGHvuvpars_dome[2][j];
      pars_ini[3] = fGHvuvpars_dome[3][j];
      s1 = interpolate(fborder_corr_angulo_dome, fborder_corr_dome[0], theta, true);
      s2 = interpolate(fborder_corr_angulo_dome, fborder_corr_dome[1], theta, true);
      s3 = interpolate(fborder_corr_angulo_dome, fborder_corr_dome[2], theta, true);
    }
    else
      std::cout << "Error: Invalid optical detector type. 0 = rectangular, 1 = dome, 2 = disk. Or "
                   "corrections for chosen optical detector type missing."
                << std::endl;

    // add border correction
    pars_ini[0] = pars_ini[0] + s1 * r;
    pars_ini[1] = pars_ini[1] + s2 * r;
    pars_ini[2] = pars_ini[2] + s3 * r;
    pars_ini[3] = pars_ini[3];

    // calculate correction
    double GH_correction = Gaisser_Hillas(distance, pars_ini);

    // calculate number photons
    double hits_rec = GH_correction * hits_geo / cosine;
    int hits_vuv = std::round(G4Poisson(hits_rec));

    return hits_vuv;
  }

  // VIS hits semi-analytic model calculation
  int OpFastScintillation::VISHits(geo::Point_t const& ScintPoint_v,
                                   OpticalDetector const& opDet,
                                   const double cathode_hits_rec,
                                   const std::array<double, 3> hotspot) const
  {
    // 1). calculate total number of hits of VUV photons on reflective foils via solid
    //     angle + Gaisser-Hillas corrections.  Done outside as it doesn't depend on
    //     OpDetPoint

    // set plane_depth for correct TPC:
    double plane_depth;
    if (ScintPoint_v.X() < 0.) { plane_depth = -fplane_depth; }
    else {
      plane_depth = fplane_depth;
    }

    // 2). calculate number of these hits which reach the optical detector from the
    //     hotspot via solid angle

    // the interface has been converted into geo::Point_t, the implementation not yet
    std::array<double, 3U> ScintPoint;
    std::array<double, 3U> OpDetPoint;
    geo::vect::fillCoords(ScintPoint, ScintPoint_v);
    geo::vect::fillCoords(OpDetPoint, opDet.OpDetPoint);

    // calculate distances and angles for application of corrections; distance from
    // hotspot to optical detector
    double distance_vis = dist(&hotspot[0], &OpDetPoint[0], 3);
    //  angle between hotspot and optical detector
    double cosine_vis = std::abs(hotspot[0] - OpDetPoint[0]) / distance_vis;
    double theta_vis = fast_acos(cosine_vis) * 180. / CLHEP::pi;

    // calculate solid angle of optical detector
    double solid_angle_detector = 0;
    // rectangular aperture
    if (opDet.type == 0) {
      // get hotspot coordinates relative to opDet
      std::array<double, 3> emission_relative = {std::abs(hotspot[0] - OpDetPoint[0]),
                                                 std::abs(hotspot[1] - OpDetPoint[1]),
                                                 std::abs(hotspot[2] - OpDetPoint[2])};
      // calculate solid angle
      solid_angle_detector = Rectangle_SolidAngle(Dims{opDet.h, opDet.w}, emission_relative);
    }
    // dome aperture
    else if (opDet.type == 1) {
      solid_angle_detector = Omega_Dome_Model(distance_vis, theta_vis);
    }
    // disk aperture
    else if (opDet.type == 2) {
      // offset in z-y plane
      double d = dist(&hotspot[1], &OpDetPoint[1], 2);
      // drift distance (in x)
      double h = std::abs(hotspot[0] - OpDetPoint[0]);
      // calculate solid angle
      solid_angle_detector = Disk_SolidAngle(d, h, fradius);
    }
    else {
      std::cout << "Error: Invalid optical detector type. 0 = rectangular, 1 = dome, 2 = disk"
                << std::endl;
    }

    // calculate number of hits via geometeric acceptance
    double hits_geo = (solid_angle_detector / (2. * CLHEP::pi)) *
                      cathode_hits_rec; // 2*pi due to presence of reflective foils

    // determine correction factor, depending on PD type
    const size_t k = (theta_vis / fdelta_angulo_vis); // off-set angle bin
    // radial distance from centre of detector (Y-Z)
    // FIXME (KJK): Original implementation only accounts for the Y component
    double r = std::abs(ScintPoint[1] - fcathode_centre.Y());
    double d_c = std::abs(ScintPoint[0] - plane_depth); // distance to cathode
    double border_correction = 0;
    // flat PDs
    if ((opDet.type == 0 || opDet.type == 2) && fIsFlatPDCorr) {
      // interpolate in d_c for each r bin
      const size_t nbins_r = fvispars_flat[k].size();
      std::vector<double> interp_vals(nbins_r, 0.0);
      {
        size_t idx = 0;
        size_t size = fvis_distances_x_flat.size();
        if (d_c >= fvis_distances_x_flat[size - 2])
          idx = size - 2;
        else {
          while (d_c > fvis_distances_x_flat[idx + 1])
            idx++;
        }
        for (size_t i = 0; i < nbins_r; ++i) {
          interp_vals[i] = interpolate(fvis_distances_x_flat, fvispars_flat[k][i], d_c, false, idx);
        }
      }
      // interpolate in r
      border_correction = interpolate(fvis_distances_r_flat, interp_vals, r, false);
    }
    // dome PDs
    else if (opDet.type == 1 && fIsDomePDCorr) {
      // interpolate in d_c for each r bin
      const size_t nbins_r = fvispars_dome[k].size();
      std::vector<double> interp_vals(nbins_r, 0.0);
      {
        size_t idx = 0;
        size_t size = fvis_distances_x_dome.size();
        if (d_c >= fvis_distances_x_dome[size - 2])
          idx = size - 2;
        else {
          while (d_c > fvis_distances_x_dome[idx + 1])
            idx++;
        }
        for (size_t i = 0; i < nbins_r; ++i) {
          interp_vals[i] = interpolate(fvis_distances_x_dome, fvispars_dome[k][i], d_c, false, idx);
        }
      }
      // interpolate in r
      border_correction = interpolate(fvis_distances_r_dome, interp_vals, r, false);
    }
    else {
      std::cout << "Error: Invalid optical detector type. 0 = rectangular, 1 = dome, 2 = disk. Or "
                   "corrections for chosen optical detector type missing."
                << std::endl;
    }

    // apply correction factor
    double hits_rec = border_correction * hits_geo / cosine_vis;
    double hits_vis = std::round(G4Poisson(hits_rec));

    return hits_vis;
  }

  bool OpFastScintillation::isOpDetInSameTPC(geo::Point_t const& ScintPoint,
                                             geo::Point_t const& OpDetPoint) const
  {
    // check optical channel is in same TPC as scintillation light, if not return 0 hits;
    // temporary method working for SBND, uBooNE, DUNE 1x2x6; to be replaced to work in
    // full DUNE geometry; check x coordinate has same sign or is close to zero, otherwise
    // return 0 hits
    if (((ScintPoint.X() < 0.) != (OpDetPoint.X() < 0.)) &&
        std::abs(OpDetPoint.X()) > 10.) { // TODO: unhardcode
      return false;
    }
    return true;
  }

  bool OpFastScintillation::isScintInActiveVolume(geo::Point_t const& ScintPoint)
  {
    //semi-analytic approach only works in the active volume
    return fActiveVolumes[0].ContainsPosition(ScintPoint);
  }

  G4double OpFastScintillation::single_exp(const G4double t, const G4double tau2) const
  {
    return std::exp(-1.0 * t / tau2) / tau2;
  }

  G4double OpFastScintillation::bi_exp(const G4double t,
                                       const G4double tau1,
                                       const G4double tau2) const
  { // TODO: what's up with this? ... / tau2 / tau2 ...
    return std::exp(-1.0 * t / tau2) * (1 - std::exp(-1.0 * t / tau1)) / tau2 / tau2 *
           (tau1 + tau2);
  }

  G4double OpFastScintillation::Gaisser_Hillas(const double x, const double* par) const
  {
    double X_mu_0 = par[3];
    double Normalization = par[0];
    double Diff = par[1] - X_mu_0;
    double Term = std::pow((x - X_mu_0) / Diff, Diff / par[2]);
    double Exponential = std::exp((par[1] - x) / par[2]);
    return (Normalization * Term * Exponential);
  }
  //======================================================================
  // Returns interpolated value at x from parallel arrays ( xData, yData ) Assumes that
  // xData has at least two elements, is sorted and is strictly monotonic increasing
  // boolean argument extrapolate determines behaviour beyond ends of array (if needed)
  double OpFastScintillation::interpolate(const std::vector<double>& xData,
                                          const std::vector<double>& yData,
                                          double x,
                                          bool extrapolate,
                                          size_t i) const
  {
    if (i == 0) {
      size_t size = xData.size();
      if (x >= xData[size - 2]) { // special case: beyond right end
        i = size - 2;
      }
      else {
        while (x > xData[i + 1])
          i++;
      }
    }
    double xL = xData[i];
    double xR = xData[i + 1];
    double yL = yData[i];
    double yR = yData[i + 1]; // points on either side (unless beyond ends)
    if (!extrapolate) {       // if beyond ends of array and not extrapolating
      if (x < xL) return yL;
      if (x > xR) return yL;
    }
    const double dydx = (yR - yL) / (xR - xL); // gradient
    return yL + dydx * (x - xL);               // linear interpolation
  }

  void OpFastScintillation::interpolate3(std::array<double, 3>& inter,
                                         const std::vector<double>& xData,
                                         const std::vector<double>& yData1,
                                         const std::vector<double>& yData2,
                                         const std::vector<double>& yData3,
                                         double x,
                                         bool extrapolate) const
  {
    size_t size = xData.size();
    size_t i = 0;               // find left end of interval for interpolation
    if (x >= xData[size - 2]) { // special case: beyond right end
      i = size - 2;
    }
    else {
      while (x > xData[i + 1])
        i++;
    }
    double xL = xData[i];
    double xR = xData[i + 1]; // points on either side (unless beyond ends)
    double yL1 = yData1[i];
    double yR1 = yData1[i + 1];
    double yL2 = yData2[i];
    double yR2 = yData2[i + 1];
    double yL3 = yData3[i];
    double yR3 = yData3[i + 1];

    if (!extrapolate) { // if beyond ends of array and not extrapolating
      if (x < xL) {
        inter[0] = yL1;
        inter[1] = yL2;
        inter[2] = yL3;
        return;
      }
      if (x > xR) {
        inter[0] = yL1;
        inter[1] = yL2;
        inter[2] = yL3;
        return;
      }
    }
    const double m = (x - xL) / (xR - xL);
    inter[0] = m * (yR1 - yL1) + yL1;
    inter[1] = m * (yR2 - yL2) + yL2;
    inter[2] = m * (yR3 - yL3) + yL3;
  }

  // solid angle of circular aperture
  // TODO: allow greater tolerance in comparisons, by default it's using:
  // std::numeric_limits<double>::epsilon(): 2.22045e-16 that's an unrealistic small
  // number, better setting constexpr double tolerance = 0.0000001; // 1 nm
  double OpFastScintillation::Disk_SolidAngle(const double d, const double h, const double b) const
  {
    if (b <= 0. || d < 0. || h <= 0.) return 0.;
    const double leg2 = (b + d) * (b + d);
    const double aa = std::sqrt(h * h / (h * h + leg2));
    if (isApproximatelyZero(d)) { return 2. * CLHEP::pi * (1. - aa); }
    double bb = 2. * std::sqrt(b * d / (h * h + leg2));
    double cc = 4. * b * d / leg2;

    if (isDefinitelyGreaterThan(d, b)) {
      try {
        return 2. * aa *
               (std::sqrt(1. - cc) * boost::math::ellint_3(bb, cc, noLDoublePromote()) -
                boost::math::ellint_1(bb, noLDoublePromote()));
      }
      catch (std::domain_error& e) {
        if (isApproximatelyEqual(d, b, 1e-9)) {
          mf::LogWarning("OpFastScintillation")
            << "Elliptic Integral in Disk_SolidAngle() given parameters outside domain."
            << "\nbb: " << bb << "\ncc: " << cc << "\nException message: " << e.what()
            << "\nRelax condition and carry on.";
          return CLHEP::pi - 2. * aa * boost::math::ellint_1(bb, noLDoublePromote());
        }
        else {
          mf::LogError("OpFastScintillation")
            << "Elliptic Integral inside Disk_SolidAngle() given parameters outside domain.\n"
            << "\nbb: " << bb << "\ncc: " << cc << "Exception message: " << e.what();
          return 0.;
        }
      }
    }
    if (isDefinitelyLessThan(d, b)) {
      try {
        return 2. * CLHEP::pi -
               2. * aa *
                 (boost::math::ellint_1(bb, noLDoublePromote()) +
                  std::sqrt(1. - cc) * boost::math::ellint_3(bb, cc, noLDoublePromote()));
      }
      catch (std::domain_error& e) {
        if (isApproximatelyEqual(d, b, 1e-9)) {
          mf::LogWarning("OpFastScintillation")
            << "Elliptic Integral in Disk_SolidAngle() given parameters outside domain."
            << "\nbb: " << bb << "\ncc: " << cc << "\nException message: " << e.what()
            << "\nRelax condition and carry on.";
          return CLHEP::pi - 2. * aa * boost::math::ellint_1(bb, noLDoublePromote());
        }
        else {
          mf::LogError("OpFastScintillation")
            << "Elliptic Integral inside Disk_SolidAngle() given parameters outside domain.\n"
            << "\nbb: " << bb << "\ncc: " << cc << "Exception message: " << e.what();
          return 0.;
        }
      }
    }
    if (isApproximatelyEqual(d, b)) {
      return CLHEP::pi - 2. * aa * boost::math::ellint_1(bb, noLDoublePromote());
    }
    return 0.;
  }

  // solid angle of rectangular aperture
  double OpFastScintillation::Rectangle_SolidAngle(const double a,
                                                   const double b,
                                                   const double d) const
  {
    double aa = a / (2. * d);
    double bb = b / (2. * d);
    double aux = (1. + aa * aa + bb * bb) / ((1. + aa * aa) * (1. + bb * bb));
    return 4. * fast_acos(std::sqrt(aux));
  }

  // TODO: allow greater tolerance in comparisons, see note above on Disk_SolidAngle()
  double OpFastScintillation::Rectangle_SolidAngle(Dims const& o,
                                                   const std::array<double, 3> v) const
  {
    // v is the position of the track segment with respect to the center position of the
    // arapuca window

    // arapuca plane fixed in x direction
    if (isApproximatelyZero(v[1]) && isApproximatelyZero(v[2])) {
      return Rectangle_SolidAngle(o.h, o.w, v[0]);
    }
    if (isDefinitelyGreaterThan(v[1], o.h * .5) && isDefinitelyGreaterThan(v[2], o.w * .5)) {
      double A = v[1] - o.h * .5;
      double B = v[2] - o.w * .5;
      double to_return = (Rectangle_SolidAngle(2. * (A + o.h), 2. * (B + o.w), v[0]) -
                          Rectangle_SolidAngle(2. * A, 2. * (B + o.w), v[0]) -
                          Rectangle_SolidAngle(2. * (A + o.h), 2. * B, v[0]) +
                          Rectangle_SolidAngle(2. * A, 2. * B, v[0])) *
                         .25;
      return to_return;
    }
    if ((v[1] <= o.h * .5) && (v[2] <= o.w * .5)) {
      double A = -v[1] + o.h * .5;
      double B = -v[2] + o.w * .5;
      double to_return = (Rectangle_SolidAngle(2. * (o.h - A), 2. * (o.w - B), v[0]) +
                          Rectangle_SolidAngle(2. * A, 2. * (o.w - B), v[0]) +
                          Rectangle_SolidAngle(2. * (o.h - A), 2. * B, v[0]) +
                          Rectangle_SolidAngle(2. * A, 2. * B, v[0])) *
                         .25;
      return to_return;
    }
    if (isDefinitelyGreaterThan(v[1], o.h * .5) && (v[2] <= o.w * .5)) {
      double A = v[1] - o.h * .5;
      double B = -v[2] + o.w * .5;
      double to_return = (Rectangle_SolidAngle(2. * (A + o.h), 2. * (o.w - B), v[0]) -
                          Rectangle_SolidAngle(2. * A, 2. * (o.w - B), v[0]) +
                          Rectangle_SolidAngle(2. * (A + o.h), 2. * B, v[0]) -
                          Rectangle_SolidAngle(2. * A, 2. * B, v[0])) *
                         .25;
      return to_return;
    }
    if ((v[1] <= o.h * .5) && isDefinitelyGreaterThan(v[2], o.w * .5)) {
      double A = -v[1] + o.h * .5;
      double B = v[2] - o.w * .5;
      double to_return = (Rectangle_SolidAngle(2. * (o.h - A), 2. * (B + o.w), v[0]) -
                          Rectangle_SolidAngle(2. * (o.h - A), 2. * B, v[0]) +
                          Rectangle_SolidAngle(2. * A, 2. * (B + o.w), v[0]) -
                          Rectangle_SolidAngle(2. * A, 2. * B, v[0])) *
                         .25;
      return to_return;
    }
    // error message if none of these cases, i.e. something has gone wrong!
    return 0.;
  }

  double OpFastScintillation::Omega_Dome_Model(const double distance, const double theta) const
  {
    // this function calculates the solid angle of a semi-sphere of radius b, as a
    // correction to the analytic formula of the on-axix solid angle, as we move off-axis
    // an angle theta. We have used 9-angular bins with delta_theta width.

    // par0 = Radius correction close
    // par1 = Radius correction far
    // par2 = breaking distance betwween "close" and "far"

    double par0[9] = {0., 0., 0., 0., 0., 0.597542, 1.00872, 1.46993, 2.04221};
    double par1[9] = {0, 0, 0.19569, 0.300449, 0.555598, 0.854939, 1.39166, 2.19141, 2.57732};
    const double delta_theta = 10.;
    int j = int(theta / delta_theta);
    // PMT radius
    const double b = fradius; // cm
    // distance form which the model parameters break (empirical value)
    const double d_break = 5 * b; //par2

    if (distance >= d_break) {
      double R_apparent_far = b - par1[j];
      return (2 * CLHEP::pi * (1 - std::sqrt(1 - std::pow(R_apparent_far / distance, 2))));
    }
    else {
      double R_apparent_close = b - par0[j];
      return (2 * CLHEP::pi * (1 - std::sqrt(1 - std::pow(R_apparent_close / distance, 2))));
    }
  }

  // ---------------------------------------------------------------------------
  std::vector<geo::BoxBoundedGeo> OpFastScintillation::extractActiveVolumes(
    geo::GeometryCore const& geom)
  {
    std::vector<geo::BoxBoundedGeo> activeVolumes;
    activeVolumes.reserve(geom.Ncryostats());

    for (geo::CryostatGeo const& cryo : geom.Iterate<geo::CryostatGeo>()) {

      // can't use it default-constructed since it would always include origin
      geo::BoxBoundedGeo box{cryo.TPC(0).ActiveBoundingBox()};

      for (geo::TPCGeo const& TPC : cryo.IterateTPCs())
        box.ExtendToInclude(TPC.ActiveBoundingBox());

      activeVolumes.push_back(std::move(box));

    } // for cryostats

    return activeVolumes;
  }

  // ---------------------------------------------------------------------------

}