#include "ISCalcCorrelated.h"

Inheritance diagram for larg4::ISCalcCorrelated:

Public Member Functions
	ISCalcCorrelated (detinfo::DetectorPropertiesData const &detProp, CLHEP::HepRandomEngine &Engine)

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

double	AngleToEFieldAtStep (double efield, sim::SimEnergyDeposit const &edep)

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

double	GetScintYield (sim::SimEnergyDeposit const &edep, bool prescale)

double	GetScintYieldRatio (sim::SimEnergyDeposit const &edep)

Private Member Functions
double	EscapingEFraction (double const dEdx) const

double	FieldCorrection (double const EF, double const dEdx) const

Private Attributes
ISTPC	fISTPC

const spacecharge::SpaceCharge *	fSCE

CLHEP::RandBinomial	fBinomialGen

double	fGeVToElectrons
	from LArG4Parameters service More...

double	fWion
	W_ion (23.6 eV) == 1/fGeVToElectrons. More...

double	fWph
	from LArG4Parameters service More...

double	fScintPreScale
	scintillation pre-scaling factor from LArProperties service More...

double	fRecombA
	from LArG4Parameters service More...

double	fRecombk
	from LArG4Parameters service More...

double	fModBoxA
	from LArG4Parameters service More...

double	fModBoxB
	from LArG4Parameters service More...

double	fEllipsModBoxA
	from LArG4Parameters service More...

double	fEllipsModBoxB
	from LArG4Parameters service More...

double	fEllipsModBoxR
	from LArG4Parameters service More...

double	fLarqlChi0A
	from LArG4Parameters service More...

double	fLarqlChi0B
	from LArG4Parameters service More...

double	fLarqlChi0C
	from LArG4Parameters service More...

double	fLarqlChi0D
	from LArG4Parameters service More...

double	fLarqlAlpha
	from LArG4Parameters service More...

double	fLarqlBeta
	from LArG4Parameters service More...

double	fQAlpha
	from LArG4Parameters service More...

bool	fUseModBoxRecomb
	from LArG4Parameters service More...

bool	fUseEllipsModBoxRecomb
	from LArG4Parameters service More...

bool	fUseModLarqlRecomb
	from LArG4Parameters service More...

bool	fUseBinomialFlucts
	from LArG4Parameters service More...

Detailed Description

Definition at line 33 of file ISCalcCorrelated.h.

Constructor & Destructor Documentation

larg4::ISCalcCorrelated::ISCalcCorrelated	(	detinfo::DetectorPropertiesData const &	detProp,
		CLHEP::HepRandomEngine &	Engine
	)

Definition at line 34 of file ISCalcCorrelated.cxx.

     : fISTPC{*(lar::providerFrom<geo::Geometry>())}
     , fSCE(lar::providerFrom<spacecharge::SpaceChargeService>())
     , fBinomialGen{CLHEP::RandBinomial(Engine)}
   {
     MF_LOG_INFO("ISCalcCorrelated") << "IonizationAndScintillation/ISCalcCorrelated Initialize.";
 
     fScintPreScale = lar::providerFrom<detinfo::LArPropertiesService>()->ScintPreScale();
 
     art::ServiceHandle<sim::LArG4Parameters const> LArG4PropHandle;
 
     // The recombination coefficient is in g/(MeVcm^2), but we report
     // energy depositions in MeV/cm, need to divide Recombk from the
     // LArG4Parameters service by the density of the argon we got
     // above.
     fRecombA = LArG4PropHandle->RecombA();
     fRecombk = LArG4PropHandle->Recombk() / detProp.Density(detProp.Temperature());
     fModBoxA = LArG4PropHandle->ModBoxA();
     fModBoxB = LArG4PropHandle->ModBoxB() / detProp.Density(detProp.Temperature());
     fEllipsModBoxA = LArG4PropHandle->EllipsModBoxA();
     fEllipsModBoxB = LArG4PropHandle->EllipsModBoxB() / detProp.Density(detProp.Temperature());
     fEllipsModBoxR = LArG4PropHandle->EllipsModBoxR();
     fUseModBoxRecomb = (bool)LArG4PropHandle->UseModBoxRecomb();
     fUseEllipsModBoxRecomb = (bool)LArG4PropHandle->UseEllipsModBoxRecomb();
     fUseModLarqlRecomb = (bool)LArG4PropHandle->UseModLarqlRecomb();
     fUseBinomialFlucts = (bool)LArG4PropHandle->UseBinomialFlucts();
     fLarqlChi0A = LArG4PropHandle->LarqlChi0A();
     fLarqlChi0B = LArG4PropHandle->LarqlChi0B();
     fLarqlChi0C = LArG4PropHandle->LarqlChi0C();
     fLarqlChi0D = LArG4PropHandle->LarqlChi0D();
     fLarqlAlpha = LArG4PropHandle->LarqlAlpha();
     fLarqlBeta = LArG4PropHandle->LarqlBeta();
     fQAlpha = LArG4PropHandle->QAlpha();
     fGeVToElectrons = LArG4PropHandle->GeVToElectrons();
 
     // ionization work function
     fWion = 1. / fGeVToElectrons * 1e3; // MeV
 
     // ion+excitation work function
     fWph = LArG4PropHandle->Wph() * 1e-6; // MeV
   }

