#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

geo::Point_t	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().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().GetCathodeCenter();
         fcathode_centre.SetY(fActiveVolumes[0].CenterY());
         fcathode_centre.SetZ(
           fActiveVolumes[0].CenterZ()); // to get full cathode dimension rather than just single tpc
 
         if (fPVS->StoreReflected()) {
           fStoreReflected = true;
           // Load corrections for VIS semi-anlytic hits
           std::cout << "Loading visible light corrections" << std::endl;
           fPVS->LoadVisSemiAnalyticProperties(fdelta_angulo_vis, fradius);
           if (fIsFlatPDCorr) {
             fPVS->LoadVisParsFlat(fvis_distances_x_flat, fvis_distances_r_flat, fvispars_flat);
           }
           if (fIsDomePDCorr) {
             fPVS->LoadVisParsDome(fvis_distances_x_dome, fvis_distances_r_dome, fvispars_dome);
           }
 
           // cathode dimensions
           fcathode_ydimension = fActiveVolumes[0].SizeY();
           fcathode_zdimension = fActiveVolumes[0].SizeZ();
 
           // set cathode plane struct for solid angle function
           fcathode_plane.h = fcathode_ydimension;
           fcathode_plane.w = fcathode_zdimension;
           fplane_depth = std::abs(fcathode_centre.X());
         }
         else
           fStoreReflected = false;
       }
     }
     tpbemission = lar::providerFrom<detinfo::LArPropertiesService>()->TpbEm();
     std::vector<double> parent;
     parent.reserve(tpbemission.size());
     for (auto iter = tpbemission.begin(); iter != tpbemission.end(); ++iter) {
       parent.push_back(iter->second);
     }
     fTPBEm = std::make_unique<CLHEP::RandGeneral>(parent.data(), parent.size());
   }

larg4::OpFastScintillation::~OpFastScintillation ( )

Definition at line 386 of file OpFastScintillation.cxx.

References theFastIntegralTable, and theSlowIntegralTable.

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

Member Function Documentation

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

inline

Definition at line 215 of file OpFastScintillation.hh.

215 { emSaturation = sat; }

larg4::OpFastScintillation::emSaturation

G4EmSaturation * emSaturation

Definition: OpFastScintillation.hh:325

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

virtual

Definition at line 399 of file OpFastScintillation.cxx.

References PostStepDoIt().

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

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

private

Definition at line 1080 of file OpFastScintillation.cxx.

