TPC: add UseGeant4Edep flag with ionisation charge fluctuation for Kr-83m calibration

Adds a boolean flag TPCDetParam.UseGeant4Edep (default false) that switches
ProcessHits' ionisation calculation from the standard Bethe-Bloch/NA49
collision-sampling model to the direct Geant4 energy deposit
(fMC->Edep()), with Fano-limited Gamma-distributed smearing via
ScaleFactorG4/FanoFactorG4 (TPCGasParam), matching the O2 Geant4-specific
ionisation model used before commit 7860b87.

A SetSpecialPhysicsCuts() override sets 1 keV VMC cuts on the TPC drift gas
mediums before simcuts.dat is read, ensuring Auger electrons (1.921 keV
minimum) are tracked and not killed early.

The flag is off by default so the normal simulation chain is unaffected.
Activate for Kr calibration via:
  --configKeyValues "TPCDetParam.UseGeant4Edep=1"

Author: Ankur Yadav

Detectors/TPC/base/include/TPCBase/ParameterDetector.h

@@ -33,6 +33,7 @@ struct ParameterDetector : public o2::conf::ConfigurableParamHelper<ParameterDet
   TimeBin TmaxTriggered = 550;   ///< Maximum time bin in case of triggered readout mode
   float DriftTimeOffset = 7.3;   ///< drift time offset in time bins (we observe ~2.4\mus before October 2023 and ~1.45 \mus after)
   bool ExcludeFCGap = true;      ///< exclude electrons created in the gap between the IFC vessel and OFC vessel and FC strips
+  bool UseGeant4Edep = false;    ///< use Geant4 energy deposit directly for ionisation (for Kr-83m calibration runs)
 
   O2ParamDef(ParameterDetector, "TPCDetParam");
 };

Detectors/TPC/simulation/include/TPCSimulation/Detector.h

@@ -146,6 +146,8 @@ class Detector : public o2::base::DetImpl<Detector>
   void PostTrack() override { ; }
   void PreTrack() override { ; }
 
+  void SetSpecialPhysicsCuts() override;
+
   void SetGeoFileName(const TString file) { mGeoFileName = file; }
   const TString& GetGeoFileName() const { return mGeoFileName; }
 

Detectors/TPC/simulation/src/Detector.cxx

@@ -193,41 +193,52 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
   Int_t numberOfElectrons = 0;
   // I.H. - the type expected in addHit is short
 
