cpu/material/stopping_power.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 // Project include(s).
0010 #include "detray/material/interaction.hpp"
0011 #include "detray/material/material.hpp"
0012 #include "detray/material/predefined_materials.hpp"
0013
0014 // Detray test include(s)
0015 #include "detray/test/framework/types.hpp"
0016
0017 // GTest include(s).
0018 #include <gtest/gtest.h>
0019
0020 using namespace detray;
0021
0022 // Test class for the stopping power
0023 // Input tuple: < material, particle type, kinetic energy, expected output >
0024 class StoppingPowerValidation
0025     : public ::testing::TestWithParam<
0026           std::tuple<material<test::scalar>, pdg_particle<test::scalar>,
0027                      test::scalar, test::scalar>> {};
0028
0029 TEST_P(StoppingPowerValidation, stopping_power) {
0030   // Interaction object
0031   interaction<test::scalar> I;
0032
0033   // Material
0034   material<test::scalar> mat = std::get<0>(GetParam());
0035
0036   // Particle
0037   pdg_particle<test::scalar> ptc = std::get<1>(GetParam());
0038
0039   // Kinetic energy
0040   const test::scalar T = std::get<2>(GetParam());
0041
0042   // Total energy
0043   const test::scalar E = T + ptc.mass();
0044
0045   // Momentum
0046   const test::scalar p = math::sqrt(E * E - ptc.mass() * ptc.mass());
0047
0048   // qoverp
0049   const test::scalar qop{ptc.charge() / p};
0050
0051   // Stopping power in MeV * cm^2 / g
0052   const test::scalar dEdx{I.compute_stopping_power(mat, ptc, {ptc, qop}) /
0053                           mat.mass_density() /
0054                           (unit<test::scalar>::MeV * unit<test::scalar>::cm2 /
0055                            unit<test::scalar>::g)};
0056
0057   const test::scalar expected_dEdx = std::get<3>(GetParam());
0058
0059   // Check if difference is within 8% error
0060   EXPECT_NEAR((expected_dEdx - dEdx) / dEdx, 0.f, 0.08f);
0061 }
0062
0063 /******************
0064  *   Muon tests
0065  ******************/
0066
0067 // From https://pdg.lbl.gov/2024/AtomicNuclearProperties/index.html
0068 // Note 1: that we took the PDG value only from Ionization loss (Radiative loss
0069 // is ignored) Note 2: assumes that the stopping powers of muon and antimuon are
0070 // the same Note 3: Test fails with He Gas and 1 GeV muons (18 % difference)
0071 INSTANTIATE_TEST_SUITE_P(
0072     muon_stopping_power_He, StoppingPowerValidation,
0073     ::testing::Values(
0074         std::make_tuple(helium_gas<test::scalar>(), muon<test::scalar>(),
0075                         100.0f * unit<test::scalar>::MeV, 2.165f),
0076         // std::make_tuple(helium_gas<test::scalar>(), muon<test::scalar>(),
0077         //                 1.f * unit<test::scalar>::GeV, 2.133f),
0078         std::make_tuple(helium_gas<test::scalar>(), muon<test::scalar>(),
0079                         10.0f * unit<test::scalar>::GeV, 2.768f),
0080         std::make_tuple(helium_gas<test::scalar>(), muon<test::scalar>(),
0081                         100.0f * unit<test::scalar>::GeV, 3.188f)));
0082
0083 INSTANTIATE_TEST_SUITE_P(
0084     muon_stopping_power_Si, StoppingPowerValidation,
0085     ::testing::Values(
0086         std::make_tuple(silicon<test::scalar>(), muon<test::scalar>(),
0087                         100.0f * unit<test::scalar>::MeV, 1.849f),
0088         std::make_tuple(silicon<test::scalar>(), muon<test::scalar>(),
0089                         1.f * unit<test::scalar>::GeV, 1.803f),
0090         std::make_tuple(silicon<test::scalar>(), muon<test::scalar>(),
0091                         10.0f * unit<test::scalar>::GeV, 2.177f),
0092         std::make_tuple(silicon<test::scalar>(), muon<test::scalar>(),
0093                         100.0f * unit<test::scalar>::GeV, 2.451f)));
0094
0095 INSTANTIATE_TEST_SUITE_P(
0096     anti_muon_stopping_power_He, StoppingPowerValidation,
0097     ::testing::Values(
0098         std::make_tuple(helium_gas<test::scalar>(), antimuon<test::scalar>(),
0099                         100.0f * unit<test::scalar>::MeV, 2.165f),
0100         // std::make_tuple(helium_gas<test::scalar>(), antimuon<test::scalar>(),
0101         //                 1.f * unit<test::scalar>::GeV, 2.133f),
0102         std::make_tuple(helium_gas<test::scalar>(), antimuon<test::scalar>(),
0103                         10.0f * unit<test::scalar>::GeV, 2.768f),
0104         std::make_tuple(helium_gas<test::scalar>(), antimuon<test::scalar>(),
0105                         100.0f * unit<test::scalar>::GeV, 3.188f)));
0106
0107 INSTANTIATE_TEST_SUITE_P(
0108     anti_muon_stopping_power_Si, StoppingPowerValidation,
0109     ::testing::Values(
0110         std::make_tuple(silicon<test::scalar>(), antimuon<test::scalar>(),
0111                         100.0f * unit<test::scalar>::MeV, 1.849f),
0112         std::make_tuple(silicon<test::scalar>(), antimuon<test::scalar>(),
0113                         1.f * unit<test::scalar>::GeV, 1.803f),
0114         std::make_tuple(silicon<test::scalar>(), antimuon<test::scalar>(),
0115                         10.0f * unit<test::scalar>::GeV, 2.177f),
0116         std::make_tuple(silicon<test::scalar>(), antimuon<test::scalar>(),
0117                         100.0f * unit<test::scalar>::GeV, 2.451f)));
0118
0119 /*********************
0120  *   Electron tests
0121  *********************/
0122
0123 // From https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html
0124 // Assumes that the stopping powers of electron and positron are the same
0125 INSTANTIATE_TEST_SUITE_P(
0126     electron_stopping_power_He, StoppingPowerValidation,
0127     ::testing::Values(std::make_tuple(helium_gas<test::scalar>(),
0128                                       electron<test::scalar>(),
0129                                       100.0f * unit<test::scalar>::MeV, 3.532f),
0130                       std::make_tuple(helium_gas<test::scalar>(),
0131                                       electron<test::scalar>(),
0132                                       1.f * unit<test::scalar>::GeV, 13.14f)));
0133
0134 INSTANTIATE_TEST_SUITE_P(
0135     electron_stopping_power_Si, StoppingPowerValidation,
0136     ::testing::Values(std::make_tuple(silicon<test::scalar>(),
0137                                       electron<test::scalar>(),
0138                                       100.0f * unit<test::scalar>::MeV, 6.017f),
0139                       std::make_tuple(silicon<test::scalar>(),
0140                                       electron<test::scalar>(),
0141                                       1.f * unit<test::scalar>::GeV, 46.69f)));
0142
0143 INSTANTIATE_TEST_SUITE_P(
0144     positron_stopping_power_He, StoppingPowerValidation,
0145     ::testing::Values(std::make_tuple(helium_gas<test::scalar>(),
0146                                       positron<test::scalar>(),
0147                                       100.0f * unit<test::scalar>::MeV, 3.532f),
0148                       std::make_tuple(helium_gas<test::scalar>(),
0149                                       positron<test::scalar>(),
0150                                       1.f * unit<test::scalar>::GeV, 13.14f)));
0151
0152 INSTANTIATE_TEST_SUITE_P(
0153     positron_stopping_power_Si, StoppingPowerValidation,
0154     ::testing::Values(std::make_tuple(silicon<test::scalar>(),
0155                                       positron<test::scalar>(),
0156                                       100.0f * unit<test::scalar>::MeV, 6.017f),
0157                       std::make_tuple(silicon<test::scalar>(),
0158                                       positron<test::scalar>(),
0159                                       1.f * unit<test::scalar>::GeV, 46.69f)));