Member Function Documentation

double larg4::ISCalcCorrelated::AngleToEFieldAtStep	(	double	efield,
		sim::SimEnergyDeposit const &	edep
	)

Definition at line 219 of file ISCalcCorrelated.cxx.

References util::abs(), geo::TPCGeo::DriftAxisWithSign(), spacecharge::SpaceCharge::EnableSimEfieldSCE(), sim::SimEnergyDeposit::End(), fISTPC, fSCE, spacecharge::SpaceCharge::GetEfieldOffsets(), larg4::ISTPC::isScintInActiveVolume(), sim::SimEnergyDeposit::MidPoint(), geo::Negative, geo::GeometryCore::PositionToTPCID(), pmtana::sign(), sim::SimEnergyDeposit::Start(), geo::to_int(), geo::GeometryCore::TPC(), geo::X, geo::Y, and geo::Z.

Referenced by CalcIonAndScint().

   {
 
     // electric field outside active volume set to zero
     if (!fISTPC.isScintInActiveVolume(edep.MidPoint())) return 0.;
 
     geo::Vector_t stepvec = edep.Start() - edep.End();
     geo::Vector_t elecvec{};
 
     art::ServiceHandle<geo::Geometry const> fGeometry;
     geo::TPCID tpcid = fGeometry->PositionToTPCID(edep.MidPoint());
     if (!bool(tpcid)) return 0.;
     const geo::TPCGeo& tpcGeo = fGeometry->TPC(tpcid);
 
     using geo::to_int;
     auto const [axis, sign] = tpcGeo.DriftAxisWithSign();
     switch (axis) {
     case geo::Coordinate::X: {
       elecvec.SetX(to_int(sign));
       break;
     }
     case geo::Coordinate::Y: {
       elecvec.SetY(to_int(sign));
       break;
     }
     case geo::Coordinate::Z: {
       elecvec.SetZ(to_int(sign));
     }
     }
 
     elecvec *= efield;
 
     // electric field inside active volume
     if (fSCE->EnableSimEfieldSCE()) {
       auto const eFieldOffsets = fSCE->GetEfieldOffsets(edep.MidPoint());
       geo::Vector_t scevec = eFieldOffsets * efield;
       if (sign == geo::DriftSign::Negative) {
         switch (axis) {
         case geo::Coordinate::X: {
           scevec.SetX(scevec.X() * -1.);
           break;
         }
         case geo::Coordinate::Y: {
           scevec.SetY(scevec.Y() * -1.);
           break;
         }
         case geo::Coordinate::Z: {
           scevec.SetZ(scevec.Z() * -1.);
         }
         }
       }
       elecvec += scevec;
     }
 
     double angle = std::acos(stepvec.Dot(elecvec) / (stepvec.R() * elecvec.R()));
 
     if (angle > TMath::PiOver2()) { angle = abs(TMath::Pi() - angle); }
 
     return angle;
   }

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

overridevirtual

Implements larg4::ISCalc.

Definition at line 78 of file ISCalcCorrelated.cxx.

