LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
larg4::OpFastScintillation Class Reference

#include "OpFastScintillation.hh"

Inheritance diagram for larg4::OpFastScintillation:

Classes

struct  Dims
 
struct  OpticalDetector
 

Public Member Functions

 OpFastScintillation (const G4String &processName="Scintillation", G4ProcessType type=fElectromagnetic)
 
 ~OpFastScintillation ()
 
virtual G4bool IsApplicable (const G4ParticleDefinition &aParticleType)
 
G4double GetMeanFreePath (const G4Track &aTrack, G4double, G4ForceCondition *)
 
G4double GetMeanLifeTime (const G4Track &aTrack, G4ForceCondition *)
 
virtual G4VParticleChange * PostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
 
virtual G4VParticleChange * AtRestDoIt (const G4Track &aTrack, const G4Step &aStep)
 
void SetTrackSecondariesFirst (const G4bool state)
 
void SetFiniteRiseTime (const G4bool state)
 
G4bool GetTrackSecondariesFirst () const
 
G4bool GetFiniteRiseTime () const
 
void SetScintillationYieldFactor (const G4double yieldfactor)
 
G4double GetScintillationYieldFactor () const
 
void SetScintillationExcitationRatio (const G4double excitationratio)
 
G4double GetScintillationExcitationRatio () const
 
G4PhysicsTable * GetFastIntegralTable () const
 
G4PhysicsTable * GetSlowIntegralTable () const
 
void AddSaturation (G4EmSaturation *sat)
 
void RemoveSaturation ()
 
G4EmSaturation * GetSaturation () const
 
void SetScintillationByParticleType (const G4bool)
 
G4bool GetScintillationByParticleType () const
 
void DumpPhysicsTable () const
 
void getVUVTimes (std::vector< double > &arrivalTimes, const double distance_in_cm, const size_t angle_bin)
 
void generateParam (const size_t index, const size_t angle_bin)
 
void getVISTimes (std::vector< double > &arrivalTimes, const TVector3 &ScintPoint, const TVector3 &OpDetPoint)
 
void detectedDirectHits (std::map< size_t, int > &DetectedNum, const double Num, geo::Point_t const &ScintPoint) const
 
void detectedReflecHits (std::map< size_t, int > &ReflDetectedNum, const double Num, geo::Point_t const &ScintPoint) const
 

Protected Member Functions

void BuildThePhysicsTable ()
 
bool RecordPhotonsProduced (const G4Step &aStep, double N)
 

Protected Attributes

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

Private Member Functions

bool usesSemiAnalyticModel () const
 Returns whether the semi-analytic visibility parametrization is being used. More...
 
int VUVHits (const double Nphotons_created, geo::Point_t const &ScintPoint, OpticalDetector const &opDet) const
 
int VISHits (geo::Point_t const &ScintPoint, OpticalDetector const &opDet, const double cathode_hits_rec, const std::array< double, 3 > hotspot) const
 
G4double single_exp (const G4double t, const G4double tau2) const
 
G4double bi_exp (const G4double t, const G4double tau1, const G4double tau2) const
 
G4double scint_time (const G4Step &aStep, G4double ScintillationTime, G4double ScintillationRiseTime) const
 
void propagationTime (std::vector< double > &arrival_time_dist, G4ThreeVector x0, const size_t OpChannel, bool Reflected=false)
 
G4double sample_time (const G4double tau1, const G4double tau2) const
 
double reemission_energy () const
 
void average_position (G4Step const &aStep, double *xzyPos) const
 
double Rectangle_SolidAngle (const double a, const double b, const double d) const
 
double Rectangle_SolidAngle (Dims const &o, const std::array< double, 3 > v) const
 
double Disk_SolidAngle (const double d, const double h, const double b) const
 
double Omega_Dome_Model (const double distance, const double theta) const
 
G4double Gaisser_Hillas (const double x, const double *par) const
 
void ProcessStep (const G4Step &step)
 
bool isOpDetInSameTPC (geo::Point_t const &ScintPoint, geo::Point_t const &OpDetPoint) const
 
bool isScintInActiveVolume (geo::Point_t const &ScintPoint)
 