Referenced by RecordPhotonsProduced().

   {
     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 1833 of file OpFastScintillation.cxx.

Referenced by sample_time().

   { // 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 853 of file OpFastScintillation.cxx.

References theFastIntegralTable, and theSlowIntegralTable.

Referenced by OpFastScintillation().

   {
     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 1478 of file OpFastScintillation.cxx.

References util::counter(), fOpDetCenter, fOpDetHeight, fOpDetLength, fOpDetType, fPVS, isOpDetInSameTPC(), phot::PhotonVisibilityService::NOpChannels(), and VUVHits().

Referenced by RecordPhotonsProduced().

   {
     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 1500 of file OpFastScintillation.cxx.

References util::abs(), util::counter(), fborder_corr_angulo_flat, fborder_corr_flat, fcathode_centre, fcathode_plane, fGHvuvpars_flat, fIsFlatPDCorr, fL_abs_vuv, fOpDetCenter, fOpDetHeight, fOpDetLength, fOpDetType, fplane_depth, fPVS, Gaisser_Hillas(), interpolate(), isOpDetInSameTPC(), phot::PhotonVisibilityService::NOpChannels(), util::pi(), r, Rectangle_SolidAngle(), and VISHits().

Referenced by RecordPhotonsProduced().

   {
     // 1). calculate total number of hits of VUV photons on
     // reflective foils via solid angle + Gaisser-Hillas
     // corrections:
 
     // set plane_depth for correct TPC:
     double plane_depth;
     if (ScintPoint.X() < 0.) { plane_depth = -fplane_depth; }
     else {
       plane_depth = fplane_depth;
     }
 
     // get scintpoint coords relative to centre of cathode plane
     std::array<double, 3> ScintPoint_relative = {std::abs(ScintPoint.X() - plane_depth),
                                                  std::abs(ScintPoint.Y() - fcathode_centre.Y()),
                                                  std::abs(ScintPoint.Z() - fcathode_centre.Z())};
     // calculate solid angle of cathode from the scintillation point
     double solid_angle_cathode = Rectangle_SolidAngle(fcathode_plane, ScintPoint_relative);
 
     // calculate distance and angle between ScintPoint and hotspot
     // vast majority of hits in hotspot region directly infront of scintpoint,
     // therefore consider attenuation for this distance and on axis GH instead of for the centre coordinate
     double distance_cathode = std::abs(plane_depth - ScintPoint.X());
     // calculate hits on cathode plane via geometric acceptance
     double cathode_hits_geo = std::exp(-1. * distance_cathode / fL_abs_vuv) *
                               (solid_angle_cathode / (4. * CLHEP::pi)) * Num;
     // determine Gaisser-Hillas correction including border effects
     // use flat correction
     double r = std::sqrt(std::pow(ScintPoint.Y() - fcathode_centre.Y(), 2) +
                          std::pow(ScintPoint.Z() - fcathode_centre.Z(), 2));
     double pars_ini[4] = {0, 0, 0, 0};
     double s1 = 0;
     double s2 = 0;
     double s3 = 0;
     if (fIsFlatPDCorr) {
       pars_ini[0] = fGHvuvpars_flat[0][0];
       pars_ini[1] = fGHvuvpars_flat[1][0];
       pars_ini[2] = fGHvuvpars_flat[2][0];
       pars_ini[3] = fGHvuvpars_flat[3][0];
       s1 = interpolate(fborder_corr_angulo_flat, fborder_corr_flat[0], 0, true);
       s2 = interpolate(fborder_corr_angulo_flat, fborder_corr_flat[1], 0, true);
       s3 = interpolate(fborder_corr_angulo_flat, fborder_corr_flat[2], 0, true);
     }
     else
       std::cout
         << "Error: flat optical detector VUV correction required for reflected semi-analytic hits."
         << std::endl;
 
     // add border correction
     pars_ini[0] = pars_ini[0] + s1 * r;
     pars_ini[1] = pars_ini[1] + s2 * r;
     pars_ini[2] = pars_ini[2] + s3 * r;
     pars_ini[3] = pars_ini[3];
 
     // calculate corrected number of hits
     double GH_correction = Gaisser_Hillas(distance_cathode, pars_ini);
     const double cathode_hits_rec = GH_correction * cathode_hits_geo;
 
     // detemine hits on each PD
     const std::array<double, 3> hotspot = {plane_depth, ScintPoint.Y(), ScintPoint.Z()};
     for (size_t const OpDet : util::counter(fPVS->NOpChannels())) {
       if (!isOpDetInSameTPC(ScintPoint, fOpDetCenter.at(OpDet))) continue;
 
       // set detector struct for solid angle function
       const OpFastScintillation::OpticalDetector op{fOpDetHeight.at(OpDet),
                                                     fOpDetLength.at(OpDet),
                                                     fOpDetCenter.at(OpDet),
                                                     fOpDetType.at(OpDet)};
 
       int const ReflDetThis = VISHits(ScintPoint, op, cathode_hits_rec, hotspot);
       if (ReflDetThis > 0) { ReflDetectedNum[OpDet] = ReflDetThis; }
     }
   }

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

private

Definition at line 1990 of file OpFastScintillation.cxx.

References d, e, larg4::isApproximatelyEqual(), larg4::isApproximatelyZero(), larg4::isDefinitelyGreaterThan(), larg4::isDefinitelyLessThan(), and util::pi().

Referenced by VISHits(), and VUVHits().

   {
     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 516 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 2149 of file OpFastScintillation.cxx.

References geo::BoxBoundedGeo::ExtendToInclude(), geo::GeometryCore::Iterate(), and geo::GeometryCore::Ncryostats().

Referenced by OpFastScintillation().

   {
     std::vector<geo::BoxBoundedGeo> activeVolumes;
     activeVolumes.reserve(geom.Ncryostats());
 
     for (geo::CryostatGeo const& cryo : geom.Iterate<geo::CryostatGeo>()) {
 
       // can't use it default-constructed since it would always include origin
       geo::BoxBoundedGeo box{cryo.TPC(0).ActiveBoundingBox()};
 
       for (geo::TPCGeo const& TPC : cryo.IterateTPCs())
         box.ExtendToInclude(TPC.ActiveBoundingBox());
 
       activeVolumes.push_back(std::move(box));
 
     } // for cryostats
 
     return activeVolumes;
   } // OpFastScintillation::extractActiveVolumes()

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

private

Definition at line 1841 of file OpFastScintillation.cxx.

Referenced by detectedReflecHits(), and VUVHits().

   {
     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 1229 of file OpFastScintillation.cxx.

References finflexion_point_distance, larg4::finter_d(), fmin_d, fparameters, fstep_size, fvuv_vgroup_max, fvuv_vgroup_mean, interpolate(), interpolate3(), larg4::model_close(), larg4::model_far(), VUV_max, VUV_min, and VUV_timing.

Referenced by getVUVTimes().

   {
     // 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 511 of file OpFastScintillation.hh.

   {
     return theFastIntegralTable.get();
   }

G4bool larg4::OpFastScintillation::GetFiniteRiseTime ( ) const

inline

Definition at line 481 of file OpFastScintillation.hh.

   {
     return fFiniteRiseTime;
   }

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

Definition at line 978 of file OpFastScintillation.cxx.

   {
     *condition = StronglyForced;
     return DBL_MAX;
   }

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

Definition at line 986 of file OpFastScintillation.cxx.

   {
     *condition = Forced;
     return DBL_MAX;
   }

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

inline

Definition at line 221 of file OpFastScintillation.hh.

221 { return emSaturation; }

larg4::OpFastScintillation::emSaturation

G4EmSaturation * emSaturation

Definition: OpFastScintillation.hh:325

G4bool larg4::OpFastScintillation::GetScintillationByParticleType ( ) const

inline

Definition at line 228 of file OpFastScintillation.hh.

228 { return scintillationByParticleType; }

larg4::OpFastScintillation::scintillationByParticleType

G4bool scintillationByParticleType

Definition: OpFastScintillation.hh:280

G4double larg4::OpFastScintillation::GetScintillationExcitationRatio ( ) const

inline

Definition at line 501 of file OpFastScintillation.hh.

   {
     return ExcitationRatio;
   }

G4double larg4::OpFastScintillation::GetScintillationYieldFactor ( ) const

inline

Definition at line 491 of file OpFastScintillation.hh.

   {
     return YieldFactor;
   }

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

inline

Definition at line 506 of file OpFastScintillation.hh.

   {
     return theSlowIntegralTable.get();
   }

G4bool larg4::OpFastScintillation::GetTrackSecondariesFirst ( ) const

inline

Definition at line 476 of file OpFastScintillation.hh.

   {
     return fTrackSecondariesFirst;
   }

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

Definition at line 1354 of file OpFastScintillation.cxx.

References util::abs(), util::counter(), fangle_bin_timing_vis, larg4::fast_acos(), fcathode_centre, fcut_off_pars, fdistances_refl, fmin_d, fplane_depth, fradial_distances_refl, fstep_size, ftau_pars, fvis_vmean, fvuv_vgroup_max, getVUVTimes(), interpolate(), util::pi(), r, util::size(), VUV_min, and x.

Referenced by propagationTime().

   {
     // *************************************************************************************************
     //     Calculation of earliest arrival times and corresponding unsmeared distribution
     // *************************************************************************************************
 
     // set plane_depth for correct TPC:
     double plane_depth;
     if (ScintPoint[0] < 0) { plane_depth = -fplane_depth; }
     else {
       plane_depth = fplane_depth;
     }
 
     // calculate point of reflection for shortest path
     TVector3 bounce_point(plane_depth, ScintPoint[1], ScintPoint[2]);
 
     // calculate distance travelled by VUV light and by vis light
     double VUVdist = (bounce_point - ScintPoint).Mag();
     double Visdist = (OpDetPoint - bounce_point).Mag();
 
     // calculate times taken by VUV part of path
     int angle_bin_vuv = 0; // on-axis by definition
     getVUVTimes(arrivalTimes, VUVdist, angle_bin_vuv);
 
     // add visible direct path transport time
     for (size_t i = 0; i < arrivalTimes.size(); ++i) {
       arrivalTimes[i] += Visdist / fvis_vmean;
     }
 
     // *************************************************************************************************
     //      Smearing of arrival time distribution
     // *************************************************************************************************
     // calculate fastest time possible
     // vis part
     double vis_time = Visdist / fvis_vmean;
     // vuv part
     double vuv_time;
     if (VUVdist < fmin_d) { vuv_time = VUVdist / fvuv_vgroup_max; }
     else {
       // find index of required parameterisation
       const size_t index = std::round((VUVdist - fmin_d) / fstep_size);
       // find shortest time
       vuv_time = VUV_min[angle_bin_vuv][index];
     }
     // sum
     double fastest_time = vis_time + vuv_time;
 
     // calculate angle theta between bound_point and optical detector
     double cosine_theta = std::abs(OpDetPoint[0] - bounce_point[0]) / Visdist;
     double theta = fast_acos(cosine_theta) * 180. / CLHEP::pi;
 
     // determine smearing parameters using interpolation of generated points:
     // 1). tau = exponential smearing factor, varies with distance and angle
     // 2). cutoff = largest smeared time allowed, preventing excessively large
     //     times caused by exponential distance to cathode
     double distance_cathode_plane = std::abs(plane_depth - ScintPoint[0]);
     // angular bin
     size_t theta_bin = theta / fangle_bin_timing_vis;
     // radial distance from centre of TPC (y,z plane)
     double r = std::sqrt(std::pow(ScintPoint.Y() - fcathode_centre.Z(), 2) +
                          std::pow(ScintPoint.Y() - fcathode_centre.Z(), 2));
 
     // cut-off and tau
     // cut-off
     // interpolate in d_c for each r bin
     std::vector<double> interp_vals(fcut_off_pars[theta_bin].size(), 0.0);
     for (size_t i = 0; i < fcut_off_pars[theta_bin].size(); i++) {
       interp_vals[i] =
         interpolate(fdistances_refl, fcut_off_pars[theta_bin][i], distance_cathode_plane, true);
     }
     // interpolate in r
     double cutoff = interpolate(fradial_distances_refl, interp_vals, r, true);
 
     // tau
     // interpolate in x for each r bin
     std::vector<double> interp_vals_tau(ftau_pars[theta_bin].size(), 0.0);
     for (size_t i = 0; i < ftau_pars[theta_bin].size(); i++) {
       interp_vals_tau[i] =
         interpolate(fdistances_refl, ftau_pars[theta_bin][i], distance_cathode_plane, true);
     }
     // interpolate in r
     double tau = interpolate(fradial_distances_refl, interp_vals_tau, r, true);
 
     if (tau < 0) { tau = 0; } // failsafe if tau extrapolate goes wrong
 
     // apply smearing:
     for (size_t i = 0; i < arrivalTimes.size(); ++i) {
       double arrival_time_smeared;
       // if time is already greater than cutoff, do not apply smearing
       if (arrivalTimes[i] >= cutoff) { continue; }
       // otherwise smear
       else {
         unsigned int counter = 0;
         // loop until time generated is within cutoff limit
         // most are within single attempt, very few take more than two
         do {
           // don't attempt smearings too many times
           if (counter >= 10) {                      // TODO: unhardcode
             arrival_time_smeared = arrivalTimes[i]; // don't smear
             break;
           }
           else {
             // generate random number in appropriate range
             double x = gRandom->Uniform(0.5, 1.0); // TODO: unhardcode
             // apply the exponential smearing
             arrival_time_smeared =
               arrivalTimes[i] + (arrivalTimes[i] - fastest_time) * (std::pow(x, -tau) - 1);
           }
           counter++;
         } while (arrival_time_smeared > cutoff);
       }
       arrivalTimes[i] = arrival_time_smeared;
     }
   }

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

Definition at line 1329 of file OpFastScintillation.cxx.

References fmin_d, fstep_size, fvuv_vgroup_mean, generateParam(), VUV_max, VUV_min, and VUV_timing.

Referenced by getVISTimes(), and propagationTime().

   {
     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 1911 of file OpFastScintillation.cxx.

References util::size().

Referenced by detectedReflecHits(), generateParam(), getVISTimes(), OpFastScintillation(), VISHits(), and VUVHits().

   {
     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 1939 of file OpFastScintillation.cxx.

References util::size().

Referenced by generateParam().

   {
     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 458 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 1809 of file OpFastScintillation.cxx.

References util::abs().

Referenced by detectedDirectHits(), detectedReflecHits(), and RecordPhotonsProduced().

   {
     // 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 1822 of file OpFastScintillation.cxx.

References fActiveVolumes.

Referenced by RecordPhotonsProduced().

   {
     //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 2118 of file OpFastScintillation.cxx.

References fradius, and util::pi().

Referenced by VISHits(), and VUVHits().

   {
     // 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 409 of file OpFastScintillation.cxx.

References larg4::IonizationAndScintillation::Instance(), larg4::IonizationAndScintillation::NumberScintillationPhotons(), RecordPhotonsProduced(), and larg4::IonizationAndScintillation::Reset().

Referenced by AtRestDoIt().

   {
     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 505 of file OpFastScintillation.cxx.

References larg4::OpDetPhotonTable::AddEnergyDeposit(), larg4::ParticleListAction::GetCurrentOrigTrackID(), larg4::ParticleListAction::GetCurrentTrackID(), and larg4::OpDetPhotonTable::Instance().

Referenced by RecordPhotonsProduced().

   {
     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(),
       ParticleListAction::GetCurrentOrigTrackID(),
       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 1009 of file OpFastScintillation.cxx.

References util::abs(), fangle_bin_timing_vuv, larg4::fast_acos(), fOpDetCenter, fPVS, getVISTimes(), getVUVTimes(), phot::PhotonVisibilityService::IncludeParPropTime(), phot::PhotonVisibilityService::IncludePropTime(), ParPropTimeTF1, util::pi(), and geo::vect::toTVector3().

Referenced by RecordPhotonsProduced().

   {
     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 529 of file OpFastScintillation.cxx.

References larg4::OpDetPhotonTable::AddLitePhoton(), larg4::OpDetPhotonTable::AddOpDetBacktrackerRecord(), larg4::OpDetPhotonTable::AddPhoton(), sim::OpDetBacktrackerRecord::AddScintillationPhotons(), average_position(), bPropagate, util::counter(), detectedDirectHits(), detectedReflecHits(), emSaturation, sim::OnePhoton::Energy, ExcitationRatio, fFiniteRiseTime, sim::LArG4Parameters::FillSimEnergyDeposits(), fOnlyActiveVolume, fOpaqueCathode, fOpDetCenter, fPVS, phot::PhotonVisibilityService::GetAllVisibilities(), larg4::ParticleListAction::GetCurrentTrackID(), phot::PhotonVisibilityService::GetReflT0s(), phot::PhotonVisibilityService::GetTimingTF1(), phot::PhotonVisibilityService::IncludeParPropTime(), sim::OnePhoton::InitialPosition, larg4::IonizationAndScintillation::Instance(), larg4::OpDetPhotonTable::Instance(), isOpDetInSameTPC(), isScintInActiveVolume(), sim::OnePhoton::MotherTrackID, phot::PhotonVisibilityService::NOpChannels(), ParPropTimeTF1, ProcessStep(), propagationTime(), reemission_energy(), ReflT0s, scint_time(), scintillationByParticleType, sim::OnePhoton::SetInSD, phot::PhotonVisibilityService::StoreReflected(), phot::PhotonVisibilityService::StoreReflT0(), theFastIntegralTable, theSlowIntegralTable, sim::OnePhoton::Time, sim::LArG4Parameters::UseLitePhotons(), usesSemiAnalyticModel(), and larg4::IonizationAndScintillation::VisibleEnergyDeposit().

Referenced by PostStepDoIt().

   {
     // 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 2051 of file OpFastScintillation.cxx.

References d, and larg4::fast_acos().

Referenced by detectedReflecHits(), Rectangle_SolidAngle(), VISHits(), and VUVHits().

   {
     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 2063 of file OpFastScintillation.cxx.

References larg4::OpFastScintillation::Dims::h, larg4::isApproximatelyZero(), larg4::isDefinitelyGreaterThan(), Rectangle_SolidAngle(), and larg4::OpFastScintillation::Dims::w.

   {
     // 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 1074 of file OpFastScintillation.cxx.

References fTPBEm, and tpbemission.

Referenced by RecordPhotonsProduced().

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

void larg4::OpFastScintillation::RemoveSaturation ( )

inline

Definition at line 218 of file OpFastScintillation.hh.

Referenced by SetScintillationByParticleType().

218 { emSaturation = NULL; }

larg4::OpFastScintillation::emSaturation

G4EmSaturation * emSaturation

Definition: OpFastScintillation.hh:325

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

private

Definition at line 1054 of file OpFastScintillation.cxx.

References bi_exp(), d, and single_exp().

Referenced by scint_time().

   {
     // 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 992 of file OpFastScintillation.cxx.

References sample_time().

Referenced by RecordPhotonsProduced().

   {
     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 471 of file OpFastScintillation.hh.

   {
     fFiniteRiseTime = state;
   }

void larg4::OpFastScintillation::SetScintillationByParticleType ( const G4bool scintType )

Definition at line 966 of file OpFastScintillation.cxx.

References emSaturation, RemoveSaturation(), and scintillationByParticleType.

   {
     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 496 of file OpFastScintillation.hh.

   {
     ExcitationRatio = excitationratio;
   }

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

inline

Definition at line 486 of file OpFastScintillation.hh.

   {
     YieldFactor = yieldfactor;
   }

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

inline

Definition at line 466 of file OpFastScintillation.hh.

   {
     fTrackSecondariesFirst = state;
   }

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

private

Definition at line 1828 of file OpFastScintillation.cxx.

Referenced by sample_time().

   {
     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 1472 of file OpFastScintillation.cxx.

References fUseNhitsModel.

Referenced by OpFastScintillation(), and RecordPhotonsProduced().

   {
     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 1681 of file OpFastScintillation.cxx.

References util::abs(), d, Disk_SolidAngle(), larg4::dist(), larg4::fast_acos(), fcathode_centre, fdelta_angulo_vis, geo::vect::fillCoords(), fIsDomePDCorr, fIsFlatPDCorr, fplane_depth, fradius, fvis_distances_r_dome, fvis_distances_r_flat, fvis_distances_x_dome, fvis_distances_x_flat, fvispars_dome, fvispars_flat, larg4::OpFastScintillation::OpticalDetector::h, interpolate(), Omega_Dome_Model(), larg4::OpFastScintillation::OpticalDetector::OpDetPoint, util::pi(), r, Rectangle_SolidAngle(), util::size(), larg4::OpFastScintillation::OpticalDetector::type, and larg4::OpFastScintillation::OpticalDetector::w.

Referenced by detectedReflecHits().

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

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

private

Definition at line 1578 of file OpFastScintillation.cxx.

References util::abs(), d, Disk_SolidAngle(), larg4::dist(), larg4::fast_acos(), fborder_corr_angulo_dome, fborder_corr_angulo_flat, fborder_corr_dome, fborder_corr_flat, fcathode_centre, fdelta_angulo_vuv, fGHvuvpars_dome, fGHvuvpars_flat, geo::vect::fillCoords(), fIsDomePDCorr, fIsFlatPDCorr, fL_abs_vuv, fradius, Gaisser_Hillas(), larg4::OpFastScintillation::OpticalDetector::h, interpolate(), Omega_Dome_Model(), larg4::OpFastScintillation::OpticalDetector::OpDetPoint, util::pi(), r, Rectangle_SolidAngle(), larg4::OpFastScintillation::OpticalDetector::type, and larg4::OpFastScintillation::OpticalDetector::w.

Referenced by detectedDirectHits().

   {
     // 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)
     // FIXME (KJK): Original implementation only accounts for the Y component
     double r = std::abs(ScintPoint[1] - fcathode_centre.Y());
     double pars_ini[4] = {0, 0, 0, 0};
     double s1 = 0;
     double s2 = 0;
     double s3 = 0;
     // flat PDs
     if ((opDet.type == 0 || opDet.type == 2) && fIsFlatPDCorr) {
       pars_ini[0] = fGHvuvpars_flat[0][j];
       pars_ini[1] = fGHvuvpars_flat[1][j];
       pars_ini[2] = fGHvuvpars_flat[2][j];
       pars_ini[3] = fGHvuvpars_flat[3][j];
       s1 = interpolate(fborder_corr_angulo_flat, fborder_corr_flat[0], theta, true);
       s2 = interpolate(fborder_corr_angulo_flat, fborder_corr_flat[1], theta, true);
       s3 = interpolate(fborder_corr_angulo_flat, fborder_corr_flat[2], theta, true);
     }
     // dome PDs
     else if (opDet.type == 1 && fIsDomePDCorr) {
       pars_ini[0] = fGHvuvpars_dome[0][j];
       pars_ini[1] = fGHvuvpars_dome[1][j];
       pars_ini[2] = fGHvuvpars_dome[2][j];
       pars_ini[3] = fGHvuvpars_dome[3][j];
       s1 = interpolate(fborder_corr_angulo_dome, fborder_corr_dome[0], theta, true);
       s2 = interpolate(fborder_corr_angulo_dome, fborder_corr_dome[1], theta, true);
       s3 = interpolate(fborder_corr_angulo_dome, fborder_corr_dome[2], theta, true);
     }
     else
       std::cout << "Error: Invalid optical detector type. 0 = rectangular, 1 = dome, 2 = disk. Or "
                    "corrections for chosen optical detector type missing."
                 << std::endl;
 
     // add border correction
     pars_ini[0] = pars_ini[0] + s1 * r;
     pars_ini[1] = pars_ini[1] + s2 * r;
     pars_ini[2] = pars_ini[2] + s3 * r;
     pars_ini[3] = pars_ini[3];
 
     // calculate correction
     double GH_correction = Gaisser_Hillas(distance, pars_ini);
 
     // calculate number photons
     double hits_rec = GH_correction * hits_geo / cosine;
     int hits_vuv = std::round(G4Poisson(hits_rec));
 
     return hits_vuv;
   }

Member Data Documentation

bool const larg4::OpFastScintillation::bPropagate

private

Whether propagation of photons is enabled.

Definition at line 412 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and RecordPhotonsProduced().

G4EmSaturation* larg4::OpFastScintillation::emSaturation

private

Definition at line 325 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), RecordPhotonsProduced(), and SetScintillationByParticleType().

G4double larg4::OpFastScintillation::ExcitationRatio

protected

Definition at line 278 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and RecordPhotonsProduced().

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

private

Definition at line 397 of file OpFastScintillation.hh.

Referenced by isScintInActiveVolume(), and OpFastScintillation().

double larg4::OpFastScintillation::fangle_bin_timing_vis

private

Definition at line 348 of file OpFastScintillation.hh.

Referenced by getVISTimes(), and OpFastScintillation().

double larg4::OpFastScintillation::fangle_bin_timing_vuv

private

Definition at line 338 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and propagationTime().

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

private

Definition at line 378 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VUVHits().

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

private

Definition at line 373 of file OpFastScintillation.hh.

Referenced by detectedReflecHits(), OpFastScintillation(), and VUVHits().

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

private

Definition at line 379 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VUVHits().

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

private

Definition at line 374 of file OpFastScintillation.hh.

Referenced by detectedReflecHits(), OpFastScintillation(), and VUVHits().

geo::Point_t larg4::OpFastScintillation::fcathode_centre

private

Definition at line 396 of file OpFastScintillation.hh.

Referenced by detectedReflecHits(), getVISTimes(), OpFastScintillation(), VISHits(), and VUVHits().

Dims larg4::OpFastScintillation::fcathode_plane

private

Definition at line 402 of file OpFastScintillation.hh.

Referenced by detectedReflecHits(), and OpFastScintillation().

double larg4::OpFastScintillation::fcathode_ydimension

private

Definition at line 395 of file OpFastScintillation.hh.

Referenced by OpFastScintillation().

double larg4::OpFastScintillation::fcathode_zdimension

private

Definition at line 395 of file OpFastScintillation.hh.

Referenced by OpFastScintillation().

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

private

Definition at line 351 of file OpFastScintillation.hh.

Referenced by getVISTimes(), and OpFastScintillation().

double larg4::OpFastScintillation::fdelta_angulo_vis

private

Definition at line 384 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

double larg4::OpFastScintillation::fdelta_angulo_vuv

private

Definition at line 369 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VUVHits().

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

private

Definition at line 349 of file OpFastScintillation.hh.

Referenced by getVISTimes(), and OpFastScintillation().

G4bool larg4::OpFastScintillation::fFiniteRiseTime

protected

Definition at line 274 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and RecordPhotonsProduced().

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

private

Definition at line 377 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VUVHits().

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

private

Definition at line 372 of file OpFastScintillation.hh.

Referenced by detectedReflecHits(), OpFastScintillation(), and VUVHits().

double larg4::OpFastScintillation::finflexion_point_distance

private

Definition at line 338 of file OpFastScintillation.hh.

Referenced by generateParam(), and OpFastScintillation().

bool larg4::OpFastScintillation::fIsDomePDCorr

private

Definition at line 376 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), VISHits(), and VUVHits().

bool larg4::OpFastScintillation::fIsFlatPDCorr

private

Definition at line 371 of file OpFastScintillation.hh.

Referenced by detectedReflecHits(), OpFastScintillation(), VISHits(), and VUVHits().

int larg4::OpFastScintillation::fL_abs_vuv

private

Definition at line 403 of file OpFastScintillation.hh.

Referenced by detectedReflecHits(), OpFastScintillation(), and VUVHits().

double larg4::OpFastScintillation::fmax_d

private

Definition at line 338 of file OpFastScintillation.hh.

Referenced by OpFastScintillation().

double larg4::OpFastScintillation::fmin_d

private

Definition at line 338 of file OpFastScintillation.hh.

Referenced by generateParam(), getVISTimes(), getVUVTimes(), and OpFastScintillation().

bool const larg4::OpFastScintillation::fOnlyActiveVolume = false

private

Whether photon propagation is performed only from active volumes.

Definition at line 420 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and RecordPhotonsProduced().

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 423 of file OpFastScintillation.hh.

Referenced by OpFastScintillation().

bool const larg4::OpFastScintillation::fOpaqueCathode = false

private

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

Definition at line 425 of file OpFastScintillation.hh.

Referenced by RecordPhotonsProduced().

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

private

Definition at line 404 of file OpFastScintillation.hh.

Referenced by detectedDirectHits(), detectedReflecHits(), OpFastScintillation(), propagationTime(), and RecordPhotonsProduced().

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

private

Definition at line 407 of file OpFastScintillation.hh.

Referenced by detectedDirectHits(), detectedReflecHits(), and OpFastScintillation().

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

private

Definition at line 406 of file OpFastScintillation.hh.

Referenced by detectedDirectHits(), detectedReflecHits(), and OpFastScintillation().

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

private

Definition at line 405 of file OpFastScintillation.hh.

Referenced by detectedDirectHits(), detectedReflecHits(), and OpFastScintillation().

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

private

Definition at line 340 of file OpFastScintillation.hh.

Referenced by generateParam(), and OpFastScintillation().

double larg4::OpFastScintillation::fplane_depth

private

Definition at line 395 of file OpFastScintillation.hh.

Referenced by detectedReflecHits(), getVISTimes(), OpFastScintillation(), and VISHits().

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

private

Photon visibility service instance.

Definition at line 415 of file OpFastScintillation.hh.

Referenced by detectedDirectHits(), detectedReflecHits(), OpFastScintillation(), propagationTime(), and RecordPhotonsProduced().

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

private

Definition at line 350 of file OpFastScintillation.hh.

Referenced by getVISTimes(), and OpFastScintillation().

double larg4::OpFastScintillation::fradius

private

Definition at line 401 of file OpFastScintillation.hh.

Referenced by Omega_Dome_Model(), OpFastScintillation(), VISHits(), and VUVHits().

double larg4::OpFastScintillation::fstep_size

private

Definition at line 338 of file OpFastScintillation.hh.

Referenced by generateParam(), getVISTimes(), getVUVTimes(), and OpFastScintillation().

bool larg4::OpFastScintillation::fStoreReflected

private

Definition at line 382 of file OpFastScintillation.hh.

Referenced by OpFastScintillation().

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

private

Definition at line 352 of file OpFastScintillation.hh.

Referenced by getVISTimes(), and OpFastScintillation().

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

private

Definition at line 321 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and reemission_energy().

G4bool larg4::OpFastScintillation::fTrackSecondariesFirst

protected

Definition at line 273 of file OpFastScintillation.hh.

Referenced by OpFastScintillation().

bool const larg4::OpFastScintillation::fUseNhitsModel = false

private

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

Definition at line 418 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and usesSemiAnalyticModel().

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

private

Definition at line 391 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

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

private

Definition at line 387 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

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

private

Definition at line 390 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

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

private

Definition at line 386 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

double larg4::OpFastScintillation::fvis_vmean

private

Definition at line 348 of file OpFastScintillation.hh.

Referenced by getVISTimes(), and OpFastScintillation().

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

private

Definition at line 392 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

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

private

Definition at line 388 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

double larg4::OpFastScintillation::fvuv_vgroup_max

private

Definition at line 338 of file OpFastScintillation.hh.

Referenced by generateParam(), getVISTimes(), and OpFastScintillation().

double larg4::OpFastScintillation::fvuv_vgroup_mean

private

Definition at line 338 of file OpFastScintillation.hh.

Referenced by generateParam(), getVUVTimes(), and OpFastScintillation().

phot::MappedFunctions_t larg4::OpFastScintillation::ParPropTimeTF1

private

Definition at line 327 of file OpFastScintillation.hh.

Referenced by propagationTime(), and RecordPhotonsProduced().

phot::MappedT0s_t larg4::OpFastScintillation::ReflT0s

private

Definition at line 328 of file OpFastScintillation.hh.

Referenced by RecordPhotonsProduced().

G4bool larg4::OpFastScintillation::scintillationByParticleType

protected

Definition at line 280 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), RecordPhotonsProduced(), and SetScintillationByParticleType().

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

protected

Definition at line 271 of file OpFastScintillation.hh.

Referenced by BuildThePhysicsTable(), RecordPhotonsProduced(), and ~OpFastScintillation().

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

protected

Definition at line 270 of file OpFastScintillation.hh.

Referenced by BuildThePhysicsTable(), RecordPhotonsProduced(), and ~OpFastScintillation().

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

private

Definition at line 320 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and reemission_energy().

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

private

Definition at line 344 of file OpFastScintillation.hh.

Referenced by generateParam(), getVUVTimes(), and OpFastScintillation().

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

private

Definition at line 345 of file OpFastScintillation.hh.

Referenced by generateParam(), getVISTimes(), getVUVTimes(), and OpFastScintillation().

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

private

Definition at line 342 of file OpFastScintillation.hh.

Referenced by generateParam(), getVUVTimes(), and OpFastScintillation().

G4double larg4::OpFastScintillation::YieldFactor

protected

Definition at line 276 of file OpFastScintillation.hh.

Referenced by OpFastScintillation().

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

larsim/v09_43_00/source/larsim/LegacyLArG4/OpFastScintillation.hh
larsim/v09_43_00/source/larsim/LegacyLArG4/OpFastScintillation.cxx

Classes

Public Member Functions

Protected Member Functions

Protected Attributes

Private Member Functions

Static Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation