OpFastScintillation.hh

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.

#ifndef OpFastScintillation_h
#define OpFastScintillation_h 1

#include "larcorealg/Geometry/BoxBoundedGeo.h"
#include "larcorealg/Geometry/fwd.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h" // geo::Point_t
#include "larsim/PhotonPropagation/PhotonVisibilityTypes.h" // phot::MappedT0s_t

#include "Geant4/G4ForceCondition.hh"
#include "Geant4/G4ParticleDefinition.hh"
#include "Geant4/G4PhysicsOrderedFreeVector.hh"
#include "Geant4/G4PhysicsTable.hh"
#include "Geant4/G4ProcessType.hh"
#include "Geant4/G4String.hh"
#include "Geant4/G4ThreeVector.hh"
#include "Geant4/G4Types.hh"
#include "Geant4/G4VRestDiscreteProcess.hh"

#include "TF1.h"
#include "TVector3.h"

#include <memory> // std::unique_ptr

class G4EmSaturation;
class G4Step;
class G4Track;
class G4VParticleChange;
namespace CLHEP {
  class RandGeneral;
}
namespace phot {
  class PhotonVisibilityService;
}

namespace larg4 {

  class OpFastScintillation : public G4VRestDiscreteProcess {
  public:
    OpFastScintillation(const G4String& processName = "Scintillation",
                        G4ProcessType type = fElectromagnetic);

    ~OpFastScintillation();

    // OpFastScintillation Process has both PostStepDoIt (for energy deposition of
    // particles in flight) and AtRestDoIt (for energy given to the medium by particles at
    // rest)

    virtual G4bool IsApplicable(const G4ParticleDefinition& aParticleType);
    // Returns true -> 'is applicable', for any particle type except for an
    // 'opticalphoton' and for short-lived particles

    G4double GetMeanFreePath(const G4Track& aTrack, G4double, G4ForceCondition*);
    // Returns infinity; i. e. the process does not limit the step, but sets the
    // 'StronglyForced' condition for the DoIt to be invoked at every step.

    G4double GetMeanLifeTime(const G4Track& aTrack, G4ForceCondition*);
    // Returns infinity; i. e. the process does not limit the time, but sets the
    // 'StronglyForced' condition for the DoIt to be invoked at every step.

    virtual G4VParticleChange* PostStepDoIt(const G4Track& aTrack, const G4Step& aStep);
    virtual G4VParticleChange* AtRestDoIt(const G4Track& aTrack, const G4Step& aStep);

    // These are the methods implementing the scintillation process.

    void SetTrackSecondariesFirst(const G4bool state);
    // If set, the primary particle tracking is interrupted and any produced scintillation
    // photons are tracked next. When all have been tracked, the tracking of the primary
    // resumes.

    void SetFiniteRiseTime(const G4bool state);
    // If set, the OpFastScintillation process expects the user to have set the constant
    // material property FAST/SLOWSCINTILLATIONRISETIME.

    G4bool GetTrackSecondariesFirst() const;
    // Returns the boolean flag for tracking secondaries first.

    G4bool GetFiniteRiseTime() const;
    // Returns the boolean flag for a finite scintillation rise time.

    void SetScintillationYieldFactor(const G4double yieldfactor);
    // Called to set the scintillation photon yield factor, needed when the yield is
    // different for different types of particles. This scales the yield obtained from the
    // G4MaterialPropertiesTable.

    G4double GetScintillationYieldFactor() const;
    // Returns the photon yield factor.

    void SetScintillationExcitationRatio(const G4double excitationratio);
    // Called to set the scintillation exciation ratio, needed when the scintillation
    // level excitation is different for different types of particles. This overwrites the
    // YieldRatio obtained from the G4MaterialPropertiesTable.

    G4double GetScintillationExcitationRatio() const;
    // Returns the scintillation level excitation ratio.

    G4PhysicsTable* GetFastIntegralTable() const;
    // Returns the address of the fast scintillation integral table.

    G4PhysicsTable* GetSlowIntegralTable() const;
    // Returns the address of the slow scintillation integral table.

    void AddSaturation(G4EmSaturation* sat) { emSaturation = sat; }
    // Adds Birks Saturation to the process.

    void RemoveSaturation() { emSaturation = NULL; }
    // Removes the Birks Saturation from the process.

    G4EmSaturation* GetSaturation() const { return emSaturation; }
    // Returns the Birks Saturation.