double interpolate (const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
 
void interpolate3 (std::array< double, 3 > &inter, const std::vector< double > &xData, const std::vector< double > &yData1, const std::vector< double > &yData2, const std::vector< double > &yData3, double x, bool extrapolate) const
 

Static Private Member Functions

static std::vector< geo::BoxBoundedGeoextractActiveVolumes (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_tfOpDetCenter
 
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 56 of file OpFastScintillation.hh.

Constructor & Destructor Documentation

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

Definition at line 143 of file OpFastScintillation.cxx.

References util::abs(), bPropagate, BuildThePhysicsTable(), art::errors::Configuration, util::counter(), emSaturation, util::enumerate(), ExcitationRatio, extractActiveVolumes(), fActiveVolumes, fangle_bin_timing_vis, fangle_bin_timing_vuv, fborder_corr_angulo_dome, fborder_corr_angulo_flat, fborder_corr_dome, fborder_corr_flat, fcathode_centre, fcathode_plane, fcathode_ydimension, fcathode_zdimension, fcut_off_pars, fdelta_angulo_vis, fdelta_angulo_vuv, fdistances_refl, fFiniteRiseTime, fGHvuvpars_dome, fGHvuvpars_flat, finflexion_point_distance, fIsDomePDCorr, fIsFlatPDCorr, fL_abs_vuv, fmax_d, fmin_d, fOnlyActiveVolume, fOnlyOneCryostat, fOpDetCenter, fOpDetHeight, fOpDetLength, fOpDetType, fparameters, fplane_depth, fPVS, fradial_distances_refl, fradius, fstep_size, fStoreReflected, ftau_pars, fTPBEm, fTrackSecondariesFirst, fUseNhitsModel, fvis_distances_r_dome, fvis_distances_r_flat, fvis_distances_x_dome, fvis_distances_x_flat, fvis_vmean, fvispars_dome, fvispars_flat, fvuv_vgroup_max, fvuv_vgroup_mean, geo::TPCGeo::GetCathodeCenter(), geo::OpDetGeo::GetCenter(), larg4::OpFastScintillation::Dims::h, geo::OpDetGeo::Height(), phot::PhotonVisibilityService::IncludePropTime(), interpolate(), geo::OpDetGeo::isBar(), geo::OpDetGeo::Length(), phot::PhotonVisibilityService::LoadGHDome(), phot::PhotonVisibilityService::LoadGHFlat(), phot::PhotonVisibilityService::LoadTimingsForVISPar(), phot::PhotonVisibilityService::LoadTimingsForVUVPar(), phot::PhotonVisibilityService::LoadVisParsDome(), phot::PhotonVisibilityService::LoadVisParsFlat(), phot::PhotonVisibilityService::LoadVisSemiAnalyticProperties(), phot::PhotonVisibilityService::LoadVUVSemiAnalyticProperties(), geo::GeometryCore::Ncryostats(), phot::PhotonVisibilityService::NOpChannels(), geo::GeometryCore::OpDetGeoFromOpDet(), detinfo::LArProperties::ScintByParticleType(), scintillationByParticleType, geo::OpDetGeo::Shape(), phot::PhotonVisibilityService::StoreReflected(), tpbemission, geo::GeometryCore::TPC(), phot::PhotonVisibilityService::UseNhitsModel(), usesSemiAnalyticModel(), lar::dump::vector(), VUV_max, VUV_min, VUV_timing, larg4::OpFastScintillation::Dims::w, and YieldFactor.

144  : G4VRestDiscreteProcess(processName, type)
145  , fActiveVolumes{extractActiveVolumes(*(lar::providerFrom<geo::Geometry>()))}
146  , bPropagate(!(art::ServiceHandle<sim::LArG4Parameters const>()->NoPhotonPropagation()))
149  // for now, limit to the active volume only if semi-analytic model is used
151  {
152  SetProcessSubType(25); // TODO: unhardcode
153  fTrackSecondariesFirst = false;
154  fFiniteRiseTime = false;
155  YieldFactor = 1.0;
156  ExcitationRatio = 1.0;
157 
158  const detinfo::LArProperties* larp = lar::providerFrom<detinfo::LArPropertiesService>();
159 
161 
162  if (verboseLevel > 0) { G4cout << GetProcessName() << " is created " << G4endl; }
163 
165  emSaturation = NULL;
166 
167  if (bPropagate) {
168  assert(fPVS);
169 
170  // Loading the position of each optical channel, neccessary for the
171  // parametrizatiuons of Nhits and prop-time
172  geo::GeometryCore const& geom = *(lar::providerFrom<geo::Geometry>());
173  {
174  auto log = mf::LogTrace("OpFastScintillation")
175  << "OpFastScintillation: active volume boundaries from " << fActiveVolumes.size()
176  << " volumes:";
177  for (auto const& [iCryo, box] : fActiveVolumes | ranges::views::enumerate) {
178  log << "\n - C:" << iCryo << ": " << box.Min() << " -- " << box.Max() << " cm";
179  } // for
180  log << "\n (scintillation photons are propagated "
181  << (fOnlyActiveVolume ? "only from active volumes" : "from anywhere") << ")";
182  } // local scope
183 
184  if (usesSemiAnalyticModel() && (geom.Ncryostats() > 1U)) {
185  if (fOnlyOneCryostat) {
186  mf::LogWarning("OpFastScintillation")
187  << std::string(80, '=') << "\nA detector with " << geom.Ncryostats()
188  << " cryostats is configured"
189  << " , and semi-analytic model is requested for scintillation photon propagation."
190  << " THIS CONFIGURATION IS NOT SUPPORTED and it is open to bugs"
191  << " (e.g. scintillation may be detected only in cryostat #0)."
192  << "\nThis would be normally a fatal error, but it has been forcibly overridden."
193  << "\n"
194  << std::string(80, '=');
195  }
196  else {
198  << "Photon propagation via semi-analytic model is not supported yet"
199  << " on detectors with more than one cryostat.";
200  }
201  } // if
202 
203  geo::Point_t const Cathode_centre{geom.TPC().GetCathodeCenter().X(),
204  fActiveVolumes[0].CenterY(),
205  fActiveVolumes[0].CenterZ()};
206  mf::LogTrace("OpFastScintillation") << "Cathode_centre: " << Cathode_centre << " cm";
207 
208  for (size_t const i : util::counter(fPVS->NOpChannels())) {
209  geo::OpDetGeo const& opDet = geom.OpDetGeoFromOpDet(i);
210  fOpDetCenter.push_back(opDet.GetCenter());
211 
212  if (dynamic_cast<TGeoSphere const*>(opDet.Shape()) != nullptr) { // sphere/dome
213  fOpDetType.push_back(1); // Dome PMTs
214  fOpDetLength.push_back(-1);
215  fOpDetHeight.push_back(-1);
216  }
217  else if (opDet.isBar()) { // box
218  fOpDetType.push_back(0); //Arapucas
219  fOpDetLength.push_back(opDet.Length());
220  fOpDetHeight.push_back(opDet.Height());
221  }
222  else { // disk
223  fOpDetType.push_back(2); // Disk PMTs
224  fOpDetLength.push_back(-1);
225  fOpDetHeight.push_back(-1);
226  }
227  }
228 
229  if (fPVS->IncludePropTime()) {
230  std::cout << "Using parameterisation of timings." << std::endl;
231  //New VUV time parapetrization
233  fstep_size,
234  fmax_d,
235  fmin_d,
240 
241  // create vector of empty TF1s that will be replaces with the parameterisations
242  // that are generated as they are required default TF1() constructor gives
243  // function with 0 dimensions, can then check numDim to qucikly see if a
244  // parameterisation has been generated
245  const size_t num_params =
246  (fmax_d - fmin_d) /
247  fstep_size; // for d < fmin_d, no parameterisaton, a delta function is used instead
248  size_t num_angles = std::round(90 / fangle_bin_timing_vuv);
249  VUV_timing = std::vector(num_angles, std::vector(num_params, TF1()));
250 
251  // initialise vectors to contain range parameterisations sampled to in each case
252  // when using TF1->GetRandom(xmin,xmax), must be in same range otherwise sampling
253  // is regenerated, this is the slow part!
254  VUV_max = std::vector(num_angles, std::vector(num_params, 0.0));
255  VUV_min = std::vector(num_angles, std::vector(num_params, 0.0));
256 
257  // VIS time parameterisation
258  if (fPVS->StoreReflected()) {
259  // load parameters
263  ftau_pars,
264  fvis_vmean,
266  }
267  }
268  if (usesSemiAnalyticModel()) {
269  mf::LogVerbatim("OpFastScintillation")
270  << "OpFastScintillation: using semi-analytic model for number of hits";
271 
272  // LAr absorption length in cm
273  std::map<double, double> abs_length_spectrum =
274  lar::providerFrom<detinfo::LArPropertiesService>()->AbsLengthSpectrum();
275  std::vector<double> x_v, y_v;
276  for (auto elem : abs_length_spectrum) {
277  x_v.push_back(elem.first);
278  y_v.push_back(elem.second);
279  }
280  fL_abs_vuv =
281  std::round(interpolate(x_v, y_v, 9.7, false)); // 9.7 eV: peak of VUV emission spectrum
282 
283  // Load Gaisser-Hillas corrections for VUV semi-analytic hits
284  std::cout << "Loading the GH corrections" << std::endl;
287  if (!fIsFlatPDCorr && !fIsDomePDCorr) {
288  throw cet::exception("OpFastScintillation")
289  << "Both isFlatPDCorr and isDomePDCorr parameters are false, at least one type of "
290  "parameterisation is required for the semi-analytic light simulation."
291  << "\n";
292  }
293  if (fIsFlatPDCorr) {
295  }
296  if (fIsDomePDCorr) {
298  }
299  // cathode center coordinates required for corrections
301  fcathode_centre.SetY(fActiveVolumes[0].CenterY());
302  fcathode_centre.SetZ(
303  fActiveVolumes[0].CenterZ()); // to get full cathode dimension rather than just single tpc
304 
305  if (fPVS->StoreReflected()) {
306  fStoreReflected = true;
307  // Load corrections for VIS semi-anlytic hits
308  std::cout << "Loading visible light corrections" << std::endl;
310  if (fIsFlatPDCorr) {
312  }
313  if (fIsDomePDCorr) {
315  }
316 
317  // cathode dimensions
318  fcathode_ydimension = fActiveVolumes[0].SizeY();
319  fcathode_zdimension = fActiveVolumes[0].SizeZ();
320 
321  // set cathode plane struct for solid angle function
325  }
326  else
327  fStoreReflected = false;
328  }
329  }
330  tpbemission = lar::providerFrom<detinfo::LArPropertiesService>()->TpbEm();
331  std::vector<double> parent;
332  parent.reserve(tpbemission.size());
333  for (auto iter = tpbemission.begin(); iter != tpbemission.end(); ++iter) {
334  parent.push_back(iter->second);
335  }
336  fTPBEm = std::make_unique<CLHEP::RandGeneral>(parent.data(), parent.size());
337  }
void LoadTimingsForVISPar(std::vector< double > &distances, std::vector< double > &radial_distances, std::vector< std::vector< std::vector< double >>> &cut_off, std::vector< std::vector< std::vector< double >>> &tau, double &vis_vmean, double &angle_bin_timing_vis) const
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Point_t GetCathodeCenter() const
Returns the expected drift direction based on geometry.
Definition: TPCGeo.cxx:139
std::vector< std::vector< std::vector< double > > > fvispars_dome
std::vector< geo::Point_t > fOpDetCenter
void LoadTimingsForVUVPar(std::vector< std::vector< double >>(&v)[7], double &step_size, double &max_d, double &min_d, double &vuv_vgroup_mean, double &vuv_vgroup_max, double &inflexion_point_distance, double &angle_bin_timing_vuv) const
std::vector< double > fOpDetHeight
static std::vector< geo::BoxBoundedGeo > extractActiveVolumes(geo::GeometryCore const &geom)
bool isBar() const
Returns whether the detector shape is a bar (TGeoBBox).
Definition: OpDetGeo.h:167
bool const fOnlyActiveVolume
Whether photon propagation is performed only from active volumes.
Point_t const & GetCenter() const
Definition: OpDetGeo.h:71
std::map< double, double > tpbemission
std::vector< std::vector< std::vector< double > > > fvispars_flat
constexpr auto abs(T v)
Returns the absolute value of the argument.
void LoadVisParsFlat(std::vector< double > &vis_distances_x_flat, std::vector< double > &vis_distances_r_flat, std::vector< std::vector< std::vector< double >>> &vispars_flat) const
std::vector< double > fvis_distances_r_flat
std::vector< double > fdistances_refl
std::vector< std::vector< TF1 > > VUV_timing
std::vector< std::vector< double > > fparameters[7]
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:303
void LoadGHDome(std::vector< std::vector< double >> &GHvuvpars_dome, std::vector< double > &border_corr_angulo_dome, std::vector< std::vector< double >> &border_corr_dome) const
std::vector< std::vector< double > > fGHvuvpars_flat
std::vector< std::vector< double > > VUV_min
std::vector< std::vector< double > > fGHvuvpars_dome
void LoadVUVSemiAnalyticProperties(bool &isFlatPDCorr, bool &isDomePDCorr, double &delta_angulo_vuv, double &radius) const
std::vector< double > fborder_corr_angulo_dome
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
auto enumerate(Iterables &&...iterables)
Range-for loop helper tracking the number of iteration.
Definition: enumerate.h:49
bool const bPropagate
Whether propagation of photons is enabled.
std::vector< std::vector< double > > fborder_corr_dome
std::vector< std::vector< double > > fborder_corr_flat
std::vector< double > fvis_distances_x_flat
void LoadGHFlat(std::vector< std::vector< double >> &GHvuvpars_flat, std::vector< double > &border_corr_angulo_flat, std::vector< std::vector< double >> &border_corr_flat) const
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:292
bool usesSemiAnalyticModel() const
Returns whether the semi-analytic visibility parametrization is being used.
void LoadVisParsDome(std::vector< double > &vis_distances_x_dome, std::vector< double > &vis_distances_r_dome, std::vector< std::vector< std::vector< double >>> &vispars_dome) const
std::vector< double > fborder_corr_angulo_flat
std::vector< double > fOpDetLength
std::unique_ptr< CLHEP::RandGeneral > fTPBEm
std::vector< geo::BoxBoundedGeo > const fActiveVolumes
std::vector< std::vector< std::vector< double > > > ftau_pars
double Length() const
Definition: OpDetGeo.h:77
Description of the physical geometry of one entire detector.
Definition: GeometryCore.h:91
std::vector< double > fradial_distances_refl
std::vector< double > fvis_distances_r_dome
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
void LoadVisSemiAnalyticProperties(double &delta_angulo_vis, double &radius) const
std::vector< std::vector< std::vector< double > > > fcut_off_pars
bool const fUseNhitsModel
Whether the semi-analytic model is being used for photon visibility.
std::vector< double > fvis_distances_x_dome
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
phot::PhotonVisibilityService const *const fPVS
Photon visibility service instance.
std::vector< std::vector< double > > VUV_max
MaybeLogger_< ELseverityLevel::ELsev_success, true > LogTrace
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:448
virtual bool ScintByParticleType() const =0
double Height() const
Definition: OpDetGeo.h:79
TGeoShape const * Shape() const
Returns the geometry object as TGeoShape.
Definition: OpDetGeo.h:126
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
OpDetGeo const & OpDetGeoFromOpDet(unsigned int OpDet) const
Returns the geo::OpDetGeo object for the given detector number.
larg4::OpFastScintillation::~OpFastScintillation ( )

Definition at line 339 of file OpFastScintillation.cxx.

References theFastIntegralTable, and theSlowIntegralTable.

340  {
341  if (theFastIntegralTable) theFastIntegralTable->clearAndDestroy();
342  if (theSlowIntegralTable) theSlowIntegralTable->clearAndDestroy();
343  }
std::unique_ptr< G4PhysicsTable > theSlowIntegralTable
std::unique_ptr< G4PhysicsTable > theFastIntegralTable

Member Function Documentation

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

Definition at line 121 of file OpFastScintillation.hh.

121 { emSaturation = sat; }
G4VParticleChange * larg4::OpFastScintillation::AtRestDoIt ( const G4Track &  aTrack,
const G4Step &  aStep 
)
virtual

Definition at line 348 of file OpFastScintillation.cxx.

References PostStepDoIt().

351  {
352  return OpFastScintillation::PostStepDoIt(aTrack, aStep);
353  }
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
void larg4::OpFastScintillation::average_position ( G4Step const &  aStep,
double *  xzyPos 
) const
private

Definition at line 935 of file OpFastScintillation.cxx.

Referenced by RecordPhotonsProduced().

936  {
937  CLHEP::Hep3Vector prePoint = (aStep.GetPreStepPoint())->GetPosition();
938  CLHEP::Hep3Vector postPoint = (aStep.GetPostStepPoint())->GetPosition();
939  xyzPos[0] = (((prePoint.getX() + postPoint.getX()) / 2.0) / CLHEP::cm);
940  xyzPos[1] = (((prePoint.getY() + postPoint.getY()) / 2.0) / CLHEP::cm);
941  xyzPos[2] = (((prePoint.getZ() + postPoint.getZ()) / 2.0) / CLHEP::cm);
942  }
G4double larg4::OpFastScintillation::bi_exp ( const G4double  t,
const G4double  tau1,
const G4double  tau2 
) const
private

Definition at line 1544 of file OpFastScintillation.cxx.

Referenced by sample_time().

1547  { // TODO: what's up with this? ... / tau2 / tau2 ...
1548  return std::exp(-1.0 * t / tau2) * (1 - std::exp(-1.0 * t / tau1)) / tau2 / tau2 *
1549  (tau1 + tau2);
1550  }
void larg4::OpFastScintillation::BuildThePhysicsTable ( )
protected

Definition at line 714 of file OpFastScintillation.cxx.

References theFastIntegralTable, and theSlowIntegralTable.

Referenced by OpFastScintillation().

715  {
717 
718  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
719  G4int numOfMaterials = G4Material::GetNumberOfMaterials();
720 
721  // create new physics table
723  theFastIntegralTable = std::make_unique<G4PhysicsTable>(numOfMaterials);
725  theSlowIntegralTable = std::make_unique<G4PhysicsTable>(numOfMaterials);
726 
727  // loop for materials
728  for (G4int i = 0; i < numOfMaterials; i++) {
729  G4PhysicsOrderedFreeVector* aPhysicsOrderedFreeVector = new G4PhysicsOrderedFreeVector();
730  G4PhysicsOrderedFreeVector* bPhysicsOrderedFreeVector = new G4PhysicsOrderedFreeVector();
731 
732  // Retrieve vector of scintillation wavelength intensity for the material from the
733  // material's optical properties table.
734  G4Material* aMaterial = (*theMaterialTable)[i];
735 
736  G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
737 
738  if (aMaterialPropertiesTable) {
739 
740  G4MaterialPropertyVector* theFastLightVector =
741  aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
742 
743  if (theFastLightVector) {
744  // Retrieve the first intensity point in vector of (photon energy, intensity)
745  // pairs
746  G4double currentIN = (*theFastLightVector)[0];
747  if (currentIN >= 0.0) {
748  // Create first (photon energy, Scintillation Integral) pair
749  G4double currentPM = theFastLightVector->Energy(0);
750  G4double currentCII = 0.0;
751 
752  aPhysicsOrderedFreeVector->InsertValues(currentPM, currentCII);
753 
754  // Set previous values to current ones prior to loop
755  G4double prevPM = currentPM;
756  G4double prevCII = currentCII;
757  G4double prevIN = currentIN;
758 
759  // loop over all (photon energy, intensity) pairs stored for this material
760  for (size_t i = 1; i < theFastLightVector->GetVectorLength(); i++) {
761  currentPM = theFastLightVector->Energy(i);
762  currentIN = (*theFastLightVector)[i];
763 
764  currentCII = 0.5 * (prevIN + currentIN);
765 
766  currentCII = prevCII + (currentPM - prevPM) * currentCII;
767 
768  aPhysicsOrderedFreeVector->InsertValues(currentPM, currentCII);
769 
770  prevPM = currentPM;
771  prevCII = currentCII;
772  prevIN = currentIN;
773  }
774  }
775  }
776 
777  G4MaterialPropertyVector* theSlowLightVector =
778  aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
779  if (theSlowLightVector) {
780  // Retrieve the first intensity point in vector of (photon energy, intensity)
781  // pairs
782  G4double currentIN = (*theSlowLightVector)[0];
783  if (currentIN >= 0.0) {
784  // Create first (photon energy, Scintillation Integral) pair
785  G4double currentPM = theSlowLightVector->Energy(0);
786  G4double currentCII = 0.0;
787 
788  bPhysicsOrderedFreeVector->InsertValues(currentPM, currentCII);
789 
790  // Set previous values to current ones prior to loop
791  G4double prevPM = currentPM;
792  G4double prevCII = currentCII;
793  G4double prevIN = currentIN;
794 
795  // loop over all (photon energy, intensity) pairs stored for this material
796  for (size_t i = 1; i < theSlowLightVector->GetVectorLength(); i++) {
797  currentPM = theSlowLightVector->Energy(i);
798  currentIN = (*theSlowLightVector)[i];
799 
800  currentCII = 0.5 * (prevIN + currentIN);
801 
802  currentCII = prevCII + (currentPM - prevPM) * currentCII;
803 
804  bPhysicsOrderedFreeVector->InsertValues(currentPM, currentCII);
805 
806  prevPM = currentPM;
807  prevCII = currentCII;
808  prevIN = currentIN;
809  }
810  }
811  }
812  }
813  // The scintillation integral(s) for a given material will be inserted in the
814  // table(s) according to the position of the material in the material table.
815  theFastIntegralTable->insertAt(i, aPhysicsOrderedFreeVector);
816  theSlowIntegralTable->insertAt(i, bPhysicsOrderedFreeVector);
817  }
818  }
std::unique_ptr< G4PhysicsTable > theSlowIntegralTable
std::unique_ptr< G4PhysicsTable > theFastIntegralTable
void larg4::OpFastScintillation::detectedDirectHits ( std::map< size_t, int > &  DetectedNum,
const double  Num,
geo::Point_t const &  ScintPoint 
) const

Definition at line 1194 of file OpFastScintillation.cxx.

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

Referenced by RecordPhotonsProduced().

1197  {
1198  for (size_t const OpDet : util::counter(fPVS->NOpChannels())) {
1199  if (!isOpDetInSameTPC(ScintPoint, fOpDetCenter.at(OpDet))) continue;
1200 
1201  // set detector struct for solid angle function
1202  const OpFastScintillation::OpticalDetector op{fOpDetHeight.at(OpDet),
1203  fOpDetLength.at(OpDet),
1204  fOpDetCenter.at(OpDet),
1205  fOpDetType.at(OpDet)};
1206  const int DetThis = VUVHits(Num, ScintPoint, op);
1207  if (DetThis > 0) { DetectedNum[OpDet] = DetThis; }
1208  }
1209  }
std::vector< geo::Point_t > fOpDetCenter
std::vector< double > fOpDetHeight
bool isOpDetInSameTPC(geo::Point_t const &ScintPoint, geo::Point_t const &OpDetPoint) const
int VUVHits(const double Nphotons_created, geo::Point_t const &ScintPoint, OpticalDetector const &opDet) const
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:292
std::vector< double > fOpDetLength
phot::PhotonVisibilityService const *const fPVS
Photon visibility service instance.
void larg4::OpFastScintillation::detectedReflecHits ( std::map< size_t, int > &  ReflDetectedNum,
const double  Num,
geo::Point_t const &  ScintPoint 
) const

Definition at line 1211 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().

1214  {
1215  // 1). calculate total number of hits of VUV photons on reflective foils via solid
1216  // angle + Gaisser-Hillas corrections:
1217 
1218  // set plane_depth for correct TPC:
1219  double plane_depth;
1220  if (ScintPoint.X() < 0.) { plane_depth = -fplane_depth; }
1221  else {
1222  plane_depth = fplane_depth;
1223  }
1224 
1225  // get scintpoint coords relative to centre of cathode plane
1226  std::array<double, 3> ScintPoint_relative = {std::abs(ScintPoint.X() - plane_depth),
1227  std::abs(ScintPoint.Y() - fcathode_centre.Y()),
1228  std::abs(ScintPoint.Z() - fcathode_centre.Z())};
1229  // calculate solid angle of cathode from the scintillation point
1230  double solid_angle_cathode = Rectangle_SolidAngle(fcathode_plane, ScintPoint_relative);
1231 
1232  // calculate distance and angle between ScintPoint and hotspot; vast majority of hits
1233  // in hotspot region directly infront of scintpoint, therefore consider attenuation
1234  // for this distance and on axis GH instead of for the centre coordinate
1235  double distance_cathode = std::abs(plane_depth - ScintPoint.X());
1236  // calculate hits on cathode plane via geometric acceptance
1237  double cathode_hits_geo = std::exp(-1. * distance_cathode / fL_abs_vuv) *
1238  (solid_angle_cathode / (4. * CLHEP::pi)) * Num;
1239  // determine Gaisser-Hillas correction including border effects use flat correction
1240  double r = std::sqrt(std::pow(ScintPoint.Y() - fcathode_centre.Y(), 2) +
1241  std::pow(ScintPoint.Z() - fcathode_centre.Z(), 2));
1242  double pars_ini[4] = {0, 0, 0, 0};
1243  double s1 = 0;
1244  double s2 = 0;
1245  double s3 = 0;
1246  if (fIsFlatPDCorr) {
1247  pars_ini[0] = fGHvuvpars_flat[0][0];
1248  pars_ini[1] = fGHvuvpars_flat[1][0];
1249  pars_ini[2] = fGHvuvpars_flat[2][0];
1250  pars_ini[3] = fGHvuvpars_flat[3][0];
1254  }
1255  else
1256  std::cout
1257  << "Error: flat optical detector VUV correction required for reflected semi-analytic hits."
1258  << std::endl;
1259 
1260  // add border correction
1261  pars_ini[0] = pars_ini[0] + s1 * r;
1262  pars_ini[1] = pars_ini[1] + s2 * r;
1263  pars_ini[2] = pars_ini[2] + s3 * r;
1264  pars_ini[3] = pars_ini[3];
1265 
1266  // calculate corrected number of hits
1267  double GH_correction = Gaisser_Hillas(distance_cathode, pars_ini);
1268  const double cathode_hits_rec = GH_correction * cathode_hits_geo;
1269 
1270  // detemine hits on each PD
1271  const std::array<double, 3> hotspot = {plane_depth, ScintPoint.Y(), ScintPoint.Z()};
1272  for (size_t const OpDet : util::counter(fPVS->NOpChannels())) {
1273  if (!isOpDetInSameTPC(ScintPoint, fOpDetCenter.at(OpDet))) continue;
1274 
1275  // set detector struct for solid angle function
1276  const OpFastScintillation::OpticalDetector op{fOpDetHeight.at(OpDet),
1277  fOpDetLength.at(OpDet),
1278  fOpDetCenter.at(OpDet),
1279  fOpDetType.at(OpDet)};
1280 
1281  int const ReflDetThis = VISHits(ScintPoint, op, cathode_hits_rec, hotspot);
1282  if (ReflDetThis > 0) { ReflDetectedNum[OpDet] = ReflDetThis; }
1283  }
1284  }
TRandom r
Definition: spectrum.C:23
std::vector< geo::Point_t > fOpDetCenter
G4double Gaisser_Hillas(const double x, const double *par) const
std::vector< double > fOpDetHeight
int VISHits(geo::Point_t const &ScintPoint, OpticalDetector const &opDet, const double cathode_hits_rec, const std::array< double, 3 > hotspot) const
constexpr auto abs(T v)
Returns the absolute value of the argument.
bool isOpDetInSameTPC(geo::Point_t const &ScintPoint, geo::Point_t const &OpDetPoint) const
std::vector< std::vector< double > > fGHvuvpars_flat
std::vector< std::vector< double > > fborder_corr_flat
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:292
std::vector< double > fborder_corr_angulo_flat
std::vector< double > fOpDetLength
double Rectangle_SolidAngle(const double a, const double b, const double d) const
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
phot::PhotonVisibilityService const *const fPVS
Photon visibility service instance.
double larg4::OpFastScintillation::Disk_SolidAngle ( const double  d,
const double  h,
const double  b 
) const
private

Definition at line 1643 of file OpFastScintillation.cxx.

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

Referenced by VISHits(), and VUVHits().

1644  {
1645  if (b <= 0. || d < 0. || h <= 0.) return 0.;
1646  const double leg2 = (b + d) * (b + d);
1647  const double aa = std::sqrt(h * h / (h * h + leg2));
1648  if (isApproximatelyZero(d)) { return 2. * CLHEP::pi * (1. - aa); }
1649  double bb = 2. * std::sqrt(b * d / (h * h + leg2));
1650  double cc = 4. * b * d / leg2;
1651 
1652  if (isDefinitelyGreaterThan(d, b)) {
1653  try {
1654  return 2. * aa *
1655  (std::sqrt(1. - cc) * boost::math::ellint_3(bb, cc, noLDoublePromote()) -
1656  boost::math::ellint_1(bb, noLDoublePromote()));
1657  }
1658  catch (std::domain_error& e) {
1659  if (isApproximatelyEqual(d, b, 1e-9)) {
1660  mf::LogWarning("OpFastScintillation")
1661  << "Elliptic Integral in Disk_SolidAngle() given parameters outside domain."
1662  << "\nbb: " << bb << "\ncc: " << cc << "\nException message: " << e.what()
1663  << "\nRelax condition and carry on.";
1664  return CLHEP::pi - 2. * aa * boost::math::ellint_1(bb, noLDoublePromote());
1665  }
1666  else {
1667  mf::LogError("OpFastScintillation")
1668  << "Elliptic Integral inside Disk_SolidAngle() given parameters outside domain.\n"
1669  << "\nbb: " << bb << "\ncc: " << cc << "Exception message: " << e.what();
1670  return 0.;
1671  }
1672  }
1673  }
1674  if (isDefinitelyLessThan(d, b)) {
1675  try {
1676  return 2. * CLHEP::pi -
1677  2. * aa *
1678  (boost::math::ellint_1(bb, noLDoublePromote()) +
1679  std::sqrt(1. - cc) * boost::math::ellint_3(bb, cc, noLDoublePromote()));
1680  }
1681  catch (std::domain_error& e) {
1682  if (isApproximatelyEqual(d, b, 1e-9)) {
1683  mf::LogWarning("OpFastScintillation")
1684  << "Elliptic Integral in Disk_SolidAngle() given parameters outside domain."
1685  << "\nbb: " << bb << "\ncc: " << cc << "\nException message: " << e.what()
1686  << "\nRelax condition and carry on.";
1687  return CLHEP::pi - 2. * aa * boost::math::ellint_1(bb, noLDoublePromote());
1688  }
1689  else {
1690  mf::LogError("OpFastScintillation")
1691  << "Elliptic Integral inside Disk_SolidAngle() given parameters outside domain.\n"
1692  << "\nbb: " << bb << "\ncc: " << cc << "Exception message: " << e.what();
1693  return 0.;
1694  }
1695  }
1696  }
1697  if (isApproximatelyEqual(d, b)) {
1698  return CLHEP::pi - 2. * aa * boost::math::ellint_1(bb, noLDoublePromote());
1699  }
1700  return 0.;
1701  }
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Float_t d
Definition: plot.C:235
static constexpr bool isDefinitelyGreaterThan(TReal a, TReal b, TReal tolerance=std::numeric_limits< TReal >::epsilon())
static constexpr bool isApproximatelyEqual(TReal a, TReal b, TReal tolerance=std::numeric_limits< TReal >::epsilon())
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
boost::math::policies::policy< boost::math::policies::promote_double< false > > noLDoublePromote
static constexpr bool isDefinitelyLessThan(TReal a, TReal b, TReal tolerance=std::numeric_limits< TReal >::epsilon())
Float_t e
Definition: plot.C:35
static constexpr bool isApproximatelyZero(TReal a, TReal tolerance=std::numeric_limits< TReal >::epsilon())
void larg4::OpFastScintillation::DumpPhysicsTable ( ) const
inline

Definition at line 393 of file OpFastScintillation.hh.

394  {
395  if (theFastIntegralTable) {
396  G4int PhysicsTableSize = theFastIntegralTable->entries();
397  G4PhysicsOrderedFreeVector* v;
398  for (G4int i = 0; i < PhysicsTableSize; i++) {
399  v = (G4PhysicsOrderedFreeVector*)(*theFastIntegralTable)[i];
400  v->DumpValues();
401  }
402  }
403  if (theSlowIntegralTable) {
404  G4int PhysicsTableSize = theSlowIntegralTable->entries();
405  G4PhysicsOrderedFreeVector* v;
406  for (G4int i = 0; i < PhysicsTableSize; i++) {
407  v = (G4PhysicsOrderedFreeVector*)(*theSlowIntegralTable)[i];
408  v->DumpValues();
409  }
410  }
411  }
std::unique_ptr< G4PhysicsTable > theSlowIntegralTable
std::unique_ptr< G4PhysicsTable > theFastIntegralTable
std::vector< geo::BoxBoundedGeo > larg4::OpFastScintillation::extractActiveVolumes ( geo::GeometryCore const &  geom)
staticprivate

Definition at line 1799 of file OpFastScintillation.cxx.

References geo::BoxBoundedGeo::ExtendToInclude(), geo::Iterable< IterationPolicy, Transform >::Iterate(), and geo::GeometryCore::Ncryostats().

Referenced by OpFastScintillation().

1801  {
1802  std::vector<geo::BoxBoundedGeo> activeVolumes;
1803  activeVolumes.reserve(geom.Ncryostats());
1804 
1805  for (geo::CryostatGeo const& cryo : geom.Iterate<geo::CryostatGeo>()) {
1806 
1807  // can't use it default-constructed since it would always include origin
1808  geo::BoxBoundedGeo box{cryo.TPC(0).ActiveBoundingBox()};
1809 
1810  for (geo::TPCGeo const& TPC : cryo.IterateTPCs())
1811  box.ExtendToInclude(TPC.ActiveBoundingBox());
1812 
1813  activeVolumes.push_back(std::move(box));
1814 
1815  } // for cryostats
1816 
1817  return activeVolumes;
1818  }
Geometry information for a single TPC.
Definition: TPCGeo.h:33
Geometry information for a single cryostat.
Definition: CryostatGeo.h:42
A base class aware of world box coordinatesAn object describing a simple shape can inherit from this ...
Definition: BoxBoundedGeo.h:31
void ExtendToInclude(Coord_t x, Coord_t y, Coord_t z)
Extends the current box to also include the specified point.
G4double larg4::OpFastScintillation::Gaisser_Hillas ( const double  x,
const double *  par 
) const
private

Definition at line 1552 of file OpFastScintillation.cxx.

Referenced by detectedReflecHits(), and VUVHits().

1553  {
1554  double X_mu_0 = par[3];
1555  double Normalization = par[0];
1556  double Diff = par[1] - X_mu_0;
1557  double Term = std::pow((x - X_mu_0) / Diff, Diff / par[2]);
1558  double Exponential = std::exp((par[1] - x) / par[2]);
1559  return (Normalization * Term * Exponential);
1560  }
Float_t x
Definition: compare.C:6
void larg4::OpFastScintillation::generateParam ( const size_t  index,
const size_t  angle_bin 
)

Definition at line 945 of file OpFastScintillation.cxx.

References finflexion_point_distance, fmin_d, fparameters, fstep_size, fvuv_vgroup_max, fvuv_vgroup_mean, interpolate(), interpolate3(), VUV_max, VUV_min, and VUV_timing.

Referenced by getVUVTimes().

946  {
947  // get distance
948  double distance_in_cm = (index * fstep_size) + fmin_d;
949 
950  // time range
951  const double signal_t_range = 5000.; // TODO: unhardcode
952 
953  // parameterisation TF1
954  TF1 fVUVTiming;
955 
956  // For very short distances the time correction is just a shift
957  double t_direct_mean = distance_in_cm / fvuv_vgroup_mean;
958  double t_direct_min = distance_in_cm / fvuv_vgroup_max;
959 
960  // Defining the model function(s) describing the photon transportation timing vs
961  // distance Getting the landau parameters from the time parametrization
962  std::array<double, 3> pars_landau;
963  interpolate3(pars_landau,
964  fparameters[0][0],
965  fparameters[2][angle_bin],
966  fparameters[3][angle_bin],
967  fparameters[1][angle_bin],
968  distance_in_cm,
969  true);
970  // Deciding which time model to use (depends on the distance) defining useful times
971  // for the VUV arrival time shapes
972  if (distance_in_cm >= finflexion_point_distance) {
973  double pars_far[4] = {t_direct_min, pars_landau[0], pars_landau[1], pars_landau[2]};
974  // Set model: Landau
975  fVUVTiming = TF1("fVUVTiming", model_far, 0, signal_t_range, 4);
976  fVUVTiming.SetParameters(pars_far);
977  }
978  else {
979  // Set model: Landau + Exponential
980  fVUVTiming = TF1("fVUVTiming", model_close, 0, signal_t_range, 7);
981  // Exponential parameters
982  double pars_expo[2];
983  // Getting the exponential parameters from the time parametrization
984  pars_expo[1] =
985  interpolate(fparameters[4][0], fparameters[5][angle_bin], distance_in_cm, true);
986  pars_expo[0] =
987  interpolate(fparameters[4][0], fparameters[6][angle_bin], distance_in_cm, true);
988  pars_expo[0] *= pars_landau[2];
989  pars_expo[0] = std::log(pars_expo[0]);
990  // this is to find the intersection point between the two functions:
991  TF1 fint = TF1("fint", finter_d, pars_landau[0], 4 * t_direct_mean, 5);
992  double parsInt[5] = {
993  pars_landau[0], pars_landau[1], pars_landau[2], pars_expo[0], pars_expo[1]};
994  fint.SetParameters(parsInt);
995  double t_int = fint.GetMinimumX();
996  double minVal = fint.Eval(t_int);
997  // the functions must intersect - output warning if they don't
998  if (minVal > 0.015) {
999  std::cout << "WARNING: Parametrization of VUV light discontinuous for distance = "
1000  << distance_in_cm << std::endl;
1001  std::cout << "WARNING: This shouldn't be happening " << std::endl;
1002  }
1003  double parsfinal[7] = {t_int,
1004  pars_landau[0],
1005  pars_landau[1],
1006  pars_landau[2],
1007  pars_expo[0],
1008  pars_expo[1],
1009  t_direct_min};
1010  fVUVTiming.SetParameters(parsfinal);
1011  }
1012 
1013  // set the number of points used to sample parameterisation for shorter distances,
1014  // peak is sharper so more sensitive sampling required
1015  int fsampling; // TODO: unhardcode
1016  if (distance_in_cm < 50) { fsampling = 10000; }
1017  else if (distance_in_cm < 100) {
1018  fsampling = 5000;
1019  }
1020  else {
1021  fsampling = 1000;
1022  }
1023  fVUVTiming.SetNpx(fsampling);
1024 
1025  // calculate max and min distance relevant to sample parameterisation max
1026  const size_t nq_max = 1;
1027  double xq_max[nq_max];
1028  double yq_max[nq_max];
1029  xq_max[0] = 0.995; // include 99.5%
1030  fVUVTiming.GetQuantiles(nq_max, yq_max, xq_max);
1031  double max = yq_max[0];
1032  // min
1033  double min = t_direct_min;
1034 
1035  // store TF1 and min/max, this allows identical TF1 to be used every time sampling the
1036  // first call of GetRandom generates the timing sampling and stores it in the TF1
1037  // object, this is the slow part all subsequent calls check if it has been generated
1038  // previously and are ~100+ times quicker
1039  VUV_timing[angle_bin][index] = fVUVTiming;
1040  VUV_max[angle_bin][index] = max;
1041  VUV_min[angle_bin][index] = min;
1042  }
std::vector< std::vector< TF1 > > VUV_timing
std::vector< std::vector< double > > fparameters[7]
std::vector< std::vector< double > > VUV_min
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
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
std::vector< std::vector< double > > VUV_max
G4PhysicsTable * larg4::OpFastScintillation::GetFastIntegralTable ( ) const
inline

Definition at line 388 of file OpFastScintillation.hh.

389  {
390  return theFastIntegralTable.get();
391  }
std::unique_ptr< G4PhysicsTable > theFastIntegralTable
G4bool larg4::OpFastScintillation::GetFiniteRiseTime ( ) const
inline

Definition at line 358 of file OpFastScintillation.hh.

359  {
360  return fFiniteRiseTime;
361  }
G4double larg4::OpFastScintillation::GetMeanFreePath ( const G4Track &  aTrack,
G4double  ,
G4ForceCondition *  condition 
)

Definition at line 834 of file OpFastScintillation.cxx.

837  {
838  *condition = StronglyForced;
839  return DBL_MAX;
840  }
G4double larg4::OpFastScintillation::GetMeanLifeTime ( const G4Track &  aTrack,
G4ForceCondition *  condition 
)

Definition at line 842 of file OpFastScintillation.cxx.

843  {
844  *condition = Forced;
845  return DBL_MAX;
846  }
G4EmSaturation* larg4::OpFastScintillation::GetSaturation ( ) const
inline

Definition at line 127 of file OpFastScintillation.hh.

127 { return emSaturation; }
G4bool larg4::OpFastScintillation::GetScintillationByParticleType ( ) const
inline

Definition at line 134 of file OpFastScintillation.hh.

G4double larg4::OpFastScintillation::GetScintillationExcitationRatio ( ) const
inline

Definition at line 378 of file OpFastScintillation.hh.

379  {
380  return ExcitationRatio;
381  }
G4double larg4::OpFastScintillation::GetScintillationYieldFactor ( ) const
inline

Definition at line 368 of file OpFastScintillation.hh.

369  {
370  return YieldFactor;
371  }
G4PhysicsTable * larg4::OpFastScintillation::GetSlowIntegralTable ( ) const
inline

Definition at line 383 of file OpFastScintillation.hh.

384  {
385  return theSlowIntegralTable.get();
386  }
std::unique_ptr< G4PhysicsTable > theSlowIntegralTable
G4bool larg4::OpFastScintillation::GetTrackSecondariesFirst ( ) const
inline

Definition at line 353 of file OpFastScintillation.hh.

354  {
355  return fTrackSecondariesFirst;
356  }
void larg4::OpFastScintillation::getVISTimes ( std::vector< double > &  arrivalTimes,
const TVector3 &  ScintPoint,
const TVector3 &  OpDetPoint 
)

Definition at line 1070 of file OpFastScintillation.cxx.

References util::abs(), util::counter(), fangle_bin_timing_vis, phot::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().

1073  {
1074  // *************************************************************************************************
1075  // Calculation of earliest arrival times and corresponding unsmeared distribution
1076  // *************************************************************************************************
1077 
1078  // set plane_depth for correct TPC:
1079  double plane_depth;
1080  if (ScintPoint[0] < 0) { plane_depth = -fplane_depth; }
1081  else {
1082  plane_depth = fplane_depth;
1083  }
1084 
1085  // calculate point of reflection for shortest path
1086  TVector3 bounce_point(plane_depth, ScintPoint[1], ScintPoint[2]);
1087 
1088  // calculate distance travelled by VUV light and by vis light
1089  double VUVdist = (bounce_point - ScintPoint).Mag();
1090  double Visdist = (OpDetPoint - bounce_point).Mag();
1091 
1092  // calculate times taken by VUV part of path
1093  int angle_bin_vuv = 0; // on-axis by definition
1094  getVUVTimes(arrivalTimes, VUVdist, angle_bin_vuv);
1095 
1096  // add visible direct path transport time
1097  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
1098  arrivalTimes[i] += Visdist / fvis_vmean;
1099  }
1100 
1101  // *************************************************************************************************
1102  // Smearing of arrival time distribution
1103  // *************************************************************************************************
1104  // calculate fastest time possible
1105  // vis part
1106  double vis_time = Visdist / fvis_vmean;
1107  // vuv part
1108  double vuv_time;
1109  if (VUVdist < fmin_d) { vuv_time = VUVdist / fvuv_vgroup_max; }
1110  else {
1111  // find index of required parameterisation
1112  const size_t index = std::round((VUVdist - fmin_d) / fstep_size);
1113  // find shortest time
1114  vuv_time = VUV_min[angle_bin_vuv][index];
1115  }
1116  // sum
1117  double fastest_time = vis_time + vuv_time;
1118 
1119  // calculate angle theta between bound_point and optical detector
1120  double cosine_theta = std::abs(OpDetPoint[0] - bounce_point[0]) / Visdist;
1121  double theta = fast_acos(cosine_theta) * 180. / CLHEP::pi;
1122 
1123  // determine smearing parameters using interpolation of generated points:
1124  // 1). tau = exponential smearing factor, varies with distance and angle
1125  // 2). cutoff = largest smeared time allowed, preventing excessively large times
1126  // caused by exponential distance to cathode
1127  double distance_cathode_plane = std::abs(plane_depth - ScintPoint[0]);
1128  // angular bin
1129  size_t theta_bin = theta / fangle_bin_timing_vis;
1130  // radial distance from centre of TPC (y,z plane)
1131  double r = std::sqrt(std::pow(ScintPoint.Y() - fcathode_centre.Z(), 2) +
1132  std::pow(ScintPoint.Y() - fcathode_centre.Z(), 2));
1133 
1134  // cut-off and tau
1135  // cut-off
1136  // interpolate in d_c for each r bin
1137  std::vector<double> interp_vals(fcut_off_pars[theta_bin].size(), 0.0);
1138  for (size_t i = 0; i < fcut_off_pars[theta_bin].size(); i++) {
1139  interp_vals[i] =
1140  interpolate(fdistances_refl, fcut_off_pars[theta_bin][i], distance_cathode_plane, true);
1141  }
1142  // interpolate in r
1143  double cutoff = interpolate(fradial_distances_refl, interp_vals, r, true);
1144 
1145  // tau
1146  // interpolate in x for each r bin
1147  std::vector<double> interp_vals_tau(ftau_pars[theta_bin].size(), 0.0);
1148  for (size_t i = 0; i < ftau_pars[theta_bin].size(); i++) {
1149  interp_vals_tau[i] =
1150  interpolate(fdistances_refl, ftau_pars[theta_bin][i], distance_cathode_plane, true);
1151  }
1152  // interpolate in r
1153  double tau = interpolate(fradial_distances_refl, interp_vals_tau, r, true);
1154 
1155  if (tau < 0) { tau = 0; } // failsafe if tau extrapolate goes wrong
1156 
1157  // apply smearing:
1158  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
1159  double arrival_time_smeared;
1160  // if time is already greater than cutoff, do not apply smearing
1161  if (arrivalTimes[i] >= cutoff) { continue; }
1162  // otherwise smear
1163  else {
1164  unsigned int counter = 0;
1165  // loop until time generated is within cutoff limit; most are within single
1166  // attempt, very few take more than two
1167  do {
1168  // don't attempt smearings too many times
1169  if (counter >= 10) { // TODO: unhardcode
1170  arrival_time_smeared = arrivalTimes[i]; // don't smear
1171  break;
1172  }
1173  else {
1174  // generate random number in appropriate range
1175  double x = gRandom->Uniform(0.5, 1.0); // TODO: unhardcode
1176  // apply the exponential smearing
1177  arrival_time_smeared =
1178  arrivalTimes[i] + (arrivalTimes[i] - fastest_time) * (std::pow(x, -tau) - 1);
1179  }
1180  counter++;
1181  } while (arrival_time_smeared > cutoff);
1182  }
1183  arrivalTimes[i] = arrival_time_smeared;
1184  }
1185  }
Float_t x
Definition: compare.C:6
TRandom r
Definition: spectrum.C:23
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::vector< double > fdistances_refl
std::vector< std::vector< double > > VUV_min
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:292
double fast_acos(double xin)
std::vector< std::vector< std::vector< double > > > ftau_pars
std::vector< double > fradial_distances_refl
std::vector< std::vector< std::vector< double > > > fcut_off_pars
void getVUVTimes(std::vector< double > &arrivalTimes, const double distance_in_cm, const size_t angle_bin)
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
void larg4::OpFastScintillation::getVUVTimes ( std::vector< double > &  arrivalTimes,
const double  distance_in_cm,
const size_t  angle_bin 
)

