cpu/material/bremsstrahlung.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/definitions/pdg_particle.hpp"
0011 #include "detray/material/interaction.hpp"
0012 #include "detray/material/material.hpp"
0013 #include "detray/material/predefined_materials.hpp"
0014
0015 // Detray test include(s)
0016 #include "detray/test/framework/types.hpp"
0017
0018 // GTest include(s).
0019 #include <gtest/gtest.h>
0020
0021 using namespace detray;
0022
0023 using scalar = test::scalar;
0024
0025 // Test class for energy loss with Bremsstrahlung
0026 // Input tuple: < material, particle type, kinetic energy, expected output >
0027 class EnergyLossBremsValidation
0028     : public ::testing::TestWithParam<
0029           std::tuple<material<scalar>, pdg_particle<scalar>, scalar, scalar>> {
0030 };
0031
0032 // This tests the material functionalities
0033 TEST_P(EnergyLossBremsValidation, bremsstrahlung) {
0034   // Interaction object
0035   interaction<scalar> I;
0036
0037   // Material
0038   material<scalar> mat = std::get<0>(GetParam());
0039
0040   // Particle
0041   pdg_particle<scalar> ptc = std::get<1>(GetParam());
0042
0043   // Kinetic energy
0044   const scalar T = std::get<2>(GetParam());
0045
0046   // Total energy
0047   const scalar E = T + ptc.mass();
0048
0049   // Momentum
0050   const scalar p = math::sqrt(E * E - ptc.mass() * ptc.mass());
0051
0052   // qoverp
0053   const scalar qop{ptc.charge() / p};
0054
0055   // Bremsstrahlung stopping power in MeV * cm^2 / g
0056   const scalar dEdx{I.compute_bremsstrahlung(mat, ptc, {ptc, qop}) /
0057                     mat.mass_density() /
0058                     (unit<scalar>::MeV * unit<scalar>::cm2 / unit<scalar>::g)};
0059
0060   const scalar expected_dEdx = std::get<3>(GetParam());
0061
0062   // We have not implemented the bremsstrahlung for the heavier charged
0063   // particles, which is negligible
0064   if ((ptc.pdg_num() == electron<scalar>().pdg_num()) ||
0065       (ptc.pdg_num() == positron<scalar>().pdg_num())) {
0066     // Check if difference is within 14% error
0067     EXPECT_NEAR((expected_dEdx - dEdx) / dEdx, 0.f, 0.14f);
0068   } else {
0069     EXPECT_FLOAT_EQ(static_cast<float>(dEdx), 0.f);
0070   }
0071 }
0072
0073 // From https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html
0074 // Assumes that the stopping powers of electron and positron are the same
0075
0076 // Electrons
0077 INSTANTIATE_TEST_SUITE_P(
0078     electron_Bremsstrahlung_100MeV_He, EnergyLossBremsValidation,
0079     ::testing::Values(std::make_tuple(helium_gas<scalar>(), electron<scalar>(),
0080                                       100.0f * unit<scalar>::MeV, 0.9229f)));
0081
0082 INSTANTIATE_TEST_SUITE_P(
0083     electron_Bremsstrahlung_100MeV_Al, EnergyLossBremsValidation,
0084     ::testing::Values(std::make_tuple(aluminium<scalar>(), electron<scalar>(),
0085                                       100.0f * unit<scalar>::MeV, 3.714f)));
0086
0087 INSTANTIATE_TEST_SUITE_P(
0088     electron_Bremsstrahlung_100MeV_Si, EnergyLossBremsValidation,
0089     ::testing::Values(std::make_tuple(silicon<scalar>(), electron<scalar>(),
0090                                       100.0f * unit<scalar>::MeV, 4.099f)));
0091
0092 INSTANTIATE_TEST_SUITE_P(
0093     electron_Bremsstrahlung_100MeV_Cu, EnergyLossBremsValidation,
0094     ::testing::Values(std::make_tuple(copper<scalar>(), electron<scalar>(),
0095                                       100.0f * unit<scalar>::MeV, 7.079f)));
0096
0097 // Positrons
0098 INSTANTIATE_TEST_SUITE_P(
0099     positron_Bremsstrahlung_100MeV_He, EnergyLossBremsValidation,
0100     ::testing::Values(std::make_tuple(helium_gas<scalar>(), positron<scalar>(),
0101                                       100.0f * unit<scalar>::MeV, 0.9229f)));
0102
0103 INSTANTIATE_TEST_SUITE_P(
0104     positron_Bremsstrahlung_100MeV_Al, EnergyLossBremsValidation,
0105     ::testing::Values(std::make_tuple(aluminium<scalar>(), positron<scalar>(),
0106                                       100.0f * unit<scalar>::MeV, 3.714f)));
0107
0108 INSTANTIATE_TEST_SUITE_P(
0109     positron_Bremsstrahlung_100MeV_Si, EnergyLossBremsValidation,
0110     ::testing::Values(std::make_tuple(silicon<scalar>(), positron<scalar>(),
0111                                       100.0f * unit<scalar>::MeV, 4.099f)));
0112
0113 INSTANTIATE_TEST_SUITE_P(
0114     positron_Bremsstrahlung_100MeV_Cu, EnergyLossBremsValidation,
0115     ::testing::Values(std::make_tuple(copper<scalar>(), positron<scalar>(),
0116                                       100.0f * unit<scalar>::MeV, 7.079f)));
0117
0118 // We have not implemented the bremsstrahlung for muons
0119 INSTANTIATE_TEST_SUITE_P(
0120     muon_Bremsstrahlung_100MeV_Cu_muon, EnergyLossBremsValidation,
0121     ::testing::Values(std::make_tuple(copper<scalar>(), muon<scalar>(),
0122                                       100.0f * unit<scalar>::MeV, 0.f)));
0123
0124 INSTANTIATE_TEST_SUITE_P(
0125     antimuon_Bremsstrahlung_100MeV_Cu_muon, EnergyLossBremsValidation,
0126     ::testing::Values(std::make_tuple(copper<scalar>(), antimuon<scalar>(),
0127                                       100.0f * unit<scalar>::MeV, 0.f)));