Skip to content

Commit

Permalink
Merge pull request #827 from beomki-yeo/positron-brems
Browse files Browse the repository at this point in the history
Update electron/positron bremss
  • Loading branch information
beomki-yeo authored Dec 2, 2024
2 parents 7275508 + 4504708 commit b4fd39a
Show file tree
Hide file tree
Showing 5 changed files with 222 additions and 20 deletions.
20 changes: 16 additions & 4 deletions core/include/detray/materials/interaction.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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)
Expand Down
1 change: 1 addition & 0 deletions tests/unit_tests/cpu/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
2 changes: 1 addition & 1 deletion tests/unit_tests/cpu/material/bethe_equation.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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>(),
Expand Down
60 changes: 45 additions & 15 deletions tests/unit_tests/cpu/material/bremsstrahlung.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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)));
159 changes: 159 additions & 0 deletions tests/unit_tests/cpu/material/stopping_power.cpp
Original file line number Diff line number Diff line change
@@ -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)));

0 comments on commit b4fd39a

Please sign in to comment.