Definition at line 1045 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().

1048  {
1049  if (distance < fmin_d) {
1050  // times are fixed shift i.e. direct path only
1051  double t_prop_correction = distance / fvuv_vgroup_mean;
1052  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
1053  arrivalTimes[i] = t_prop_correction;
1054  }
1055  }
1056  else { // distance >= fmin_d
1057  // determine nearest parameterisation in discretisation
1058  int index = std::round((distance - fmin_d) / fstep_size);
1059  // check whether required parameterisation has been generated, generating if not
1060  if (VUV_timing[angle_bin][index].GetNdim() == 0) { generateParam(index, angle_bin); }
1061  // randomly sample parameterisation for each photon
1062  for (size_t i = 0; i < arrivalTimes.size(); ++i) {
1063  arrivalTimes[i] = VUV_timing[angle_bin][index].GetRandom(VUV_min[angle_bin][index],
1064  VUV_max[angle_bin][index]);
1065  }
1066  }
1067  }
std::vector< std::vector< TF1 > > VUV_timing
std::vector< std::vector< double > > VUV_min
void generateParam(const size_t index, const size_t angle_bin)
std::vector< std::vector< double > > VUV_max
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 1565 of file OpFastScintillation.cxx.

References util::size().

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

1570  {
1571  if (i == 0) {
1572  size_t size = xData.size();
1573  if (x >= xData[size - 2]) { // special case: beyond right end
1574  i = size - 2;
1575  }
1576  else {
1577  while (x > xData[i + 1])
1578  i++;
1579  }
1580  }
1581  double xL = xData[i];
1582  double xR = xData[i + 1];
1583  double yL = yData[i];
1584  double yR = yData[i + 1]; // points on either side (unless beyond ends)
1585  if (!extrapolate) { // if beyond ends of array and not extrapolating
1586  if (x < xL) return yL;
1587  if (x > xR) return yL;
1588  }
1589  const double dydx = (yR - yL) / (xR - xL); // gradient
1590  return yL + dydx * (x - xL); // linear interpolation
1591  }
Float_t x
Definition: compare.C:6
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
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 1593 of file OpFastScintillation.cxx.

