larg4::OpFastScintillation Class Reference

#include <OpFastScintillation.hh>

Inheritance diagram for larg4::OpFastScintillation:

Classes
struct	Dims
struct	OpticalDetector

Public Member Functions
	OpFastScintillation (const G4String &processName="Scintillation", G4ProcessType type=fElectromagnetic)
	~OpFastScintillation ()
virtual G4bool	IsApplicable (const G4ParticleDefinition &aParticleType)
G4double	GetMeanFreePath (const G4Track &aTrack, G4double, G4ForceCondition *)
G4double	GetMeanLifeTime (const G4Track &aTrack, G4ForceCondition *)
virtual G4VParticleChange *	PostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
virtual G4VParticleChange *	AtRestDoIt (const G4Track &aTrack, const G4Step &aStep)
void	SetTrackSecondariesFirst (const G4bool state)
void	SetFiniteRiseTime (const G4bool state)
G4bool	GetTrackSecondariesFirst () const
G4bool	GetFiniteRiseTime () const
void	SetScintillationYieldFactor (const G4double yieldfactor)
G4double	GetScintillationYieldFactor () const
void	SetScintillationExcitationRatio (const G4double excitationratio)
G4double	GetScintillationExcitationRatio () const
G4PhysicsTable *	GetFastIntegralTable () const
G4PhysicsTable *	GetSlowIntegralTable () const
void	AddSaturation (G4EmSaturation *sat)
void	RemoveSaturation ()
G4EmSaturation *	GetSaturation () const
void	SetScintillationByParticleType (const G4bool)
G4bool	GetScintillationByParticleType () const
void	DumpPhysicsTable () const
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)
void	getVISTimes (std::vector< double > &arrivalTimes, const TVector3 &ScintPoint, const TVector3 &OpDetPoint)
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 Member Functions
void	BuildThePhysicsTable ()
bool	RecordPhotonsProduced (const G4Step &aStep, double N)

Protected Attributes
std::unique_ptr< G4PhysicsTable >	theSlowIntegralTable
std::unique_ptr< G4PhysicsTable >	theFastIntegralTable
G4bool	fTrackSecondariesFirst
G4bool	fFiniteRiseTime
G4double	YieldFactor
G4double	ExcitationRatio
G4bool	scintillationByParticleType

Private Member Functions
bool	usesSemiAnalyticModel () const
	Returns whether the semi-analytic visibility parametrization is being used. More...
int	VUVHits (const double Nphotons_created, geo::Point_t const &ScintPoint, OpticalDetector const &opDet) const
int	VISHits (geo::Point_t const &ScintPoint, OpticalDetector const &opDet, const double cathode_hits_rec, const std::array< double, 3 > hotspot) const
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)
G4double	sample_time (const G4double tau1, const G4double tau2) const
double	reemission_energy () const
void	average_position (G4Step const &aStep, double *xzyPos) const
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
double	Disk_SolidAngle (const double d, const double h, const double b) const
double	Omega_Dome_Model (const double distance, const double theta) const
G4double	Gaisser_Hillas (const double x, const double *par) const
void	ProcessStep (const G4Step &step)
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 Private Member Functions
static std::vector< geo::BoxBoundedGeo >	extractActiveVolumes (geo::GeometryCore const &geom)

