ISCalculationSeparate.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file  ISCalculationSeparate.cxx
/// \brief Interface to algorithm class for calculating ionization electrons
///        and scintillation photons using separate algorithms for each
///
/// \author  brebel@fnal.gov
////////////////////////////////////////////////////////////////////////

#include "Geant4/G4EmSaturation.hh"
#include "Geant4/G4LossTableManager.hh"
#include "Geant4/G4ParticleTypes.hh"

#include "larcore/CoreUtils/ServiceUtil.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "lardata/DetectorInfoServices/LArPropertiesService.h"
#include "larsim/LegacyLArG4/ISCalculationSeparate.h"
#include "larsim/Simulation/LArG4Parameters.h"
#include "larsim/Simulation/LArVoxelCalculator.h"

#include "cetlib_except/exception.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

namespace larg4 {

  //----------------------------------------------------------------------------
  ISCalculationSeparate::ISCalculationSeparate()
  {
    art::ServiceHandle<sim::LArG4Parameters const> lgpHandle;
    const detinfo::LArProperties* larp = lar::providerFrom<detinfo::LArPropertiesService>();
    auto const detProp =
      art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();

    double density = detProp.Density(detProp.Temperature());
    fEfield = detProp.Efield();
    fScintByParticleType = larp->ScintByParticleType();
    fGeVToElectrons = lgpHandle->GeVToElectrons();

    // \todo get scintillation yield from LArG4Parameters or LArProperties
    fScintYieldFactor = 1.;

    // 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 = lgpHandle->RecombA();
    fRecombk = lgpHandle->Recombk() / density;
    fModBoxA = lgpHandle->ModBoxA();
    fModBoxB = lgpHandle->ModBoxB() / density;
    fUseModBoxRecomb = lgpHandle->UseModBoxRecomb();

    // Use Birks Correction in the Scintillation process
    fEMSaturation = G4LossTableManager::Instance()->EmSaturation();

    // determine the step size using the voxel sizes
    art::ServiceHandle<sim::LArVoxelCalculator const> lvc;
    double maxsize =
      std::max(lvc->VoxelSizeX(), std::max(lvc->VoxelSizeY(), lvc->VoxelSizeZ())) * CLHEP::cm;

    fStepSize = 0.1 * maxsize;
  }

  //----------------------------------------------------------------------------
  // fNumIonElectrons returns a value that is not corrected for life time effects
  void
  ISCalculationSeparate::Reset()
  {
    fEnergyDeposit = 0.;
    fNumScintPhotons = 0.;
    fNumIonElectrons = 0.;
  }