References util::size().

Referenced by generateParam().

1600  {
1601  size_t size = xData.size();
1602  size_t i = 0; // find left end of interval for interpolation
1603  if (x >= xData[size - 2]) { // special case: beyond right end
1604  i = size - 2;
1605  }
1606  else {
1607  while (x > xData[i + 1])
1608  i++;
1609  }
1610  double xL = xData[i];
1611  double xR = xData[i + 1]; // points on either side (unless beyond ends)
1612  double yL1 = yData1[i];
1613  double yR1 = yData1[i + 1];
1614  double yL2 = yData2[i];
1615  double yR2 = yData2[i + 1];
1616  double yL3 = yData3[i];
1617  double yR3 = yData3[i + 1];
1618 
1619  if (!extrapolate) { // if beyond ends of array and not extrapolating
1620  if (x < xL) {
1621  inter[0] = yL1;
1622  inter[1] = yL2;
1623  inter[2] = yL3;
1624  return;
1625  }
1626  if (x > xR) {
1627  inter[0] = yL1;
1628  inter[1] = yL2;
1629  inter[2] = yL3;
1630  return;
1631  }
1632  }
1633  const double m = (x - xL) / (xR - xL);
1634  inter[0] = m * (yR1 - yL1) + yL1;
1635  inter[1] = m * (yR2 - yL2) + yL2;
1636  inter[2] = m * (yR3 - yL3) + yL3;
1637  }
Float_t x
Definition: compare.C:6
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
G4bool larg4::OpFastScintillation::IsApplicable ( const G4ParticleDefinition &  aParticleType)
inlinevirtual

