larg4::ISCalcNESTLAr Class Reference

#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
double	GetScintYield (sim::SimEnergyDeposit const &edep, bool prescale)
double	GetScintYieldRatio (sim::SimEnergyDeposit const &edep)

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

Private Attributes
CLHEP::HepRandomEngine &	fEngine
const spacecharge::SpaceCharge *	fSCE

Detailed Description

Definition at line 23 of file ISCalcNESTLAr.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 33 of file ISCalcNESTLAr.cxx.

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

Member Function Documentation

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

private

Definition at line 174 of file ISCalcNESTLAr.cxx.

References fEngine, and pmtana::mean().

Referenced by CalcIonAndScint().

  {
    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 199 of file ISCalcNESTLAr.cxx.

Referenced by CalcIonAndScint().

  {
    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 40 of file ISCalcNESTLAr.cxx.

References util::abs(), BinomFluct(), CalcElectronLET(), detinfo::DetectorPropertiesData::Density(), Density_LAr, edep, detinfo::DetectorPropertiesData::Efield(), EFieldAtStep(), sim::SimEnergyDeposit::EndX(), sim::SimEnergyDeposit::EndY(), sim::SimEnergyDeposit::EndZ(), sim::SimEnergyDeposit::Energy(), fEngine, larg4::ISCalc::GetScintYieldRatio(), sim::SimEnergyDeposit::PdgCode(), sim::SimEnergyDeposit::StartX(), sim::SimEnergyDeposit::StartY(), and sim::SimEnergyDeposit::StartZ().

  {
    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 218 of file ISCalcNESTLAr.cxx.

References spacecharge::SpaceCharge::EnableSimEfieldSCE(), fSCE, spacecharge::SpaceCharge::GetEfieldOffsets(), and sim::SimEnergyDeposit::MidPoint().

Referenced by CalcIonAndScint().

  {
    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::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 CalcIonAndScint(), larg4::ISCalcCorrelated::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();
    }
  }

Member Data Documentation

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

private

Definition at line 34 of file ISCalcNESTLAr.h.

Referenced by BinomFluct(), and CalcIonAndScint().

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

private

Definition at line 35 of file ISCalcNESTLAr.h.

Referenced by EFieldAtStep().

---

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

• larsim/v09_43_00/source/larsim/IonizationScintillation/ISCalcNESTLAr.h
• larsim/v09_43_00/source/larsim/IonizationScintillation/ISCalcNESTLAr.cxx