26 #include "CLHEP/Random/RandBinomial.h" 35 CLHEP::HepRandomEngine& Engine)
36 : fISTPC{*(lar::providerFrom<geo::Geometry>())}
37 ,
fSCE(lar::providerFrom<spacecharge::SpaceChargeService>())
40 MF_LOG_INFO(
"ISCalcCorrelated") <<
"IonizationAndScintillation/ISCalcCorrelated Initialize.";
42 fScintPreScale = lar::providerFrom<detinfo::LArPropertiesService>()->ScintPreScale();
51 fRecombk = LArG4PropHandle->
Recombk() / detProp.Density(detProp.Temperature());
53 fModBoxB = LArG4PropHandle->
ModBoxB() / detProp.Density(detProp.Temperature());
74 fWph = LArG4PropHandle->
Wph() * 1
e-6;
82 double const energy_deposit = edep.
Energy();
85 double num_ions = 0.0;
86 if (energy_deposit >=
fWion) num_ions = energy_deposit /
fWion;
87 double num_quanta = energy_deposit /
fWph;
90 double dEdx = (ds <= 0.0) ? 0.0 : energy_deposit /
ds;
91 dEdx = (dEdx < 1.) ? 1. :
dEdx;
93 double recomb = 0., num_electrons = 0.;
96 if (EFieldStep > 0.) {
100 double Xi =
fModBoxB * dEdx / EFieldStep;
101 recomb = std::log(
fModBoxA + Xi) / Xi;
107 if (std::isnan(phi)) {
108 double Xi =
fModBoxB * dEdx / EFieldStep;
109 recomb = std::log(
fModBoxA + Xi) / Xi;
114 (EFieldStep * std::hypot(std::sin(phi), std::cos(phi) /
fEllipsModBoxR));
126 edep.
PdgCode() != 1000020040) {
133 <<
"Recombination survival fraction is lower than 0.: " << recomb <<
", fixing it to 0.";
136 else if (recomb > 1.) {
138 <<
"Recombination survival fraction is higher than 1.: " << recomb <<
", fixing it to 1.";
148 double num_photons = (num_quanta - num_electrons) *
fScintPreScale;
150 if (edep.
PdgCode() == 1000020040) {
151 num_electrons = num_electrons *
fQAlpha;
152 num_photons = num_photons *
fQAlpha;
156 <<
"With " << energy_deposit <<
" MeV of deposited energy, " 157 <<
"and a recombination of " << recomb <<
", \nthere are " << num_electrons
158 <<
" electrons, and " << num_photons <<
" photons.";
173 if (!
bool(tpcid))
return 0.;
201 scevec.SetX(scevec.X() * -1.);
205 scevec.SetY(scevec.Y() * -1.);
209 scevec.SetZ(scevec.Z() * -1.);
230 if (!
bool(tpcid))
return 0.;
258 scevec.SetX(scevec.X() * -1.);
262 scevec.SetY(scevec.Y() * -1.);
266 scevec.SetZ(scevec.Z() * -1.);
273 double angle = std::acos(stepvec.Dot(elecvec) / (stepvec.R() * elecvec.R()));
275 if (angle > TMath::PiOver2()) { angle =
abs(TMath::Pi() - angle); }
virtual bool EnableSimEfieldSCE() const =0
Store parameters for running LArG4.
geo::Length_t StepLength() const
Utilities related to art service access.
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
double LarqlChi0B() const
Geometry information for a single TPC.
double LarqlAlpha() const
constexpr auto abs(T v)
Returns the absolute value of the argument.
bool UseModBoxRecomb() const
double LarqlChi0D() const
double EllipsModBoxA() const
DriftAxis DriftAxisWithSign() const
Returns the expected drift direction based on geometry.
bool UseBinomialFlucts() const
double Efield(unsigned int planegap=0) const
kV/cm
double EllipsModBoxR() const
double LarqlChi0C() const
geo::Point_t Start() const
bool UseModLarqlRecomb() const
Drift towards negative values.
constexpr int to_int(Coordinate const coord) noexcept
Enumerate the possible plane projections.
geo::Point_t MidPoint() const
The data type to uniquely identify a TPC.
double LarqlChi0A() const
#define MF_LOG_INFO(category)
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
virtual geo::Vector_t GetEfieldOffsets(geo::Point_t const &point) const =0
double EllipsModBoxB() const
contains information for a single step in the detector simulation
bool UseEllipsModBoxRecomb() const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Energy deposition in the active material.
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
bool isScintInActiveVolume(geo::Point_t const &ScintPoint)
double GetScintYieldRatio(sim::SimEnergyDeposit const &edep)
TPCID PositionToTPCID(Point_t const &point) const
Returns the ID of the TPC at specified location.
double GeVToElectrons() const
art framework interface to geometry description