    void SetScintillationByParticleType(const G4bool);
    // Called by the user to set the scintillation yield as a function of energy deposited
    // by particle type

    G4bool GetScintillationByParticleType() const { return scintillationByParticleType; }
    // Return the boolean that determines the method of scintillation production

    void DumpPhysicsTable() const;
    // Prints the fast and slow scintillation integral tables.

    void getVUVTimes(std::vector<double>& arrivalTimes,
                     const double distance_in_cm,
                     const size_t angle_bin);
    void generateParam(const size_t index, const size_t angle_bin);
    // Functions for vuv component Landau + Exponential timing parameterisation, updated
    // method

    void getVISTimes(std::vector<double>& arrivalTimes,
                     const TVector3& ScintPoint,
                     const TVector3& OpDetPoint);
    // Visible component timing parameterisation

    void detectedDirectHits(std::map<size_t, int>& DetectedNum,
                            const double Num,
                            geo::Point_t const& ScintPoint) const;
    void detectedReflecHits(std::map<size_t, int>& ReflDetectedNum,
                            const double Num,
                            geo::Point_t const& ScintPoint) const;

  protected:
    void BuildThePhysicsTable();
    // It builds either the fast or slow scintillation integral table; or both.

    bool RecordPhotonsProduced(const G4Step& aStep, double N);
    // Note the production of N photons in at point xyz.  pass on to generate detector
    // response, etc.

    std::unique_ptr<G4PhysicsTable> theSlowIntegralTable;
    std::unique_ptr<G4PhysicsTable> theFastIntegralTable;

    G4bool fTrackSecondariesFirst;
    G4bool fFiniteRiseTime;

    G4double YieldFactor;

    G4double ExcitationRatio;

    G4bool scintillationByParticleType;

  private:
    struct OpticalDetector {
      double h; // height
      double w; // width
      geo::Point_t OpDetPoint;
      int type;
    };

    bool usesSemiAnalyticModel() const;

    int VUVHits(const double Nphotons_created,
                geo::Point_t const& ScintPoint,
                OpticalDetector const& opDet) const;
    // Calculates semi-analytic model number of hits for vuv component

    int VISHits(geo::Point_t const& ScintPoint,
                OpticalDetector const& opDet,
                const double cathode_hits_rec,
                const std::array<double, 3> hotspot) const;
    // Calculates semi-analytic model number of hits for visible component

    G4double single_exp(const G4double t, const G4double tau2) const;
    G4double bi_exp(const G4double t, const G4double tau1, const G4double tau2) const;

    G4double scint_time(const G4Step& aStep,
                        G4double ScintillationTime,
                        G4double ScintillationRiseTime) const;
    void propagationTime(std::vector<double>& arrival_time_dist,
                         G4ThreeVector x0,
                         const size_t OpChannel,
                         bool Reflected = false); //const;

    // emission time distribution when there is a finite rise time
    G4double sample_time(const G4double tau1, const G4double tau2) const;

    // Facility for TPB emission energies
    double reemission_energy() const;
    std::map<double, double> tpbemission;
    std::unique_ptr<CLHEP::RandGeneral> fTPBEm;

    void average_position(G4Step const& aStep, double* xzyPos) const;

    G4EmSaturation* emSaturation;
    // functions and parameters for the propagation time parametrization
    phot::MappedFunctions_t ParPropTimeTF1;
    phot::MappedT0s_t ReflT0s;

    //For new VUV time parametrization
    double fstep_size, fmin_d, fmax_d, fvuv_vgroup_mean, fvuv_vgroup_max, finflexion_point_distance,
      fangle_bin_timing_vuv;
    std::vector<std::vector<double>> fparameters[7];
    // vector containing generated VUV timing parameterisations
    std::vector<std::vector<TF1>> VUV_timing;
    // vector containing min and max range VUV timing parameterisations are sampled to
    std::vector<std::vector<double>> VUV_max;
    std::vector<std::vector<double>> VUV_min;

    // For new VIS time parameterisation
    double fvis_vmean, fangle_bin_timing_vis;
    std::vector<double> fdistances_refl;
    std::vector<double> fradial_distances_refl;
    std::vector<std::vector<std::vector<double>>> fcut_off_pars;
    std::vector<std::vector<std::vector<double>>> ftau_pars;

    struct Dims {
      double h, w; // height, width
    };

