Add energy and range straggeling to energy loss models (ATTPC#256)

* Add interface for straggling (range and energy)

* Add and test range straggling to Eloss models (CATIMA and SRIM)

* Fix density scaling to SRIM (dE/dx was scaled, energy loss was not)

* Remove scaling parameter in AtELossModel. Enforce scaling in SetDensity function.

* Add range test for CATIMA at different density

* Set environment variables before running tests

Author: anthoak13

.github/workflows/CI-tests.yml

@@ -8,6 +8,7 @@ on: [push, pull_request]
 env:
   # Customize the CMake build type here (Release, Debug, RelWithDebInfo, etc.)
   BUILD_TYPE: RELWITHDEBINFO
+  
 
 jobs:
   build:
@@ -16,7 +17,6 @@ jobs:
     # See: https://docs.github.com/en/free-pro-team@latest/actions/learn-github-actions/managing-complex-workflows#using-a-build-matrix
     runs-on: ubuntu-latest
     container: anthoak13/attpcroot-deps
-
     steps:
     - uses: actions/checkout@v4
 
@@ -33,5 +33,6 @@ jobs:
       # working-directory: ${{github.workspace}}/build
       # Execute tests defined by the CMake configuration.
       # See https://cmake.org/cmake/help/latest/manual/ctest.1.html for more detail
-      run: cd ${{github.workspace}}/build && ctest
+      run: cd ${{github.workspace}}/build && source config.sh && ctest
+      shell: bash
 

AtTools/AtELossCATIMA.cxx

@@ -1,37 +1,34 @@
 #include "AtELossCATIMA.h"
 
 #include <FairLogger.h>
+namespace AtTools {
 
-AtTools::AtELossCATIMA::AtELossCATIMA(double density, std::vector<std::tuple<int, int, int>> materialComponents)
+AtELossCATIMA::AtELossCATIMA(double density, std::vector<std::tuple<int, int, int>> materialComponents)
    : AtELossModel(density)
 {
-   fMaterial = std::make_unique<catima::Material>();
-
-   for (auto materialComponent : materialComponents)
-      fMaterial->add_element(std::get<0>(materialComponent), std::get<1>(materialComponent),
-                             std::get<2>(materialComponent));
+   SetMaterial(materialComponents);
 }
 
-double AtTools::AtELossCATIMA::GetdEdx(double energy) const
+double AtELossCATIMA::GetdEdx(double energy) const
 {
    if (fProjectile == nullptr) {
       LOG(warning)
          << " Warning in AtTools::AtELossCATIMA::GetdEdx : The projectile was not set! GetdEdx will return 0!";
       return 0;
    }
 
-   if (fProjectileMassUma <= 0) {
+   if (fProjectileMassAmu <= 0) {
       LOG(error) << " Error in AtTools::AtELossCATIMA::GetdEdx : The projectile's mass in umas can not be <= 0! "
                     "GetdEdx will return 0!";
       return 0;
    }
 
-   catima::Result result = catima::calculate(*fProjectile, *fMaterial, energy / fProjectileMassUma);
-   double dEdx = result.dEdxi * fDensity;
-   return dEdx;
+   catima::Result result = catima::calculate(*fProjectile, *fMaterial, energy / fProjectileMassAmu);
+   double dEdx = result.dEdxi * fDensity; // MeV/cm
+   return dEdx / 10.0;                    // convert to MeV/mm
 }
 
-double AtTools::AtELossCATIMA::GetRange(double energyIni, double energyFin) const
+double AtELossCATIMA::GetRange(double energyIni, double energyFin) const
 {
    if (energyFin < 0) {
       LOG(warning) << " Warning in AtTools::AtELossCATIMA::GetRange : The final energy was set to a negative value! "
@@ -40,14 +37,14 @@ double AtTools::AtELossCATIMA::GetRange(double energyIni, double energyFin) cons
    }
 
    if (energyFin == 0) {
-      fProjectile->T = energyIni / fProjectileMassUma;
+      fProjectile->T = energyIni / fProjectileMassAmu;
       return catima::range(*fProjectile, *fMaterial) / fDensity * 10.;
    }
 
    double remainingEnergy{energyIni};
    double range{0};
    while (remainingEnergy > energyFin) {
-      catima::Result result = catima::calculate(*fProjectile, *fMaterial, remainingEnergy / fProjectileMassUma);
+      catima::Result result = catima::calculate(*fProjectile, *fMaterial, remainingEnergy / fProjectileMassAmu);
       double dEdx = result.dEdxi * fDensity;
       double DE = dEdx * fRangeStepSize / 10.;
 
@@ -63,12 +60,12 @@ double AtTools::AtELossCATIMA::GetRange(double energyIni, double energyFin) cons
    return range;
 }
 
-double AtTools::AtELossCATIMA::GetEnergy(double energyIni, double distance) const
+double AtELossCATIMA::GetEnergy(double energyIni, double distance) const
 {
    double remainingEnergy{energyIni};
    double range{0};
    while (range < distance) {
-      catima::Result result = catima::calculate(*fProjectile, *fMaterial, remainingEnergy / fProjectileMassUma);
+      catima::Result result = catima::calculate(*fProjectile, *fMaterial, remainingEnergy / fProjectileMassAmu);
       double dEdx = result.dEdxi * fDensity;
       double DE{};
 
@@ -86,9 +83,55 @@ double AtTools::AtELossCATIMA::GetEnergy(double energyIni, double distance) cons
 
    return remainingEnergy;
 }
+double AtELossCATIMA::GetRangeVariance(double energy) const
+{
+   if (fProjectile == nullptr || fMaterial == nullptr) {
+      LOG(error) << "Projectile or material not set. Range variance is 0.";
+      return 0;
+   }
+   auto range_var =
+      catima::range_variance(*fProjectile, energy / fProjectileMassAmu, *fMaterial); // range var in (g/cm^2)^2
+   LOG(debug) << "Range variance in (g/cm^2)^2: " << range_var << " for energy: " << energy / fProjectileMassAmu
+              << " MeV/u";
+   range_var /= fDensity * fDensity; // convert to (cm)^2
+   LOG(debug) << "Range variance in (cm)^2: " << range_var << " for energy: " << energy / fProjectileMassAmu
+              << " MeV/u";
+   return range_var * 100; // convert to mm^2
+}
+
+double AtELossCATIMA::GetElossStraggling(double energyIni, double energyFin) const
+{
+   if (fProjectile == nullptr || fMaterial == nullptr) {
+      LOG(error) << "Projectile or material not set.";
+      return 0;
+   }
+   if (energyFin > energyIni) {
+      LOG(error) << "Final energy must be less than initial energy!";
+      return 0;
+   }
+   auto energy_strag = catima::energy_straggling_from_E(*fProjectile, energyIni / fProjectileMassAmu,
+                                                        energyFin / fProjectileMassAmu, *fMaterial);
+   return energy_strag;
+}
+double AtELossCATIMA::GetdEdxStraggling(double energyIni, double energyFin) const
+{
+   if (fProjectile == nullptr || fMaterial == nullptr) {
+      LOG(error) << "Projectile or material not set. dEdx straggling is 0.";
+      return 0;
+   }
+   auto dedx_min = GetdEdx(energyIni);
+   auto dedx_max = GetdEdx(energyFin);
+   if (std::abs(dedx_min - dedx_max) / dedx_min > 0.01) {
+      LOG(warning) << "From " << energyIni << " to " << energyFin
+                   << " MeV the dEdx is not constant. dEdx straggling calculation is unreliable.";
+   }
+   auto dE_st = GetElossStraggling(energyIni, energyFin);
+   auto factor = dE_st / (energyIni - energyFin);
+   return factor * dedx_min;
+}
 
 std::vector<std::pair<double, double>>
-AtTools::AtELossCATIMA::GetBraggCurve(double energy, double rangeStepSize, double totalFractionELoss) const
+AtELossCATIMA::GetBraggCurve(double energy, double rangeStepSize, double totalFractionELoss) const
 {
    if (rangeStepSize == 0)
       return GetBraggCurve(energy, fRangeStepSize, totalFractionELoss);
@@ -99,7 +142,7 @@ AtTools::AtELossCATIMA::GetBraggCurve(double energy, double rangeStepSize, doubl
    double range{};
    while (remainingEnergy / energy > totalFractionELoss) {
 
-      catima::Result result = catima::calculate(*fProjectile, *fMaterial, remainingEnergy / fProjectileMassUma);
+      catima::Result result = catima::calculate(*fProjectile, *fMaterial, remainingEnergy / fProjectileMassAmu);
       double dEdx = result.dEdxi * fDensity;
       braggCurve.push_back(std::make_pair(dEdx, range));
 
@@ -113,3 +156,14 @@ AtTools::AtELossCATIMA::GetBraggCurve(double energy, double rangeStepSize, doubl
 
    return braggCurve;
 }
+
+void AtELossCATIMA::SetMaterial(std::vector<std::tuple<int, int, int>> materialComponents)
+{
+   fMaterial = std::make_unique<catima::Material>();
+   for (const auto &materialComponent : materialComponents) {
+      fMaterial->add_element(std::get<0>(materialComponent), std::get<1>(materialComponent),
+                             std::get<2>(materialComponent));
+   }
+}
+
+} // namespace AtTools

AtTools/AtELossCATIMA.h

@@ -16,18 +16,23 @@ class AtELossCATIMA : public AtELossModel {
    std::unique_ptr<catima::Material> fMaterial{nullptr};
    std::unique_ptr<catima::Projectile> fProjectile{nullptr};
 
-   double fProjectileMassUma{-1};
+   double fProjectileMassAmu{-1}; /// Mass of the projectile in amu (atomic mass units).
 
    double fRangeStepSize{0.1}; // mm
 
 public:
    /**
     * Initializer of the CATIMA AtELossModel wrapper.
-    * @param[in] density Density of the material.
+    * @param[in] density Density of the material (g/cm^2).
     * @param[in] materialComponents Components of the material. They are passed as a vector of tuples (A, Z,
     * stoichiometry).
     */
+   AtELossCATIMA(double density) : AtELossModel(density) {}
    AtELossCATIMA(double density, std::vector<std::tuple<int, int, int>> materialComponents);
+   AtELossCATIMA(double density, const catima::Material &material)
+      : AtELossModel(density), fMaterial(std::make_unique<catima::Material>(material))
+   {
+   }
 
    virtual double GetdEdx(double energy) const override;
    virtual double GetRange(double energyIni, double energyFin = 0) const override;
@@ -37,20 +42,27 @@ class AtELossCATIMA : public AtELossModel {
    }
    virtual double GetEnergy(double energyIni, double distance) const override;
 
+   virtual double GetRangeVariance(double energy) const override;
+   virtual double GetElossStraggling(double energyIni, double energyFin) const override;
+   virtual double GetdEdxStraggling(double energyIni, double energyFin) const override;
+
    virtual std::vector<std::pair<double, double>>
    GetBraggCurve(double energy, double rangeStepSize = 0, double totalFractionELoss = 0.001) const override;
 
    /**
     * Setter of the catima projectile used for calculations.
     * @param[in] A Mass number of the projectile.
     * @param[in] Z Charge number of the projectile.
-    * @param[in] massUma Mass of the projectile in umas.
+    * @param[in] massAmu Mass of the projectile in amu.
     */
-   void SetProjectile(double A, double Z, double massUma)
+   void SetProjectile(double A, double Z, double massAmu)
    {
       fProjectile = std::make_unique<catima::Projectile>(A, Z);
-      fProjectileMassUma = massUma;
+      fProjectileMassAmu = massAmu;
    }
+   void SetMaterial(const catima::Material &material) { fMaterial = std::make_unique<catima::Material>(material); }
+   void SetMaterial(std::vector<std::tuple<int, int, int>> materialComponents);
+
    /**
     * Setter of the range step size used for calculations. By default it is set to 0.1mm.
     * @param[in] stepSize The step size used for ranges. It must be input in mm.

AtTools/AtELossCATIMATest.cxx

@@ -1,11 +1,99 @@
 #include "AtELossCATIMA.h"
 
+#include <FairLogger.h>
+
+#include <catima/catima.h>
+#include <catima/structures.h>
 #include <cmath>
 #include <gtest/gtest.h>
 #include <tuple>
 #include <utility>
 #include <vector>
 
+/**
+ * Test fixture for the AtELossCATIMA class that initializes a H2 gas w/ proton model.
+ */
+class AtELossCATIMATestFixture : public ::testing::Test {
+protected:
+   AtTools::AtELossCATIMA model;
+   double mass{1.007825031898}; // Mass of proton in amu
+
+   AtELossCATIMATestFixture() : model(6.5643e-5)
+   {
+      // Initialize the CATIMA model with H2 gas density and components and set projectile to proton.
+      model.SetMaterial(catima::Material(1, 1)); // Set material to H2
+      model.SetProjectile(1, 1, mass);           // Set projectile to proton
+   }
+};
+
+TEST_F(AtELossCATIMATestFixture, ConstructCATIMAModel)
+{
+   // Basic check: ensure model is constructed and can compute range
+   double range = model.GetRange(1.0); // 1 MeV
+   EXPECT_GT(range, 0.0);
+}
+
+TEST_F(AtELossCATIMATestFixture, TestRangeStraggling)
+{
+   // Check dEdx for a known energy
+   double range_var = model.GetRangeVariance(1.0); // 1 MeV
+   double expected = 1.99;                         // Expected value from LISE for H2 at 1 MeV
+   ASSERT_NEAR(range_var, expected * expected, 0.2 * expected * expected);
+
+   double range_straggling = model.GetRangeStraggling(1.0); // 1 MeV
+   ASSERT_NEAR(range_straggling, expected, 0.1 * expected);
+
+   // Check straggling at 10 MeV
+   range_straggling = model.GetRangeStraggling(10.0); // 10 MeV
+   expected = 95.933;                                 // mm
+   ASSERT_NEAR(range_straggling, expected, 0.1 * expected);
+
+   model.SetDensity(4e-5);
+   range_straggling = model.GetRangeStraggling(1.0); // 1 MeV
+   expected = 3.27;                                  // mm
+   ASSERT_NEAR(range_straggling, expected, 0.1 * expected);
+}
+
+TEST_F(AtELossCATIMATestFixture, TestEnergyLossStraggling)
+{
+
+   double expectedSigma = 0.0084 * mass;                                 // Expected sigma from LISE
+   double eloss_straggling = model.GetElossStraggling(1.0, 0.75 * mass); // 1 MeV to 10 MeV
+   ASSERT_NEAR(eloss_straggling, expectedSigma, 0.1 * expectedSigma);
+
+   expectedSigma = 0.0376 * mass;
+   eloss_straggling = model.GetElossStraggling(5.0, 3.58164 * mass);
+   ASSERT_NEAR(eloss_straggling, expectedSigma, 0.1 * expectedSigma);
+}
+
+TEST_F(AtELossCATIMATestFixture, TestEnergyLossStragglingDistance)
+{
+   double expectedSigma = 0.0084 * mass;                                  // Expected sigma from LISE
+   double eloss_straggling = model.GetElossStragglingDistance(1.0, 50.0); // 1 MeV over 10 mm
+   ASSERT_NEAR(eloss_straggling, expectedSigma, 0.1 * expectedSigma);
+
+   expectedSigma = 0.0376 * mass;
+   eloss_straggling = model.GetElossStragglingDistance(5.0, 1000.0); // 5 MeV over 100 mm
+   ASSERT_NEAR(eloss_straggling, expectedSigma, 0.1 * expectedSigma);
+}
+
+TEST_F(AtELossCATIMATestFixture, TestdEdxStraggling)
+{
+   double dedx = model.GetdEdx(5.0);
+   double dE = dedx * 50;
+   double expected_dE = model.GetEnergyLoss(5.0, 50.0);
+
+   ASSERT_NEAR(dE, expected_dE, 0.01 * expected_dE); // Verify linear assumption is true
+
+   double E_st = model.GetElossStragglingDistance(5.0, 50.0);
+   double dedx_straggling = model.GetdEdxStraggling(5.0, model.GetEnergy(5.0, 50.0));
+   ASSERT_NEAR(dedx_straggling * 50, E_st, 0.01 * E_st); // Check dEdx straggling
+
+   double e_min = expected_dE - E_st;
+   double e_min_dedx = dE - dedx_straggling * 50;
+   ASSERT_NEAR(e_min, e_min_dedx, 0.01 * e_min); // Check minimum energy loss
+}
+
 TEST(AtELossCATIMATest, LISE_Match)
 {
    // Create the vector with the gas components for H2.

AtTools/AtELossModel.cxx

@@ -9,11 +9,7 @@ namespace AtTools {
  */
 void AtELossModel::SetDensity(double density)
 {
-   if (fDensityIni == 0)
-      throw std::invalid_argument("Cannot set the density if the density of in model is not known");
-
    fDensity = density;
-   fdEdxScale = fDensity / fDensityIni;
 }
 
 std::vector<std::pair<double, double>>