Definition at line 335 of file OpFastScintillation.hh.

336  {
337  if (aParticleType.GetParticleName() == "opticalphoton") return false;
338  if (aParticleType.IsShortLived()) return false;
339 
340  return true;
341  }
bool larg4::OpFastScintillation::isOpDetInSameTPC ( geo::Point_t const &  ScintPoint,
geo::Point_t const &  OpDetPoint 
) const
private

Definition at line 1519 of file OpFastScintillation.cxx.

References util::abs().

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

1521  {
1522  // check optical channel is in same TPC as scintillation light, if not return 0 hits;
1523  // temporary method working for SBND, uBooNE, DUNE 1x2x6; to be replaced to work in
1524  // full DUNE geometry; check x coordinate has same sign or is close to zero, otherwise
1525  // return 0 hits
1526  if (((ScintPoint.X() < 0.) != (OpDetPoint.X() < 0.)) &&
1527  std::abs(OpDetPoint.X()) > 10.) { // TODO: unhardcode
1528  return false;
1529  }
1530  return true;
1531  }
constexpr auto abs(T v)
Returns the absolute value of the argument.
bool larg4::OpFastScintillation::isScintInActiveVolume ( geo::Point_t const &  ScintPoint)
private

Definition at line 1533 of file OpFastScintillation.cxx.

References fActiveVolumes.

Referenced by RecordPhotonsProduced().

1534  {
1535  //semi-analytic approach only works in the active volume
1536  return fActiveVolumes[0].ContainsPosition(ScintPoint);
1537  }
std::vector< geo::BoxBoundedGeo > const fActiveVolumes
double larg4::OpFastScintillation::Omega_Dome_Model ( const double  distance,
const double  theta 
) const
private

Definition at line 1769 of file OpFastScintillation.cxx.

References fradius, and util::pi().

Referenced by VISHits(), and VUVHits().

1770  {
1771  // this function calculates the solid angle of a semi-sphere of radius b, as a
1772  // correction to the analytic formula of the on-axix solid angle, as we move off-axis
1773  // an angle theta. We have used 9-angular bins with delta_theta width.
1774 
1775  // par0 = Radius correction close
1776  // par1 = Radius correction far
1777  // par2 = breaking distance betwween "close" and "far"
1778 
1779  double par0[9] = {0., 0., 0., 0., 0., 0.597542, 1.00872, 1.46993, 2.04221};
1780  double par1[9] = {0, 0, 0.19569, 0.300449, 0.555598, 0.854939, 1.39166, 2.19141, 2.57732};
1781  const double delta_theta = 10.;
1782  int j = int(theta / delta_theta);
1783  // PMT radius
1784  const double b = fradius; // cm
1785  // distance form which the model parameters break (empirical value)
1786  const double d_break = 5 * b; //par2
1787 
1788  if (distance >= d_break) {
1789  double R_apparent_far = b - par1[j];
1790  return (2 * CLHEP::pi * (1 - std::sqrt(1 - std::pow(R_apparent_far / distance, 2))));
1791  }
1792  else {
1793  double R_apparent_close = b - par0[j];
1794  return (2 * CLHEP::pi * (1 - std::sqrt(1 - std::pow(R_apparent_close / distance, 2))));
1795  }
1796  }
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
G4VParticleChange * larg4::OpFastScintillation::PostStepDoIt ( const G4Track &  aTrack,
const G4Step &  aStep 
)
virtual

Definition at line 358 of file OpFastScintillation.cxx.

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

Referenced by AtRestDoIt().

363  {
364  aParticleChange.Initialize(aTrack);
365  // Check that we are in a material with a properties table, if not just return
366  const G4Material* aMaterial = aTrack.GetMaterial();
367  G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
368  if (!aMaterialPropertiesTable) return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
369 
370  G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
371  G4ThreeVector x0 = pPreStepPoint->GetPosition();
372  G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
373 
374  // We don't want to produce any trackable G4 secondaries
375  aParticleChange.SetNumberOfSecondaries(0);
376 
377  // Retrieve the Scintillation Integral for this material new
378  // G4PhysicsOrderedFreeVector allocated to hold CII's
379 
380  // Some explanation for later improvements to scint yield code:
381  //
382  // What does G4 do here?
383  // It produces light in 2 steps, fast (scnt=1) then slow (scnt=2)
384  //
385  // The ratio of slow photons to fast photons is related by the yieldratio parameter.
386  // G4's poisson fluctuating scheme is a bit different to ours - we should check that
387  // they are equivalent.
388  //
389  // G4 poisson fluctuates the number of initial photons then divides them with a
390  // constant factor between fast + slow, whereas we poisson fluctuate separateyly the
391  // fast and slow detection numbers.
392  //
393 
394  // get the number of photons produced from the IonizationAndScintillation singleton
396  double MeanNumberOfPhotons =
398 
399  RecordPhotonsProduced(aStep, MeanNumberOfPhotons); //, stepEnergy);
400  if (verboseLevel > 0) {
401  G4cout << "\n Exiting from OpFastScintillation::DoIt -- NumberOfSecondaries = "
402  << aParticleChange.GetNumberOfSecondaries() << G4endl;
403  }
404  return G4VRestDiscreteProcess::PostStepDoIt(aTrack, aStep);
405  }
bool RecordPhotonsProduced(const G4Step &aStep, double N)
static IonizationAndScintillation * Instance()
void larg4::OpFastScintillation::ProcessStep ( const G4Step &  step)
private

Definition at line 407 of file OpFastScintillation.cxx.

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

Referenced by RecordPhotonsProduced().

408  {
409  if (step.GetTotalEnergyDeposit() <= 0) return;
410 
412  -1,
413  -1,
414  1.0, //scintillation yield
415  (double)(step.GetTotalEnergyDeposit() / CLHEP::MeV), //energy in MeV
416  (float)(step.GetPreStepPoint()->GetPosition().x() / CLHEP::cm),
417  (float)(step.GetPreStepPoint()->GetPosition().y() / CLHEP::cm),
418  (float)(step.GetPreStepPoint()->GetPosition().z() / CLHEP::cm),
419  (float)(step.GetPostStepPoint()->GetPosition().x() / CLHEP::cm),
420  (float)(step.GetPostStepPoint()->GetPosition().y() / CLHEP::cm),
421  (float)(step.GetPostStepPoint()->GetPosition().z() / CLHEP::cm),
422  (double)(step.GetPreStepPoint()->GetGlobalTime()),
423  (double)(step.GetPostStepPoint()->GetGlobalTime()),
425  step.GetTrack()->GetParticleDefinition()->GetPDGEncoding(),
427  step.GetPreStepPoint()->GetPhysicalVolume()->GetName());
428  }
static OpDetPhotonTable * Instance(bool LitePhotons=false)
void AddEnergyDeposit(int n_photon, int n_elec, double scint_yield, double energy, float start_x, float start_y, float start_z, float end_x, float end_y, float end_z, double start_time, double end_time, int trackid, int pdgcode, int g4trackid, std::string const &vol="EMPTY")
void larg4::OpFastScintillation::propagationTime ( std::vector< double > &  arrival_time_dist,
G4ThreeVector  x0,
const size_t  OpChannel,
bool  Reflected = false 
)
private

Definition at line 865 of file OpFastScintillation.cxx.

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

Referenced by RecordPhotonsProduced().

869  {
870  assert(fPVS);
872  throw cet::exception("OpFastScintillation")
873  << "Cannot have both propagation time models simultaneously.";
874  }
875  else if (fPVS->IncludeParPropTime() &&
876  !(ParPropTimeTF1 && (ParPropTimeTF1[OpChannel].GetNdim() == 1))) {
877  // Warning: TF1::GetNdim()==1 will tell us if the TF1 is really defined or it is the
878  // default one. This will fix a segfault when using timing and interpolation.
879  G4cout << "WARNING: Requested parameterized timing, but no function found. Not applying "
880  "propagation time."
881  << G4endl;
882  }
883  else if (fPVS->IncludeParPropTime()) {
884  if (Reflected)
885  throw cet::exception("OpFastScintillation")
886  << "No parameterized propagation time for reflected light";
887  for (size_t i = 0; i < arrival_time_dist.size(); ++i) {
888  arrival_time_dist[i] = ParPropTimeTF1[OpChannel].GetRandom();
889  }
890  }
891  else if (fPVS->IncludePropTime()) {
892  // Get VUV photons arrival time distribution from the parametrization
893  geo::Point_t const& opDetCenter = fOpDetCenter.at(OpChannel);
894  if (!Reflected) {
895  const G4ThreeVector OpDetPoint(
896  opDetCenter.X() * CLHEP::cm, opDetCenter.Y() * CLHEP::cm, opDetCenter.Z() * CLHEP::cm);
897  double distance_in_cm = (x0 - OpDetPoint).mag() / CLHEP::cm; // this must be in CENTIMETERS!
898  double cosine = std::abs(x0[0] * CLHEP::cm - OpDetPoint[0] * CLHEP::cm) / distance_in_cm;
899  double theta = fast_acos(cosine) * 180. / CLHEP::pi;
900  int angle_bin = theta / fangle_bin_timing_vuv;
901  getVUVTimes(arrival_time_dist, distance_in_cm, angle_bin); // in ns
902  }
903  else {
904  TVector3 const ScintPoint(x0[0] / CLHEP::cm, x0[1] / CLHEP::cm, x0[2] / CLHEP::cm); // in cm
905  getVISTimes(arrival_time_dist, ScintPoint, geo::vect::toTVector3(opDetCenter)); // in ns
906  }
907  }
908  }
std::vector< geo::Point_t > fOpDetCenter
constexpr auto abs(T v)
Returns the absolute value of the argument.
void getVISTimes(std::vector< double > &arrivalTimes, const TVector3 &ScintPoint, const TVector3 &OpDetPoint)
double fast_acos(double xin)
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
void getVUVTimes(std::vector< double > &arrivalTimes, const double distance_in_cm, const size_t angle_bin)
TVector3 toTVector3(Vector const &v)
Converts a vector into a TVector3.
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
phot::PhotonVisibilityService const *const fPVS
Photon visibility service instance.
phot::MappedFunctions_t ParPropTimeTF1
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool larg4::OpFastScintillation::RecordPhotonsProduced ( const G4Step &  aStep,
double  N 
)
protected

Definition at line 430 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().