-  // ---| Stepsize in cm |---
-  const double stepSize = fMC->TrackStep();
-
-  double betaGamma = momentum.P() / fMC->TrackMass();
-  betaGamma = TMath::Max(betaGamma, 7.e-3); // protection against too small bg
-
-  // ---| number of primary ionisations per cm |---
-  const double primaryElectronsPerCM =
-    gasParam.Nprim * BetheBlochAleph(static_cast<float>(betaGamma), gasParam.BetheBlochParam[0],
-                                     gasParam.BetheBlochParam[1], gasParam.BetheBlochParam[2],
-                                     gasParam.BetheBlochParam[3], gasParam.BetheBlochParam[4]);
-
-  // ---| mean number of collisions and random for this event |---
-  const double meanNcoll = stepSize * trackCharge * trackCharge * primaryElectronsPerCM;
-  const int nColl = static_cast<int>(fMC->GetRandom()->Poisson(meanNcoll));
-
-  // Variables needed to generate random powerlaw distributed energy loss
-  const double alpha_p1 = 1. - gasParam.Exp; // NA49/G3 value
-  const double oneOverAlpha_p1 = 1. / alpha_p1;
-  const double eMin = gasParam.Ipot;
-  const double eMax = gasParam.Eend;
-  const double kMin = TMath::Power(eMin, alpha_p1);
-  const double kMax = TMath::Power(eMax, alpha_p1);
-  const double wIon = gasParam.Wion;
-
-  for (Int_t n = 0; n < nColl; n++) {
-    // Use GEANT3 / NA49 expression:
-    // P(eDep) ~ k * edep^-gasParam.getExp()
-    // eMin(~I) < eDep < eMax(300 electrons)
-    // k fixed so that Int_Emin^EMax P(Edep) = 1.
-    const double rndm = fMC->GetRandom()->Rndm();
-    const double eDep = TMath::Power((kMax - kMin) * rndm + kMin, oneOverAlpha_p1);
-    int nel_step = static_cast<int>(((eDep - eMin) / wIon) + 1);
-    nel_step = TMath::Min(nel_step, 300); // 300 electrons corresponds to 10 keV
-    numberOfElectrons += nel_step;
+  // use Geant4 energy deposit directly for ionisation (Kr-83m calibration simulations)
+  if (detParam.UseGeant4Edep) {
+    // We have multiple collisions and add fluctuations: smear nel using
+    // gamma distr with mean = meanIon and variance = meanIon*FanoFactorG4.
+    // These parameters were tuned for GEANT4.
+    const double meanIon = fMC->Edep() / (gasParam.Wion * gasParam.ScaleFactorG4);
+    if (meanIon > 0.) {
+      numberOfElectrons = static_cast<int>(gasParam.FanoFactorG4 * Gamma(meanIon / gasParam.FanoFactorG4));
+    }
+  } else {
+    // ---| Stepsize in cm |---
+    const double stepSize = fMC->TrackStep();
+
+    double betaGamma = momentum.P() / fMC->TrackMass();
+    betaGamma = TMath::Max(betaGamma, 7.e-3); // protection against too small bg
+
+    // ---| number of primary ionisations per cm |---
+    const double primaryElectronsPerCM =
+      gasParam.Nprim * BetheBlochAleph(static_cast<float>(betaGamma), gasParam.BetheBlochParam[0],
+                                       gasParam.BetheBlochParam[1], gasParam.BetheBlochParam[2],
+                                       gasParam.BetheBlochParam[3], gasParam.BetheBlochParam[4]);
+
+    // ---| mean number of collisions and random for this event |---
+    const double meanNcoll = stepSize * trackCharge * trackCharge * primaryElectronsPerCM;
+    const int nColl = static_cast<int>(fMC->GetRandom()->Poisson(meanNcoll));
+
+    // Variables needed to generate random powerlaw distributed energy loss
+    const double alpha_p1 = 1. - gasParam.Exp; // NA49/G3 value
+    const double oneOverAlpha_p1 = 1. / alpha_p1;
+    const double eMin = gasParam.Ipot;
+    const double eMax = gasParam.Eend;
+    const double kMin = TMath::Power(eMin, alpha_p1);
+    const double kMax = TMath::Power(eMax, alpha_p1);
+    const double wIon = gasParam.Wion;
+
+    for (Int_t n = 0; n < nColl; n++) {
+      // Use GEANT3 / NA49 expression:
+      // P(eDep) ~ k * edep^-gasParam.getExp()
+      // eMin(~I) < eDep < eMax(300 electrons)
+      // k fixed so that Int_Emin^EMax P(Edep) = 1.
+      const double rndm = fMC->GetRandom()->Rndm();
+      const double eDep = TMath::Power((kMax - kMin) * rndm + kMin, oneOverAlpha_p1);
+      int nel_step = static_cast<int>(((eDep - eMin) / wIon) + 1);
+      nel_step = TMath::Min(nel_step, 300); // 300 electrons corresponds to 10 keV
+      numberOfElectrons += nel_step;
+    }
   }
 
   // LOG(info) << "tpc::AddHit" << FairLogger::endl << "Eloss: "
@@ -3240,6 +3251,23 @@ std::string Detector::getHitBranchNames(int probe) const
   return std::string();
 }
 
+void Detector::SetSpecialPhysicsCuts()
+{
+  // lower energy threshold to track low-energy electrons for Kr-83m calibration
+  auto const& detParam = ParameterDetector::Instance();
+  LOG(info) << "TPC SetSpecialPhysicsCuts: UseGeant4Edep=" << detParam.UseGeant4Edep;
+  if (detParam.UseGeant4Edep) {
+    auto& matmgr = o2::base::MaterialManager::Instance();
+    for (int med : {(int)kDriftGas1, (int)kDriftGas2, (int)kCO2}) {
+      matmgr.SpecialCut(GetName(), med, o2::base::ECut::kCUTELE, 1e-6f);
+      matmgr.SpecialCut(GetName(), med, o2::base::ECut::kCUTGAM, 1e-6f);
+      matmgr.SpecialCut(GetName(), med, o2::base::ECut::kDCUTE,  1e-6f);
+      matmgr.SpecialCut(GetName(), med, o2::base::ECut::kBCUTE,  1e-6f);
+    }
+  }
+  o2::base::Detector::SetSpecialPhysicsCuts();
+}
+
 ClassImp(o2::tpc::Detector);
 
 // Define Factory method for calling from the outside