TRD ionisation model independent of the Geant4 energy loss

Detectors/TRD/simulation/CMakeLists.txt

@@ -12,6 +12,7 @@
 o2_add_library(TRDSimulation
                TARGETVARNAME targetName
                SOURCES src/Detector.cxx
+                       src/IonisationModel.cxx
                        src/TRsim.cxx
                        src/DiffAndTimeStructEstimator.cxx
                        src/Digitizer.cxx
@@ -66,6 +67,12 @@ o2_add_test(PileupTool
             ENVIRONMENT VMCWORKDIR=${CMAKE_BINARY_DIR}/stage
             LABELS trd)
 
+o2_add_test(IonisationModel
+            SOURCES test/testTRDIonisationModel.cxx
+            COMPONENT_NAME trd
+            PUBLIC_LINK_LIBRARIES O2::TRDSimulation
+            LABELS trd)
+
 o2_add_test(DiffusionCoefficient
             SOURCES test/testTRDDiffusionCoefficient.cxx
             COMPONENT_NAME trd

Detectors/TRD/simulation/README.md

@@ -25,6 +25,17 @@ o2-sim -m PIPE MAG TRD -n 1000 -g boxgen --configKeyValues 'BoxGun.pdg=211;BoxGu
 --configKeyValues 'TRDSimParams.doTR=false;TRDSimParams.maxMCStepSize=0.1'
 ```
 
+The charge of a hit in the gas comes from one of two models, selected with `TRDSimParams.ionisationModel`:
+- `0`: the Geant4 energy deposit in the gas divided by Wion. Geant4 then uses the PAI model in the TRD gas
+  (`Detectors/gconfig/g4config_trd_pai.in`, applied by `g4Config.C`).
+- `1`: the TRD collision model (`src/IonisationModel.cxx`). Geant4 runs without PAI in the gas. The tables are made
+  with `macros/makeIonisationTables.cxx` from Geant4's PAI cross-section, for the 10 keV delta-ray cut of the gas.
+  Slow tracks keep the energy deposit.
+```
+--configKeyValues 'TRDSimParams.ionisationModel=1'
+```
+A custom `G4.configMacroFile` replaces `g4config.in` and the PAI file; it must contain the PAI lines itself if wanted.
+
 # Generate Raw data from MonteCarlo:
 
 ```

Detectors/TRD/simulation/data/simcuts.dat

@@ -4,6 +4,6 @@
 *    Med GAM   ELEC  NHAD  CHAD  MUON  EBREM MUHAB EDEL  MUDEL MUPA  ANNI  BREM  COMP  DCAY  DRAY  HADR  LOSS  MULS  PAIR  PHOT  RAYL  STRA
 * Mylar
 TRD  27  1e-5  1e-5  1e-3  1e-3  1e-5  1e-5  1e-5  1e-5  1e-5  -1.   1     1     1     1     1     1     3     1     1     1     1     0
-* Gas mixture
+* Gas mixture; its DCUTE must stay the delta-ray cut of the TRD ionisation tables (IonisationModel)
 TRD   9  1e-5  1e-5  1e-3  1e-3  1e-5  1e-5  1e-5  1e-5  1e-5  -1.   1     1     1     1     1     1     3     1     1     1     1     1
 *

Detectors/TRD/simulation/include/TRDSimulation/Detector.h

@@ -12,12 +12,14 @@
 #ifndef ALICEO2_TRD_DETECTOR_H_
 #define ALICEO2_TRD_DETECTOR_H_
 
+#include <memory>
 #include <vector>
 #include "DetectorsBase/Detector.h"
 #include "SimulationDataFormat/BaseHits.h"
 #include "CommonUtils/ShmAllocator.h"
 
 #include "DataFormatsTRD/Hit.h"
+#include "TRDSimulation/IonisationModel.h"
 
 class FairVolume;
 
@@ -90,6 +92,17 @@ class Detector : public o2::base::DetImpl<Detector>
 
   float mWion; // Ionization potential
 
+  // Energy (eV) from which the charge of the current step is made
+  double ionisationEnergy(float edep) const;
+
+  // Checks once that the gas has the delta-ray cut of the collision tables
+  void checkDeltaRayCut();
+  bool mDeltaRayCutChecked = false; //!
+
+  std::unique_ptr<IonisationModel> mIonisation; //! collision model in the gas, null when the energy deposit is used
+  double mIonMinBetaGamma = 0.;                 //!
+  double mIonMinEkin = 0.;                      //! eV
+
   Geometry* mGeom = nullptr;
 
   // Volume-id lookup tables, resolved once at initialisation so that ProcessHits does

Detectors/TRD/simulation/include/TRDSimulation/IonisationModel.h

@@ -0,0 +1,61 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#ifndef ALICEO2_TRD_IONISATIONMODEL_H_
+#define ALICEO2_TRD_IONISATIONMODEL_H_
+
+#include "TRDSimulation/SimParam.h"
+
+class TRandom;
+
+namespace o2
+{
+namespace trd
+{
+
+/// Primary ionisation in the TRD gas from tabulated collision spectra, independent of the transport engine's
+/// energy-loss model. The tables hold the number of collisions per cm with an energy transfer below the delta-ray
+/// production cut of the gas (10 keV) and their energy distribution, as functions of beta*gamma, taken from
+/// Geant4's photo-absorption ionisation (PAI) cross-section for the Xe/CO2 and Ar/CO2 mixtures.
+class IonisationModel
+{
+ public:
+  explicit IonisationModel(SimParam::GasMixture gas);
+
+  /// Energy (eV) lost in collisions along a step of length (cm) at beta*gamma, for a particle of charge squared z2
+  double stepEnergy(double length, double betaGamma, double z2, TRandom& rnd) const;
+
+  /// Mean number of collisions per cm for unit charge
+  double collisionsPerCm(double betaGamma) const;
+
+  /// Mean energy loss (eV/cm) below the delta-ray cut for unit charge, from the tables
+  double restrictedDEdx(double betaGamma) const;
+
+  /// Lowest beta*gamma of the tables
+  static double minBetaGamma();
+
+  /// Delta-ray production cut (eV) the tables were made for
+  static double deltaRayCut();
+
+ private:
+  double sampleCollision(int row, double u) const;
+
+  int mNEnergy = 0;
+  const double* mEnergy = nullptr;       ///< energy grid (eV)
+  const double* mCollisions = nullptr;   ///< collisions per cm, per beta*gamma row
+  const double* mDEdx = nullptr;         ///< restricted dE/dx (eV/cm), per beta*gamma row
+  const float* mFractionAbove = nullptr; ///< fraction of collisions above each grid energy, per row
+};
+
+} // namespace trd
+} // namespace o2
+
+#endif

Detectors/TRD/simulation/include/TRDSimulation/TRDSimParams.h

@@ -37,6 +37,10 @@ struct TRDSimParams : public o2::conf::ConfigurableParamHelper<TRDSimParams> {
   float maxMCStepSize = 0.1;                              ///< maximum size of MC steps
   bool doTR = true;                                       ///< switch for transition radiation
   SimParam::GasMixture gas = SimParam::GasMixture::Xenon; ///< the gas mixture in the TRD
+  // ionisation in the gas
+  int ionisationModel = 0;             ///< charge from 0: Geant4 energy deposit with PAI, 1: TRD collision tables, no PAI
+  float ionisationMinBetaGamma = 0.04; ///< below this beta*gamma the charge comes from the energy deposit
+  float ionisationMinEkin = 30e3;      ///< (eV) below this kinetic energy the charge comes from the energy deposit
   // TRF parameters
   int trf = 1;                  // Sampled TRF function. 0: default TRF, 1: TRF described in TRF TDR, 2: No TRF, 3: Landau dist as TRF (parameters specified below)
   float trf_landau_mu = 0.;     // Mu of the Landau distribution used to describe TRF

Detectors/TRD/simulation/macros/makeIonisationTables.cxx

@@ -0,0 +1,137 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+// Writes the TRD gas collision tables (C++) from Geant4's PAI cross-section: number of collisions per cm with an
+// energy transfer below the delta-ray cut, restricted dE/dx, and N(>E)/N on a log energy grid, per beta*gamma.
+// Standalone program, not built with O2. Build: g++ -O2 $(geant4-config --cflags) makeIonisationTables.cxx -o makeIonisationTables $(geant4-config --libs)
+#include "G4Element.hh"
+#include "G4Material.hh"
+#include "G4PAIxSection.hh"
+#include "G4SandiaTable.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4Version.hh"
+#include <cmath>
+#include <cstdio>
+#include <string>
+#include <vector>
+
+constexpr int kNBetaGamma = 31;
+constexpr double kLogBetaGammaMin = -1.4;
+constexpr double kLogBetaGammaStep = 0.2;
+
+struct Gas {
+  std::string name;
+  double z, a, nNoble, nC, nO, fNoble, dNoble;
+};
+
+static void table(const Gas& gas, double cut)
+{
+  const double dco = 0.00186; // CO2 at 20 C, as in Detector::createMaterials
+  const double dens = (gas.fNoble * gas.dNoble + (1. - gas.fNoble) * dco) * g / cm3;
+  auto* el = new G4Element(gas.name, gas.name, gas.z, gas.a * g / mole);
+  auto* c = new G4Element("C_" + gas.name, "C", 6., 12.0107 * g / mole);
+  auto* o = new G4Element("O_" + gas.name, "O", 8., 15.9994 * g / mole);
+  const double w[3] = {gas.nNoble * gas.a, gas.nC * 12.0107, gas.nO * 15.9994};
+  const double ws = w[0] + w[1] + w[2];
+  auto* mat = new G4Material("TRD_" + gas.name + "CO2", dens, 3, kStateGas, 293.15 * kelvin);
+  mat->AddElement(el, w[0] / ws);
+  mat->AddElement(c, w[1] / ws);
+  mat->AddElement(o, w[2] / ws);
+
+  std::vector<double> lbgs;
+  for (int k = 0; k < kNBetaGamma; ++k) {
+    lbgs.push_back(kLogBetaGammaMin + kLogBetaGammaStep * k);
+  }
+  std::vector<double> egrid;
+  std::vector<std::vector<double>> ccdf;
+  std::vector<double> nsub, dedx;
+  for (double lbg : lbgs) {
+    const double bg = std::pow(10., lbg);
+    G4SandiaTable sandia;
+    sandia.Initialize(mat);
+    G4PAIxSection xs;
+    xs.Initialize(mat, 1. * MeV, bg * bg, &sandia);
+    const int n = xs.GetSplineSize();
+    std::vector<double> en(n), in(n), de(n);
+    for (int i = 1; i <= n; ++i) {
+      en[i - 1] = xs.GetSplineEnergy(i);
+      in[i - 1] = xs.GetIntegralPAIxSection(i);
+      de[i - 1] = xs.GetIntegralPAIdEdx(i);
+    }
+    auto interp = [&](const std::vector<double>& v, double e) {
+      if (e <= en.front()) {
+        return v.front();
+      }
+      for (int i = 1; i < n; ++i) {
+        if (en[i] >= e) {
+          const double t = std::log(e / en[i - 1]) / std::log(en[i] / en[i - 1]);
+          return v[i - 1] + t * (v[i] - v[i - 1]);
+        }
+      }
+      return v.back();
+    };
+    if (egrid.empty()) {
+      egrid.push_back(en[0] / eV);
+      for (int j = 0; j <= 30; ++j) {
+        const double e = std::pow(10., 1. + 0.1 * j);
+        if (e > en[0] / eV * 1.001) {
+          egrid.push_back(e);
+        }
+      }
+    }
+    const double ns = (in[0] - interp(in, cut)) * 10.; // per cm
+    nsub.push_back(ns);
+    dedx.push_back((de[0] - interp(de, cut)) * 10. / eV);
+    std::vector<double> row;
+    for (double e : egrid) {
+      row.push_back(std::max((interp(in, e * eV) - interp(in, cut)) * 10. / ns, 0.));
+    }
+    row.front() = 1.;
+    row.back() = 0.;
+    ccdf.push_back(row);
+  }
+  std::printf("// %s/CO2: density %.6g g/cm3, mean excitation energy %.4g eV\n", gas.name.c_str(), dens / (g / cm3),
+              mat->GetIonisation()->GetMeanExcitationEnergy() / eV);
+  std::printf("constexpr double k%sEnergy[%zu] = {", gas.name.c_str(), egrid.size());
+  for (size_t j = 0; j < egrid.size(); ++j) {
+    std::printf("%s%.5g", j ? ", " : "", egrid[j]);
+  }
+  std::printf("};\nconstexpr double k%sCollisionsPerCm[%zu] = {", gas.name.c_str(), nsub.size());
+  for (size_t k = 0; k < nsub.size(); ++k) {
+    std::printf("%s%.6g", k ? ", " : "", nsub[k]);
+  }
+  std::printf("};\nconstexpr double k%sRestrictedDEdx[%zu] = {", gas.name.c_str(), dedx.size());
+  for (size_t k = 0; k < dedx.size(); ++k) {
+    std::printf("%s%.6g", k ? ", " : "", dedx[k]);
+  }
+  std::printf("};\nconstexpr float k%sFractionAbove[%zu][%zu] = {\n", gas.name.c_str(), ccdf.size(), egrid.size());
+  for (auto& row : ccdf) {
+    std::printf("  {");
+    for (size_t j = 0; j < row.size(); ++j) {
+      std::printf("%s%.5g", j ? ", " : "", row[j]);
+    }
+    std::printf("},\n");
+  }
+  std::printf("};\n");
+}
+
+int main()
+{
+  const double cut = 10. * keV;
+  std::printf("// Generated with Geant4 %d G4PAIxSection, delta-ray cut %.4g keV\n", G4VERSION_NUMBER, cut / keV);
+  std::printf("constexpr double kDeltaRayCut = %.6g; // eV\n", cut / eV);
+  std::printf("constexpr int kNBetaGamma = %d;\n", kNBetaGamma);
+  std::printf("constexpr double kLogBetaGammaMin = %.6g; // beta*gamma of row k is 10^(min + k * step)\n", kLogBetaGammaMin);
+  std::printf("constexpr double kLogBetaGammaStep = %.6g;\n", kLogBetaGammaStep);
+  table({"Xe", 54., 131.29, 8.5, 1.5, 3.0, 0.85, 0.00549}, cut);
+  table({"Ar", 18., 39.948, 8.2, 1.8, 3.6, 0.82, 0.00166}, cut);
+  return 0;
+}

Detectors/TRD/simulation/src/Detector.cxx

@@ -18,13 +18,15 @@
 #include "DataFormatsTRD/Hit.h"
 
 #include "CommonUtils/ShmAllocator.h"
+#include "DetectorsBase/MaterialManager.h"
 #include "DetectorsBase/Stack.h"
 #include "FairVolume.h"
 #include "FairRootManager.h"
 
 #include <TGeoManager.h>
 #include <TVirtualMC.h>
 
+#include <cmath>
 #include <cstdio>
 #include <vector>
 #include <stdexcept>
@@ -120,6 +122,56 @@ void Detector::InitializeParams()
   }
   mTR = new TRsim();
   mMaxMCStepDef = TRDSimParams::Instance().maxMCStepSize;
+
+  const auto& simParams = TRDSimParams::Instance();
+  if (simParams.ionisationModel == 1) {
+    mIonisation = std::make_unique<IonisationModel>(simParams.gas);
+    // the tables start at minBetaGamma; below it the energy deposit is used
+    mIonMinBetaGamma = std::max<double>(simParams.ionisationMinBetaGamma, IonisationModel::minBetaGamma());
+    mIonMinEkin = simParams.ionisationMinEkin;
+    LOG(info) << "TRD: ionisation from collision tables, energy deposit below beta*gamma " << mIonMinBetaGamma
+              << " or kinetic energy " << mIonMinEkin << " eV";
+  } else if (simParams.ionisationModel != 0) {
+    LOG(fatal) << "TRD: unknown ionisationModel " << simParams.ionisationModel;
+  }
+}
+
+void Detector::checkDeltaRayCut()
+{
+  // The collision tables end at the delta-ray cut they were made for. Above it Geant4 makes delta electrons
+  // with their own hits, so a different cut in the gas (medium 9, Gas-mix) would double count or lose energy.
+  auto& mgr = o2::base::MaterialManager::Instance();
+  std::vector<std::pair<o2::base::ECut, Float_t>> cuts;
+  mgr.getSpecialCuts(GetName(), 9, cuts);
+  double dcute = -1.; // eV
+  for (const auto& [cut, value] : cuts) {
+    if (cut == o2::base::ECut::kDCUTE) {
+      dcute = value * 1e9;
+    }
+  }
+  const double tableCut = IonisationModel::deltaRayCut();
+  if (!mgr.specialCutsEnabled() || std::abs(dcute - tableCut) > 1e-3 * tableCut) {
+    LOG(fatal) << "TRD: the ionisation tables need DCUTE = " << tableCut << " eV in the gas, the medium has "
+               << (mgr.specialCutsEnabled() ? dcute : -1.) << " eV; use TRDSimParams.ionisationModel=0 with other cuts";
+  }
+  mDeltaRayCutChecked = true;
+}
+
+double Detector::ionisationEnergy(float edep) const
+{
+  double px, py, pz, etot;
+  fMC->TrackMomentum(px, py, pz, etot);
+  const double mass = fMC->TrackMass();
+  const double epre = etot + fMC->Edep(); // total energy before the step
+  const double ekin = (epre - mass) * 1e9;
+  const double betaGamma = std::sqrt(std::max(epre * epre - mass * mass, 0.)) / mass;
+  // slow tracks, including the step that stops a track, keep the energy deposit
+  if (betaGamma < mIonMinBetaGamma || ekin < mIonMinEkin) {
+    return edep;
+  }
+  // bare charge, no effective charge: an ion is close to bare for beta >> Z/137, for alphas above the beta*gamma cut
+  const double z = fMC->TrackCharge();
+  return std::min(mIonisation->stepEnergy(fMC->TrackStep(), betaGamma, z * z, *o2::base::simRandom()), ekin);
 }
 
 bool Detector::ProcessHits(FairVolume* v)
@@ -221,10 +273,14 @@ bool Detector::ProcessHits(FairVolume* v)
     trkStat = 2;
   }
 
-  // Calculate the charge according to GEANT Edep
+  // Calculate the charge from the collision model, or from the transport engine's energy deposit
   // Create a new dEdx hit
-  const float enDep = std::max(fMC->Edep(), 0.0) * 1e9; // Energy in eV
-  const int totalChargeDep = (int)(enDep / mWion);      // Total charge
+  const float edep = std::max(fMC->Edep(), 0.0) * 1e9; // Energy in eV
+  if (mIonisation && !mDeltaRayCutChecked) {
+    checkDeltaRayCut();
+  }
+  const float enDep = mIonisation ? ionisationEnergy(edep) : edep;
+  const int totalChargeDep = (int)(enDep / mWion); // Total charge
 
   // Store those hits with enDep bigger than the ionization potential of the gas mixture for in-flight tracks
   // or store hits of tracks that are entering or exiting