#include <ISCalcNESTLAr.h>

Inheritance diagram for larg4::ISCalcNESTLAr:

Public Member Functions
	ISCalcNESTLAr (CLHEP::HepRandomEngine &fEngine)

double	EFieldAtStep (double efield, sim::SimEnergyDeposit const &edep) override

ISCalcData	CalcIonAndScint (detinfo::DetectorPropertiesData const &detProp, sim::SimEnergyDeposit const &edep) override

Public Member Functions inherited from larg4::ISCalc
virtual	~ISCalc ()=default

Private Member Functions
int	BinomFluct (int N0, double prob)

double	CalcElectronLET (double E)

double	GetScintYieldRatio (sim::SimEnergyDeposit const &edep)

Private Attributes
CLHEP::HepRandomEngine &	fEngine

const spacecharge::SpaceCharge *	fSCE

const detinfo::LArProperties *	fLArProp

Detailed Description

Definition at line 23 of file ISCalcNESTLAr.h.

Constructor & Destructor Documentation

larg4::ISCalcNESTLAr::ISCalcNESTLAr ( CLHEP::HepRandomEngine & fEngine )

explicit

Definition at line 31 of file ISCalcNESTLAr.cxx.

     : fEngine(Engine)
     , fSCE{lar::providerFrom<spacecharge::SpaceChargeService>()}
     , fLArProp{lar::providerFrom<detinfo::LArPropertiesService>()}
   {
     std::cout << "ISCalcNESTLAr Initialize." << std::endl;
   }

Member Function Documentation

int larg4::ISCalcNESTLAr::BinomFluct	(	int	N0,
		double	prob
	)

private

Definition at line 176 of file ISCalcNESTLAr.cxx.

   {
     CLHEP::RandGauss GaussGen(fEngine);
     CLHEP::RandFlat UniformGen(fEngine);
 
     double mean = N0 * prob;
     double sigma = sqrt(N0 * prob * (1 - prob));
     int N1 = 0;
     if (prob == 0.00) { return N1; }
     if (prob == 1.00) { return N0; }
 
     if (N0 < 10) {
       for (int i = 0; i < N0; i++) {
         if (UniformGen.fire() < prob) { N1++; }
       }
     }
     else {
       N1 = int(floor(GaussGen.fire(mean, sigma) + 0.5));
     }
     if (N1 > N0) { N1 = N0; }
     if (N1 < 0) { N1 = 0; }
     return N1;
   }

double larg4::ISCalcNESTLAr::CalcElectronLET ( double E )

private

Definition at line 202 of file ISCalcNESTLAr.cxx.

   {
     double LET;
 
     if (E >= 1) {
       LET = 116.70 - 162.97 * log10(E) + 99.361 * pow(log10(E), 2) - 33.405 * pow(log10(E), 3) +
             6.5069 * pow(log10(E), 4) - 0.69334 * pow(log10(E), 5) + .031563 * pow(log10(E), 6);
     }
     else if (E > 0 && E < 1) {
       LET = 100;
     }
     else {
       LET = 0;
     }
 
     return LET;
   }

ISCalcData larg4::ISCalcNESTLAr::CalcIonAndScint	(	detinfo::DetectorPropertiesData const &	detProp,
		sim::SimEnergyDeposit const &	edep
	)

overridevirtual

Implements larg4::ISCalc.