Private Attributes
std::map< double, double >	tpbemission
std::unique_ptr< CLHEP::RandGeneral >	fTPBEm
G4EmSaturation *	emSaturation
phot::MappedFunctions_t	ParPropTimeTF1
phot::MappedT0s_t	ReflT0s
double	fstep_size
double	fmin_d
double	fmax_d
double	fvuv_vgroup_mean
double	fvuv_vgroup_max
double	finflexion_point_distance
double	fangle_bin_timing_vuv
std::vector< std::vector< double > >	fparameters [7]
std::vector< std::vector< TF1 > >	VUV_timing
std::vector< std::vector< double > >	VUV_max
std::vector< std::vector< double > >	VUV_min
double	fvis_vmean
double	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
double	fdelta_angulo_vuv
bool	fIsFlatPDCorr
std::vector< std::vector< double > >	fGHvuvpars_flat
std::vector< double >	fborder_corr_angulo_flat
std::vector< std::vector< double > >	fborder_corr_flat
bool	fIsDomePDCorr
std::vector< std::vector< double > >	fGHvuvpars_dome
std::vector< double >	fborder_corr_angulo_dome
std::vector< std::vector< double > >	fborder_corr_dome
bool	fStoreReflected
double	fdelta_angulo_vis
std::vector< double >	fvis_distances_x_flat
std::vector< double >	fvis_distances_r_flat
std::vector< std::vector< std::vector< double > > >	fvispars_flat
std::vector< double >	fvis_distances_x_dome
std::vector< double >	fvis_distances_r_dome
std::vector< std::vector< std::vector< double > > >	fvispars_dome
double	fplane_depth
double	fcathode_zdimension
double	fcathode_ydimension
TVector3	fcathode_centre
std::vector< geo::BoxBoundedGeo > const	fActiveVolumes
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
bool const	bPropagate
	Whether propagation of photons is enabled. More...
phot::PhotonVisibilityService const *const	fPVS
	Photon visibility service instance. More...
bool const	fUseNhitsModel = false
	Whether the semi-analytic model is being used for photon visibility. More...
bool const	fOnlyActiveVolume = false
	Whether photon propagation is performed only from active volumes. More...
bool const	fOnlyOneCryostat = false
bool const	fOpaqueCathode = false
	Whether the cathodes are fully opaque; currently hard coded "no". More...

Detailed Description

Definition at line 130 of file OpFastScintillation.hh.

Constructor & Destructor Documentation

larg4::OpFastScintillation::OpFastScintillation	(	const G4String &	processName = "Scintillation",
		G4ProcessType	type = fElectromagnetic
	)

Definition at line 180 of file OpFastScintillation.cxx.

    : 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] : util::enumerate(fActiveVolumes)) {
          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(0, 0).GetCathodeCenter().X(),
                                        fActiveVolumes[0].CenterY(),
                                        fActiveVolumes[0].CenterZ()};
      mf::LogTrace("OpFastScintillation") << "Cathode_centre: " << Cathode_centre << " cm";

      // std::cout << "\nInitialize acos_arr with " << acos_bins+1
      //           << " hence with a resolution of " << 1./acos_bins << std::endl;
      // for(size_t i=0; i<=acos_bins; ++i){
      //   acos_arr[i] = std::acos(i/double(acos_bins));
      // }

      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
          //    std::cout<<"Radio: "<<geom.OpDetGeoFromOpDet(i).RMax()<<std::endl;
          fOpDetLength.push_back(-1);
          fOpDetHeight.push_back(-1);
        }
      }

      if (fPVS->IncludePropTime()) {
        std::cout << "Using parameterisation of timings." << std::endl;
        //OLD VUV time parapetrization (to be removed soon)
        //fPVS->SetDirectLightPropFunctions(functions_vuv, fd_break, fd_max, ftf1_sampling_factor);
        //fPVS->SetReflectedCOLightPropFunctions(functions_vis, ft0_max, ft0_break_point);
        //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(0, 0).GetCathodeCenter();
        fcathode_centre[1] = fActiveVolumes[0].CenterY();
        fcathode_centre[2] = 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[0]);
        }
        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());
  }

larg4::OpFastScintillation::~OpFastScintillation

(

)

Definition at line 400 of file OpFastScintillation.cxx.

  {
    if (theFastIntegralTable) theFastIntegralTable->clearAndDestroy();
    if (theSlowIntegralTable) theSlowIntegralTable->clearAndDestroy();
  }

Member Function Documentation

void larg4::OpFastScintillation::AddSaturation ( G4EmSaturation *  sat )

inline

Definition at line 216 of file OpFastScintillation.hh.

    {
      emSaturation = sat;
    }

G4VParticleChange * larg4::OpFastScintillation::AtRestDoIt ( const G4Track &  aTrack, const G4Step &  aStep  )

virtual

Definition at line 414 of file OpFastScintillation.cxx.

  {
    return OpFastScintillation::PostStepDoIt(aTrack, aStep);
  }

void larg4::OpFastScintillation::average_position ( G4Step const &  aStep, double *  xzyPos  ) const

private