431  {
432  // make sure that whatever happens afterwards, the energy deposition is stored
434  if (lgp->FillSimEnergyDeposits()) ProcessStep(aStep);
435 
436  // Get the pointer to the fast scintillation table
437  OpDetPhotonTable* fst = OpDetPhotonTable::Instance();
438 
439  const G4Track* aTrack = aStep.GetTrack();
440 
441  G4StepPoint const* pPreStepPoint = aStep.GetPreStepPoint();
442 
443  const G4DynamicParticle* aParticle = aTrack->GetDynamicParticle();
444  const G4Material* aMaterial = aTrack->GetMaterial();
445 
446  G4int materialIndex = aMaterial->GetIndex();
447  G4int tracknumber = aStep.GetTrack()->GetTrackID();
448 
449  G4ThreeVector x0 = pPreStepPoint->GetPosition();
450  G4ThreeVector p0 = aStep.GetDeltaPosition().unit();
451  G4double t0 = pPreStepPoint->GetGlobalTime();
452 
453  G4MaterialPropertiesTable* aMaterialPropertiesTable = aMaterial->GetMaterialPropertiesTable();
454 
455  G4MaterialPropertyVector* Fast_Intensity =
456  aMaterialPropertiesTable->GetProperty("FASTCOMPONENT");
457  G4MaterialPropertyVector* Slow_Intensity =
458  aMaterialPropertiesTable->GetProperty("SLOWCOMPONENT");
459 
460  if (!Fast_Intensity && !Slow_Intensity) return 1;
461 
462  if (!bPropagate) return 0;
463 
464  // Get the visibility vector for this point
465  assert(fPVS);
466 
467  G4int nscnt = 1;
468  if (Fast_Intensity && Slow_Intensity) nscnt = 2;
469 
470  double Num = 0;
471  double YieldRatio = 0;
472 
474  // The scintillation response is a function of the energy deposited by particle
475  // types.
476 
477  // Get the definition of the current particle
478  G4ParticleDefinition* pDef = aParticle->GetDefinition();
479 
480  // Obtain the G4MaterialPropertyVectory containing the scintillation light yield as
481  // a function of the deposited energy for the current particle type
482 
483  // Protons
484  if (pDef == G4Proton::ProtonDefinition()) {
485  YieldRatio = aMaterialPropertiesTable->GetConstProperty("PROTONYIELDRATIO");
486  }
487  // Muons
488  else if (pDef == G4MuonPlus::MuonPlusDefinition() ||
489  pDef == G4MuonMinus::MuonMinusDefinition()) {
490  YieldRatio = aMaterialPropertiesTable->GetConstProperty("MUONYIELDRATIO");
491  }
492  // Pions
493  else if (pDef == G4PionPlus::PionPlusDefinition() ||
494  pDef == G4PionMinus::PionMinusDefinition()) {
495  YieldRatio = aMaterialPropertiesTable->GetConstProperty("PIONYIELDRATIO");
496  }
497  // Kaons
498  else if (pDef == G4KaonPlus::KaonPlusDefinition() ||
499  pDef == G4KaonMinus::KaonMinusDefinition()) {
500  YieldRatio = aMaterialPropertiesTable->GetConstProperty("KAONYIELDRATIO");
501  }
502  // Alphas
503  else if (pDef == G4Alpha::AlphaDefinition()) {
504  YieldRatio = aMaterialPropertiesTable->GetConstProperty("ALPHAYIELDRATIO");
505  }
506  // Electrons (must also account for shell-binding energy attributed to gamma from
507  // standard PhotoElectricEffect)
508  else if (pDef == G4Electron::ElectronDefinition() || pDef == G4Gamma::GammaDefinition()) {
509  YieldRatio = aMaterialPropertiesTable->GetConstProperty("ELECTRONYIELDRATIO");
510  }
511  // Default for particles not enumerated/listed above
512  else {
513  YieldRatio = aMaterialPropertiesTable->GetConstProperty("ELECTRONYIELDRATIO");
514  }
515  // If the user has not specified yields for (p,d,t,a,carbon) then these unspecified
516  // particles will default to the electron's scintillation yield
517  if (YieldRatio == 0) {
518  YieldRatio = aMaterialPropertiesTable->GetConstProperty("ELECTRONYIELDRATIO");
519  }
520  }
521 
522  geo::Point_t const ScintPoint = {x0[0] / CLHEP::cm, x0[1] / CLHEP::cm, x0[2] / CLHEP::cm};
523  if (fOnlyActiveVolume && !isScintInActiveVolume(ScintPoint)) return 0;
524  const phot::MappedCounts_t& Visibilities = fPVS->GetAllVisibilities(ScintPoint);
525 
526  phot::MappedCounts_t ReflVisibilities;
527 
528  // Store timing information in the object for use in propagationTime method
529  if (fPVS->StoreReflected()) {
530  ReflVisibilities = fPVS->GetAllVisibilities(ScintPoint, true);
531  if (fPVS->StoreReflT0()) ReflT0s = fPVS->GetReflT0s(ScintPoint);
532  }
533  if (fPVS->IncludeParPropTime()) { ParPropTimeTF1 = fPVS->GetTimingTF1(ScintPoint); }
534 
535  for (G4int scnt = 1; scnt <= nscnt; scnt++) {
536  G4double ScintillationTime = 0. * CLHEP::ns;
537  G4double ScintillationRiseTime = 0. * CLHEP::ns;
538  G4PhysicsOrderedFreeVector* ScintillationIntegral = NULL;
539  if (scnt == 1) {
540  if (nscnt == 1) {
541  if (Fast_Intensity) {
542  ScintillationTime = aMaterialPropertiesTable->GetConstProperty("FASTTIMECONSTANT");
543  if (fFiniteRiseTime) {
544  ScintillationRiseTime =
545  aMaterialPropertiesTable->GetConstProperty("FASTSCINTILLATIONRISETIME");
546  }
547  ScintillationIntegral =
548  (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
549  }
550  if (Slow_Intensity) {
551  ScintillationTime = aMaterialPropertiesTable->GetConstProperty("SLOWTIMECONSTANT");
552  if (fFiniteRiseTime) {
553  ScintillationRiseTime =
554  aMaterialPropertiesTable->GetConstProperty("SLOWSCINTILLATIONRISETIME");
555  }
556  ScintillationIntegral =
557  (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
558  }
559  } //endif nscnt=1
560  else {
561  if (YieldRatio == 0)
562  YieldRatio = aMaterialPropertiesTable->GetConstProperty("YIELDRATIO");
563 
564  if (ExcitationRatio == 1.0) { Num = std::min(YieldRatio, 1.0) * MeanNumberOfPhotons; }
565  else {
566  Num = std::min(ExcitationRatio, 1.0) * MeanNumberOfPhotons;
567  }
568  ScintillationTime = aMaterialPropertiesTable->GetConstProperty("FASTTIMECONSTANT");
569  if (fFiniteRiseTime) {
570  ScintillationRiseTime =
571  aMaterialPropertiesTable->GetConstProperty("FASTSCINTILLATIONRISETIME");
572  }
573  ScintillationIntegral =
574  (G4PhysicsOrderedFreeVector*)((*theFastIntegralTable)(materialIndex));
575  } //endif nscnt!=1
576  } //end scnt=1
577 
578  else {
579  Num = MeanNumberOfPhotons - Num;
580  ScintillationTime = aMaterialPropertiesTable->GetConstProperty("SLOWTIMECONSTANT");
581  if (fFiniteRiseTime) {
582  ScintillationRiseTime =
583  aMaterialPropertiesTable->GetConstProperty("SLOWSCINTILLATIONRISETIME");
584  }
585  ScintillationIntegral =
586  (G4PhysicsOrderedFreeVector*)((*theSlowIntegralTable)(materialIndex));
587  }
588 
589  if (!ScintillationIntegral) continue;
590 
591  if (!Visibilities && !usesSemiAnalyticModel()) continue;
592 
593  // detected photons from direct light
594  std::map<size_t, int> DetectedNum;
595  if (Visibilities && !usesSemiAnalyticModel()) {
596  for (size_t const OpDet : util::counter(fPVS->NOpChannels())) {
597  if (fOpaqueCathode && !isOpDetInSameTPC(ScintPoint, fOpDetCenter.at(OpDet))) continue;
598  int const DetThis = std::round(G4Poisson(Visibilities[OpDet] * Num));
599  if (DetThis > 0) DetectedNum[OpDet] = DetThis;
600  }
601  }
602  else {
603  detectedDirectHits(DetectedNum, Num, ScintPoint);
604  }
605 
606  // detected photons from reflected light
607  std::map<size_t, int> ReflDetectedNum;
608  if (fPVS->StoreReflected()) {
609  if (!usesSemiAnalyticModel()) {
610  for (size_t const OpDet : util::counter(fPVS->NOpChannels())) {
611  if (fOpaqueCathode && !isOpDetInSameTPC(ScintPoint, fOpDetCenter.at(OpDet))) continue;
612  int const ReflDetThis = std::round(G4Poisson(ReflVisibilities[OpDet] * Num));
613  if (ReflDetThis > 0) ReflDetectedNum[OpDet] = ReflDetThis;
614  }
615  }
616  else {
617  detectedReflecHits(ReflDetectedNum, Num, ScintPoint);
618  }
619  }
620 
621  std::vector<double> arrival_time_dist;
622  // Now we run through each PMT figuring out num of detected photons
623  for (size_t Reflected = 0; Reflected <= 1; ++Reflected) {
624  // Only do the reflected loop if we have reflected visibilities
625  if (Reflected && !fPVS->StoreReflected()) continue;
626 
629  if (Reflected) {
630  itstart = ReflDetectedNum.begin();
631  itend = ReflDetectedNum.end();
632  }
633  else {
634  itstart = DetectedNum.begin();
635  itend = DetectedNum.end();
636  }
637  for (auto itdetphot = itstart; itdetphot != itend; ++itdetphot) {
638  const size_t OpChannel = itdetphot->first;
639  const int NPhotons = itdetphot->second;
640 
641  // Set up the OpDetBTR information
642  sim::OpDetBacktrackerRecord tmpOpDetBTRecord(OpChannel);
643  int thisG4TrackID = ParticleListAction::GetCurrentTrackID();
644  double xyzPos[3];
645  average_position(aStep, xyzPos);
646  double Edeposited = 0;
648  // We use this when it is the only sensical information. It may be of limited
649  // use to end users.
650  Edeposited = aStep.GetTotalEnergyDeposit();
651  }
652  else if (emSaturation) {
653  // If Birk Coefficient used, log VisibleEnergies.
654  Edeposited =
656  }
657  else {
658  // We use this when it is the only sensical information. It may be of limited
659  // use to end users.
660  Edeposited = aStep.GetTotalEnergyDeposit();
661  }
662 
663  // Get the transport time distribution
664  arrival_time_dist.resize(NPhotons);
665  propagationTime(arrival_time_dist, x0, OpChannel, Reflected);
666 
667  // We need to split the energy up by the number of photons so that we never try
668  // to write a 0 energy.
669  Edeposited = Edeposited / double(NPhotons);
670 
671  // Loop through the photons
672  for (G4int i = 0; i < NPhotons; ++i) {
673  G4double Time = t0 + scint_time(aStep, ScintillationTime, ScintillationRiseTime) +
674  arrival_time_dist[i] * CLHEP::ns;
675 
676  // Always store the BTR
677  tmpOpDetBTRecord.AddScintillationPhotons(thisG4TrackID, Time, 1, xyzPos, Edeposited);
678 
679  // Store as lite photon or as OnePhoton
680  if (lgp->UseLitePhotons()) {
681  fst->AddLitePhoton(OpChannel, static_cast<int>(Time), 1, Reflected);
682  }
683  else {
684  // The sim photon in this case stores its production point and time
685  TVector3 PhotonPosition(x0[0], x0[1], x0[2]);
686 
687  float PhotonEnergy = 0;
688  if (Reflected)
689  PhotonEnergy = reemission_energy() * CLHEP::eV;
690  else
691  PhotonEnergy = 9.7 * CLHEP::eV; // 9.7 eV peak of VUV emission spectrum
692 
693  // Make a photon object for the collection
694  sim::OnePhoton PhotToAdd;
695  PhotToAdd.InitialPosition = PhotonPosition;
696  PhotToAdd.Energy = PhotonEnergy;
697  PhotToAdd.Time = Time;
698  PhotToAdd.SetInSD = false;
699  PhotToAdd.MotherTrackID = tracknumber;
700 
701  fst->AddPhoton(OpChannel, std::move(PhotToAdd), Reflected);
702  }
703  }
704  fst->AddOpDetBacktrackerRecord(tmpOpDetBTRecord, Reflected);
705  }
706  }
707  }
708  return 0;
709  }
code to link reconstructed objects back to the MC truth information
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
void average_position(G4Step const &aStep, double *xzyPos) const
bool isScintInActiveVolume(geo::Point_t const &ScintPoint)
std::vector< geo::Point_t > fOpDetCenter
void propagationTime(std::vector< double > &arrival_time_dist, G4ThreeVector x0, const size_t OpChannel, bool Reflected=false)
bool const fOnlyActiveVolume
Whether photon propagation is performed only from active volumes.
All information of a photon entering the sensitive optical detector volume.
Definition: SimPhotons.h:60
bool isOpDetInSameTPC(geo::Point_t const &ScintPoint, geo::Point_t const &OpDetPoint) const
MappedFunctions_t GetTimingTF1(Point const &p) const
void ProcessStep(const G4Step &step)
std::unique_ptr< G4PhysicsTable > theSlowIntegralTable
intermediate_table::const_iterator const_iterator
G4double scint_time(const G4Step &aStep, G4double ScintillationTime, G4double ScintillationRiseTime) const
Energy deposited on a readout Optical Detector by simulated tracks.
geo::Point_t InitialPosition
Scintillation position in world coordinates [cm].
Definition: SimPhotons.h:63
bool const bPropagate
Whether propagation of photons is enabled.
auto counter(T begin, T end)
Returns an object to iterate values from begin to end in a range-for loop.
Definition: counter.h:292
bool usesSemiAnalyticModel() const
Returns whether the semi-analytic visibility parametrization is being used.
MappedT0s_t GetReflT0s(Point const &p) const
static IonizationAndScintillation * Instance()
bool FillSimEnergyDeposits() const
MappedCounts_t GetAllVisibilities(Point const &p, bool wantReflected=false) const
static OpDetPhotonTable * Instance(bool LitePhotons=false)
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
A container for photon visibility mapping data.
bool SetInSD
Whether the photon reaches the sensitive detector.
Definition: SimPhotons.h:84
phot::PhotonVisibilityService const *const fPVS
Photon visibility service instance.
bool const fOpaqueCathode
Whether the cathodes are fully opaque; currently hard coded "no".
float Energy
Scintillation photon energy [GeV].
Definition: SimPhotons.h:78
phot::MappedFunctions_t ParPropTimeTF1
bool UseLitePhotons() const
std::unique_ptr< G4PhysicsTable > theFastIntegralTable
int MotherTrackID
ID of the GEANT4 track causing the scintillation.
Definition: SimPhotons.h:81
double larg4::OpFastScintillation::Rectangle_SolidAngle ( const double  a,
const double  b,
const double  d 
) const
private

Definition at line 1704 of file OpFastScintillation.cxx.

References d, and phot::fast_acos().

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

1707  {
1708  double aa = a / (2. * d);
1709  double bb = b / (2. * d);
1710  double aux = (1. + aa * aa + bb * bb) / ((1. + aa * aa) * (1. + bb * bb));
1711  return 4. * fast_acos(std::sqrt(aux));
1712  }
double fast_acos(double xin)
Float_t d
Definition: plot.C:235
double larg4::OpFastScintillation::Rectangle_SolidAngle ( Dims const &  o,
const std::array< double, 3 >  v 
) const
private

Definition at line 1715 of file OpFastScintillation.cxx.

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