Definition at line 41 of file ISCalcNESTLAr.cxx.

   {
     CLHEP::RandGauss GaussGen(fEngine);
     CLHEP::RandFlat UniformGen(fEngine);
 
     double yieldFactor = 1.0; // default quenching factor, for electronic recoils
     double excitationRatio =
       0.21; // ratio for light particle in LAr, such as e-, mu-, Aprile et. al book
 
     double const energyDeposit = edep.Energy();
     if (energyDeposit < 1 * CLHEP::eV) // too small energy deposition
     {
       return {0., 0., 0., 0.};
     }
 
     int pdgcode = edep.PdgCode();
 
     geo::Length_t startx = edep.StartX();
     geo::Length_t starty = edep.StartY();
     geo::Length_t startz = edep.StartZ();
     geo::Length_t endx = edep.EndX();
     geo::Length_t endy = edep.EndY();
     geo::Length_t endz = edep.EndZ();
 
     double DokeBirks[3];
 
     double eField = EFieldAtStep(detProp.Efield(), edep);
     if (eField) {
       DokeBirks[0] = 0.07 * pow((eField / 1.0e3), -0.85);
       DokeBirks[2] = 0.00;
     }
     else {
       DokeBirks[0] = 0.0003;
       DokeBirks[2] = 0.75;
     }
 
     double Density = detProp.Density() /
                      (CLHEP::g / CLHEP::cm3); // argon density at the temperature from Temperature()
 
     // nuclear recoil quenching "L" factor: total yield is
     // reduced for nuclear recoil as per Lindhard theory
     double epsilon = 11.5 * (energyDeposit / CLHEP::keV) * pow(LAr_Z, (-7. / 3.));
 
     if (pdgcode == 2112 || pdgcode == -2112) //nuclear recoil
     {
       yieldFactor = 0.23 * (1 + exp(-5 * epsilon)); //liquid argon L_eff
       excitationRatio = 0.69337 + 0.3065 * exp(-0.008806 * pow(eField, 0.76313));
     }
 
     // determine ultimate number of quanta from current E-deposition (ph+e-) total mean number of exc/ions
     //the total number of either quanta produced is equal to product of the
     //work function, the energy deposited, and yield reduction, for NR
     double MeanNumQuanta = scint_yield * energyDeposit;
     double sigma = sqrt(resolution_scale * MeanNumQuanta); //Fano
     int NumQuanta = int(floor(GaussGen.fire(MeanNumQuanta, sigma) + 0.5));
     double LeffVar = GaussGen.fire(yieldFactor, 0.25 * yieldFactor);
     LeffVar = std::clamp(LeffVar, 0., 1.);
 
     if (yieldFactor < 1) //nuclear reocils
     {
       NumQuanta = BinomFluct(NumQuanta, LeffVar);
     }
 
     //if Edep below work function, can't make any quanta, and if NumQuanta
     //less than zero because Gaussian fluctuated low, update to zero
     if (energyDeposit < 1 / scint_yield || NumQuanta < 0) { NumQuanta = 0; }
 
     // next section binomially assigns quanta to excitons and ions
     int NumExcitons = BinomFluct(NumQuanta, excitationRatio / (1 + excitationRatio));
     int NumIons = NumQuanta - NumExcitons;
 
     // this section calculates recombination following the modified Birks'Law of Doke, deposition by deposition,
     // may be overridden later in code if a low enough energy necessitates switching to the
     // Thomas-Imel box model for recombination instead (determined by site)
     double dE = energyDeposit / CLHEP::MeV;
     double dx = 0.0;
     double LET = 0.0;
     double recombProb;
 
     if (pdgcode != 11 && pdgcode != -11 && pdgcode != 13 &&
         pdgcode != -13) //e-: 11, e+: -11, mu-: 13, mu+: -13
     {
       //in other words, if it's a gamma,ion,proton,alpha,pion,et al. do not
       //use the step length provided by Geant4 because it's not relevant,
       //instead calculate an estimated LET and range of the electrons that
       //would have been produced if Geant4 could track them
       LET = CalcElectronLET(1000 * dE);
 
       if (LET) {
         dx = dE / (Density * LET); //find the range based on the LET
       }
 
       if (abs(pdgcode) == 2112) //nuclear recoils
       {
         dx = 0;
       }
     }
     else //normal case of an e-/+ energy deposition recorded by Geant
     {
       dx = std::hypot(startx - endx, starty - endy, startz - endz) / CLHEP::cm;
       if (dx) {
         LET = (dE / dx) * (1 / Density); //lin. energy xfer (prop. to dE/dx)
       }
       if (LET > 0 && dE > 0 && dx > 0) {
         double ratio = CalcElectronLET(dE * 1e3) / LET;
         if (ratio < 0.7 && pdgcode == 11) {
           dx /= ratio;
           LET *= ratio;
         }
       }
     }
 
     DokeBirks[1] = DokeBirks[0] / (1 - DokeBirks[2]); //B=A/(1-C) (see paper)    r
     recombProb = (DokeBirks[0] * LET) / (1 + DokeBirks[1] * LET) +
                  DokeBirks[2]; //Doke/Birks' Law as spelled out in the NEST pape
     recombProb *= (Density / Density_LAr);
 
     //check against unphysicality resulting from rounding errors
     recombProb = std::clamp(recombProb, 0., 1.);
 
     //use binomial distribution to assign photons, electrons, where photons
     //are excitons plus recombined ionization electrons, while final
     //collected electrons are the "escape" (non-recombined) electrons
     int const NumPhotons = NumExcitons + BinomFluct(NumIons, recombProb);
     int const NumElectrons = NumQuanta - NumPhotons;
 
     return {energyDeposit,
             static_cast<double>(NumElectrons),
             static_cast<double>(NumPhotons),
             GetScintYieldRatio(edep)};
   }