References AngleToEFieldAtStep(), tca::dEdx(), edep, detinfo::DetectorPropertiesData::Efield(), EFieldAtStep(), sim::SimEnergyDeposit::Energy(), EscapingEFraction(), fBinomialGen, fEllipsModBoxA, fEllipsModBoxB, fEllipsModBoxR, FieldCorrection(), fModBoxA, fModBoxB, fQAlpha, fRecombA, fRecombk, fScintPreScale, fUseBinomialFlucts, fUseEllipsModBoxRecomb, fUseModBoxRecomb, fUseModLarqlRecomb, fWion, fWph, larg4::ISCalc::GetScintYieldRatio(), MF_LOG_DEBUG, sim::SimEnergyDeposit::PdgCode(), and sim::SimEnergyDeposit::StepLength().

   {
 
     double const energy_deposit = edep.Energy();
 
     // calculate total quanta (ions + excitons)
     double num_ions = 0.0; //check if the deposited energy is above ionization threshold
     if (energy_deposit >= fWion) num_ions = energy_deposit / fWion;
     double num_quanta = energy_deposit / fWph;
 
     double ds = edep.StepLength();
     double dEdx = (ds <= 0.0) ? 0.0 : energy_deposit / ds;
     dEdx = (dEdx < 1.) ? 1. : dEdx;
     double EFieldStep = EFieldAtStep(detProp.Efield(), edep);
     double recomb = 0., num_electrons = 0.;
 
     //calculate recombination survival fraction value inside, otherwise zero
     if (EFieldStep > 0.) {
       // calculate recombination survival fraction
       // ...using Modified Box model
       if (fUseModBoxRecomb) {
         double Xi = fModBoxB * dEdx / EFieldStep;
         recomb = std::log(fModBoxA + Xi) / Xi;
       }
       else if (fUseEllipsModBoxRecomb) {
 
         double phi = AngleToEFieldAtStep(detProp.Efield(), edep);
 
         if (std::isnan(phi)) {
           double Xi = fModBoxB * dEdx / EFieldStep;
           recomb = std::log(fModBoxA + Xi) / Xi;
         }
         else {
           double B_ellips =
             fEllipsModBoxB * dEdx /
             (EFieldStep * std::hypot(std::sin(phi), std::cos(phi) / fEllipsModBoxR));
 
           recomb = std::log(fEllipsModBoxA + B_ellips) / B_ellips;
         }
       }
       // ... or using Birks/Doke
       else {
         recomb = fRecombA / (1. + dEdx * fRecombk / EFieldStep);
       }
     }
 
     if (fUseModLarqlRecomb &&
         edep.PdgCode() != 1000020040) { //Use corrections from LArQL model (except for alpha)
       recomb += EscapingEFraction(dEdx) * FieldCorrection(EFieldStep, dEdx); //Correction for low EF
     }
 
     // Guard against unphysical recombination values
     if (recomb < 0.) {
       mf::LogWarning("ISCalcCorrelated")
         << "Recombination survival fraction is lower than 0.: " << recomb << ", fixing it to 0.";
       recomb = 0.;
     }
     else if (recomb > 1.) {
       mf::LogWarning("ISCalcCorrelated")
         << "Recombination survival fraction is higher than 1.: " << recomb << ", fixing it to 1.";
       recomb = 1.;
     }
 
     // using this recombination, calculate number energy_deposit of ionization electrons
     if (num_ions > 0.)
       num_electrons =
         (fUseBinomialFlucts) ? fBinomialGen.fire(num_ions, recomb) : (num_ions * recomb);
 
     // calculate scintillation photons
     double num_photons = (num_quanta - num_electrons) * fScintPreScale;
 
     if (edep.PdgCode() == 1000020040) {
       num_electrons = num_electrons * fQAlpha;
       num_photons = num_photons * fQAlpha;
     }
 
     MF_LOG_DEBUG("ISCalcCorrelated")
       << "With " << energy_deposit << " MeV of deposited energy, "
       << "and a recombination of " << recomb << ", \nthere are " << num_electrons
       << " electrons, and " << num_photons << " photons.";
 
     return {energy_deposit, num_electrons, num_photons, GetScintYieldRatio(edep)};
   }

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

overridevirtual

Implements larg4::ISCalc.

Definition at line 164 of file ISCalcCorrelated.cxx.

References geo::TPCGeo::DriftAxisWithSign(), spacecharge::SpaceCharge::EnableSimEfieldSCE(), fISTPC, fSCE, spacecharge::SpaceCharge::GetEfieldOffsets(), larg4::ISTPC::isScintInActiveVolume(), sim::SimEnergyDeposit::MidPoint(), geo::Negative, geo::GeometryCore::PositionToTPCID(), pmtana::sign(), geo::to_int(), geo::GeometryCore::TPC(), geo::X, geo::Y, and geo::Z.

