Update electron/positron bremss

core/include/detray/materials/interaction.hpp

@@ -23,7 +23,7 @@ struct interaction {
     using relativistic_quantities =
         detail::relativistic_quantities<scalar_type>;
 
-    // @returns the total stopping power
+    // @returns the total stopping power (<-dE/dx>)
     DETRAY_HOST_DEVICE scalar_type
     compute_stopping_power(const detray::material<scalar_type>& mat,
                            const pdg_particle<scalar_type>& ptc,
@@ -68,9 +68,20 @@ struct interaction {
 
         scalar_type stopping_power{0.f};
 
-        // Only consider electrons at the moment
+        // Only consider electrons and positrons at the moment
         // For middle-heavy particles muons, the bremss is negligibe
-        if (ptc.pdg_num() == electron<scalar_type>().pdg_num()) {
+        //
+        // Also, over 10 MeV, positron and electron bremss might be similar.
+        // In ICRU 37, it is written that "In our tabulations, the radiative
+        // stopping power for positrons has been assumed to be the same as that
+        // for electrons, which is a good approximation at energies above, say,
+        // 10 MeV. However, it should be mentioned that exploratory calculations
+        // by Feng et at. (1981), employing the same method as that previously
+        // used by them for electrons, indicate significant differences between
+        // positrons and electrons in regard to the differential bremsstrahlung
+        // cross sections in oxygen and uranium at 500, 50, and 10 keV"
+        if (ptc.pdg_num() == electron<scalar_type>().pdg_num() ||
+            ptc.pdg_num() == positron<scalar_type>().pdg_num()) {
             // Stopping power ~ E/X (B. B. Rossi, High-energy particles,1952)
             // This approximation gets poor in low energy below 10 MeV
             stopping_power = rq.m_gamma * ptc.mass() / mat.X0();
@@ -148,7 +159,8 @@ struct interaction {
 
         scalar_type derivative{0.f};
 
-        if (ptc.pdg_num() == electron<scalar_type>().pdg_num()) {
+        if (ptc.pdg_num() == electron<scalar_type>().pdg_num() ||
+            ptc.pdg_num() == positron<scalar_type>().pdg_num()) {
             // Stopping power ~ E/X = gamma * m/X
             // Derivative = dgamma/dqop * m/X = -beta^2 gamma/qop * m/X
             // (Eq (D.5) of arXiv:2403.16720)

tests/unit_tests/cpu/CMakeLists.txt

@@ -63,6 +63,7 @@ macro(detray_add_cpu_test algebra)
        "material/material_maps.cpp"
        "material/materials.cpp"
        "material/stopping_power_derivative.cpp"
+       "material/stopping_power.cpp"
        "navigation/intersection/cuboid_intersector.cpp"
        "navigation/intersection/cylinder_intersector.cpp"
        "navigation/intersection/helix_intersector.cpp"

tests/unit_tests/cpu/material/bethe_equation.cpp

@@ -95,7 +95,7 @@ INSTANTIATE_TEST_SUITE_P(
                                       0.1003f * unit<scalar>::GeV, 3.082f)));
 
 /*
-//@fixme: Test fails with He Gas and 10 GeV muons (18 % difference)
+//@fixme: Test fails with He Gas and 1 GeV muons (18 % difference)
 INSTANTIATE_TEST_SUITE_P(
     detray_material_Bethe_1GeV_HeGas, EnergyLossBetheValidation,
     ::testing::Values(std::make_tuple(helium_gas<scalar>(), muon<scalar>(),

tests/unit_tests/cpu/material/bremsstrahlung.cpp

@@ -60,37 +60,67 @@ TEST_P(EnergyLossBremsValidation, bremsstrahlung) {
 
     // We have not implemented the bremsstrahlung for the heavier charged
     // particles, which is negligible
-    if (ptc.pdg_num() == electron<scalar>().pdg_num()) {
-        // Check if difference is within 11% error
-        EXPECT_NEAR((expected_dEdx - dEdx) / dEdx, 0.f, 0.11f);
+    if ((ptc.pdg_num() == electron<scalar>().pdg_num()) ||
+        (ptc.pdg_num() == positron<scalar>().pdg_num())) {
+        // Check if difference is within 14% error
+        EXPECT_NEAR((expected_dEdx - dEdx) / dEdx, 0.f, 0.14f);
     } else {
         EXPECT_FLOAT_EQ(static_cast<float>(dEdx), 0.f);
     }
 }
 
-// REFERENCE
-//
-// Atomic Data and Nuclear Data Tables Volume 4, March 1972, Pages 1-27,
-//"Energy loss, range, and bremsstrahlung yield for 10-keV to 100-MeV electrons
-// in various elements and chemical compounds"
+// From https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html
+// Assumes that the stopping powers of electron and positron are the same
 
+// Electrons
 INSTANTIATE_TEST_SUITE_P(
-    detray_material_Bremsstrahlung_100MeV_He, EnergyLossBremsValidation,
+    electron_Bremsstrahlung_100MeV_He, EnergyLossBremsValidation,
     ::testing::Values(std::make_tuple(helium_gas<scalar>(), electron<scalar>(),
-                                      100.0f * unit<scalar>::MeV, 0.95886f)));
+                                      100.0f * unit<scalar>::MeV, 0.9229f)));
 
 INSTANTIATE_TEST_SUITE_P(
-    detray_material_Bremsstrahlung_100MeV_Al, EnergyLossBremsValidation,
+    electron_Bremsstrahlung_100MeV_Al, EnergyLossBremsValidation,
     ::testing::Values(std::make_tuple(aluminium<scalar>(), electron<scalar>(),
-                                      100.0f * unit<scalar>::MeV, 3.8172f)));
+                                      100.0f * unit<scalar>::MeV, 3.714f)));
 
 INSTANTIATE_TEST_SUITE_P(
-    detray_material_Bremsstrahlung_100MeV_Cu, EnergyLossBremsValidation,
+    electron_Bremsstrahlung_100MeV_Si, EnergyLossBremsValidation,
+    ::testing::Values(std::make_tuple(silicon<scalar>(), electron<scalar>(),
+                                      100.0f * unit<scalar>::MeV, 4.099f)));
+
+INSTANTIATE_TEST_SUITE_P(
+    electron_Bremsstrahlung_100MeV_Cu, EnergyLossBremsValidation,
     ::testing::Values(std::make_tuple(copper<scalar>(), electron<scalar>(),
-                                      100.0f * unit<scalar>::MeV, 7.2365f)));
+                                      100.0f * unit<scalar>::MeV, 7.079f)));
+
+// Positrons
+INSTANTIATE_TEST_SUITE_P(
+    positron_Bremsstrahlung_100MeV_He, EnergyLossBremsValidation,
+    ::testing::Values(std::make_tuple(helium_gas<scalar>(), positron<scalar>(),
+                                      100.0f * unit<scalar>::MeV, 0.9229f)));
+
+INSTANTIATE_TEST_SUITE_P(
+    positron_Bremsstrahlung_100MeV_Al, EnergyLossBremsValidation,
+    ::testing::Values(std::make_tuple(aluminium<scalar>(), positron<scalar>(),
+                                      100.0f * unit<scalar>::MeV, 3.714f)));
+
+INSTANTIATE_TEST_SUITE_P(
+    positron_Bremsstrahlung_100MeV_Si, EnergyLossBremsValidation,
+    ::testing::Values(std::make_tuple(silicon<scalar>(), positron<scalar>(),
+                                      100.0f * unit<scalar>::MeV, 4.099f)));
+
+INSTANTIATE_TEST_SUITE_P(
+    positron_Bremsstrahlung_100MeV_Cu, EnergyLossBremsValidation,
+    ::testing::Values(std::make_tuple(copper<scalar>(), positron<scalar>(),
+                                      100.0f * unit<scalar>::MeV, 7.079f)));
 
 // We have not implemented the bremsstrahlung for muons
 INSTANTIATE_TEST_SUITE_P(
-    detray_material_Bremsstrahlung_100MeV_Cu_muon, EnergyLossBremsValidation,
+    muon_Bremsstrahlung_100MeV_Cu_muon, EnergyLossBremsValidation,
     ::testing::Values(std::make_tuple(copper<scalar>(), muon<scalar>(),
                                       100.0f * unit<scalar>::MeV, 0.f)));
+
+INSTANTIATE_TEST_SUITE_P(
+    antimuon_Bremsstrahlung_100MeV_Cu_muon, EnergyLossBremsValidation,
+    ::testing::Values(std::make_tuple(copper<scalar>(), antimuon<scalar>(),
+                                      100.0f * unit<scalar>::MeV, 0.f)));

tests/unit_tests/cpu/material/stopping_power.cpp

@@ -0,0 +1,159 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2024 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s).
+#include "detray/materials/interaction.hpp"
+#include "detray/materials/material.hpp"
+#include "detray/materials/predefined_materials.hpp"
+
+// Detray test include(s)
+#include "detray/test/utils/types.hpp"
+
+// GTest include(s).
+#include <gtest/gtest.h>
+
+using namespace detray;
+
+// Test class for the stopping power
+// Input tuple: < material, particle type, kinetic energy, expected output >
+class StoppingPowerValidation
+    : public ::testing::TestWithParam<
+          std::tuple<material<test::scalar>, pdg_particle<test::scalar>,
+                     test::scalar, test::scalar>> {};
+
+TEST_P(StoppingPowerValidation, stopping_power) {
+
+    // Interaction object
+    interaction<test::scalar> I;
+
+    // Material
+    material<test::scalar> mat = std::get<0>(GetParam());
+
+    // Particle
+    pdg_particle<test::scalar> ptc = std::get<1>(GetParam());
+
+    // Kinetic energy
+    const test::scalar T = std::get<2>(GetParam());
+
+    // Total energy
+    const test::scalar E = T + ptc.mass();
+
+    // Momentum
+    const test::scalar p = math::sqrt(E * E - ptc.mass() * ptc.mass());
+
+    // qoverp
+    const test::scalar qop{ptc.charge() / p};
+
+    // Stopping power in MeV * cm^2 / g
+    const test::scalar dEdx{I.compute_stopping_power(mat, ptc, {ptc, qop}) /
+                            mat.mass_density() /
+                            (unit<test::scalar>::MeV * unit<test::scalar>::cm2 /
+                             unit<test::scalar>::g)};
+
+    const test::scalar expected_dEdx = std::get<3>(GetParam());
+
+    // Check if difference is within 8% error
+    EXPECT_NEAR((expected_dEdx - dEdx) / dEdx, 0.f, 0.08f);
+}
+
+/******************
+ *   Muon tests
+ ******************/
+
+// From https://pdg.lbl.gov/2024/AtomicNuclearProperties/index.html
+// Note 1: that we took the PDG value only from Ionization loss (Radiative loss
+// is ignored) Note 2: assumes that the stopping powers of muon and antimuon are
+// the same Note 3: Test fails with He Gas and 1 GeV muons (18 % difference)
+INSTANTIATE_TEST_SUITE_P(
+    muon_stopping_power_He, StoppingPowerValidation,
+    ::testing::Values(
+        std::make_tuple(helium_gas<test::scalar>(), muon<test::scalar>(),
+                        100.0f * unit<test::scalar>::MeV, 2.165f),
+        // std::make_tuple(helium_gas<test::scalar>(), muon<test::scalar>(),
+        //                 1.f * unit<test::scalar>::GeV, 2.133f),
+        std::make_tuple(helium_gas<test::scalar>(), muon<test::scalar>(),
+                        10.0f * unit<test::scalar>::GeV, 2.768f),
+        std::make_tuple(helium_gas<test::scalar>(), muon<test::scalar>(),
+                        100.0f * unit<test::scalar>::GeV, 3.188f)));
+
+INSTANTIATE_TEST_SUITE_P(
+    muon_stopping_power_Si, StoppingPowerValidation,
+    ::testing::Values(
+        std::make_tuple(silicon<test::scalar>(), muon<test::scalar>(),
+                        100.0f * unit<test::scalar>::MeV, 1.849f),
+        std::make_tuple(silicon<test::scalar>(), muon<test::scalar>(),
+                        1.f * unit<test::scalar>::GeV, 1.803f),
+        std::make_tuple(silicon<test::scalar>(), muon<test::scalar>(),
+                        10.0f * unit<test::scalar>::GeV, 2.177f),
+        std::make_tuple(silicon<test::scalar>(), muon<test::scalar>(),
+                        100.0f * unit<test::scalar>::GeV, 2.451f)));
+
+INSTANTIATE_TEST_SUITE_P(
+    anti_muon_stopping_power_He, StoppingPowerValidation,
+    ::testing::Values(
+        std::make_tuple(helium_gas<test::scalar>(), antimuon<test::scalar>(),
+                        100.0f * unit<test::scalar>::MeV, 2.165f),
+        // std::make_tuple(helium_gas<test::scalar>(), antimuon<test::scalar>(),
+        //                 1.f * unit<test::scalar>::GeV, 2.133f),
+        std::make_tuple(helium_gas<test::scalar>(), antimuon<test::scalar>(),
+                        10.0f * unit<test::scalar>::GeV, 2.768f),
+        std::make_tuple(helium_gas<test::scalar>(), antimuon<test::scalar>(),
+                        100.0f * unit<test::scalar>::GeV, 3.188f)));
+
+INSTANTIATE_TEST_SUITE_P(
+    anti_muon_stopping_power_Si, StoppingPowerValidation,
+    ::testing::Values(
+        std::make_tuple(silicon<test::scalar>(), antimuon<test::scalar>(),
+                        100.0f * unit<test::scalar>::MeV, 1.849f),
+        std::make_tuple(silicon<test::scalar>(), antimuon<test::scalar>(),
+                        1.f * unit<test::scalar>::GeV, 1.803f),
+        std::make_tuple(silicon<test::scalar>(), antimuon<test::scalar>(),
+                        10.0f * unit<test::scalar>::GeV, 2.177f),
+        std::make_tuple(silicon<test::scalar>(), antimuon<test::scalar>(),
+                        100.0f * unit<test::scalar>::GeV, 2.451f)));
+
+/*********************
+ *   Electron tests
+ *********************/
+
+// From https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html
+// Assumes that the stopping powers of electron and positron are the same
+INSTANTIATE_TEST_SUITE_P(
+    electron_stopping_power_He, StoppingPowerValidation,
+    ::testing::Values(std::make_tuple(helium_gas<test::scalar>(),
+                                      electron<test::scalar>(),
+                                      100.0f * unit<test::scalar>::MeV, 3.532f),
+                      std::make_tuple(helium_gas<test::scalar>(),
+                                      electron<test::scalar>(),
+                                      1.f * unit<test::scalar>::GeV, 13.14f)));
+
+INSTANTIATE_TEST_SUITE_P(
+    electron_stopping_power_Si, StoppingPowerValidation,
+    ::testing::Values(std::make_tuple(silicon<test::scalar>(),
+                                      electron<test::scalar>(),
+                                      100.0f * unit<test::scalar>::MeV, 6.017f),
+                      std::make_tuple(silicon<test::scalar>(),
+                                      electron<test::scalar>(),
+                                      1.f * unit<test::scalar>::GeV, 46.69f)));
+
+INSTANTIATE_TEST_SUITE_P(
+    positron_stopping_power_He, StoppingPowerValidation,
+    ::testing::Values(std::make_tuple(helium_gas<test::scalar>(),
+                                      positron<test::scalar>(),
+                                      100.0f * unit<test::scalar>::MeV, 3.532f),
+                      std::make_tuple(helium_gas<test::scalar>(),
+                                      positron<test::scalar>(),
+                                      1.f * unit<test::scalar>::GeV, 13.14f)));
+
+INSTANTIATE_TEST_SUITE_P(
+    positron_stopping_power_Si, StoppingPowerValidation,
+    ::testing::Values(std::make_tuple(silicon<test::scalar>(),
+                                      positron<test::scalar>(),
+                                      100.0f * unit<test::scalar>::MeV, 6.017f),
+                      std::make_tuple(silicon<test::scalar>(),
+                                      positron<test::scalar>(),
+                                      1.f * unit<test::scalar>::GeV, 46.69f)));