Definition at line 1102 of file OpFastScintillation.cxx.

  {
    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);
  }

G4double larg4::OpFastScintillation::bi_exp ( const G4double  t, const G4double  tau1, const G4double  tau2  ) const

private

Definition at line 1849 of file OpFastScintillation.cxx.

  { // 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);
  }

void larg4::OpFastScintillation::BuildThePhysicsTable ( )

protected

Definition at line 870 of file OpFastScintillation.cxx.

  {
    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);
    }
  }

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

Definition at line 1499 of file OpFastScintillation.cxx.

  {
    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;
        //   mf::LogInfo("OpFastScintillation") << "FastScint: " <<
        //   //   it->second<<" " << Num << " " << DetThisPMT;
        //det_photon_ctr += DetThisPMT; // CASE-DEBUG DO NOT REMOVE THIS COMMENT
      }
    }
  }

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

Definition at line 1521 of file OpFastScintillation.cxx.

  {
    // 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[1]),
                                                 std::abs(ScintPoint.Z() - fcathode_centre[2])};
    // 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[1], 2) + std::pow(ScintPoint.Z() - fcathode_centre[2], 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; }
    }
  }

double larg4::OpFastScintillation::Disk_SolidAngle ( const double  d, const double  h, const double  b  ) const

private

Definition at line 2013 of file OpFastScintillation.cxx.

  {
    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.;
  }

void larg4::OpFastScintillation::DumpPhysicsTable ( ) const

inline

Definition at line 541 of file OpFastScintillation.hh.

  {
    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();
      }
    }
  }

std::vector< geo::BoxBoundedGeo > larg4::OpFastScintillation::extractActiveVolumes ( geo::GeometryCore const &  geom )

staticprivate

Definition at line 2172 of file OpFastScintillation.cxx.

  {
    std::vector<geo::BoxBoundedGeo> activeVolumes;
    activeVolumes.reserve(geom.Ncryostats());

    for (geo::CryostatGeo const& cryo : geom.IterateCryostats()) {

      // 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;
  } // OpFastScintillation::extractActiveVolumes()

G4double larg4::OpFastScintillation::Gaisser_Hillas ( const double  x, const double *  par  ) const

private

Definition at line 1856 of file OpFastScintillation.cxx.

  {
    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);
  }

void larg4::OpFastScintillation::generateParam	(	const size_t	index,
		const size_t	angle_bin
	)

Definition at line 1252 of file OpFastScintillation.cxx.

  {
    // 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;
  }

G4PhysicsTable * larg4::OpFastScintillation::GetFastIntegralTable ( ) const

inline

Definition at line 535 of file OpFastScintillation.hh.

  {
    return theFastIntegralTable.get();
  }

G4bool larg4::OpFastScintillation::GetFiniteRiseTime ( ) const

inline

Definition at line 499 of file OpFastScintillation.hh.

  {
    return fFiniteRiseTime;
  }

G4double larg4::OpFastScintillation::GetMeanFreePath	(	const G4Track &	aTrack,
		G4double	,
		G4ForceCondition *	condition
	)

Definition at line 997 of file OpFastScintillation.cxx.

  {
    *condition = StronglyForced;
    return DBL_MAX;
  }

G4double larg4::OpFastScintillation::GetMeanLifeTime	(	const G4Track &	aTrack,
		G4ForceCondition *	condition
	)

Definition at line 1004 of file OpFastScintillation.cxx.

  {
    *condition = Forced;
    return DBL_MAX;
  }

G4EmSaturation* larg4::OpFastScintillation::GetSaturation ( ) const

inline

Definition at line 230 of file OpFastScintillation.hh.

    {
      return emSaturation;
    }

G4bool larg4::OpFastScintillation::GetScintillationByParticleType ( ) const

inline

Definition at line 241 of file OpFastScintillation.hh.

    {
      return scintillationByParticleType;
    }

G4double larg4::OpFastScintillation::GetScintillationExcitationRatio ( ) const

inline

Definition at line 523 of file OpFastScintillation.hh.

  {
    return ExcitationRatio;
  }

G4double larg4::OpFastScintillation::GetScintillationYieldFactor ( ) const

inline

Definition at line 511 of file OpFastScintillation.hh.

  {
    return YieldFactor;
  }