double larg4::ISCalcNESTLAr::EFieldAtStep	(	double	efield,
		sim::SimEnergyDeposit const &	edep
	)

overridevirtual

Implements larg4::ISCalc.

Definition at line 222 of file ISCalcNESTLAr.cxx.

   {
     geo::Point_t pos = edep.MidPoint();
     double EField = efield;
 
     geo::Vector_t eFieldOffsets;
 
     if (fSCE->EnableSimEfieldSCE()) {
       eFieldOffsets = fSCE->GetEfieldOffsets(pos);
       EField =
         std::sqrt((efield + efield * eFieldOffsets.X()) * (efield + efield * eFieldOffsets.X()) +
                   (efield * eFieldOffsets.Y() * efield * eFieldOffsets.Y()) +
                   (efield * eFieldOffsets.Z() * efield * eFieldOffsets.Z()));
     }
 
     return EField;
   }

double larg4::ISCalcNESTLAr::GetScintYieldRatio ( sim::SimEnergyDeposit const & edep )

private

Definition at line 242 of file ISCalcNESTLAr.cxx.

   {
     if (!fLArProp->ScintByParticleType()) { return fLArProp->ScintYieldRatio(); }
     switch (edep.PdgCode()) {
     case 2212: return fLArProp->ProtonScintYieldRatio();
     case 13:
     case -13: return fLArProp->MuonScintYieldRatio();
     case 211:
     case -211: return fLArProp->PionScintYieldRatio();
     case 321:
     case -321: return fLArProp->KaonScintYieldRatio();
     case 1000020040: return fLArProp->AlphaScintYieldRatio();
     case 11:
     case -11:
     case 22: return fLArProp->ElectronScintYieldRatio();
     default: return fLArProp->ElectronScintYieldRatio();
     }
   }

Member Data Documentation

CLHEP::HepRandomEngine& larg4::ISCalcNESTLAr::fEngine

private

Definition at line 34 of file ISCalcNESTLAr.h.

const detinfo::LArProperties* larg4::ISCalcNESTLAr::fLArProp

private

Definition at line 36 of file ISCalcNESTLAr.h.

const spacecharge::SpaceCharge* larg4::ISCalcNESTLAr::fSCE

private

Definition at line 35 of file ISCalcNESTLAr.h.

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

larsim/larsim/IonizationScintillation/ISCalcNESTLAr.h
larsim/larsim/IonizationScintillation/ISCalcNESTLAr.cxx

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation