diff --git a/core/include/detray/materials/interaction.hpp b/core/include/detray/materials/interaction.hpp index 9c8961305..458bcb60e 100644 --- a/core/include/detray/materials/interaction.hpp +++ b/core/include/detray/materials/interaction.hpp @@ -23,7 +23,7 @@ struct interaction { using relativistic_quantities = detail::relativistic_quantities; - // @returns the total stopping power + // @returns the total stopping power (<-dE/dx>) DETRAY_HOST_DEVICE scalar_type compute_stopping_power(const detray::material& mat, const pdg_particle& 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().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().pdg_num() || + ptc.pdg_num() == positron().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().pdg_num()) { + if (ptc.pdg_num() == electron().pdg_num() || + ptc.pdg_num() == positron().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) diff --git a/tests/unit_tests/cpu/CMakeLists.txt b/tests/unit_tests/cpu/CMakeLists.txt index 191828168..65efe611e 100644 --- a/tests/unit_tests/cpu/CMakeLists.txt +++ b/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" diff --git a/tests/unit_tests/cpu/material/bethe_equation.cpp b/tests/unit_tests/cpu/material/bethe_equation.cpp index b1ea5fb0c..9916e09de 100644 --- a/tests/unit_tests/cpu/material/bethe_equation.cpp +++ b/tests/unit_tests/cpu/material/bethe_equation.cpp @@ -95,7 +95,7 @@ INSTANTIATE_TEST_SUITE_P( 0.1003f * unit::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(), muon(), diff --git a/tests/unit_tests/cpu/material/bremsstrahlung.cpp b/tests/unit_tests/cpu/material/bremsstrahlung.cpp index af0942351..64098a3e8 100644 --- a/tests/unit_tests/cpu/material/bremsstrahlung.cpp +++ b/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().pdg_num()) { - // Check if difference is within 11% error - EXPECT_NEAR((expected_dEdx - dEdx) / dEdx, 0.f, 0.11f); + if ((ptc.pdg_num() == electron().pdg_num()) || + (ptc.pdg_num() == positron().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(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(), electron(), - 100.0f * unit::MeV, 0.95886f))); + 100.0f * unit::MeV, 0.9229f))); INSTANTIATE_TEST_SUITE_P( - detray_material_Bremsstrahlung_100MeV_Al, EnergyLossBremsValidation, + electron_Bremsstrahlung_100MeV_Al, EnergyLossBremsValidation, ::testing::Values(std::make_tuple(aluminium(), electron(), - 100.0f * unit::MeV, 3.8172f))); + 100.0f * unit::MeV, 3.714f))); INSTANTIATE_TEST_SUITE_P( - detray_material_Bremsstrahlung_100MeV_Cu, EnergyLossBremsValidation, + electron_Bremsstrahlung_100MeV_Si, EnergyLossBremsValidation, + ::testing::Values(std::make_tuple(silicon(), electron(), + 100.0f * unit::MeV, 4.099f))); + +INSTANTIATE_TEST_SUITE_P( + electron_Bremsstrahlung_100MeV_Cu, EnergyLossBremsValidation, ::testing::Values(std::make_tuple(copper(), electron(), - 100.0f * unit::MeV, 7.2365f))); + 100.0f * unit::MeV, 7.079f))); + +// Positrons +INSTANTIATE_TEST_SUITE_P( + positron_Bremsstrahlung_100MeV_He, EnergyLossBremsValidation, + ::testing::Values(std::make_tuple(helium_gas(), positron(), + 100.0f * unit::MeV, 0.9229f))); + +INSTANTIATE_TEST_SUITE_P( + positron_Bremsstrahlung_100MeV_Al, EnergyLossBremsValidation, + ::testing::Values(std::make_tuple(aluminium(), positron(), + 100.0f * unit::MeV, 3.714f))); + +INSTANTIATE_TEST_SUITE_P( + positron_Bremsstrahlung_100MeV_Si, EnergyLossBremsValidation, + ::testing::Values(std::make_tuple(silicon(), positron(), + 100.0f * unit::MeV, 4.099f))); + +INSTANTIATE_TEST_SUITE_P( + positron_Bremsstrahlung_100MeV_Cu, EnergyLossBremsValidation, + ::testing::Values(std::make_tuple(copper(), positron(), + 100.0f * unit::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(), muon(), 100.0f * unit::MeV, 0.f))); + +INSTANTIATE_TEST_SUITE_P( + antimuon_Bremsstrahlung_100MeV_Cu_muon, EnergyLossBremsValidation, + ::testing::Values(std::make_tuple(copper(), antimuon(), + 100.0f * unit::MeV, 0.f))); diff --git a/tests/unit_tests/cpu/material/stopping_power.cpp b/tests/unit_tests/cpu/material/stopping_power.cpp new file mode 100644 index 000000000..28ce0b545 --- /dev/null +++ b/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 + +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, pdg_particle, + test::scalar, test::scalar>> {}; + +TEST_P(StoppingPowerValidation, stopping_power) { + + // Interaction object + interaction I; + + // Material + material mat = std::get<0>(GetParam()); + + // Particle + pdg_particle 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::MeV * unit::cm2 / + unit::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(), muon(), + 100.0f * unit::MeV, 2.165f), + // std::make_tuple(helium_gas(), muon(), + // 1.f * unit::GeV, 2.133f), + std::make_tuple(helium_gas(), muon(), + 10.0f * unit::GeV, 2.768f), + std::make_tuple(helium_gas(), muon(), + 100.0f * unit::GeV, 3.188f))); + +INSTANTIATE_TEST_SUITE_P( + muon_stopping_power_Si, StoppingPowerValidation, + ::testing::Values( + std::make_tuple(silicon(), muon(), + 100.0f * unit::MeV, 1.849f), + std::make_tuple(silicon(), muon(), + 1.f * unit::GeV, 1.803f), + std::make_tuple(silicon(), muon(), + 10.0f * unit::GeV, 2.177f), + std::make_tuple(silicon(), muon(), + 100.0f * unit::GeV, 2.451f))); + +INSTANTIATE_TEST_SUITE_P( + anti_muon_stopping_power_He, StoppingPowerValidation, + ::testing::Values( + std::make_tuple(helium_gas(), antimuon(), + 100.0f * unit::MeV, 2.165f), + // std::make_tuple(helium_gas(), antimuon(), + // 1.f * unit::GeV, 2.133f), + std::make_tuple(helium_gas(), antimuon(), + 10.0f * unit::GeV, 2.768f), + std::make_tuple(helium_gas(), antimuon(), + 100.0f * unit::GeV, 3.188f))); + +INSTANTIATE_TEST_SUITE_P( + anti_muon_stopping_power_Si, StoppingPowerValidation, + ::testing::Values( + std::make_tuple(silicon(), antimuon(), + 100.0f * unit::MeV, 1.849f), + std::make_tuple(silicon(), antimuon(), + 1.f * unit::GeV, 1.803f), + std::make_tuple(silicon(), antimuon(), + 10.0f * unit::GeV, 2.177f), + std::make_tuple(silicon(), antimuon(), + 100.0f * unit::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(), + electron(), + 100.0f * unit::MeV, 3.532f), + std::make_tuple(helium_gas(), + electron(), + 1.f * unit::GeV, 13.14f))); + +INSTANTIATE_TEST_SUITE_P( + electron_stopping_power_Si, StoppingPowerValidation, + ::testing::Values(std::make_tuple(silicon(), + electron(), + 100.0f * unit::MeV, 6.017f), + std::make_tuple(silicon(), + electron(), + 1.f * unit::GeV, 46.69f))); + +INSTANTIATE_TEST_SUITE_P( + positron_stopping_power_He, StoppingPowerValidation, + ::testing::Values(std::make_tuple(helium_gas(), + positron(), + 100.0f * unit::MeV, 3.532f), + std::make_tuple(helium_gas(), + positron(), + 1.f * unit::GeV, 13.14f))); + +INSTANTIATE_TEST_SUITE_P( + positron_stopping_power_Si, StoppingPowerValidation, + ::testing::Values(std::make_tuple(silicon(), + positron(), + 100.0f * unit::MeV, 6.017f), + std::make_tuple(silicon(), + positron(), + 1.f * unit::GeV, 46.69f)));