1717  {
1718  // v is the position of the track segment with respect to the center position of the
1719  // arapuca window
1720 
1721  // arapuca plane fixed in x direction
1722  if (isApproximatelyZero(v[1]) && isApproximatelyZero(v[2])) {
1723  return Rectangle_SolidAngle(o.h, o.w, v[0]);
1724  }
1725  if (isDefinitelyGreaterThan(v[1], o.h * .5) && isDefinitelyGreaterThan(v[2], o.w * .5)) {
1726  double A = v[1] - o.h * .5;
1727  double B = v[2] - o.w * .5;
1728  double to_return = (Rectangle_SolidAngle(2. * (A + o.h), 2. * (B + o.w), v[0]) -
1729  Rectangle_SolidAngle(2. * A, 2. * (B + o.w), v[0]) -
1730  Rectangle_SolidAngle(2. * (A + o.h), 2. * B, v[0]) +
1731  Rectangle_SolidAngle(2. * A, 2. * B, v[0])) *
1732  .25;
1733  return to_return;
1734  }
1735  if ((v[1] <= o.h * .5) && (v[2] <= o.w * .5)) {
1736  double A = -v[1] + o.h * .5;
1737  double B = -v[2] + o.w * .5;
1738  double to_return = (Rectangle_SolidAngle(2. * (o.h - A), 2. * (o.w - B), v[0]) +
1739  Rectangle_SolidAngle(2. * A, 2. * (o.w - B), v[0]) +
1740  Rectangle_SolidAngle(2. * (o.h - A), 2. * B, v[0]) +
1741  Rectangle_SolidAngle(2. * A, 2. * B, v[0])) *
1742  .25;
1743  return to_return;
1744  }
1745  if (isDefinitelyGreaterThan(v[1], o.h * .5) && (v[2] <= o.w * .5)) {
1746  double A = v[1] - o.h * .5;
1747  double B = -v[2] + o.w * .5;
1748  double to_return = (Rectangle_SolidAngle(2. * (A + o.h), 2. * (o.w - B), v[0]) -
1749  Rectangle_SolidAngle(2. * A, 2. * (o.w - B), v[0]) +
1750  Rectangle_SolidAngle(2. * (A + o.h), 2. * B, v[0]) -
1751  Rectangle_SolidAngle(2. * A, 2. * B, v[0])) *
1752  .25;
1753  return to_return;
1754  }
1755  if ((v[1] <= o.h * .5) && isDefinitelyGreaterThan(v[2], o.w * .5)) {
1756  double A = -v[1] + o.h * .5;
1757  double B = v[2] - o.w * .5;
1758  double to_return = (Rectangle_SolidAngle(2. * (o.h - A), 2. * (B + o.w), v[0]) -
1759  Rectangle_SolidAngle(2. * (o.h - A), 2. * B, v[0]) +
1760  Rectangle_SolidAngle(2. * A, 2. * (B + o.w), v[0]) -
1761  Rectangle_SolidAngle(2. * A, 2. * B, v[0])) *
1762  .25;
1763  return to_return;
1764  }
1765  // error message if none of these cases, i.e. something has gone wrong!
1766  return 0.;
1767  }
static constexpr bool isDefinitelyGreaterThan(TReal a, TReal b, TReal tolerance=std::numeric_limits< TReal >::epsilon())
double Rectangle_SolidAngle(const double a, const double b, const double d) const
static constexpr bool isApproximatelyZero(TReal a, TReal tolerance=std::numeric_limits< TReal >::epsilon())
double larg4::OpFastScintillation::reemission_energy ( ) const
private

Definition at line 929 of file OpFastScintillation.cxx.

References fTPBEm, and tpbemission.

Referenced by RecordPhotonsProduced().

930  {
931  return fTPBEm->fire() * ((*(--tpbemission.end())).first - (*tpbemission.begin()).first) +
932  (*tpbemission.begin()).first;
933  }
std::map< double, double > tpbemission
std::unique_ptr< CLHEP::RandGeneral > fTPBEm
void larg4::OpFastScintillation::RemoveSaturation ( )
inline

Definition at line 124 of file OpFastScintillation.hh.

Referenced by SetScintillationByParticleType().

124 { emSaturation = NULL; }
G4double larg4::OpFastScintillation::sample_time ( const G4double  tau1,
const G4double  tau2 
) const
private

Definition at line 910 of file OpFastScintillation.cxx.

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

Referenced by scint_time().

911  {
912  // tau1: rise time and tau2: decay time
913  while (1) {
914  // two random numbers
915  G4double ran1 = G4UniformRand();
916  G4double ran2 = G4UniformRand();
917  //
918  // exponential distribution as envelope function: very efficient
919  //
920  G4double d = (tau1 + tau2) / tau2;
921  // make sure the envelope function is always larger than the bi-exponential
922  G4double t = -1.0 * tau2 * std::log(1 - ran1);
923  G4double g = d * single_exp(t, tau2);
924  if (ran2 <= bi_exp(t, tau1, tau2) / g) return t;
925  }
926  return -1.0;
927  }
G4double bi_exp(const G4double t, const G4double tau1, const G4double tau2) const
Float_t d
Definition: plot.C:235
G4double single_exp(const G4double t, const G4double tau2) const
G4double larg4::OpFastScintillation::scint_time ( const G4Step &  aStep,
G4double  ScintillationTime,
G4double  ScintillationRiseTime 
) const
private

Definition at line 848 of file OpFastScintillation.cxx.

References sample_time().

Referenced by RecordPhotonsProduced().

851  {
852  G4StepPoint const* pPreStepPoint = aStep.GetPreStepPoint();
853  G4StepPoint const* pPostStepPoint = aStep.GetPostStepPoint();
854  G4double avgVelocity = (pPreStepPoint->GetVelocity() + pPostStepPoint->GetVelocity()) / 2.;
855  G4double deltaTime = aStep.GetStepLength() / avgVelocity;
856  if (ScintillationRiseTime == 0.0) {
857  deltaTime = deltaTime - ScintillationTime * std::log(G4UniformRand());
858  }
859  else {
860  deltaTime = deltaTime + sample_time(ScintillationRiseTime, ScintillationTime);
861  }
862  return deltaTime;
863  }
G4double sample_time(const G4double tau1, const G4double tau2) const
void larg4::OpFastScintillation::SetFiniteRiseTime ( const G4bool  state)
inline

Definition at line 348 of file OpFastScintillation.hh.

349  {
350  fFiniteRiseTime = state;
351  }
void larg4::OpFastScintillation::SetScintillationByParticleType ( const G4bool  scintType)

Definition at line 822 of file OpFastScintillation.cxx.

References emSaturation, RemoveSaturation(), and scintillationByParticleType.

823  {
824  if (emSaturation) {
825  G4Exception("OpFastScintillation::SetScintillationByParticleType",
826  "Scint02",
827  JustWarning,
828  "Redefinition: Birks Saturation is replaced by ScintillationByParticleType!");
830  }
831  scintillationByParticleType = scintType;
832  }
void larg4::OpFastScintillation::SetScintillationExcitationRatio ( const G4double  excitationratio)
inline

Definition at line 373 of file OpFastScintillation.hh.

374  {
375  ExcitationRatio = excitationratio;
376  }
void larg4::OpFastScintillation::SetScintillationYieldFactor ( const G4double  yieldfactor)
inline

Definition at line 363 of file OpFastScintillation.hh.

364  {
365  YieldFactor = yieldfactor;
366  }
void larg4::OpFastScintillation::SetTrackSecondariesFirst ( const G4bool  state)
inline

Definition at line 343 of file OpFastScintillation.hh.

344  {
345  fTrackSecondariesFirst = state;
346  }
G4double larg4::OpFastScintillation::single_exp ( const G4double  t,
const G4double  tau2 
) const
private

Definition at line 1539 of file OpFastScintillation.cxx.

Referenced by sample_time().

1540  {
1541  return std::exp(-1.0 * t / tau2) / tau2;
1542  }
bool larg4::OpFastScintillation::usesSemiAnalyticModel ( ) const
private

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

Definition at line 1188 of file OpFastScintillation.cxx.

References fUseNhitsModel.

Referenced by OpFastScintillation(), and RecordPhotonsProduced().

1189  {
1190  return fUseNhitsModel;
1191  }
bool const fUseNhitsModel
Whether the semi-analytic model is being used for photon visibility.
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 1392 of file OpFastScintillation.cxx.

References util::abs(), d, Disk_SolidAngle(), larg4::dist(), phot::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().

1396  {
1397  // 1). calculate total number of hits of VUV photons on reflective foils via solid
1398  // angle + Gaisser-Hillas corrections. Done outside as it doesn't depend on
1399  // OpDetPoint
1400 
1401  // set plane_depth for correct TPC:
1402  double plane_depth;
1403  if (ScintPoint_v.X() < 0.) { plane_depth = -fplane_depth; }
1404  else {
1405  plane_depth = fplane_depth;
1406  }
1407 
1408  // 2). calculate number of these hits which reach the optical detector from the
1409  // hotspot via solid angle
1410 
1411  // the interface has been converted into geo::Point_t, the implementation not yet
1412  std::array<double, 3U> ScintPoint;
1413  std::array<double, 3U> OpDetPoint;
1414  geo::vect::fillCoords(ScintPoint, ScintPoint_v);
1415  geo::vect::fillCoords(OpDetPoint, opDet.OpDetPoint);
1416 
1417  // calculate distances and angles for application of corrections; distance from
1418  // hotspot to optical detector
1419  double distance_vis = dist(&hotspot[0], &OpDetPoint[0], 3);
1420  // angle between hotspot and optical detector
1421  double cosine_vis = std::abs(hotspot[0] - OpDetPoint[0]) / distance_vis;
1422  double theta_vis = fast_acos(cosine_vis) * 180. / CLHEP::pi;
1423 
1424  // calculate solid angle of optical detector
1425  double solid_angle_detector = 0;
1426  // rectangular aperture
1427  if (opDet.type == 0) {
1428  // get hotspot coordinates relative to opDet
1429  std::array<double, 3> emission_relative = {std::abs(hotspot[0] - OpDetPoint[0]),
1430  std::abs(hotspot[1] - OpDetPoint[1]),
1431  std::abs(hotspot[2] - OpDetPoint[2])};
1432  // calculate solid angle
1433  solid_angle_detector = Rectangle_SolidAngle(Dims{opDet.h, opDet.w}, emission_relative);
1434  }
1435  // dome aperture
1436  else if (opDet.type == 1) {
1437  solid_angle_detector = Omega_Dome_Model(distance_vis, theta_vis);
1438  }
1439  // disk aperture
1440  else if (opDet.type == 2) {
1441  // offset in z-y plane
1442  double d = dist(&hotspot[1], &OpDetPoint[1], 2);
1443  // drift distance (in x)
1444  double h = std::abs(hotspot[0] - OpDetPoint[0]);
1445  // calculate solid angle
1446  solid_angle_detector = Disk_SolidAngle(d, h, fradius);
1447  }
1448  else {
1449  std::cout << "Error: Invalid optical detector type. 0 = rectangular, 1 = dome, 2 = disk"
1450  << std::endl;
1451  }
1452 
1453  // calculate number of hits via geometeric acceptance
1454  double hits_geo = (solid_angle_detector / (2. * CLHEP::pi)) *
1455  cathode_hits_rec; // 2*pi due to presence of reflective foils
1456 
1457  // determine correction factor, depending on PD type
1458  const size_t k = (theta_vis / fdelta_angulo_vis); // off-set angle bin
1459  // radial distance from centre of detector (Y-Z)
1460  // FIXME (KJK): Original implementation only accounts for the Y component
1461  double r = std::abs(ScintPoint[1] - fcathode_centre.Y());
1462  double d_c = std::abs(ScintPoint[0] - plane_depth); // distance to cathode
1463  double border_correction = 0;
1464  // flat PDs
1465  if ((opDet.type == 0 || opDet.type == 2) && fIsFlatPDCorr) {
1466  // interpolate in d_c for each r bin
1467  const size_t nbins_r = fvispars_flat[k].size();
1468  std::vector<double> interp_vals(nbins_r, 0.0);
1469  {
1470  size_t idx = 0;
1471  size_t size = fvis_distances_x_flat.size();
1472  if (d_c >= fvis_distances_x_flat[size - 2])
1473  idx = size - 2;
1474  else {
1475  while (d_c > fvis_distances_x_flat[idx + 1])
1476  idx++;
1477  }
1478  for (size_t i = 0; i < nbins_r; ++i) {
1479  interp_vals[i] = interpolate(fvis_distances_x_flat, fvispars_flat[k][i], d_c, false, idx);
1480  }
1481  }
1482  // interpolate in r
1483  border_correction = interpolate(fvis_distances_r_flat, interp_vals, r, false);
1484  }
1485  // dome PDs
1486  else if (opDet.type == 1 && fIsDomePDCorr) {
1487  // interpolate in d_c for each r bin
1488  const size_t nbins_r = fvispars_dome[k].size();
1489  std::vector<double> interp_vals(nbins_r, 0.0);
1490  {
1491  size_t idx = 0;
1492  size_t size = fvis_distances_x_dome.size();
1493  if (d_c >= fvis_distances_x_dome[size - 2])
1494  idx = size - 2;
1495  else {
1496  while (d_c > fvis_distances_x_dome[idx + 1])
1497  idx++;
1498  }
1499  for (size_t i = 0; i < nbins_r; ++i) {
1500  interp_vals[i] = interpolate(fvis_distances_x_dome, fvispars_dome[k][i], d_c, false, idx);
1501  }
1502  }
1503  // interpolate in r
1504  border_correction = interpolate(fvis_distances_r_dome, interp_vals, r, false);
1505  }
1506  else {
1507  std::cout << "Error: Invalid optical detector type. 0 = rectangular, 1 = dome, 2 = disk. Or "
1508  "corrections for chosen optical detector type missing."
1509  << std::endl;
1510  }
1511 
1512  // apply correction factor
1513  double hits_rec = border_correction * hits_geo / cosine_vis;
1514  double hits_vis = std::round(G4Poisson(hits_rec));
1515 
1516  return hits_vis;
1517  }
TRandom r
Definition: spectrum.C:23
std::vector< std::vector< std::vector< double > > > fvispars_dome
std::vector< std::vector< std::vector< double > > > fvispars_flat
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::vector< double > fvis_distances_r_flat
unsigned int fillCoords(Coords &dest, Vector const &src)
Fills a coordinate array with the coordinates of a vector.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
std::vector< double > fvis_distances_x_flat
double fast_acos(double xin)
Float_t d
Definition: plot.C:235
std::vector< double > fvis_distances_r_dome
double dist(const TReal *x, const TReal *y, const unsigned int dimension)
std::vector< double > fvis_distances_x_dome
double Rectangle_SolidAngle(const double a, const double b, const double d) const
double Omega_Dome_Model(const double distance, const double theta) const
double Disk_SolidAngle(const double d, const double h, const double b) const
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
int larg4::OpFastScintillation::VUVHits ( const double  Nphotons_created,
geo::Point_t const &  ScintPoint,
OpticalDetector const &  opDet 
) const
private