G4PhysicsTable * larg4::OpFastScintillation::GetSlowIntegralTable ( ) const

inline

Definition at line 529 of file OpFastScintillation.hh.

  {
    return theSlowIntegralTable.get();
  }

G4bool larg4::OpFastScintillation::GetTrackSecondariesFirst ( ) const

inline

Definition at line 493 of file OpFastScintillation.hh.

  {
    return fTrackSecondariesFirst;
  }

void larg4::OpFastScintillation::getVISTimes	(	std::vector< double > &	arrivalTimes,
		const TVector3 &	ScintPoint,
		const TVector3 &	OpDetPoint
	)

Definition at line 1374 of file OpFastScintillation.cxx.

  {
    // *************************************************************************************************
    //     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[1] - fcathode_centre[1], 2) + std::pow(ScintPoint[2] - fcathode_centre[2], 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;
    }
  }

void larg4::OpFastScintillation::getVUVTimes	(	std::vector< double > &	arrivalTimes,
		const double	distance_in_cm,
		const size_t	angle_bin
	)

Definition at line 1351 of file OpFastScintillation.cxx.

  {
    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]);
      }
    }
  }

double larg4::OpFastScintillation::interpolate ( const std::vector< double > &  xData, const std::vector< double > &  yData, double  x, bool  extrapolate, size_t  i = 0  ) const

private

Definition at line 1932 of file OpFastScintillation.cxx.

  {
    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 larg4::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

private

Definition at line 1961 of file OpFastScintillation.cxx.

  {
    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;
  }

G4bool larg4::OpFastScintillation::IsApplicable ( const G4ParticleDefinition &  aParticleType )

inlinevirtual

Definition at line 472 of file OpFastScintillation.hh.

  {
    if (aParticleType.GetParticleName() == "opticalphoton") return false;
    if (aParticleType.IsShortLived()) return false;

    return true;
  }

bool larg4::OpFastScintillation::isOpDetInSameTPC ( geo::Point_t const &  ScintPoint, geo::Point_t const &  OpDetPoint  ) const

private

Definition at line 1822 of file OpFastScintillation.cxx.

  {
    // 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 larg4::OpFastScintillation::isScintInActiveVolume ( geo::Point_t const &  ScintPoint )

private

Definition at line 1836 of file OpFastScintillation.cxx.

  {
    //semi-analytic approach only works in the active volume
    return fActiveVolumes[0].ContainsPosition(ScintPoint);
  }

double larg4::OpFastScintillation::Omega_Dome_Model ( const double  distance, const double  theta  ) const

private

Definition at line 2140 of file OpFastScintillation.cxx.

                                                                                              {
  // 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))));
  }
}

G4VParticleChange * larg4::OpFastScintillation::PostStepDoIt ( const G4Track &  aTrack, const G4Step &  aStep  )

virtual

Definition at line 425 of file OpFastScintillation.cxx.

  {
    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();

    ///////////////////////////////////////////////////////////////////////////////////
    //   This is the old G4 way - but we do things differently - Ben J, Oct Nov 2012.
    ///////////////////////////////////////////////////////////////////////////////////
    //
    //     if (MeanNumberOfPhotons > 10.)
    //      {
    //        G4double sigma = ResolutionScale * std::sqrt(MeanNumberOfPhotons);
    //        NumPhotons = G4int(G4RandGauss::shoot(MeanNumberOfPhotons,sigma)+0.5);
    //      }
    //     else
    //      {
    //        NumPhotons = G4int(G4Poisson(MeanNumberOfPhotons));
    //      }
    //
    //
    //
    //        if (NumPhotons <= 0)
    //        {
    //  // return unchanged particle and no secondaries
    //
    //           aParticleChange.SetNumberOfSecondaries(0);
    //
    //           return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
    //        }
    //
    //
    //       aParticleChange.SetNumberOfSecondaries(NumPhotons);
    //
    //
    //        if (fTrackSecondariesFirst) {
    //           if (aTrack.GetTrackStatus() == fAlive )
    //                  aParticleChange.ProposeTrackStatus(fSuspend);
    //        }
    //
    //
    //
    //
    ////////////////////////////////////////////////////////////////////////////////////
    //

    ////////////////////////////////////////////////////////////////////////////////////
    //  The fast sim way - Ben J, Nov 2012
    ////////////////////////////////////////////////////////////////////////////////////
    //
    //
    // 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();
    // double stepEnergy          = larg4::IonizationAndScintillation::Instance()->VisibleEnergyDeposit()/CLHEP::MeV;
    RecordPhotonsProduced(aStep, MeanNumberOfPhotons); //, stepEnergy);
    if (verboseLevel > 0) {
      G4cout << "\n Exiting from OpFastScintillation::DoIt -- NumberOfSecondaries = "
             << aParticleChange.GetNumberOfSecondaries() << G4endl;
    }
    return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
  }

void larg4::OpFastScintillation::ProcessStep ( const G4Step &  step )

private

Definition at line 522 of file OpFastScintillation.cxx.

  {
    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()),
      //step.GetTrack()->GetTrackID(),
      ParticleListAction::GetCurrentTrackID(),
      step.GetTrack()->GetParticleDefinition()->GetPDGEncoding(),
      step.GetPreStepPoint()->GetPhysicalVolume()->GetName());
  }