    // solid angle of rectangular aperture calculation functions
    double Rectangle_SolidAngle(const double a, const double b, const double d) const;
    double Rectangle_SolidAngle(Dims const& o, const std::array<double, 3> v) const;
    // solid angle of circular aperture calculation functions
    double Disk_SolidAngle(const double d, const double h, const double b) const;
    // solid angle of a dome aperture calculation functions
    double Omega_Dome_Model(const double distance, const double theta) const;

    // For VUV semi-analytic hits
    // Gaisser-Hillas correction parameters for VUV Nhits estimation
    G4double Gaisser_Hillas(const double x, const double* par) const;
    double fdelta_angulo_vuv;
    // flat PDs
    bool fIsFlatPDCorr;
    std::vector<std::vector<double>> fGHvuvpars_flat;
    std::vector<double> fborder_corr_angulo_flat;
    std::vector<std::vector<double>> fborder_corr_flat;
    // dome PDs
    bool fIsDomePDCorr;
    std::vector<std::vector<double>> fGHvuvpars_dome;
    std::vector<double> fborder_corr_angulo_dome;
    std::vector<std::vector<double>> fborder_corr_dome;

    // For VIS semi-analytic hits
    bool fStoreReflected;
    // correction parameters for VIS Nhits estimation
    double fdelta_angulo_vis;
    // flat PDs
    std::vector<double> fvis_distances_x_flat;
    std::vector<double> fvis_distances_r_flat;
    std::vector<std::vector<std::vector<double>>> fvispars_flat;
    // dome PDs
    std::vector<double> fvis_distances_x_dome;
    std::vector<double> fvis_distances_r_dome;
    std::vector<std::vector<std::vector<double>>> fvispars_dome;

    // geometry properties
    double fplane_depth, fcathode_zdimension, fcathode_ydimension;
    geo::Point_t fcathode_centre;
    std::vector<geo::BoxBoundedGeo> const fActiveVolumes;

    // Optical detector properties for semi-analytic hits
    double fradius;
    Dims fcathode_plane;
    int fL_abs_vuv;
    std::vector<geo::Point_t> fOpDetCenter;
    std::vector<int> fOpDetType;
    std::vector<double> fOpDetLength;
    std::vector<double> fOpDetHeight;
    //double fGlobalTimeOffset;

    void ProcessStep(const G4Step& step);

    bool const bPropagate; 

    phot::PhotonVisibilityService const* const fPVS;

    bool const fUseNhitsModel = false;
    bool const fOnlyActiveVolume = false;
    bool const fOnlyOneCryostat = false;
    bool const fOpaqueCathode = false;

    bool isOpDetInSameTPC(geo::Point_t const& ScintPoint, geo::Point_t const& OpDetPoint) const;
    bool isScintInActiveVolume(geo::Point_t const& ScintPoint);
    double interpolate(const std::vector<double>& xData,
                       const std::vector<double>& yData,
                       double x,
                       bool extrapolate,
                       size_t i = 0) const;
    void 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;

    static std::vector<geo::BoxBoundedGeo> extractActiveVolumes(geo::GeometryCore const& geom);

  }; // class OpFastScintillation

  inline G4bool OpFastScintillation::IsApplicable(const G4ParticleDefinition& aParticleType)
  {
    if (aParticleType.GetParticleName() == "opticalphoton") return false;
    if (aParticleType.IsShortLived()) return false;

    return true;
  }

  inline void OpFastScintillation::SetTrackSecondariesFirst(const G4bool state)
  {
    fTrackSecondariesFirst = state;
  }

  inline void OpFastScintillation::SetFiniteRiseTime(const G4bool state)
  {
    fFiniteRiseTime = state;
  }

  inline G4bool OpFastScintillation::GetTrackSecondariesFirst() const
  {
    return fTrackSecondariesFirst;
  }

  inline G4bool OpFastScintillation::GetFiniteRiseTime() const
  {
    return fFiniteRiseTime;
  }

  inline void OpFastScintillation::SetScintillationYieldFactor(const G4double yieldfactor)
  {
    YieldFactor = yieldfactor;
  }

  inline G4double OpFastScintillation::GetScintillationYieldFactor() const
  {
    return YieldFactor;
  }

  inline void OpFastScintillation::SetScintillationExcitationRatio(const G4double excitationratio)
  {
    ExcitationRatio = excitationratio;
  }

  inline G4double OpFastScintillation::GetScintillationExcitationRatio() const
  {
    return ExcitationRatio;
  }