Definition at line 1287 of file OpFastScintillation.cxx.

References util::abs(), d, Disk_SolidAngle(), larg4::dist(), phot::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().

1290  {
1291  // the interface has been converted into geo::Point_t, the implementation not yet
1292  std::array<double, 3U> ScintPoint;
1293  std::array<double, 3U> OpDetPoint;
1294  geo::vect::fillCoords(ScintPoint, ScintPoint_v);
1295  geo::vect::fillCoords(OpDetPoint, opDet.OpDetPoint);
1296 
1297  // distance and angle between ScintPoint and OpDetPoint
1298  double distance = dist(&ScintPoint[0], &OpDetPoint[0], 3);
1299  double cosine = std::abs(ScintPoint[0] - OpDetPoint[0]) / distance;
1300  // double theta = std::acos(cosine) * 180. / CLHEP::pi;
1301  double theta = fast_acos(cosine) * 180. / CLHEP::pi;
1302 
1303  // calculate solid angle:
1304  double solid_angle = 0;
1305  // Arapucas/Bars (rectangle)
1306  if (opDet.type == 0) {
1307  // get scintillation point coordinates relative to arapuca window centre
1308  std::array<double, 3> ScintPoint_rel = {std::abs(ScintPoint[0] - OpDetPoint[0]),
1309  std::abs(ScintPoint[1] - OpDetPoint[1]),
1310  std::abs(ScintPoint[2] - OpDetPoint[2])};
1311  // calculate solid angle
1312  solid_angle = Rectangle_SolidAngle(Dims{opDet.h, opDet.w}, ScintPoint_rel);
1313  }
1314  // PMTs (dome)
1315  else if (opDet.type == 1) {
1316  solid_angle = Omega_Dome_Model(distance, theta);
1317  }
1318  // PMTs (disk approximation)
1319  else if (opDet.type == 2) {
1320  // offset in z-y plane
1321  double d = dist(&ScintPoint[1], &OpDetPoint[1], 2);
1322  // drift distance (in x)
1323  double h = std::abs(ScintPoint[0] - OpDetPoint[0]);
1324  // Solid angle of a disk
1325  solid_angle = Disk_SolidAngle(d, h, fradius);
1326  }
1327  else {
1328  std::cout << "Error: Invalid optical detector type. 0 = rectangular, 1 = dome, 2 = disk"
1329  << std::endl;
1330  }
1331 
1332  // calculate number of photons hits by geometric acceptance: accounting for solid
1333  // angle and LAr absorption length
1334  double hits_geo =
1335  std::exp(-1. * distance / fL_abs_vuv) * (solid_angle / (4 * CLHEP::pi)) * Nphotons_created;
1336 
1337  // apply Gaisser-Hillas correction for Rayleigh scattering distance and angular
1338  // dependence offset angle bin
1339  const size_t j = (theta / fdelta_angulo_vuv);
1340 
1341  // determine GH parameters, accounting for border effects radial distance from centre
1342  // of detector (Y-Z)
1343 
1344  // FIXME (KJK): Original implementation only accounts for the Y component
1345  double r = std::abs(ScintPoint[1] - fcathode_centre.Y());
1346  double pars_ini[4] = {0, 0, 0, 0};
1347  double s1 = 0;
1348  double s2 = 0;
1349  double s3 = 0;
1350  // flat PDs
1351  if ((opDet.type == 0 || opDet.type == 2) && fIsFlatPDCorr) {
1352  pars_ini[0] = fGHvuvpars_flat[0][j];
1353  pars_ini[1] = fGHvuvpars_flat[1][j];
1354  pars_ini[2] = fGHvuvpars_flat[2][j];
1355  pars_ini[3] = fGHvuvpars_flat[3][j];
1359  }
1360  // dome PDs
1361  else if (opDet.type == 1 && fIsDomePDCorr) {
1362  pars_ini[0] = fGHvuvpars_dome[0][j];
1363  pars_ini[1] = fGHvuvpars_dome[1][j];
1364  pars_ini[2] = fGHvuvpars_dome[2][j];
1365  pars_ini[3] = fGHvuvpars_dome[3][j];
1369  }
1370  else
1371  std::cout << "Error: Invalid optical detector type. 0 = rectangular, 1 = dome, 2 = disk. Or "
1372  "corrections for chosen optical detector type missing."
1373  << std::endl;
1374 
1375  // add border correction
1376  pars_ini[0] = pars_ini[0] + s1 * r;
1377  pars_ini[1] = pars_ini[1] + s2 * r;
1378  pars_ini[2] = pars_ini[2] + s3 * r;
1379  pars_ini[3] = pars_ini[3];
1380 
1381  // calculate correction
1382  double GH_correction = Gaisser_Hillas(distance, pars_ini);
1383 
1384  // calculate number photons
1385  double hits_rec = GH_correction * hits_geo / cosine;
1386  int hits_vuv = std::round(G4Poisson(hits_rec));
1387 
1388  return hits_vuv;
1389  }
TRandom r
Definition: spectrum.C:23
G4double Gaisser_Hillas(const double x, const double *par) const
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::vector< std::vector< double > > fGHvuvpars_flat
std::vector< std::vector< double > > fGHvuvpars_dome
unsigned int fillCoords(Coords &dest, Vector const &src)
Fills a coordinate array with the coordinates of a vector.
std::vector< double > fborder_corr_angulo_dome
std::vector< std::vector< double > > fborder_corr_dome
std::vector< std::vector< double > > fborder_corr_flat
std::vector< double > fborder_corr_angulo_flat
double fast_acos(double xin)
Float_t d
Definition: plot.C:235
double dist(const TReal *x, const TReal *y, const unsigned int dimension)
double Rectangle_SolidAngle(const double a, const double b, const double d) const
double Omega_Dome_Model(const double distance, const double theta) const
double Disk_SolidAngle(const double d, const double h, const double b) const
double interpolate(const std::vector< double > &xData, const std::vector< double > &yData, double x, bool extrapolate, size_t i=0) const
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)

Member Data Documentation

bool const larg4::OpFastScintillation::bPropagate
private

Whether propagation of photons is enabled.

Definition at line 301 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and RecordPhotonsProduced().

G4EmSaturation* larg4::OpFastScintillation::emSaturation
private
G4double larg4::OpFastScintillation::ExcitationRatio
protected

Definition at line 175 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and RecordPhotonsProduced().

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

Definition at line 287 of file OpFastScintillation.hh.

Referenced by isScintInActiveVolume(), and OpFastScintillation().

double larg4::OpFastScintillation::fangle_bin_timing_vis
private

Definition at line 238 of file OpFastScintillation.hh.

Referenced by getVISTimes(), and OpFastScintillation().

double larg4::OpFastScintillation::fangle_bin_timing_vuv
private

Definition at line 228 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and propagationTime().

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

Definition at line 268 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VUVHits().

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

Definition at line 263 of file OpFastScintillation.hh.

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

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

Definition at line 269 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VUVHits().

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

Definition at line 264 of file OpFastScintillation.hh.

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

geo::Point_t larg4::OpFastScintillation::fcathode_centre
private
Dims larg4::OpFastScintillation::fcathode_plane
private

Definition at line 291 of file OpFastScintillation.hh.

Referenced by detectedReflecHits(), and OpFastScintillation().

double larg4::OpFastScintillation::fcathode_ydimension
private

Definition at line 285 of file OpFastScintillation.hh.

Referenced by OpFastScintillation().

double larg4::OpFastScintillation::fcathode_zdimension
private

Definition at line 285 of file OpFastScintillation.hh.

Referenced by OpFastScintillation().

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

Definition at line 241 of file OpFastScintillation.hh.

Referenced by getVISTimes(), and OpFastScintillation().

double larg4::OpFastScintillation::fdelta_angulo_vis
private

Definition at line 274 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

double larg4::OpFastScintillation::fdelta_angulo_vuv
private

Definition at line 259 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VUVHits().

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

Definition at line 239 of file OpFastScintillation.hh.

Referenced by getVISTimes(), and OpFastScintillation().

G4bool larg4::OpFastScintillation::fFiniteRiseTime
protected

Definition at line 171 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and RecordPhotonsProduced().

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

Definition at line 267 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VUVHits().

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

Definition at line 262 of file OpFastScintillation.hh.

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

double larg4::OpFastScintillation::finflexion_point_distance
private

Definition at line 228 of file OpFastScintillation.hh.

Referenced by generateParam(), and OpFastScintillation().

bool larg4::OpFastScintillation::fIsDomePDCorr
private

Definition at line 266 of file OpFastScintillation.hh.

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

bool larg4::OpFastScintillation::fIsFlatPDCorr
private

Definition at line 261 of file OpFastScintillation.hh.

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

int larg4::OpFastScintillation::fL_abs_vuv
private

Definition at line 292 of file OpFastScintillation.hh.

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

double larg4::OpFastScintillation::fmax_d
private

Definition at line 228 of file OpFastScintillation.hh.

Referenced by OpFastScintillation().

double larg4::OpFastScintillation::fmin_d
private
bool const larg4::OpFastScintillation::fOnlyActiveVolume = false
private

Whether photon propagation is performed only from active volumes.

Definition at line 309 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 312 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 314 of file OpFastScintillation.hh.

Referenced by RecordPhotonsProduced().

std::vector<geo::Point_t> larg4::OpFastScintillation::fOpDetCenter
private
std::vector<double> larg4::OpFastScintillation::fOpDetHeight
private
std::vector<double> larg4::OpFastScintillation::fOpDetLength
private
std::vector<int> larg4::OpFastScintillation::fOpDetType
private
std::vector<std::vector<double> > larg4::OpFastScintillation::fparameters[7]
private

Definition at line 230 of file OpFastScintillation.hh.

Referenced by generateParam(), and OpFastScintillation().

double larg4::OpFastScintillation::fplane_depth
private
phot::PhotonVisibilityService const* const larg4::OpFastScintillation::fPVS
private

Photon visibility service instance.

Definition at line 304 of file OpFastScintillation.hh.

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

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

Definition at line 240 of file OpFastScintillation.hh.

Referenced by getVISTimes(), and OpFastScintillation().

double larg4::OpFastScintillation::fradius
private

Definition at line 290 of file OpFastScintillation.hh.

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

double larg4::OpFastScintillation::fstep_size
private
bool larg4::OpFastScintillation::fStoreReflected
private

Definition at line 272 of file OpFastScintillation.hh.

Referenced by OpFastScintillation().

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

Definition at line 242 of file OpFastScintillation.hh.

Referenced by getVISTimes(), and OpFastScintillation().

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

Definition at line 218 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and reemission_energy().

G4bool larg4::OpFastScintillation::fTrackSecondariesFirst
protected

Definition at line 170 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 307 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and usesSemiAnalyticModel().

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

Definition at line 281 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

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

Definition at line 277 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

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

Definition at line 280 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

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

Definition at line 276 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

double larg4::OpFastScintillation::fvis_vmean
private

Definition at line 238 of file OpFastScintillation.hh.

Referenced by getVISTimes(), and OpFastScintillation().

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

Definition at line 282 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

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

Definition at line 278 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and VISHits().

double larg4::OpFastScintillation::fvuv_vgroup_max
private

Definition at line 228 of file OpFastScintillation.hh.

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

double larg4::OpFastScintillation::fvuv_vgroup_mean
private

Definition at line 228 of file OpFastScintillation.hh.

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

phot::MappedFunctions_t larg4::OpFastScintillation::ParPropTimeTF1
private

Definition at line 224 of file OpFastScintillation.hh.

Referenced by propagationTime(), and RecordPhotonsProduced().

phot::MappedT0s_t larg4::OpFastScintillation::ReflT0s
private

Definition at line 225 of file OpFastScintillation.hh.

Referenced by RecordPhotonsProduced().

G4bool larg4::OpFastScintillation::scintillationByParticleType
protected
std::unique_ptr<G4PhysicsTable> larg4::OpFastScintillation::theFastIntegralTable
protected
std::unique_ptr<G4PhysicsTable> larg4::OpFastScintillation::theSlowIntegralTable
protected
std::map<double, double> larg4::OpFastScintillation::tpbemission
private

Definition at line 217 of file OpFastScintillation.hh.

Referenced by OpFastScintillation(), and reemission_energy().

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

Definition at line 234 of file OpFastScintillation.hh.

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

std::vector<std::vector<double> > larg4::OpFastScintillation::VUV_min
private
std::vector<std::vector<TF1> > larg4::OpFastScintillation::VUV_timing
private

Definition at line 232 of file OpFastScintillation.hh.

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

G4double larg4::OpFastScintillation::YieldFactor
protected

Definition at line 173 of file OpFastScintillation.hh.

Referenced by OpFastScintillation().


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