void larg4::OpFastScintillation::propagationTime ( std::vector< double > &  arrival_time_dist, G4ThreeVector  x0, const size_t  OpChannel, bool  Reflected = false  )

private

Definition at line 1028 of file OpFastScintillation.cxx.

  {
    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
      }
    }
  }

bool larg4::OpFastScintillation::RecordPhotonsProduced ( const G4Step &  aStep, double  N  )

protected

Definition at line 545 of file OpFastScintillation.cxx.

  {
    // 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();
    // unused G4StepPoint const* pPostStepPoint = aStep.GetPostStepPoint();

    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() - fGlobalTimeOffset;
    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 Kazu to debug # photons generated using csv file, by default should be commented out
    // but don't remove as it's useful. Separated portion of code relevant to this block
    // is labeled as "CASE-DEBUG DO NOT REMOVE THIS COMMENT"
    // (2017 May 2nd ... "kazuhiro at nevis dot columbia dot edu")
    //
    static bool first_time=false;
    if(!first_time) {
    std::cout<<"LOGMEid,pdg,mass,ke,te,de,x,y,z,t,dr,mean_npe,gen_npe,det_npe"<<std::endl;
    first_time=true;
    }

    std::cout<<"LOGME"
    <<aStep.GetTrack()->GetTrackID()<<","
    <<aParticle->GetDefinition()->GetPDGEncoding()<<","
    <<aParticle->GetDefinition()->GetPDGMass()/CLHEP::MeV<<","
    <<pPreStepPoint->GetKineticEnergy()<<","
    <<pPreStepPoint->GetTotalEnergy()<<","
    <<aStep.GetTotalEnergyDeposit()<<","
    <<ScintPoint<<","<<t0<<","
    <<aStep.GetDeltaPosition().mag()<<","
    <<MeanNumberOfPhotons<<","<<std::flush;

    double gen_photon_ctr=0;
    double det_photon_ctr=0;
    */
    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;
      //gen_photon_ctr += Num; // CASE-DEBUG DO NOT REMOVE THIS COMMENT
      // Max Scintillation Integral
      //            G4double CIImax = ScintillationIntegral->GetMaxValue();
      //std::cout << "++++++++++++" << Num << "++++++++++" << std::endl;

      // here we go: now if visibilities are invalid, we are in trouble
      //if (!Visibilities && (fPVS->NOpChannels() > 0)) {
      //  throw cet::exception("OpFastScintillator")
      //    << "Photon library does not cover point " << ScintPoint << " cm.\n";
      //}

      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) {
            //std::cout<<"VUV time correction: "<<arrival_time_dist[i]<<std::endl;
            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);
        }
      }
    }
    //std::cout<<gen_photon_ctr<<","<<det_photon_ctr<<std::endl; // CASE-DEBUG DO NOT REMOVE THIS COMMENT
    return 0;
  }

double larg4::OpFastScintillation::Rectangle_SolidAngle ( const double  a, const double  b, const double  d  ) const

private

