cpu/material/stopping_power_derivative.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 using scalar = test::scalar;
0023
0024 GTEST_TEST(detray_material, derivative_test_beta2) {
0025   const pdg_particle<scalar> ptc = muon<scalar>();
0026
0027   // mass
0028   const scalar m = ptc.mass();
0029
0030   // charge
0031   const scalar q = ptc.charge();
0032
0033   // Displacement for numerical differentiaion
0034   const scalar h = 0.001f;
0035
0036   // Iterate from 1 GeV to 10 GeV
0037   for (unsigned int i = 1u; i < 10; i++) {
0038     const scalar p = static_cast<scalar>(i) * detray::unit<scalar>::GeV;
0039     const scalar qop = q / p;
0040
0041     detail::relativistic_quantities rq(m, qop, q);
0042     detail::relativistic_quantities rq1(m, qop + h, q);
0043     detail::relativistic_quantities rq2(m, qop - h, q);
0044
0045     const scalar numerical = (rq1.m_beta2 - rq2.m_beta2) / (2.f * h);
0046
0047     const scalar evaluated = rq.derive_beta2();
0048
0049     EXPECT_NEAR((numerical - evaluated) / numerical, 0, 0.012);
0050   }
0051 }
0052
0053 // Test class for the derivative of bethe stopping power
0054 class DerivativeOfBetheEquationValidation
0055     : public ::testing::TestWithParam<
0056           std::tuple<pdg_particle<scalar>, material<scalar>>> {};
0057
0058 TEST_P(DerivativeOfBetheEquationValidation,
0059        derivative_of_bethe_stopping_power) {
0060   // Interaction object
0061   interaction<scalar> Interactor;
0062
0063   const pdg_particle<scalar> ptc = std::get<0>(GetParam());
0064
0065   // Material
0066   const auto mat = std::get<1>(GetParam());
0067
0068   // mass
0069   const scalar m = ptc.mass();
0070
0071   // charge
0072   const scalar q = ptc.charge();
0073
0074   // Displacement for numerical differentiaion
0075   const scalar h = 1e-3f;
0076
0077   // Iterate from 2 GeV to 100 GeV
0078   for (unsigned int i = 2u; i < 100; i++) {
0079     const scalar p = static_cast<scalar>(i) * detray::unit<scalar>::GeV;
0080     const scalar qop = q / p;
0081
0082     const detail::relativistic_quantities<scalar> rq(m, qop, q);
0083     const detail::relativistic_quantities<scalar> rq1(m, qop + h, q);
0084     const detail::relativistic_quantities<scalar> rq2(m, qop - h, q);
0085
0086     // Log term
0087     const scalar log_term1 = rq1.compute_bethe_bloch_log_term(mat);
0088     const scalar log_term2 = rq2.compute_bethe_bloch_log_term(mat);
0089
0090     const scalar numerical_log_term = (log_term1 - log_term2) / (2.f * h);
0091     const scalar evaluated_log_term = rq.derive_bethe_bloch_log_term();
0092
0093     EXPECT_NEAR(numerical_log_term, evaluated_log_term,
0094                 numerical_log_term * 0.01f);
0095
0096     // delta half
0097     const scalar dhalf1 = rq1.compute_delta_half(mat);
0098     const scalar dhalf2 = rq2.compute_delta_half(mat);
0099
0100     const scalar numerical_dhalf = (dhalf1 - dhalf2) / (2.f * h);
0101     const scalar evaluated_dhalf = rq.derive_delta_half(mat);
0102
0103     EXPECT_NEAR(numerical_dhalf, evaluated_dhalf, numerical_dhalf * 0.01f);
0104
0105     // Bethe equation
0106     const scalar bethe1 = Interactor.compute_bethe_bloch(mat, ptc, rq1);
0107     const scalar bethe2 = Interactor.compute_bethe_bloch(mat, ptc, rq2);
0108
0109     const scalar numerical_bethe = (bethe1 - bethe2) / (2.f * h);
0110
0111     const scalar evaluated_bethe = Interactor.derive_bethe_bloch(mat, ptc, rq);
0112
0113     EXPECT_NEAR(numerical_bethe, evaluated_bethe, numerical_bethe * 0.01f);
0114   }
0115 }
0116
0117 INSTANTIATE_TEST_SUITE_P(
0118     detray_material_BetheEquationDerivative,
0119     DerivativeOfBetheEquationValidation,
0120     ::testing::Values(
0121         std::make_tuple(muon<scalar>(), hydrogen_gas<scalar>()),
0122         std::make_tuple(muon<scalar>(), helium_gas<scalar>()),
0123         std::make_tuple(muon<scalar>(), isobutane<scalar>()),
0124         std::make_tuple(muon<scalar>(), aluminium<scalar>()),
0125         std::make_tuple(muon<scalar>(), silicon<scalar>()),
0126         std::make_tuple(muon<scalar>(), tungsten<scalar>()),
0127         std::make_tuple(muon<scalar>(), gold<scalar>()),
0128         std::make_tuple(muon<scalar>(), cesium_iodide<scalar>()),
0129         std::make_tuple(muon<scalar>(), silicon_with_ded<scalar>()),
0130         std::make_tuple(muon<scalar>(), cesium_iodide_with_ded<scalar>())));
0131
0132 // Test class for the derivative of bremsstrahlung stopping power
0133 class DerivativeOfBremsstrahlungValidation
0134     : public ::testing::TestWithParam<
0135           std::tuple<pdg_particle<scalar>, material<scalar>>> {};
0136
0137 TEST_P(DerivativeOfBremsstrahlungValidation,
0138        derivative_of_brems_stopping_power) {
0139   // Interaction object
0140   interaction<scalar> Interactor;
0141
0142   const pdg_particle<scalar> ptc = std::get<0>(GetParam());
0143
0144   // Material
0145   const auto mat = std::get<1>(GetParam());
0146
0147   // mass
0148   const scalar m = ptc.mass();
0149
0150   // charge
0151   const scalar q = ptc.charge();
0152
0153   // Displacement for numerical differentiaion
0154   const scalar h = 1e-3f;
0155
0156   // Iterate from 1 GeV to 100 GeV
0157   for (unsigned int i = 1u; i < 100; i++) {
0158     const scalar p = static_cast<scalar>(i) * detray::unit<scalar>::GeV;
0159     const scalar qop = q / p;
0160
0161     const detail::relativistic_quantities<scalar> rq(m, qop, q);
0162     const detail::relativistic_quantities<scalar> rq1(m, qop + h, q);
0163     const detail::relativistic_quantities<scalar> rq2(m, qop - h, q);
0164
0165     // Brems equation
0166     const scalar brems1 = Interactor.compute_bremsstrahlung(mat, ptc, rq1);
0167     const scalar brems2 = Interactor.compute_bremsstrahlung(mat, ptc, rq2);
0168
0169     const scalar numerical_brems = (brems1 - brems2) / (2.f * h);
0170
0171     const scalar evaluated_brems =
0172         Interactor.derive_bremsstrahlung(mat, ptc, rq);
0173
0174     EXPECT_NEAR(numerical_brems, evaluated_brems, numerical_brems * 0.01f);
0175   }
0176 }
0177
0178 INSTANTIATE_TEST_SUITE_P(
0179     detray_material_BremsstrahlungDerivative,
0180     DerivativeOfBremsstrahlungValidation,
0181     ::testing::Values(
0182         std::make_tuple(electron<scalar>(), hydrogen_gas<scalar>()),
0183         std::make_tuple(electron<scalar>(), helium_gas<scalar>()),
0184         std::make_tuple(electron<scalar>(), isobutane<scalar>()),
0185         std::make_tuple(electron<scalar>(), aluminium<scalar>()),
0186         std::make_tuple(electron<scalar>(), silicon<scalar>()),
0187         std::make_tuple(electron<scalar>(), tungsten<scalar>()),
0188         std::make_tuple(electron<scalar>(), gold<scalar>()),
0189         std::make_tuple(electron<scalar>(), cesium_iodide<scalar>()),
0190         std::make_tuple(electron<scalar>(), silicon_with_ded<scalar>()),
0191         std::make_tuple(electron<scalar>(), cesium_iodide_with_ded<scalar>())));