  inline G4PhysicsTable* OpFastScintillation::GetSlowIntegralTable() const
  {
    return theSlowIntegralTable.get();
  }

  inline G4PhysicsTable* OpFastScintillation::GetFastIntegralTable() const
  {
    return theFastIntegralTable.get();
  }

  inline void OpFastScintillation::DumpPhysicsTable() const
  {
    if (theFastIntegralTable) {
      G4int PhysicsTableSize = theFastIntegralTable->entries();
      G4PhysicsOrderedFreeVector* v;
      for (G4int i = 0; i < PhysicsTableSize; i++) {
        v = (G4PhysicsOrderedFreeVector*)(*theFastIntegralTable)[i];
        v->DumpValues();
      }
    }
    if (theSlowIntegralTable) {
      G4int PhysicsTableSize = theSlowIntegralTable->entries();
      G4PhysicsOrderedFreeVector* v;
      for (G4int i = 0; i < PhysicsTableSize; i++) {
        v = (G4PhysicsOrderedFreeVector*)(*theSlowIntegralTable)[i];
        v->DumpValues();
      }
    }
  }

  template <typename TReal>
  inline double dist(const TReal* x, const TReal* y, const unsigned int dimension)
  {
    double d = 0.;
    for (unsigned int p = 0; p < dimension; ++p) {
      d += (*(x + p) - *(y + p)) * (*(x + p) - *(y + p));
    }
    return std::sqrt(d);
  }

  template <typename TVector3>
  inline double dist(const std::array<double, 3> x,
                     const TVector3 y,
                     const unsigned int dimension,
                     const unsigned int start)
  {
    double d = 0.;
    for (unsigned int p = start; p < dimension; ++p) {
      d += (x[p] - y[p]) * (x[p] - y[p]);
    }
    return std::sqrt(d);
  }

  // implements relative method - do not use for comparing with zero
  // use this most of the time, tolerance needs to be meaningful in your context
  template <typename TReal>
  inline constexpr static bool
  isApproximatelyEqual(TReal a, TReal b, TReal tolerance = std::numeric_limits<TReal>::epsilon())
  {
    TReal diff = std::fabs(a - b);
    if (diff <= tolerance) return true;
    if (diff < std::fmax(std::fabs(a), std::fabs(b)) * tolerance) return true;
    return false;
  }

  // supply tolerance that is meaningful in your context
  // for example, default tolerance may not work if you are comparing double with float
  template <typename TReal>
  inline constexpr static bool isApproximatelyZero(
    TReal a,
    TReal tolerance = std::numeric_limits<TReal>::epsilon())
  {
    if (std::fabs(a) <= tolerance) return true;
    return false;
  }

  // use this when you want to be on safe side
  // for example, don't start rover unless signal is above 1
  template <typename TReal>
  inline constexpr static bool
  isDefinitelyLessThan(TReal a, TReal b, TReal tolerance = std::numeric_limits<TReal>::epsilon())
  {
    TReal diff = a - b;
    if (diff < tolerance) return true;
    if (diff < std::fmax(std::fabs(a), std::fabs(b)) * tolerance) return true;
    return false;
  }

  template <typename TReal>
  inline constexpr static bool
  isDefinitelyGreaterThan(TReal a, TReal b, TReal tolerance = std::numeric_limits<TReal>::epsilon())
  {
    TReal diff = a - b;
    if (diff > tolerance) return true;
    if (diff > std::fmax(std::fabs(a), std::fabs(b)) * tolerance) return true;
    return false;
  }

  // template<typename Function, typename... Args>
  // auto OpFastScintillation::invoke_memoized(Function function, Args... args)
  // {
  //   using key_type   = std::tuple<Args...>;
  //   using value_type = std::invoke_result_t<Function, Args...>;
  //   static_assert(! std::is_same_v<Function, std::function<value_type(Args...)>>,
  //                 "cannot memoize on std::function (use a lambda instead)");
  //   static_assert(! std::is_same_v<Function, value_type(*)(Args...)>,
  //                 "cannot memoize on function pointer (use a lambda instead)");
  //   static std::mutex mutex;
  //   static std::map<key_type, value_type> cache;
  //   auto key  = std::tuple(args...);
  //   auto lock = std::lock_guard<std::mutex>(mutex);
  //   if (cache.count(key)){
  //     return cache[key];
  //   }
  //   return cache[key] = std::apply(function, key);
  // }

} // namespace larg4

#endif /* OpFastScintillation_h */