  //----------------------------------------------------------------------------
  // fNumIonElectrons returns a value that is not corrected for life time effects
  void
  ISCalculationSeparate::CalculateIonizationAndScintillation(const G4Step* step)
  {
    fEnergyDeposit = step->GetTotalEnergyDeposit() / CLHEP::MeV;

    // Get the recombination factor for this voxel - Nucl.Instrum.Meth.A523:275-286,2004
    // R = A/(1 + (dE/dx)*k)
    // dE/dx is given by the voxel energy deposition, but have to convert it to MeV/cm
    // from GeV/voxel width
    // A = 0.800 +/- 0.003
    // k = (0.097+/-0.001) g/(MeVcm^2)
    // the dx depends on the trajectory of the step
    // k should be divided by the density as our dE/dx is in MeV/cm,
    // the division is handled in the constructor when we set fRecombk
    // B.Baller: Add Modified Box recombination - ArgoNeuT result submitted to JINST

    G4ThreeVector totstep = step->GetPostStepPoint()->GetPosition();
    totstep -= step->GetPreStepPoint()->GetPosition();

    double dx = totstep.mag() / CLHEP::cm;
    double recomb = 0.;
    double dEdx = (dx == 0.0) ? 0.0 : fEnergyDeposit / dx;
    double EFieldStep = EFieldAtStep(fEfield, step);

    // Guard against spurious values of dE/dx. Note: assumes density of LAr
    if (dEdx < 1.) dEdx = 1.;

    if (fUseModBoxRecomb) {
      if (dx) {
        double Xi = fModBoxB * dEdx / EFieldStep;
        recomb = log(fModBoxA + Xi) / Xi;
      }
      else
        recomb = 0;
    }
    else {
      recomb = fRecombA / (1. + dEdx * fRecombk / EFieldStep);
    }

    // 1.e-3 converts fEnergyDeposit to GeV
    fNumIonElectrons = fGeVToElectrons * 1.e-3 * fEnergyDeposit * recomb;

    MF_LOG_DEBUG("ISCalculationSeparate")
      << " Electrons produced for " << fEnergyDeposit << " MeV deposited with " << recomb
      << " recombination: " << fNumIonElectrons;

    // Now do the scintillation
    G4MaterialPropertiesTable* mpt = step->GetTrack()->GetMaterial()->GetMaterialPropertiesTable();
    if (!mpt)
      throw cet::exception("ISCalculationSeparate")
        << "Cannot find materials property table"
        << " for this step! " << step->GetTrack()->GetMaterial() << "\n";

    // if not doing the scintillation by particle type, use the saturation
    double scintYield = mpt->GetConstProperty("SCINTILLATIONYIELD");

    if (fScintByParticleType) {

      MF_LOG_DEBUG("ISCalculationSeparate") << "scintillating by particle type";

      // Get the definition of the current particle
      G4ParticleDefinition* pDef = step->GetTrack()->GetDynamicParticle()->GetDefinition();

      // Obtain the G4MaterialPropertyVectory containing the
      // scintillation light yield as a function of the deposited
      // energy for the current particle type

      // Protons
      if (pDef == G4Proton::ProtonDefinition()) {
        scintYield = mpt->GetConstProperty("PROTONSCINTILLATIONYIELD");
      }
      // Muons
      else if (pDef == G4MuonPlus::MuonPlusDefinition() ||
               pDef == G4MuonMinus::MuonMinusDefinition()) {
        scintYield = mpt->GetConstProperty("MUONSCINTILLATIONYIELD");
      }
      // Pions
      else if (pDef == G4PionPlus::PionPlusDefinition() ||
               pDef == G4PionMinus::PionMinusDefinition()) {
        scintYield = mpt->GetConstProperty("PIONSCINTILLATIONYIELD");
      }
      // Kaons
      else if (pDef == G4KaonPlus::KaonPlusDefinition() ||
               pDef == G4KaonMinus::KaonMinusDefinition()) {
        scintYield = mpt->GetConstProperty("KAONSCINTILLATIONYIELD");
      }
      // Alphas
      else if (pDef == G4Alpha::AlphaDefinition()) {
        scintYield = mpt->GetConstProperty("ALPHASCINTILLATIONYIELD");
      }
      // Electrons (must also account for shell-binding energy
      // attributed to gamma from standard PhotoElectricEffect)
      else if (pDef == G4Electron::ElectronDefinition() || pDef == G4Gamma::GammaDefinition()) {
        scintYield = mpt->GetConstProperty("ELECTRONSCINTILLATIONYIELD");
      }
      // Default for particles not enumerated/listed above
      else {
        scintYield = mpt->GetConstProperty("ELECTRONSCINTILLATIONYIELD");
      }

      // If the user has not specified yields for (p,d,t,a,carbon)
      // then these unspecified particles will default to the
      // electron's scintillation yield

      // Throw an exception if no scintillation yield is found
      if (!scintYield)
        throw cet::exception("ISCalculationSeparate")
          << "Request for scintillation yield for energy "
          << "deposit and particle type without correct "
          << "entry in MaterialPropertiesTable\n"
          << "ScintillationByParticleType requires at "
          << "minimum that ELECTRONSCINTILLATIONYIELD is "
          << "set by the user\n";

      fNumScintPhotons = scintYield * fEnergyDeposit;
    }
    else if (fEMSaturation) {
      // The default linear scintillation process
      fVisibleEnergyDeposition = fEMSaturation->VisibleEnergyDepositionAtAStep(step);
      fNumScintPhotons = fScintYieldFactor * scintYield * fVisibleEnergyDeposition;
    }
    else {
      fNumScintPhotons = fScintYieldFactor * scintYield * fEnergyDeposit;
    }

    MF_LOG_DEBUG("ISCalculationSeparate")
      << "number photons: " << fNumScintPhotons << " energy: " << fEnergyDeposit / CLHEP::MeV
      << " saturation: " << fEMSaturation->VisibleEnergyDepositionAtAStep(step)
      << " step length: " << step->GetStepLength() / CLHEP::cm;
  }

} // namespace