Skip to content

Commit

Permalink
Add tests
Browse files Browse the repository at this point in the history
  • Loading branch information
beomki-yeo committed Nov 22, 2024
1 parent 7ec2bfb commit c66745a
Show file tree
Hide file tree
Showing 4 changed files with 193 additions and 16 deletions.
2 changes: 1 addition & 1 deletion 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
1 change: 1 addition & 0 deletions tests/unit_tests/cpu/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,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
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)));
146 changes: 146 additions & 0 deletions tests/unit_tests/cpu/material/stopping_power.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,146 @@
/** 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<scalar>, pdg_particle<scalar>, scalar, scalar>> {
};

TEST_P(StoppingPowerValidation, stopping_power) {

// Interaction object
interaction<scalar> I;

// Material
material<scalar> mat = std::get<0>(GetParam());

// Particle
pdg_particle<scalar> ptc = std::get<1>(GetParam());

// Kinetic energy
const scalar T = std::get<2>(GetParam());

// Total energy
const scalar E = T + ptc.mass();

// Momentum
const scalar p = math::sqrt(E * E - ptc.mass() * ptc.mass());

// qoverp
const scalar qop{ptc.charge() / p};

// Stopping power in MeV * cm^2 / g
const scalar dEdx{
I.compute_stopping_power(mat, ptc, {ptc, qop}) / mat.mass_density() /
(unit<scalar>::MeV * unit<scalar>::cm2 / unit<scalar>::g)};

const 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 that we took the PDG value only from Ionization loss (Radiative loss is
// ignored)
// Also assumes that the stopping powers of muon and antimuon are the same
INSTANTIATE_TEST_SUITE_P(
muon_stopping_power_He, StoppingPowerValidation,
::testing::Values(std::make_tuple(helium_gas<scalar>(), muon<scalar>(),
100.0f * unit<scalar>::MeV, 2.165f),
std::make_tuple(helium_gas<scalar>(), muon<scalar>(),
1.f * unit<scalar>::GeV, 2.133f),
std::make_tuple(helium_gas<scalar>(), muon<scalar>(),
10.0f * unit<scalar>::GeV, 2.768f),
std::make_tuple(helium_gas<scalar>(), muon<scalar>(),
100.0f * unit<scalar>::GeV, 3.188f)));

INSTANTIATE_TEST_SUITE_P(
muon_stopping_power_Si, StoppingPowerValidation,
::testing::Values(std::make_tuple(silicon<scalar>(), muon<scalar>(),
100.0f * unit<scalar>::MeV, 1.849f),
std::make_tuple(silicon<scalar>(), muon<scalar>(),
1.f * unit<scalar>::GeV, 1.803f),
std::make_tuple(silicon<scalar>(), muon<scalar>(),
10.0f * unit<scalar>::GeV, 2.177f),
std::make_tuple(silicon<scalar>(), muon<scalar>(),
100.0f * unit<scalar>::GeV, 2.451f)));

INSTANTIATE_TEST_SUITE_P(
anti_muon_stopping_power_He, StoppingPowerValidation,
::testing::Values(std::make_tuple(helium_gas<scalar>(), antimuon<scalar>(),
100.0f * unit<scalar>::MeV, 2.165f),
std::make_tuple(helium_gas<scalar>(), antimuon<scalar>(),
1.f * unit<scalar>::GeV, 2.133f),
std::make_tuple(helium_gas<scalar>(), antimuon<scalar>(),
10.0f * unit<scalar>::GeV, 2.768f),
std::make_tuple(helium_gas<scalar>(), antimuon<scalar>(),
100.0f * unit<scalar>::GeV, 3.188f)));

INSTANTIATE_TEST_SUITE_P(
anti_muon_stopping_power_Si, StoppingPowerValidation,
::testing::Values(std::make_tuple(silicon<scalar>(), antimuon<scalar>(),
100.0f * unit<scalar>::MeV, 1.849f),
std::make_tuple(silicon<scalar>(), antimuon<scalar>(),
1.f * unit<scalar>::GeV, 1.803f),
std::make_tuple(silicon<scalar>(), antimuon<scalar>(),
10.0f * unit<scalar>::GeV, 2.177f),
std::make_tuple(silicon<scalar>(), antimuon<scalar>(),
100.0f * unit<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<scalar>(), electron<scalar>(),
100.0f * unit<scalar>::MeV, 3.532f),
std::make_tuple(helium_gas<scalar>(), electron<scalar>(),
1.f * unit<scalar>::GeV, 13.14f)));

INSTANTIATE_TEST_SUITE_P(
electron_stopping_power_Si, StoppingPowerValidation,
::testing::Values(std::make_tuple(silicon<scalar>(), electron<scalar>(),
100.0f * unit<scalar>::MeV, 6.017f),
std::make_tuple(silicon<scalar>(), electron<scalar>(),
1.f * unit<scalar>::GeV, 46.69f)));

INSTANTIATE_TEST_SUITE_P(
positron_stopping_power_He, StoppingPowerValidation,
::testing::Values(std::make_tuple(helium_gas<scalar>(), positron<scalar>(),
100.0f * unit<scalar>::MeV, 3.532f),
std::make_tuple(helium_gas<scalar>(), positron<scalar>(),
1.f * unit<scalar>::GeV, 13.14f)));

INSTANTIATE_TEST_SUITE_P(
positron_stopping_power_Si, StoppingPowerValidation,
::testing::Values(std::make_tuple(silicon<scalar>(), positron<scalar>(),
100.0f * unit<scalar>::MeV, 6.017f),
std::make_tuple(silicon<scalar>(), positron<scalar>(),
1.f * unit<scalar>::GeV, 46.69f)));

0 comments on commit c66745a

Please sign in to comment.