LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
ISCalcCorrelated.cxx
Go to the documentation of this file.
1 // Class: ISCalcCorrelated
3 // Plugin Type: algorithm
4 // File: ISCalcCorrelated.h and ISCalcCorrelated.cxx
5 // Description: Interface to algorithm class for a specific calculation of
6 // ionization electrons and scintillation photons, based on
7 // simple microphysics arguments to establish an anticorrelation
8 // between these two quantities.
9 // Input: 'sim::SimEnergyDeposit'
10 // Output: num of Photons and Electrons
11 // May 2020 by W Foreman
12 // Modified: Adding corrections for low electric field (LArQL model)
13 // Jun 2020 by L. Paulucci and F. Marinho
15 
17 
25 
26 #include "CLHEP/Random/RandBinomial.h"
29 
30 #include <iostream>
31 
32 namespace larg4 {
33  //----------------------------------------------------------------------------
35  CLHEP::HepRandomEngine& Engine)
36  : fISTPC{*(lar::providerFrom<geo::Geometry>())}
37  , fSCE(lar::providerFrom<spacecharge::SpaceChargeService>())
38  , fBinomialGen{CLHEP::RandBinomial(Engine)}
39  {
40  MF_LOG_INFO("ISCalcCorrelated") << "IonizationAndScintillation/ISCalcCorrelated Initialize.";
41 
42  fScintPreScale = lar::providerFrom<detinfo::LArPropertiesService>()->ScintPreScale();
43 
45 
46  // The recombination coefficient is in g/(MeVcm^2), but we report
47  // energy depositions in MeV/cm, need to divide Recombk from the
48  // LArG4Parameters service by the density of the argon we got
49  // above.
50  fRecombA = LArG4PropHandle->RecombA();
51  fRecombk = LArG4PropHandle->Recombk() / detProp.Density(detProp.Temperature());
52  fModBoxA = LArG4PropHandle->ModBoxA();
53  fModBoxB = LArG4PropHandle->ModBoxB() / detProp.Density(detProp.Temperature());
54  fEllipsModBoxA = LArG4PropHandle->EllipsModBoxA();
55  fEllipsModBoxB = LArG4PropHandle->EllipsModBoxB() / detProp.Density(detProp.Temperature());
56  fEllipsModBoxR = LArG4PropHandle->EllipsModBoxR();
57  fUseModBoxRecomb = (bool)LArG4PropHandle->UseModBoxRecomb();
58  fUseEllipsModBoxRecomb = (bool)LArG4PropHandle->UseEllipsModBoxRecomb();
59  fUseModLarqlRecomb = (bool)LArG4PropHandle->UseModLarqlRecomb();
60  fUseBinomialFlucts = (bool)LArG4PropHandle->UseBinomialFlucts();
61  fLarqlChi0A = LArG4PropHandle->LarqlChi0A();
62  fLarqlChi0B = LArG4PropHandle->LarqlChi0B();
63  fLarqlChi0C = LArG4PropHandle->LarqlChi0C();
64  fLarqlChi0D = LArG4PropHandle->LarqlChi0D();
65  fLarqlAlpha = LArG4PropHandle->LarqlAlpha();
66  fLarqlBeta = LArG4PropHandle->LarqlBeta();
67  fQAlpha = LArG4PropHandle->QAlpha();
68  fGeVToElectrons = LArG4PropHandle->GeVToElectrons();
69 
70  // ionization work function
71  fWion = 1. / fGeVToElectrons * 1e3; // MeV
72 
73  // ion+excitation work function
74  fWph = LArG4PropHandle->Wph() * 1e-6; // MeV
75  }
76 
77  //----------------------------------------------------------------------------
80  {
81 
82  double const energy_deposit = edep.Energy();
83 
84  // calculate total quanta (ions + excitons)
85  double num_ions = 0.0; //check if the deposited energy is above ionization threshold
86  if (energy_deposit >= fWion) num_ions = energy_deposit / fWion;
87  double num_quanta = energy_deposit / fWph;
88 
89  double ds = edep.StepLength();
90  double dEdx = (ds <= 0.0) ? 0.0 : energy_deposit / ds;
91  dEdx = (dEdx < 1.) ? 1. : dEdx;
92  double EFieldStep = EFieldAtStep(detProp.Efield(), edep);
93  double recomb = 0., num_electrons = 0.;
94 
95  //calculate recombination survival fraction value inside, otherwise zero
96  if (EFieldStep > 0.) {
97  // calculate recombination survival fraction
98  // ...using Modified Box model
99  if (fUseModBoxRecomb) {
100  double Xi = fModBoxB * dEdx / EFieldStep;
101  recomb = std::log(fModBoxA + Xi) / Xi;
102  }
103  else if (fUseEllipsModBoxRecomb) {
104 
105  double phi = AngleToEFieldAtStep(detProp.Efield(), edep);
106 
107  if (std::isnan(phi)) {
108  double Xi = fModBoxB * dEdx / EFieldStep;
109  recomb = std::log(fModBoxA + Xi) / Xi;
110  }
111  else {
112  double B_ellips =
113  fEllipsModBoxB * dEdx /
114  (EFieldStep * std::hypot(std::sin(phi), std::cos(phi) / fEllipsModBoxR));
115 
116  recomb = std::log(fEllipsModBoxA + B_ellips) / B_ellips;
117  }
118  }
119  // ... or using Birks/Doke
120  else {
121  recomb = fRecombA / (1. + dEdx * fRecombk / EFieldStep);
122  }
123  }
124 
125  if (fUseModLarqlRecomb &&
126  edep.PdgCode() != 1000020040) { //Use corrections from LArQL model (except for alpha)
127  recomb += EscapingEFraction(dEdx) * FieldCorrection(EFieldStep, dEdx); //Correction for low EF
128  }
129 
130  // Guard against unphysical recombination values
131  if (recomb < 0.) {
132  mf::LogWarning("ISCalcCorrelated")
133  << "Recombination survival fraction is lower than 0.: " << recomb << ", fixing it to 0.";
134  recomb = 0.;
135  }
136  else if (recomb > 1.) {
137  mf::LogWarning("ISCalcCorrelated")
138  << "Recombination survival fraction is higher than 1.: " << recomb << ", fixing it to 1.";
139  recomb = 1.;
140  }
141 
142  // using this recombination, calculate number energy_deposit of ionization electrons
143  if (num_ions > 0.)
144  num_electrons =
145  (fUseBinomialFlucts) ? fBinomialGen.fire(num_ions, recomb) : (num_ions * recomb);
146 
147  // calculate scintillation photons
148  double num_photons = (num_quanta - num_electrons) * fScintPreScale;
149 
150  if (edep.PdgCode() == 1000020040) {
151  num_electrons = num_electrons * fQAlpha;
152  num_photons = num_photons * fQAlpha;
153  }
154 
155  MF_LOG_DEBUG("ISCalcCorrelated")
156  << "With " << energy_deposit << " MeV of deposited energy, "
157  << "and a recombination of " << recomb << ", \nthere are " << num_electrons
158  << " electrons, and " << num_photons << " photons.";
159 
160  return {energy_deposit, num_electrons, num_photons, GetScintYieldRatio(edep)};
161  }
162 
163  //----------------------------------------------------------------------------
165  {
166  // electric field outside active volume set to zero
167  if (!fISTPC.isScintInActiveVolume(edep.MidPoint())) return 0.;
168 
169  geo::Vector_t elecvec{};
170 
172  geo::TPCID tpcid = fGeometry->PositionToTPCID(edep.MidPoint());
173  if (!bool(tpcid)) return 0.;
174  const geo::TPCGeo& tpcGeo = fGeometry->TPC(tpcid);
175 
176  using geo::to_int;
177  auto const [axis, sign] = tpcGeo.DriftAxisWithSign();
178  switch (axis) {
179  case geo::Coordinate::X: {
180  elecvec.SetX(to_int(sign));
181  break;
182  }
183  case geo::Coordinate::Y: {
184  elecvec.SetY(to_int(sign));
185  break;
186  }
187  case geo::Coordinate::Z: {
188  elecvec.SetZ(to_int(sign));
189  }
190  }
191 
192  elecvec *= efield;
193 
194  // electric field inside active volume
195  if (fSCE->EnableSimEfieldSCE()) {
196  auto const eFieldOffsets = fSCE->GetEfieldOffsets(edep.MidPoint());
197  geo::Vector_t scevec = eFieldOffsets * efield;
199  switch (axis) {
200  case geo::Coordinate::X: {
201  scevec.SetX(scevec.X() * -1.);
202  break;
203  }
204  case geo::Coordinate::Y: {
205  scevec.SetY(scevec.Y() * -1.);
206  break;
207  }
208  case geo::Coordinate::Z: {
209  scevec.SetZ(scevec.Z() * -1.);
210  }
211  }
212  }
213  elecvec += scevec;
214  }
215 
216  return elecvec.R();
217  }
218  //----------------------------------------------------------------------------
220  {
221 
222  // electric field outside active volume set to zero
223  if (!fISTPC.isScintInActiveVolume(edep.MidPoint())) return 0.;
224 
225  geo::Vector_t stepvec = edep.Start() - edep.End();
226  geo::Vector_t elecvec{};
227 
229  geo::TPCID tpcid = fGeometry->PositionToTPCID(edep.MidPoint());
230  if (!bool(tpcid)) return 0.;
231  const geo::TPCGeo& tpcGeo = fGeometry->TPC(tpcid);
232 
233  using geo::to_int;
234  auto const [axis, sign] = tpcGeo.DriftAxisWithSign();
235  switch (axis) {
236  case geo::Coordinate::X: {
237  elecvec.SetX(to_int(sign));
238  break;
239  }
240  case geo::Coordinate::Y: {
241  elecvec.SetY(to_int(sign));
242  break;
243  }
244  case geo::Coordinate::Z: {
245  elecvec.SetZ(to_int(sign));
246  }
247  }
248 
249  elecvec *= efield;
250 
251  // electric field inside active volume
252  if (fSCE->EnableSimEfieldSCE()) {
253  auto const eFieldOffsets = fSCE->GetEfieldOffsets(edep.MidPoint());
254  geo::Vector_t scevec = eFieldOffsets * efield;
256  switch (axis) {
257  case geo::Coordinate::X: {
258  scevec.SetX(scevec.X() * -1.);
259  break;
260  }
261  case geo::Coordinate::Y: {
262  scevec.SetY(scevec.Y() * -1.);
263  break;
264  }
265  case geo::Coordinate::Z: {
266  scevec.SetZ(scevec.Z() * -1.);
267  }
268  }
269  }
270  elecvec += scevec;
271  }
272 
273  double angle = std::acos(stepvec.Dot(elecvec) / (stepvec.R() * elecvec.R()));
274 
275  if (angle > TMath::PiOver2()) { angle = abs(TMath::Pi() - angle); }
276 
277  return angle;
278  }
279 
280  //----------------------------------------------------------------------------
281  // LArQL chi0 function = fraction of escaping electrons
282  double ISCalcCorrelated::EscapingEFraction(double const dEdx) const
283  {
284  return fLarqlChi0A / (fLarqlChi0B + std::exp(fLarqlChi0C + fLarqlChi0D * dEdx));
285  }
286 
287  //----------------------------------------------------------------------------
288  // LArQL f_corr function = correction factor for electric field dependence
289  double ISCalcCorrelated::FieldCorrection(double const EF, double const dEdx) const
290  {
291  return std::exp(-EF / (fLarqlAlpha * std::log(dEdx) + fLarqlBeta));
292  }
293 
294 } // namespace larg4
virtual bool EnableSimEfieldSCE() const =0
Store parameters for running LArG4.
double FieldCorrection(double const EF, double const dEdx) const
double EscapingEFraction(double const dEdx) const
geo::Length_t StepLength() const
double fModBoxA
from LArG4Parameters service
double ModBoxA() const
Utilities related to art service access.
ISCalcData CalcIonAndScint(detinfo::DetectorPropertiesData const &detProp, sim::SimEnergyDeposit const &edep) override
double fRecombk
from LArG4Parameters service
double fEllipsModBoxA
from LArG4Parameters service
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:160
double fWion
W_ion (23.6 eV) == 1/fGeVToElectrons.
CLHEP::RandBinomial fBinomialGen
double fGeVToElectrons
from LArG4Parameters service
double LarqlChi0B() const
double fModBoxB
from LArG4Parameters service
Geometry information for a single TPC.
Definition: TPCGeo.h:33
double fQAlpha
from LArG4Parameters service
double fEllipsModBoxR
from LArG4Parameters service
bool fUseModLarqlRecomb
from LArG4Parameters service
double LarqlBeta() const
double LarqlAlpha() const
constexpr auto abs(T v)
Returns the absolute value of the argument.
bool UseModBoxRecomb() const
double AngleToEFieldAtStep(double efield, sim::SimEnergyDeposit const &edep)
Geant4 interface.
double LarqlChi0D() const
double fEllipsModBoxB
from LArG4Parameters service
double EllipsModBoxA() const
DriftAxis DriftAxisWithSign() const
Returns the expected drift direction based on geometry.
Definition: TPCGeo.h:78
bool UseBinomialFlucts() const
double QAlpha() const
double Efield(unsigned int planegap=0) const
kV/cm
double EllipsModBoxR() const
double fScintPreScale
scintillation pre-scaling factor from LArProperties service
double EFieldAtStep(double efield, sim::SimEnergyDeposit const &edep) override
double LarqlChi0C() const
double fWph
from LArG4Parameters service
geo::Point_t Start() const
geo::Point_t End() const
bool fUseBinomialFlucts
from LArG4Parameters service
bool UseModLarqlRecomb() const
double fLarqlChi0B
from LArG4Parameters service
Drift towards negative values.
double Wph() const
constexpr int to_int(Coordinate const coord) noexcept
Enumerate the possible plane projections.
Definition: geo_types.h:124
bool fUseEllipsModBoxRecomb
from LArG4Parameters service
double fLarqlChi0D
from LArG4Parameters service
geo::Point_t MidPoint() const
bool fUseModBoxRecomb
from LArG4Parameters service
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
double RecombA() const
double fLarqlChi0C
from LArG4Parameters service
double LarqlChi0A() const
#define MF_LOG_INFO(category)
Double_t edep
Definition: macro.C:13
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2671
double fLarqlChi0A
from LArG4Parameters service
int sign(double val)
Definition: UtilFunc.cxx:104
virtual geo::Vector_t GetEfieldOffsets(geo::Point_t const &point) const =0
double EllipsModBoxB() const
double ModBoxB() const
contains information for a single step in the detector simulation
double Recombk() const
bool UseEllipsModBoxRecomb() const
double fLarqlBeta
from LArG4Parameters service
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Energy deposition in the active material.
#define MF_LOG_DEBUG(id)
double fLarqlAlpha
from LArG4Parameters service
double Energy() const
TPCGeo const & TPC(TPCID const &tpcid=details::tpc_zero) const
Returns the specified TPC.
Definition: GeometryCore.h:448
const spacecharge::SpaceCharge * fSCE
Float_t e
Definition: plot.C:35
bool isScintInActiveVolume(geo::Point_t const &ScintPoint)
Definition: ISTPC.cxx:37
double fRecombA
from LArG4Parameters service
double GetScintYieldRatio(sim::SimEnergyDeposit const &edep)
Definition: ISCalc.cxx:41
ISCalcCorrelated(detinfo::DetectorPropertiesData const &detProp, CLHEP::HepRandomEngine &Engine)
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