From ef2b628466991ec88d1c69c4b8241598f334ed91 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Sun, 15 Nov 2026 15:01:30 +0100 Subject: [PATCH 1/2] Add a TRD ionisation model independent of the Geant4 energy loss This adds a collision model for the ionisation in the TRD gas, selected with TRDSimParams.ionisationModel, and a unit test. - With ionisationModel=0, the default, the charge of a hit stays the Geant4 energy deposit divided by Wion. - With ionisationModel=1 the charge comes from a Poisson number of collisions along the step, from the step length, beta*gamma before the step and the charge squared. - The collision rate and energy spectrum below the 10 keV delta-ray cut are tabulated from the Geant4 PAI cross-section for Xe/CO2 and Ar/CO2, with macros/makeIonisationTables.cxx. - Tracks below beta*gamma 0.04 or 30 keV kinetic energy keep the energy deposit, including the step that stops a track. - The run stops if the delta-ray cut of the gas differs from the one of the tables. - Hit class and hit creation points are unchanged. Co-Authored-By: Claude Opus 5.5 --- Detectors/TRD/simulation/CMakeLists.txt | 7 + Detectors/TRD/simulation/data/simcuts.dat | 2 +- .../include/TRDSimulation/Detector.h | 13 ++ .../include/TRDSimulation/IonisationModel.h | 61 +++++ .../include/TRDSimulation/TRDSimParams.h | 4 + .../macros/makeIonisationTables.cxx | 137 ++++++++++++ Detectors/TRD/simulation/src/Detector.cxx | 62 +++++- .../TRD/simulation/src/IonisationModel.cxx | 210 ++++++++++++++++++ .../test/testTRDIonisationModel.cxx | 105 +++++++++ 9 files changed, 597 insertions(+), 4 deletions(-) create mode 100644 Detectors/TRD/simulation/include/TRDSimulation/IonisationModel.h create mode 100644 Detectors/TRD/simulation/macros/makeIonisationTables.cxx create mode 100644 Detectors/TRD/simulation/src/IonisationModel.cxx create mode 100644 Detectors/TRD/simulation/test/testTRDIonisationModel.cxx diff --git a/Detectors/TRD/simulation/CMakeLists.txt b/Detectors/TRD/simulation/CMakeLists.txt index decfa46836f86..9dd7162be1188 100644 --- a/Detectors/TRD/simulation/CMakeLists.txt +++ b/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 diff --git a/Detectors/TRD/simulation/data/simcuts.dat b/Detectors/TRD/simulation/data/simcuts.dat index 6f85d9d9726e2..b32179d2b5955 100644 --- a/Detectors/TRD/simulation/data/simcuts.dat +++ b/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 * diff --git a/Detectors/TRD/simulation/include/TRDSimulation/Detector.h b/Detectors/TRD/simulation/include/TRDSimulation/Detector.h index 9a10af8462e1f..db1e171143fb0 100644 --- a/Detectors/TRD/simulation/include/TRDSimulation/Detector.h +++ b/Detectors/TRD/simulation/include/TRDSimulation/Detector.h @@ -12,12 +12,14 @@ #ifndef ALICEO2_TRD_DETECTOR_H_ #define ALICEO2_TRD_DETECTOR_H_ +#include #include #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 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 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 diff --git a/Detectors/TRD/simulation/include/TRDSimulation/IonisationModel.h b/Detectors/TRD/simulation/include/TRDSimulation/IonisationModel.h new file mode 100644 index 0000000000000..419a9561c6e37 --- /dev/null +++ b/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 diff --git a/Detectors/TRD/simulation/include/TRDSimulation/TRDSimParams.h b/Detectors/TRD/simulation/include/TRDSimulation/TRDSimParams.h index 37b59b6ab8071..a187d41d058a5 100644 --- a/Detectors/TRD/simulation/include/TRDSimulation/TRDSimParams.h +++ b/Detectors/TRD/simulation/include/TRDSimulation/TRDSimParams.h @@ -37,6 +37,10 @@ struct TRDSimParams : public o2::conf::ConfigurableParamHelper { 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 diff --git a/Detectors/TRD/simulation/macros/makeIonisationTables.cxx b/Detectors/TRD/simulation/macros/makeIonisationTables.cxx new file mode 100644 index 0000000000000..b2aff54199beb --- /dev/null +++ b/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 +#include +#include +#include + +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 lbgs; + for (int k = 0; k < kNBetaGamma; ++k) { + lbgs.push_back(kLogBetaGammaMin + kLogBetaGammaStep * k); + } + std::vector egrid; + std::vector> ccdf; + std::vector 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 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& 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 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; +} diff --git a/Detectors/TRD/simulation/src/Detector.cxx b/Detectors/TRD/simulation/src/Detector.cxx index 5c951e0abbc25..6760da57f2d52 100644 --- a/Detectors/TRD/simulation/src/Detector.cxx +++ b/Detectors/TRD/simulation/src/Detector.cxx @@ -18,6 +18,7 @@ #include "DataFormatsTRD/Hit.h" #include "CommonUtils/ShmAllocator.h" +#include "DetectorsBase/MaterialManager.h" #include "DetectorsBase/Stack.h" #include "FairVolume.h" #include "FairRootManager.h" @@ -25,6 +26,7 @@ #include #include +#include #include #include #include @@ -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(simParams.gas); + // the tables start at minBetaGamma; below it the energy deposit is used + mIonMinBetaGamma = std::max(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> 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 diff --git a/Detectors/TRD/simulation/src/IonisationModel.cxx b/Detectors/TRD/simulation/src/IonisationModel.cxx new file mode 100644 index 0000000000000..d68085428afef --- /dev/null +++ b/Detectors/TRD/simulation/src/IonisationModel.cxx @@ -0,0 +1,210 @@ +// 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. + +#include "TRDSimulation/IonisationModel.h" + +#include +#include + +#include +#include + +namespace +{ +// Collision tables per gas, output of macros/makeIonisationTables.cxx. N(>E)/N is the fraction of collisions with an +// energy transfer above the grid energy. +// Generated with Geant4 1120 G4PAIxSection, delta-ray cut 10 keV +constexpr double kDeltaRayCut = 10000; // eV +constexpr int kNBetaGamma = 31; +constexpr double kLogBetaGammaMin = -1.4; // beta*gamma of row k is 10^(min + k * step) +constexpr double kLogBetaGammaStep = 0.2; +// Xe/CO2: density 0.0049455 g/cm3, mean excitation energy 431 eV +constexpr double kXeEnergy[31] = {12.191, 12.589, 15.849, 19.953, 25.119, 31.623, 39.811, 50.119, 63.096, 79.433, 100, 125.89, 158.49, 199.53, 251.19, 316.23, 398.11, 501.19, 630.96, 794.33, 1000, 1258.9, 1584.9, 1995.3, 2511.9, 3162.3, 3981.1, 5011.9, 6309.6, 7943.3, 10000}; +constexpr double kXeCollisionsPerCm[31] = {3169.21, 2138.01, 1335.37, 766.94, 406.804, 206.947, 109.762, 67.1239, 50.0173, 44.2345, 43.3844, 44.6683, 46.8794, 49.5104, 52.2097, 54.4158, 56.0136, 57.159, 57.9521, 58.4512, 58.7541, 58.9436, 59.0677, 59.1368, 59.1718, 59.1889, 59.197, 59.2006, 59.2022, 59.2028, 59.2031}; +constexpr double kXeRestrictedDEdx[31] = {292507, 198636, 124651, 71843.3, 38210.5, 19478.9, 10347.8, 6336.06, 4726.44, 4184.25, 4107.61, 4232.18, 4443.12, 4690.79, 4951.45, 5208.33, 5454.25, 5685.98, 5897.58, 6085.11, 6249.98, 6391.66, 6508.87, 6591.46, 6645.77, 6680.04, 6700.21, 6710.87, 6715.9, 6718.08, 6718.98}; +constexpr float kXeFractionAbove[31][31] = { + {1, 0.97153, 0.72262, 0.52171, 0.40152, 0.33137, 0.2899, 0.26152, 0.23989, 0.21232, 0.15153, 0.093752, 0.07254, 0.060813, 0.050836, 0.0421, 0.034771, 0.028948, 0.024297, 0.018761, 0.014504, 0.01119, 0.0085655, 0.0065282, 0.0048584, 0.00352, 0.0024333, 0.0015685, 0.00095822, 0.00044041, 0}, + {1, 0.97156, 0.72276, 0.52206, 0.40274, 0.33398, 0.29405, 0.26712, 0.24668, 0.21917, 0.15457, 0.092924, 0.072114, 0.061525, 0.052307, 0.043959, 0.036831, 0.03117, 0.026624, 0.019912, 0.014901, 0.011085, 0.0082296, 0.006147, 0.0045158, 0.003264, 0.0022793, 0.0015009, 0.00089382, 0.00040144, 0}, + {1, 0.97158, 0.72286, 0.5223, 0.40363, 0.33586, 0.29706, 0.27117, 0.2516, 0.22413, 0.15676, 0.092324, 0.071805, 0.062039, 0.05337, 0.045303, 0.03832, 0.032775, 0.028305, 0.020743, 0.015187, 0.011008, 0.0079862, 0.0058711, 0.004268, 0.0030788, 0.002168, 0.0014519, 0.00084722, 0.00037325, 0}, + {1, 0.9716, 0.72293, 0.52249, 0.40429, 0.33727, 0.29932, 0.27422, 0.25529, 0.22785, 0.15841, 0.091871, 0.071571, 0.062425, 0.054169, 0.046314, 0.039441, 0.033984, 0.029571, 0.021369, 0.015401, 0.01095, 0.0078025, 0.0056632, 0.0040812, 0.0029393, 0.0020841, 0.001415, 0.00081212, 0.00035202, 0}, + {1, 0.97161, 0.72299, 0.52263, 0.40481, 0.33837, 0.30106, 0.27657, 0.25815, 0.23073, 0.15969, 0.091521, 0.071391, 0.062723, 0.054786, 0.047093, 0.040305, 0.034915, 0.030546, 0.021851, 0.015567, 0.010904, 0.0076608, 0.0055028, 0.0039372, 0.0028317, 0.0020194, 0.0013866, 0.00078504, 0.00033565, 0}, + {1, 0.97162, 0.72303, 0.52274, 0.4052, 0.33921, 0.30241, 0.27839, 0.26036, 0.23296, 0.16067, 0.091249, 0.071251, 0.062953, 0.055263, 0.047698, 0.040974, 0.035637, 0.031302, 0.022224, 0.015695, 0.010869, 0.0075507, 0.0053783, 0.0038255, 0.0027482, 0.0019692, 0.0013645, 0.00076404, 0.00032296, 0}, + {1, 0.97163, 0.72307, 0.52282, 0.40551, 0.33986, 0.30345, 0.27979, 0.26206, 0.23467, 0.16143, 0.09104, 0.071143, 0.063131, 0.055631, 0.048163, 0.04149, 0.036193, 0.031884, 0.022511, 0.015793, 0.010842, 0.0074659, 0.0052824, 0.0037395, 0.002684, 0.0019306, 0.0013475, 0.00074787, 0.00031318, 0}, + {1, 0.97163, 0.72308, 0.52288, 0.40574, 0.34038, 0.30427, 0.2809, 0.2634, 0.23602, 0.16203, 0.090875, 0.071058, 0.06327, 0.05592, 0.048529, 0.041896, 0.036631, 0.032342, 0.022738, 0.01587, 0.01082, 0.0073991, 0.0052068, 0.0036717, 0.0026334, 0.0019001, 0.0013341, 0.00073512, 0.00030547, 0}, + {1, 0.97163, 0.72308, 0.52292, 0.40595, 0.34082, 0.30498, 0.28186, 0.26457, 0.23721, 0.16255, 0.090729, 0.070982, 0.063392, 0.056173, 0.048849, 0.042251, 0.037014, 0.032744, 0.022936, 0.015938, 0.010801, 0.0073403, 0.0051404, 0.0036121, 0.0025889, 0.0018734, 0.0013223, 0.00072392, 0.00029871, 0}, + {1, 0.97162, 0.72303, 0.52293, 0.40612, 0.34124, 0.30566, 0.28279, 0.26569, 0.23834, 0.16305, 0.090585, 0.070907, 0.063507, 0.056416, 0.049157, 0.042593, 0.037383, 0.03313, 0.023127, 0.016003, 0.010782, 0.0072834, 0.0050763, 0.0035547, 0.002546, 0.0018475, 0.0013109, 0.00071312, 0.00029219, 0}, + {1, 0.97157, 0.72288, 0.52288, 0.40626, 0.34161, 0.30629, 0.28365, 0.26675, 0.23941, 0.16352, 0.090442, 0.070831, 0.063614, 0.056645, 0.049449, 0.042917, 0.037734, 0.033498, 0.023307, 0.016064, 0.010764, 0.0072286, 0.0050147, 0.0034995, 0.0025048, 0.0018227, 0.0013, 0.00070274, 0.00028593, 0}, + {1, 0.97146, 0.7225, 0.52267, 0.40629, 0.3419, 0.30682, 0.28441, 0.26768, 0.24035, 0.16392, 0.090298, 0.070752, 0.063704, 0.056848, 0.049711, 0.04321, 0.038051, 0.033832, 0.023471, 0.016118, 0.010745, 0.0071772, 0.0049573, 0.0034482, 0.0024666, 0.0017997, 0.0012897, 0.00069309, 0.00028013, 0}, + {1, 0.97115, 0.72156, 0.52212, 0.40612, 0.34202, 0.30719, 0.28498, 0.26843, 0.24112, 0.16421, 0.090136, 0.070658, 0.063765, 0.057013, 0.049931, 0.04346, 0.038325, 0.034122, 0.02361, 0.016161, 0.010725, 0.0071287, 0.0049041, 0.0034011, 0.0024314, 0.0017784, 0.0012801, 0.0006842, 0.00027482, 0}, + {1, 0.97023, 0.71935, 0.52083, 0.40553, 0.34185, 0.30729, 0.2853, 0.26891, 0.24162, 0.16433, 0.08993, 0.07053, 0.063779, 0.057124, 0.050096, 0.043656, 0.038546, 0.034359, 0.02372, 0.01619, 0.010701, 0.007081, 0.0048537, 0.0033569, 0.0023985, 0.0017583, 0.0012707, 0.00067585, 0.00026992, 0}, + {1, 0.96868, 0.71643, 0.5196, 0.40534, 0.34217, 0.30787, 0.28605, 0.26978, 0.24245, 0.16466, 0.089894, 0.070542, 0.063907, 0.057328, 0.050336, 0.043912, 0.038814, 0.034635, 0.023862, 0.016246, 0.010699, 0.0070516, 0.0048176, 0.0033236, 0.0023736, 0.0017435, 0.0012647, 0.00066964, 0.00026605, 0}, + {1, 0.96816, 0.71678, 0.52191, 0.40853, 0.34563, 0.31137, 0.28953, 0.27321, 0.24554, 0.16646, 0.09075, 0.071274, 0.064679, 0.058103, 0.051072, 0.044597, 0.039459, 0.035243, 0.024238, 0.016464, 0.010807, 0.0070977, 0.0048346, 0.0033278, 0.0023756, 0.0017482, 0.0012723, 0.00067093, 0.00026528, 0}, + {1, 0.96831, 0.7197, 0.52714, 0.4147, 0.352, 0.31766, 0.29567, 0.27915, 0.25075, 0.16957, 0.092488, 0.072744, 0.066121, 0.059476, 0.052333, 0.045739, 0.040505, 0.036208, 0.024861, 0.016854, 0.011031, 0.0072214, 0.0049055, 0.0033696, 0.0024044, 0.0017725, 0.0012939, 0.00067977, 0.00026758, 0}, + {1, 0.96864, 0.72318, 0.533, 0.42166, 0.35932, 0.32497, 0.30283, 0.28605, 0.25657, 0.173, 0.094709, 0.074682, 0.067995, 0.061243, 0.053942, 0.047187, 0.041823, 0.037414, 0.025652, 0.01736, 0.011333, 0.0073974, 0.0050127, 0.0034367, 0.0024514, 0.00181, 0.0013249, 0.00069371, 0.00027195, 0}, + {1, 0.96895, 0.72611, 0.53792, 0.42759, 0.36567, 0.33142, 0.30922, 0.29219, 0.2614, 0.17609, 0.097137, 0.076905, 0.070146, 0.063272, 0.055794, 0.048853, 0.043337, 0.038794, 0.026564, 0.01795, 0.011692, 0.0076122, 0.0051468, 0.0035226, 0.0025118, 0.0018574, 0.0013631, 0.00071145, 0.00027785, 0}, + {1, 0.96917, 0.72811, 0.54127, 0.43167, 0.3701, 0.33597, 0.31377, 0.29655, 0.26482, 0.17875, 0.099594, 0.079256, 0.07244, 0.065454, 0.0578, 0.050671, 0.044992, 0.040301, 0.027563, 0.018604, 0.012096, 0.0078585, 0.0053031, 0.0036241, 0.0025835, 0.001913, 0.0014072, 0.00073238, 0.00028501, 0}, + {1, 0.96931, 0.72937, 0.5434, 0.43428, 0.37295, 0.33892, 0.31675, 0.29944, 0.26728, 0.181, 0.10187, 0.081495, 0.074647, 0.067582, 0.05979, 0.0525, 0.046671, 0.041828, 0.028577, 0.019279, 0.012521, 0.0081222, 0.0054728, 0.0037355, 0.0026625, 0.0019741, 0.0014554, 0.00075542, 0.00029303, 0}, + {1, 0.9694, 0.73019, 0.54478, 0.43597, 0.37481, 0.34086, 0.31872, 0.30138, 0.26907, 0.18281, 0.10376, 0.083399, 0.076539, 0.069434, 0.061557, 0.054159, 0.048215, 0.043219, 0.029509, 0.019917, 0.012936, 0.0083866, 0.0056463, 0.0038511, 0.0027451, 0.002038, 0.0015056, 0.00077961, 0.0003015, 0}, + {1, 0.96946, 0.73073, 0.5457, 0.43711, 0.37606, 0.34217, 0.32006, 0.30272, 0.27038, 0.18419, 0.10525, 0.0849, 0.07804, 0.070919, 0.063004, 0.055548, 0.049536, 0.044355, 0.030253, 0.020452, 0.013303, 0.0086328, 0.0058136, 0.0039656, 0.0028281, 0.0021026, 0.0015564, 0.00080422, 0.00031018, 0}, + {1, 0.9695, 0.73104, 0.54621, 0.43774, 0.37676, 0.34291, 0.32082, 0.30349, 0.27114, 0.18501, 0.10613, 0.085798, 0.07894, 0.071815, 0.063883, 0.056402, 0.050356, 0.045071, 0.030767, 0.020846, 0.013593, 0.0088419, 0.005964, 0.0040736, 0.0029086, 0.0021661, 0.0016065, 0.00082875, 0.00031893, 0}, + {1, 0.96952, 0.7312, 0.54648, 0.43807, 0.37712, 0.34328, 0.3212, 0.30388, 0.27154, 0.18544, 0.1066, 0.08627, 0.079415, 0.072288, 0.064351, 0.05686, 0.050799, 0.045468, 0.031078, 0.0211, 0.013794, 0.0089988, 0.006085, 0.0041663, 0.0029809, 0.0022246, 0.0016531, 0.00085204, 0.00032738, 0}, + {1, 0.96953, 0.73127, 0.5466, 0.43823, 0.3773, 0.34347, 0.32139, 0.30407, 0.27173, 0.18565, 0.10683, 0.086507, 0.079652, 0.072525, 0.064586, 0.057091, 0.051024, 0.045675, 0.03125, 0.021247, 0.013919, 0.0091022, 0.0061704, 0.0042361, 0.0030384, 0.0022728, 0.0016923, 0.00087224, 0.00033495, 0}, + {1, 0.96953, 0.73131, 0.54667, 0.4383, 0.37738, 0.34356, 0.32148, 0.30416, 0.27183, 0.18575, 0.10694, 0.08662, 0.079766, 0.07264, 0.0647, 0.057204, 0.051134, 0.045777, 0.031338, 0.021326, 0.013988, 0.0091625, 0.0062227, 0.0042812, 0.0030772, 0.0023064, 0.0017202, 0.00088737, 0.00034092, 0}, + {1, 0.96953, 0.73132, 0.54669, 0.43834, 0.37742, 0.34359, 0.32153, 0.30421, 0.27187, 0.1858, 0.10699, 0.086672, 0.079819, 0.072692, 0.064752, 0.057255, 0.051185, 0.045825, 0.031381, 0.021365, 0.014022, 0.0091936, 0.0062504, 0.0043058, 0.003099, 0.0023259, 0.0017366, 0.00089673, 0.00034483, 0}, + {1, 0.96953, 0.73133, 0.5467, 0.43835, 0.37743, 0.34361, 0.32154, 0.30422, 0.27189, 0.18582, 0.10701, 0.086695, 0.079841, 0.072715, 0.064775, 0.057278, 0.051207, 0.045846, 0.0314, 0.021382, 0.014038, 0.009208, 0.0062634, 0.0043175, 0.0031096, 0.0023354, 0.0017447, 0.00090155, 0.00034693, 0}, + {1, 0.96953, 0.73133, 0.54671, 0.43836, 0.37744, 0.34362, 0.32155, 0.30423, 0.2719, 0.18583, 0.10702, 0.086704, 0.079851, 0.072724, 0.064784, 0.057287, 0.051216, 0.045855, 0.031408, 0.021389, 0.014045, 0.0092141, 0.006269, 0.0043226, 0.0031142, 0.0023396, 0.0017483, 0.00090374, 0.00034791, 0}, + {1, 0.96953, 0.73134, 0.54671, 0.43836, 0.37744, 0.34362, 0.32155, 0.30423, 0.2719, 0.18583, 0.10703, 0.086708, 0.079855, 0.072728, 0.064788, 0.057291, 0.05122, 0.045858, 0.031411, 0.021392, 0.014048, 0.0092167, 0.0062713, 0.0043247, 0.0031161, 0.0023413, 0.0017499, 0.00090466, 0.00034833, 0}, +}; +// Ar/CO2: density 0.001696 g/cm3, mean excitation energy 161.2 eV +constexpr double kArEnergy[30] = {13.688, 15.849, 19.953, 25.119, 31.623, 39.811, 50.119, 63.096, 79.433, 100, 125.89, 158.49, 199.53, 251.19, 316.23, 398.11, 501.19, 630.96, 794.33, 1000, 1258.9, 1584.9, 1995.3, 2511.9, 3162.3, 3981.1, 5011.9, 6309.6, 7943.3, 10000}; +constexpr double kArCollisionsPerCm[31] = {3724.53, 2226.97, 1197.5, 588.316, 276.252, 131.995, 68.6451, 41.8075, 31.1234, 27.5122, 26.973, 27.7608, 29.1231, 30.7457, 32.5214, 34.4369, 35.8399, 36.6344, 37.0882, 37.3632, 37.5381, 37.6497, 37.7155, 37.7489, 37.765, 37.7726, 37.7761, 37.7778, 37.7785, 37.7788, 37.7789}; +constexpr double kArRestrictedDEdx[31] = {293869, 169704, 88909.5, 42809.7, 19786.1, 9336.67, 4808.6, 2906.03, 2148.6, 1886.66, 1837.75, 1879.93, 1960.88, 2058.29, 2162.73, 2269.49, 2365.28, 2447.45, 2519.75, 2583.85, 2638.98, 2683.33, 2715.75, 2737.35, 2751.56, 2760.89, 2766.78, 2770.17, 2771.91, 2772.71, 2773.05}; +constexpr float kArFractionAbove[31][30] = { + {1, 0.98121, 0.75234, 0.51, 0.32696, 0.22817, 0.18328, 0.15014, 0.12279, 0.10098, 0.083992, 0.070972, 0.061105, 0.052883, 0.039929, 0.030061, 0.022953, 0.017544, 0.013418, 0.010354, 0.0079634, 0.0060835, 0.0045964, 0.0034072, 0.002458, 0.0017487, 0.0011717, 0.00069875, 0.00031243, 0}, + {1, 0.98109, 0.74989, 0.50468, 0.31972, 0.2213, 0.17835, 0.14688, 0.1208, 0.099974, 0.083815, 0.071531, 0.062338, 0.054525, 0.04026, 0.029495, 0.021971, 0.016357, 0.012216, 0.0092777, 0.0070637, 0.0053708, 0.0040592, 0.0030256, 0.0022092, 0.0015475, 0.001024, 0.00060526, 0.00026891, 0}, + {1, 0.981, 0.74811, 0.50081, 0.31446, 0.21632, 0.17477, 0.14451, 0.11936, 0.099244, 0.083684, 0.071935, 0.063231, 0.055714, 0.040498, 0.029082, 0.021258, 0.015496, 0.011343, 0.0084965, 0.0064107, 0.0048535, 0.0036693, 0.0027486, 0.0020286, 0.0014014, 0.00091679, 0.00053741, 0.00023732, 0}, + {1, 0.98093, 0.74677, 0.4979, 0.31049, 0.21255, 0.17206, 0.14272, 0.11827, 0.098692, 0.083585, 0.072239, 0.063904, 0.056611, 0.040677, 0.02877, 0.020718, 0.014845, 0.010684, 0.0079066, 0.0059176, 0.0044628, 0.0033748, 0.0025394, 0.0018923, 0.0012912, 0.00083585, 0.00048617, 0.00021347, 0}, + {1, 0.98088, 0.74574, 0.49564, 0.30742, 0.20964, 0.16997, 0.14134, 0.11742, 0.098265, 0.083507, 0.072473, 0.064423, 0.057304, 0.040815, 0.028528, 0.0203, 0.014342, 0.010174, 0.0074508, 0.0055366, 0.004161, 0.0031473, 0.0023777, 0.0017869, 0.001206, 0.0007733, 0.00044658, 0.00019504, 0}, + {1, 0.98084, 0.74493, 0.49389, 0.30504, 0.20738, 0.16834, 0.14026, 0.11677, 0.097933, 0.083447, 0.072654, 0.064826, 0.057842, 0.040921, 0.028339, 0.019976, 0.013951, 0.0097788, 0.0070968, 0.0052406, 0.0039266, 0.0029706, 0.0022522, 0.001705, 0.0011398, 0.00072473, 0.00041584, 0.00018073, 0}, + {1, 0.98081, 0.74431, 0.49254, 0.3032, 0.20564, 0.16709, 0.13943, 0.11626, 0.097676, 0.0834, 0.072794, 0.065137, 0.058256, 0.041003, 0.028194, 0.019725, 0.013649, 0.0094739, 0.0068239, 0.0050126, 0.0037459, 0.0028344, 0.0021554, 0.0016419, 0.0010888, 0.00068729, 0.00039214, 0.0001697, 0}, + {1, 0.98079, 0.74381, 0.49147, 0.30175, 0.20426, 0.1661, 0.13878, 0.11586, 0.097474, 0.083363, 0.072904, 0.065381, 0.058582, 0.041068, 0.028079, 0.019528, 0.013411, 0.0092333, 0.0066087, 0.0048326, 0.0036034, 0.002727, 0.0020791, 0.0015921, 0.0010485, 0.00065776, 0.00037345, 0.000161, 0}, + {1, 0.98076, 0.74337, 0.49052, 0.30047, 0.20305, 0.16523, 0.1382, 0.11551, 0.097295, 0.083329, 0.072999, 0.065595, 0.058868, 0.041124, 0.027978, 0.019354, 0.013202, 0.009022, 0.0064196, 0.0046746, 0.0034782, 0.0026326, 0.002012, 0.0015484, 0.0010132, 0.00063182, 0.00035704, 0.00015335, 0}, + {1, 0.98074, 0.74291, 0.48958, 0.29923, 0.20188, 0.16439, 0.13764, 0.11517, 0.097121, 0.083295, 0.07309, 0.065801, 0.059143, 0.041177, 0.02788, 0.019186, 0.013001, 0.0088181, 0.0062372, 0.0045222, 0.0033574, 0.0025416, 0.0019473, 0.0015062, 0.0009791, 0.0006068, 0.0003412, 0.00014598, 0}, + {1, 0.9807, 0.7424, 0.48862, 0.29801, 0.20075, 0.16358, 0.1371, 0.11484, 0.09695, 0.083259, 0.073174, 0.065994, 0.059403, 0.041226, 0.027785, 0.019025, 0.012807, 0.0086226, 0.0060623, 0.004376, 0.0032417, 0.0024543, 0.0018853, 0.0014658, 0.00094641, 0.00058282, 0.00032603, 0.00013892, 0}, + {1, 0.98065, 0.74177, 0.48761, 0.29682, 0.19968, 0.16281, 0.13659, 0.11452, 0.096783, 0.083218, 0.073244, 0.066166, 0.059638, 0.041268, 0.027695, 0.018874, 0.012627, 0.0084416, 0.0059006, 0.004241, 0.0031347, 0.0023737, 0.001828, 0.0014284, 0.00091621, 0.00056067, 0.00031202, 0.0001324, 0}, + {1, 0.98054, 0.74081, 0.48637, 0.2956, 0.19864, 0.16207, 0.13609, 0.1142, 0.096609, 0.08316, 0.07329, 0.066305, 0.059834, 0.041295, 0.027607, 0.018734, 0.012462, 0.0082767, 0.0057536, 0.0041184, 0.0030377, 0.0023006, 0.001776, 0.0013944, 0.00088882, 0.00054061, 0.00029934, 0.00012651, 0}, + {1, 0.98028, 0.73893, 0.48448, 0.29413, 0.19755, 0.16129, 0.13556, 0.11385, 0.096392, 0.083058, 0.073286, 0.066389, 0.059973, 0.041294, 0.027512, 0.018598, 0.012308, 0.0081245, 0.005619, 0.0040065, 0.0029493, 0.002234, 0.0017285, 0.0013632, 0.00086383, 0.00052237, 0.00028784, 0.00012116, 0}, + {1, 0.97927, 0.73454, 0.48092, 0.29193, 0.19616, 0.1603, 0.13484, 0.11334, 0.096043, 0.082834, 0.073167, 0.066359, 0.060002, 0.041228, 0.027384, 0.018447, 0.012151, 0.0079763, 0.0054903, 0.0039006, 0.002866, 0.0021711, 0.0016835, 0.0013333, 0.00084022, 0.00050531, 0.00027714, 0.00011622, 0}, + {1, 0.97113, 0.72438, 0.47406, 0.28827, 0.19412, 0.15885, 0.13376, 0.11254, 0.095446, 0.082393, 0.07285, 0.066142, 0.059856, 0.041051, 0.027191, 0.018261, 0.011977, 0.0078218, 0.0053602, 0.0037953, 0.0027836, 0.002109, 0.0016387, 0.0013029, 0.00081682, 0.00048869, 0.00026683, 0.00011148, 0}, + {1, 0.96445, 0.71843, 0.47112, 0.28776, 0.19466, 0.15964, 0.13463, 0.1134, 0.096277, 0.083191, 0.073626, 0.066914, 0.060598, 0.041491, 0.027416, 0.01836, 0.011996, 0.0077973, 0.0053216, 0.003756, 0.0027502, 0.0020839, 0.0016223, 0.0012946, 0.0008077, 0.00048086, 0.00026144, 0.00010884, 0}, + {1, 0.96209, 0.71747, 0.47222, 0.29025, 0.19763, 0.16259, 0.13741, 0.11595, 0.098583, 0.085285, 0.075559, 0.068737, 0.062288, 0.042584, 0.02808, 0.018759, 0.012214, 0.0079049, 0.0053746, 0.0037822, 0.002765, 0.0020954, 0.0016342, 0.0013086, 0.00081272, 0.00048158, 0.00026077, 0.00010819, 0}, + {1, 0.96132, 0.71807, 0.4744, 0.2935, 0.20125, 0.16619, 0.14087, 0.11917, 0.10152, 0.087973, 0.078041, 0.071066, 0.064431, 0.04399, 0.028961, 0.019309, 0.012533, 0.0080808, 0.0054752, 0.0038424, 0.0028049, 0.0021259, 0.0016608, 0.0013343, 0.0008251, 0.00048675, 0.00026254, 0.00010856, 0}, + {1, 0.96111, 0.71901, 0.47655, 0.29649, 0.20459, 0.16958, 0.14419, 0.12235, 0.10452, 0.090767, 0.080653, 0.07353, 0.066682, 0.045478, 0.029915, 0.019918, 0.012898, 0.0082889, 0.005599, 0.0039195, 0.0028574, 0.0021659, 0.0016949, 0.0013658, 0.0008412, 0.00049416, 0.00026553, 0.00010945, 0}, + {1, 0.9611, 0.71987, 0.47831, 0.29888, 0.20727, 0.17233, 0.14694, 0.12506, 0.10714, 0.093277, 0.083043, 0.075806, 0.068738, 0.046859, 0.030833, 0.020524, 0.013271, 0.0085089, 0.0057334, 0.0040051, 0.0029165, 0.0022111, 0.001733, 0.0014005, 0.00085932, 0.00050279, 0.0002692, 0.00011062, 0}, + {1, 0.96114, 0.72052, 0.47958, 0.30061, 0.20921, 0.17434, 0.14898, 0.12709, 0.10914, 0.095234, 0.084941, 0.077636, 0.070309, 0.047964, 0.031612, 0.021066, 0.01362, 0.0087246, 0.0058698, 0.0040943, 0.0029791, 0.0022592, 0.0017734, 0.0014372, 0.00087868, 0.00051217, 0.00027329, 0.00011196, 0}, + {1, 0.96117, 0.72094, 0.48038, 0.30169, 0.21043, 0.1756, 0.15026, 0.12839, 0.11044, 0.096516, 0.086199, 0.078862, 0.071305, 0.048715, 0.032182, 0.021489, 0.013912, 0.0089184, 0.0059994, 0.0041827, 0.0030427, 0.0023088, 0.0018151, 0.0014748, 0.00089878, 0.00052205, 0.00027766, 0.00011343, 0}, + {1, 0.9612, 0.72115, 0.48079, 0.30225, 0.21107, 0.17626, 0.15094, 0.12907, 0.11112, 0.097195, 0.086869, 0.079517, 0.071861, 0.049162, 0.032544, 0.021776, 0.014127, 0.0090736, 0.0061114, 0.0042641, 0.003104, 0.0023577, 0.0018567, 0.0015123, 0.0009191, 0.00053218, 0.00028221, 0.00011498, 0}, + {1, 0.96121, 0.72126, 0.481, 0.30253, 0.21138, 0.17658, 0.15127, 0.1294, 0.11146, 0.097531, 0.087202, 0.079845, 0.072148, 0.049405, 0.032751, 0.02195, 0.014267, 0.0091834, 0.0061975, 0.0043317, 0.0031581, 0.0024025, 0.0018956, 0.0015476, 0.00093862, 0.00054212, 0.00028674, 0.00011654, 0}, + {1, 0.96121, 0.72131, 0.4811, 0.30266, 0.21152, 0.17674, 0.15142, 0.12956, 0.11162, 0.097693, 0.087363, 0.080004, 0.072291, 0.049531, 0.032863, 0.022048, 0.01435, 0.0092537, 0.0062567, 0.0043817, 0.0032006, 0.0024394, 0.0019285, 0.0015778, 0.00095596, 0.00055126, 0.00029102, 0.00011805, 0}, + {1, 0.96122, 0.72134, 0.48114, 0.30272, 0.2116, 0.17681, 0.1515, 0.12964, 0.1117, 0.09777, 0.087441, 0.080081, 0.072361, 0.049594, 0.032921, 0.022101, 0.014397, 0.009295, 0.0062934, 0.0044142, 0.0032296, 0.0024656, 0.0019524, 0.0016001, 0.00096951, 0.00055882, 0.00029473, 0.00011941, 0}, + {1, 0.96122, 0.72135, 0.48116, 0.30275, 0.21163, 0.17684, 0.15153, 0.12967, 0.11173, 0.097807, 0.087477, 0.080117, 0.072395, 0.049625, 0.03295, 0.022127, 0.014421, 0.0093171, 0.0063135, 0.0044326, 0.0032466, 0.0024813, 0.0019671, 0.0016139, 0.00097843, 0.00056414, 0.00029751, 0.00012048, 0}, + {1, 0.96122, 0.72135, 0.48117, 0.30276, 0.21164, 0.17686, 0.15155, 0.12969, 0.11175, 0.097823, 0.087493, 0.080133, 0.07241, 0.049639, 0.032963, 0.02214, 0.014433, 0.0093278, 0.0063234, 0.0044418, 0.0032552, 0.0024894, 0.0019747, 0.0016211, 0.00098334, 0.00056724, 0.00029923, 0.00012118, 0}, + {1, 0.96122, 0.72136, 0.48118, 0.30277, 0.21165, 0.17687, 0.15156, 0.1297, 0.11176, 0.09783, 0.0875, 0.08014, 0.072416, 0.049645, 0.032969, 0.022145, 0.014438, 0.0093325, 0.0063278, 0.004446, 0.0032591, 0.0024931, 0.0019782, 0.0016245, 0.00098568, 0.00056878, 0.00030012, 0.00012156, 0}, + {1, 0.96122, 0.72136, 0.48118, 0.30277, 0.21165, 0.17687, 0.15156, 0.1297, 0.11176, 0.097832, 0.087502, 0.080143, 0.072419, 0.049648, 0.032971, 0.022148, 0.01444, 0.0093345, 0.0063297, 0.0044477, 0.0032607, 0.0024946, 0.0019797, 0.0016259, 0.0009867, 0.00056946, 0.00030053, 0.00012174, 0}, +}; + +// Position of beta*gamma in the table rows: lower row and fraction towards the next one +inline void rowOf(double betaGamma, int& row, double& frac) +{ + const double x = (std::log10(std::max(betaGamma, 1e-6)) - kLogBetaGammaMin) / kLogBetaGammaStep; + if (x <= 0.) { + row = 0; + frac = 0.; + } else if (x >= kNBetaGamma - 1) { + row = kNBetaGamma - 2; + frac = 1.; + } else { + row = static_cast(x); + frac = x - row; + } +} + +// Log-log interpolation of a per-row table +inline double interpolate(const double* table, int row, double frac) +{ + return std::exp((1. - frac) * std::log(table[row]) + frac * std::log(table[row + 1])); +} +} // namespace + +namespace o2 +{ +namespace trd +{ + +IonisationModel::IonisationModel(SimParam::GasMixture gas) +{ + if (gas == SimParam::GasMixture::Xenon) { + mNEnergy = std::size(kXeEnergy); + mEnergy = kXeEnergy; + mCollisions = kXeCollisionsPerCm; + mDEdx = kXeRestrictedDEdx; + mFractionAbove = &kXeFractionAbove[0][0]; + } else if (gas == SimParam::GasMixture::Argon) { + mNEnergy = std::size(kArEnergy); + mEnergy = kArEnergy; + mCollisions = kArCollisionsPerCm; + mDEdx = kArRestrictedDEdx; + mFractionAbove = &kArFractionAbove[0][0]; + } else { + LOG(fatal) << "TRD ionisation model: no collision tables for this gas mixture"; + } +} + +double IonisationModel::minBetaGamma() +{ + return std::pow(10., kLogBetaGammaMin); +} + +double IonisationModel::deltaRayCut() +{ + return kDeltaRayCut; +} + +double IonisationModel::collisionsPerCm(double betaGamma) const +{ + int row; + double frac; + rowOf(betaGamma, row, frac); + return interpolate(mCollisions, row, frac); +} + +double IonisationModel::restrictedDEdx(double betaGamma) const +{ + int row; + double frac; + rowOf(betaGamma, row, frac); + return interpolate(mDEdx, row, frac); +} + +double IonisationModel::sampleCollision(int row, double u) const +{ + // u is the fraction of collisions above the sampled energy; the tabulated fractions fall with energy + const float* f = mFractionAbove + row * mNEnergy; + int j = 1; + while (j < mNEnergy - 1 && f[j] > u) { + ++j; + } + const double t = (f[j - 1] - u) / std::max(f[j - 1] - f[j], 1e-12f); + return mEnergy[j - 1] * std::exp(t * std::log(mEnergy[j] / mEnergy[j - 1])); +} + +double IonisationModel::stepEnergy(double length, double betaGamma, double z2, TRandom& rnd) const +{ + int row; + double frac; + rowOf(betaGamma, row, frac); + const int n = rnd.Poisson(length * z2 * interpolate(mCollisions, row, frac)); + if (n == 0) { + return 0.; + } + // one of the two neighbouring spectra, with the interpolation weight as probability + if (rnd.Rndm() < frac) { + ++row; + } + double energy = 0.; + for (int i = 0; i < n; ++i) { + energy += sampleCollision(row, rnd.Rndm()); + } + return energy; +} + +} // namespace trd +} // namespace o2 diff --git a/Detectors/TRD/simulation/test/testTRDIonisationModel.cxx b/Detectors/TRD/simulation/test/testTRDIonisationModel.cxx new file mode 100644 index 0000000000000..80bace19d9bfe --- /dev/null +++ b/Detectors/TRD/simulation/test/testTRDIonisationModel.cxx @@ -0,0 +1,105 @@ +// 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. + +/// \file testTRDIonisationModel.cxx +/// \brief Checks the TRD gas collision model against its own tables + +#define BOOST_TEST_MODULE Test TRD IonisationModel +#define BOOST_TEST_MAIN +#define BOOST_TEST_DYN_LINK +#include +#include "TRDSimulation/IonisationModel.h" + +#include +#include + +namespace o2 +{ +namespace trd +{ + +/// The sampled mean energy per cm reproduces the tabulated restricted dE/dx +BOOST_AUTO_TEST_CASE(TRDIonisationModel_meanEnergy) +{ + TRandom3 rnd(4711); + for (auto gas : {SimParam::GasMixture::Xenon, SimParam::GasMixture::Argon}) { + IonisationModel model(gas); + for (double bg : {0.08, 0.5, 4., 30., 3000.}) { + const int n = 20000; + double sum = 0.; + for (int i = 0; i < n; ++i) { + sum += model.stepEnergy(1., bg, 1., rnd); + } + BOOST_TEST(sum / n == model.restrictedDEdx(bg), boost::test_tools::tolerance(0.02)); + } + } +} + +/// The fraction of steps without a collision is Poisson, and scales with the charge squared +BOOST_AUTO_TEST_CASE(TRDIonisationModel_zeroCollisions) +{ + TRandom3 rnd(815); + IonisationModel model(SimParam::GasMixture::Xenon); + const double bg = 4., length = 0.02; + for (double z2 : {1., 4.}) { + const int n = 100000; + int zero = 0; + for (int i = 0; i < n; ++i) { + zero += model.stepEnergy(length, bg, z2, rnd) == 0.; + } + const double expected = std::exp(-length * z2 * model.collisionsPerCm(bg)); + BOOST_TEST(double(zero) / n == expected, boost::test_tools::tolerance(0.02)); + } +} + +/// The collision rate has its minimum near beta*gamma 3-4, at 40-47 collisions per cm in Xe/CO2, and a plateau +BOOST_AUTO_TEST_CASE(TRDIonisationModel_minimum) +{ + IonisationModel model(SimParam::GasMixture::Xenon); + double bgMin = 0., nMin = 1e9; + for (double lbg = 0.; lbg <= 1.; lbg += 0.01) { + const double n = model.collisionsPerCm(std::pow(10., lbg)); + if (n < nMin) { + nMin = n; + bgMin = std::pow(10., lbg); + } + } + BOOST_TEST(nMin > 40.); + BOOST_TEST(nMin < 47.); + BOOST_TEST(bgMin > 2.5); + BOOST_TEST(bgMin < 5.); + BOOST_TEST(model.collisionsPerCm(1e5) == model.collisionsPerCm(1e6)); +} + +/// A single collision transfers between the first grid energy (the ionisation potential) and the delta-ray cut +BOOST_AUTO_TEST_CASE(TRDIonisationModel_collisionEnergy) +{ + TRandom3 rnd(42); + IonisationModel model(SimParam::GasMixture::Xenon); + // steps short enough that almost every step with energy has exactly one collision + const double length = 1e-4, bg = 4.; + int n = 0; + double sum = 0.; + for (int i = 0; i < 2000000; ++i) { + const double e = model.stepEnergy(length, bg, 1., rnd); + if (e > 0.) { + BOOST_TEST(e >= 12.); + BOOST_TEST(e <= IonisationModel::deltaRayCut()); + sum += e; + ++n; + } + } + BOOST_TEST(n > 5000); + BOOST_TEST(sum / n == model.restrictedDEdx(bg) / model.collisionsPerCm(bg), boost::test_tools::tolerance(0.2)); +} + +} // namespace trd +} // namespace o2 From 99013d7b9b04d210325efa9048d07e15fdb978c4 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Sun, 15 Nov 2026 15:01:30 +0100 Subject: [PATCH 2/2] Apply PAI in the TRD gas only with the energy-deposit charge This moves the PAI model for the TRD gas from g4config.in to g4config_trd_pai.in, which g4Config.C applies only with TRDSimParams.ionisationModel=0. - ionisationModel=0 keeps the current production: energy deposit with PAI. - ionisationModel=1 runs without PAI, since the TRD collision model does not use the Geant4 energy loss in the gas. - A custom G4.configMacroFile is applied as before, without PAI lines added. - The README describes the two modes. Co-Authored-By: Claude Opus 5.5 --- Detectors/TRD/simulation/README.md | 11 +++++++++++ Detectors/gconfig/g4Config.C | 17 +++++++++++++++++ Detectors/gconfig/g4config.in | 6 ------ Detectors/gconfig/g4config_trd_pai.in | 6 ++++++ 4 files changed, 34 insertions(+), 6 deletions(-) create mode 100644 Detectors/gconfig/g4config_trd_pai.in diff --git a/Detectors/TRD/simulation/README.md b/Detectors/TRD/simulation/README.md index c8512f29b28fc..563497d0c71ab 100644 --- a/Detectors/TRD/simulation/README.md +++ b/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: ``` diff --git a/Detectors/gconfig/g4Config.C b/Detectors/gconfig/g4Config.C index 6173066309d98..bbbc2539bfad3 100644 --- a/Detectors/gconfig/g4Config.C +++ b/Detectors/gconfig/g4Config.C @@ -72,6 +72,8 @@ R__LOAD_LIBRARY(libgeant4vmc) #include "DetectorsBase/VMCSeederService.h" #include "CLHEP/Random/Random.h" #include "TRandom.h" +#include "CommonUtils/ConfigurableParam.h" +#include #include "DetectorsBase/MTWorkerHooks.h" #include "G4RunManager.hh" #include "G4Threading.hh" @@ -84,6 +86,16 @@ R__LOAD_LIBRARY(libgeant4vmc) // skips particles unknown to Geant4 (e.g. some ions) with a warning instead of a crash constexpr const char* kSkipUnknownParticles = "/mcPrimaryGenerator/skipUnknownParticles true"; +// whether the TRD charge comes from the Geant4 energy deposit (TRDSimParams.ionisationModel=0), true without TRD +bool trdUsesEnergyDeposit() +{ + try { + return o2::conf::ConfigurableParam::getValueAs("TRDSimParams.ionisationModel") == 0; + } catch (...) { + return true; + } +} + // Configuration macro for Geant4 VirtualMC void Config() { @@ -233,6 +245,11 @@ void Config() /// Customise Geant4 setting /// (verbose level, global range cut, ..) geant4->ProcessGeantMacro(configm1.Data()); + + // PAI in the TRD gas only when the TRD charge comes from the energy deposit + if (trdUsesEnergyDeposit()) { + geant4->ProcessGeantMacro((configm + "/Detectors/gconfig/g4config_trd_pai.in").Data()); + } } if (!mtMode) { geant4->ProcessGeantCommand(kSkipUnknownParticles); diff --git a/Detectors/gconfig/g4config.in b/Detectors/gconfig/g4config.in index cfc10ee20918b..7fa3cca1766c0 100644 --- a/Detectors/gconfig/g4config.in +++ b/Detectors/gconfig/g4config.in @@ -28,12 +28,6 @@ /process/optical/cerenkov/setTrackSecondariesFirst false /mcMagField/stepperType NystromRK4 -# PAI for TRD -# Geant4 VMC >= v3.2 -/mcPhysics/emModel/setEmModel PAI -/mcPhysics/emModel/setRegions TRD_Gas-mix -/mcPhysics/emModel/setParticles e- e+ mu- mu+ pi- pi+ kaon- kaon+ proton anti_proton deuteron triton He3 alpha - # # Precise Msc for EMCAL # diff --git a/Detectors/gconfig/g4config_trd_pai.in b/Detectors/gconfig/g4config_trd_pai.in new file mode 100644 index 0000000000000..137a3adab438a --- /dev/null +++ b/Detectors/gconfig/g4config_trd_pai.in @@ -0,0 +1,6 @@ +# PAI for the TRD gas, applied by g4Config.C only when the TRD charge comes from the Geant4 +# energy deposit (TRDSimParams.ionisationModel=0) +# Geant4 VMC >= v3.2 +/mcPhysics/emModel/setEmModel PAI +/mcPhysics/emModel/setRegions TRD_Gas-mix +/mcPhysics/emModel/setParticles e- e+ mu- mu+ pi- pi+ kaon- kaon+ proton anti_proton deuteron triton He3 alpha