Definition at line 2075 of file OpFastScintillation.cxx.

  {
    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 * std::acos(std::sqrt(aux));
    return 4. * fast_acos(std::sqrt(aux));
  }

double larg4::OpFastScintillation::Rectangle_SolidAngle ( Dims const &  o, const std::array< double, 3 >  v  ) const

private

Definition at line 2086 of file OpFastScintillation.cxx.

  {
    // 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!
    // std::cout << "Warning: invalid solid angle call." << std::endl;
    return 0.;
  }

double larg4::OpFastScintillation::reemission_energy ( ) const

private

Definition at line 1095 of file OpFastScintillation.cxx.

  {
    return fTPBEm->fire() * ((*(--tpbemission.end())).first - (*tpbemission.begin()).first) +
           (*tpbemission.begin()).first;
  }

void larg4::OpFastScintillation::RemoveSaturation ( )

inline

Definition at line 223 of file OpFastScintillation.hh.

    {
      emSaturation = NULL;
    }

G4double larg4::OpFastScintillation::sample_time ( const G4double  tau1, const G4double  tau2  ) const

private

Definition at line 1074 of file OpFastScintillation.cxx.

  {
    // 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;
  }

G4double larg4::OpFastScintillation::scint_time ( const G4Step &  aStep, G4double  ScintillationTime, G4double  ScintillationRiseTime  ) const

private