Referenced by CalcIonAndScint().

   {
     // electric field outside active volume set to zero
     if (!fISTPC.isScintInActiveVolume(edep.MidPoint())) return 0.;
 
     geo::Vector_t elecvec{};
 
     art::ServiceHandle<geo::Geometry const> fGeometry;
     geo::TPCID tpcid = fGeometry->PositionToTPCID(edep.MidPoint());
     if (!bool(tpcid)) return 0.;
     const geo::TPCGeo& tpcGeo = fGeometry->TPC(tpcid);
 
     using geo::to_int;
     auto const [axis, sign] = tpcGeo.DriftAxisWithSign();
     switch (axis) {
     case geo::Coordinate::X: {
       elecvec.SetX(to_int(sign));
       break;
     }
     case geo::Coordinate::Y: {
       elecvec.SetY(to_int(sign));
       break;
     }
     case geo::Coordinate::Z: {
       elecvec.SetZ(to_int(sign));
     }
     }
 
     elecvec *= efield;
 
     // electric field inside active volume
     if (fSCE->EnableSimEfieldSCE()) {
       auto const eFieldOffsets = fSCE->GetEfieldOffsets(edep.MidPoint());
       geo::Vector_t scevec = eFieldOffsets * efield;
       if (sign == geo::DriftSign::Negative) {
         switch (axis) {
         case geo::Coordinate::X: {
           scevec.SetX(scevec.X() * -1.);
           break;
         }
         case geo::Coordinate::Y: {
           scevec.SetY(scevec.Y() * -1.);
           break;
         }
         case geo::Coordinate::Z: {
           scevec.SetZ(scevec.Z() * -1.);
         }
         }
       }
       elecvec += scevec;
     }
 
     return elecvec.R();
   }

double larg4::ISCalcCorrelated::EscapingEFraction ( double const dEdx ) const

private

Definition at line 282 of file ISCalcCorrelated.cxx.

References fLarqlChi0A, fLarqlChi0B, fLarqlChi0C, and fLarqlChi0D.

Referenced by CalcIonAndScint().

   {
     return fLarqlChi0A / (fLarqlChi0B + std::exp(fLarqlChi0C + fLarqlChi0D * dEdx));
   }

double larg4::ISCalcCorrelated::FieldCorrection	(	double const	EF,
		double const	dEdx
	)		const

private

Definition at line 289 of file ISCalcCorrelated.cxx.

References fLarqlAlpha, and fLarqlBeta.

Referenced by CalcIonAndScint().

   {
     return std::exp(-EF / (fLarqlAlpha * std::log(dEdx) + fLarqlBeta));
   }

double larg4::ISCalc::GetScintYield	(	sim::SimEnergyDeposit const &	edep,
		bool	prescale
	)

inherited

Definition at line 21 of file ISCalc.cxx.

References detinfo::LArProperties::AlphaScintYield(), detinfo::LArProperties::ElectronScintYield(), larg4::ISCalc::fLArProp, detinfo::LArProperties::KaonScintYield(), detinfo::LArProperties::MuonScintYield(), sim::SimEnergyDeposit::PdgCode(), detinfo::LArProperties::PionScintYield(), detinfo::LArProperties::ProtonScintYield(), detinfo::LArProperties::ScintByParticleType(), and detinfo::LArProperties::ScintYield().

Referenced by larg4::ISCalcSeparate::CalcScint().

   {
     if (!fLArProp->ScintByParticleType()) return fLArProp->ScintYield(true);
     switch (edep.PdgCode()) {
     case 2212: return fLArProp->ProtonScintYield(true);
     case 13:
     case -13: return fLArProp->MuonScintYield(true);
     case 211:
     case -211: return fLArProp->PionScintYield(true);
     case 321:
     case -321: return fLArProp->KaonScintYield(true);
     case 1000020040: return fLArProp->AlphaScintYield(true);
     case 11:
     case -11:
     case 22: return fLArProp->ElectronScintYield(true);
     default: return fLArProp->ElectronScintYield(true);
     }
   }

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

inherited

Definition at line 41 of file ISCalc.cxx.

References detinfo::LArProperties::AlphaScintYieldRatio(), detinfo::LArProperties::ElectronScintYieldRatio(), larg4::ISCalc::fLArProp, detinfo::LArProperties::KaonScintYieldRatio(), detinfo::LArProperties::MuonScintYieldRatio(), sim::SimEnergyDeposit::PdgCode(), detinfo::LArProperties::PionScintYieldRatio(), detinfo::LArProperties::ProtonScintYieldRatio(), detinfo::LArProperties::ScintByParticleType(), and detinfo::LArProperties::ScintYieldRatio().

Referenced by larg4::ISCalcNESTLAr::CalcIonAndScint(), CalcIonAndScint(), and larg4::ISCalcSeparate::CalcScint().

   {
     // ScintByParticleType option only controls the scintillation
     // yield ratio, which is the ratio of fast light (singlet
     // component) to the total light (singlet+triplet components).
 
     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();
     }
   }