Definition at line 1011 of file OpFastScintillation.cxx.

  {
    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 larg4::OpFastScintillation::SetFiniteRiseTime ( const G4bool  state )

inline

Definition at line 487 of file OpFastScintillation.hh.

  {
    fFiniteRiseTime = state;
  }

void larg4::OpFastScintillation::SetScintillationByParticleType

(

const G4bool

scintType

)

Definition at line 984 of file OpFastScintillation.cxx.

  {
    if (emSaturation) {
      G4Exception("OpFastScintillation::SetScintillationByParticleType",
                  "Scint02",
                  JustWarning,
                  "Redefinition: Birks Saturation is replaced by ScintillationByParticleType!");
      RemoveSaturation();
    }
    scintillationByParticleType = scintType;
  }

void larg4::OpFastScintillation::SetScintillationExcitationRatio ( const G4double  excitationratio )

inline

Definition at line 517 of file OpFastScintillation.hh.

  {
    ExcitationRatio = excitationratio;
  }

void larg4::OpFastScintillation::SetScintillationYieldFactor ( const G4double  yieldfactor )

inline

Definition at line 505 of file OpFastScintillation.hh.

  {
    YieldFactor = yieldfactor;
  }

void larg4::OpFastScintillation::SetTrackSecondariesFirst ( const G4bool  state )

inline

Definition at line 481 of file OpFastScintillation.hh.

  {
    fTrackSecondariesFirst = state;
  }

G4double larg4::OpFastScintillation::single_exp ( const G4double  t, const G4double  tau2  ) const

private

Definition at line 1843 of file OpFastScintillation.cxx.

  {
    return std::exp(-1.0 * t / tau2) / tau2;
  }

bool larg4::OpFastScintillation::usesSemiAnalyticModel ( ) const

private

Returns whether the semi-analytic visibility parametrization is being used.

Definition at line 1492 of file OpFastScintillation.cxx.

  {
    return fUseNhitsModel;
  } // OpFastScintillation::usesSemiAnalyticModel()

int larg4::OpFastScintillation::VISHits ( geo::Point_t const &  ScintPoint, OpticalDetector const &  opDet, const double  cathode_hits_rec, const std::array< double, 3 >  hotspot  ) const

private

Definition at line 1691 of file OpFastScintillation.cxx.

  {
    // 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 = std::acos(cosine_vis) * 180. / CLHEP::pi;
    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
    double r = dist(ScintPoint, fcathode_centre, 2, 1);      // radial distance from centre of detector (Y-Z)
    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;
  }

int larg4::OpFastScintillation::VUVHits ( const double  Nphotons_created, geo::Point_t const &  ScintPoint, OpticalDetector const &  opDet  ) const

private

Definition at line 1595 of file OpFastScintillation.cxx.

  {
    // 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 absorbtion 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)
    double r = dist(ScintPoint, fcathode_centre, 2, 1);
    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;
  }

Member Data Documentation

bool const larg4::OpFastScintillation::bPropagate

private

Whether propagation of photons is enabled.

Definition at line 425 of file OpFastScintillation.hh.

G4EmSaturation* larg4::OpFastScintillation::emSaturation

private

Definition at line 338 of file OpFastScintillation.hh.

G4double larg4::OpFastScintillation::ExcitationRatio

protected

Definition at line 290 of file OpFastScintillation.hh.

std::vector<geo::BoxBoundedGeo> const larg4::OpFastScintillation::fActiveVolumes

private

Definition at line 409 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fangle_bin_timing_vis

private

Definition at line 360 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fangle_bin_timing_vuv

private

Definition at line 351 of file OpFastScintillation.hh.

std::vector<double> larg4::OpFastScintillation::fborder_corr_angulo_dome

private

Definition at line 390 of file OpFastScintillation.hh.

std::vector<double> larg4::OpFastScintillation::fborder_corr_angulo_flat

private

Definition at line 385 of file OpFastScintillation.hh.

std::vector<std::vector<double> > larg4::OpFastScintillation::fborder_corr_dome

private

Definition at line 391 of file OpFastScintillation.hh.

std::vector<std::vector<double> > larg4::OpFastScintillation::fborder_corr_flat

private

Definition at line 386 of file OpFastScintillation.hh.

TVector3 larg4::OpFastScintillation::fcathode_centre

private

Definition at line 408 of file OpFastScintillation.hh.

Dims larg4::OpFastScintillation::fcathode_plane

private

Definition at line 415 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fcathode_ydimension

private

Definition at line 407 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fcathode_zdimension

private

Definition at line 407 of file OpFastScintillation.hh.

std::vector<std::vector<std::vector<double> > > larg4::OpFastScintillation::fcut_off_pars

private

Definition at line 363 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fdelta_angulo_vis

private

Definition at line 396 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fdelta_angulo_vuv

private

Definition at line 381 of file OpFastScintillation.hh.

std::vector<double> larg4::OpFastScintillation::fdistances_refl

private

Definition at line 361 of file OpFastScintillation.hh.

G4bool larg4::OpFastScintillation::fFiniteRiseTime

protected

Definition at line 286 of file OpFastScintillation.hh.

std::vector<std::vector<double> > larg4::OpFastScintillation::fGHvuvpars_dome

private

Definition at line 389 of file OpFastScintillation.hh.

std::vector<std::vector<double> > larg4::OpFastScintillation::fGHvuvpars_flat

private

Definition at line 384 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::finflexion_point_distance

private

Definition at line 351 of file OpFastScintillation.hh.

bool larg4::OpFastScintillation::fIsDomePDCorr

private

Definition at line 388 of file OpFastScintillation.hh.

bool larg4::OpFastScintillation::fIsFlatPDCorr

private

Definition at line 383 of file OpFastScintillation.hh.

int larg4::OpFastScintillation::fL_abs_vuv

private

Definition at line 416 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fmax_d

private

Definition at line 351 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fmin_d

private

Definition at line 351 of file OpFastScintillation.hh.

bool const larg4::OpFastScintillation::fOnlyActiveVolume = false

private

Whether photon propagation is performed only from active volumes.

Definition at line 433 of file OpFastScintillation.hh.

bool const larg4::OpFastScintillation::fOnlyOneCryostat = false

private

Allows running even if light on cryostats C:1 and higher is not supported. Currently hard coded "no"

Definition at line 436 of file OpFastScintillation.hh.

bool const larg4::OpFastScintillation::fOpaqueCathode = false

private

Whether the cathodes are fully opaque; currently hard coded "no".

Definition at line 438 of file OpFastScintillation.hh.

std::vector<geo::Point_t> larg4::OpFastScintillation::fOpDetCenter

private

Definition at line 417 of file OpFastScintillation.hh.

std::vector<double> larg4::OpFastScintillation::fOpDetHeight

private

Definition at line 420 of file OpFastScintillation.hh.

std::vector<double> larg4::OpFastScintillation::fOpDetLength

private

Definition at line 419 of file OpFastScintillation.hh.

std::vector<int> larg4::OpFastScintillation::fOpDetType

private

Definition at line 418 of file OpFastScintillation.hh.

std::vector<std::vector<double> > larg4::OpFastScintillation::fparameters[7]

private

Definition at line 352 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fplane_depth

private

Definition at line 407 of file OpFastScintillation.hh.

phot::PhotonVisibilityService const* const larg4::OpFastScintillation::fPVS

private

Photon visibility service instance.

Definition at line 428 of file OpFastScintillation.hh.

std::vector<double> larg4::OpFastScintillation::fradial_distances_refl

private

Definition at line 362 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fradius

private

Definition at line 414 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fstep_size

private

Definition at line 351 of file OpFastScintillation.hh.

bool larg4::OpFastScintillation::fStoreReflected

private

Definition at line 394 of file OpFastScintillation.hh.

std::vector<std::vector<std::vector<double> > > larg4::OpFastScintillation::ftau_pars

private

Definition at line 364 of file OpFastScintillation.hh.

std::unique_ptr<CLHEP::RandGeneral> larg4::OpFastScintillation::fTPBEm

private

Definition at line 334 of file OpFastScintillation.hh.

G4bool larg4::OpFastScintillation::fTrackSecondariesFirst

protected

Definition at line 285 of file OpFastScintillation.hh.

bool const larg4::OpFastScintillation::fUseNhitsModel = false

private

Whether the semi-analytic model is being used for photon visibility.

Definition at line 431 of file OpFastScintillation.hh.

std::vector<double> larg4::OpFastScintillation::fvis_distances_r_dome

private

Definition at line 403 of file OpFastScintillation.hh.

std::vector<double> larg4::OpFastScintillation::fvis_distances_r_flat

private

Definition at line 399 of file OpFastScintillation.hh.

std::vector<double> larg4::OpFastScintillation::fvis_distances_x_dome

private

Definition at line 402 of file OpFastScintillation.hh.

std::vector<double> larg4::OpFastScintillation::fvis_distances_x_flat

private

Definition at line 398 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fvis_vmean

private

Definition at line 360 of file OpFastScintillation.hh.

std::vector<std::vector<std::vector<double> > > larg4::OpFastScintillation::fvispars_dome

private

Definition at line 404 of file OpFastScintillation.hh.

std::vector<std::vector<std::vector<double> > > larg4::OpFastScintillation::fvispars_flat

private

Definition at line 400 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fvuv_vgroup_max

private

Definition at line 351 of file OpFastScintillation.hh.

double larg4::OpFastScintillation::fvuv_vgroup_mean

private

Definition at line 351 of file OpFastScintillation.hh.

phot::MappedFunctions_t larg4::OpFastScintillation::ParPropTimeTF1

private

Definition at line 340 of file OpFastScintillation.hh.

phot::MappedT0s_t larg4::OpFastScintillation::ReflT0s

private

Definition at line 341 of file OpFastScintillation.hh.

G4bool larg4::OpFastScintillation::scintillationByParticleType

protected

Definition at line 292 of file OpFastScintillation.hh.

std::unique_ptr<G4PhysicsTable> larg4::OpFastScintillation::theFastIntegralTable

protected

Definition at line 283 of file OpFastScintillation.hh.

std::unique_ptr<G4PhysicsTable> larg4::OpFastScintillation::theSlowIntegralTable

protected

Definition at line 282 of file OpFastScintillation.hh.

std::map<double, double> larg4::OpFastScintillation::tpbemission

private

Definition at line 333 of file OpFastScintillation.hh.

std::vector<std::vector<double> > larg4::OpFastScintillation::VUV_max

private

Definition at line 356 of file OpFastScintillation.hh.

std::vector<std::vector<double> > larg4::OpFastScintillation::VUV_min

private

Definition at line 357 of file OpFastScintillation.hh.

std::vector<std::vector<TF1> > larg4::OpFastScintillation::VUV_timing

private

Definition at line 354 of file OpFastScintillation.hh.

G4double larg4::OpFastScintillation::YieldFactor

protected

Definition at line 288 of file OpFastScintillation.hh.

---

The documentation for this class was generated from the following files:

• larsim/larsim/LegacyLArG4/OpFastScintillation.hh
• larsim/larsim/LegacyLArG4/OpFastScintillation.cxx