Core/Material/InteractionsTests.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 #include <boost/test/data/test_case.hpp>
0010 #include <boost/test/unit_test.hpp>
0011
0012 #include "Acts/Definitions/PdgParticle.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/Material/Interactions.hpp"
0015 #include "Acts/Material/Material.hpp"
0016 #include "Acts/Material/MaterialSlab.hpp"
0017 #include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp"
0018
0019 #include <utility>
0020
0021 namespace data = boost::unit_test::data;
0022 using namespace Acts::UnitLiterals;
0023
0024 // fixed material
0025 static const Acts::Material material = Acts::Test::makeSilicon();
0026 // variable values for other parameters
0027 // thickness
0028 static const double valuesThickness[] = {200_um, 1_mm};
0029 static auto thickness = data::make(valuesThickness);
0030 // particle type, mass, and charge
0031 static const Acts::PdgParticle pdg[] = {Acts::eElectron, Acts::eMuon,
0032                                         Acts::ePionPlus, Acts::eProton};
0033 static const double mass[] = {511_keV, 105.7_MeV, 139.6_MeV, 938.3_MeV};
0034 static const double charge[] = {-1_e, -1_e, 1_e, 1_e};
0035 static const auto particle =
0036     data::make(pdg) ^ data::make(mass) ^ data::make(charge);
0037 // momentum range
0038 static const auto momentum_low = data::xrange(100_MeV, 10_GeV, 100_MeV);
0039 static const auto momentum_med = data::xrange(10_GeV, 100_GeV, 10_GeV);
0040 static const auto momentum_high = data::xrange(100_GeV, 10_TeV, 100_GeV);
0041 static const auto momentum = momentum_low + momentum_med + momentum_high;
0042
0043 BOOST_AUTO_TEST_SUITE(interactions)
0044
0045 // consistency checks for the energy loss values
0046 BOOST_DATA_TEST_CASE(energy_loss_consistency, thickness* particle* momentum, x,
0047                      i, m, q, p) {
0048   const auto slab = Acts::MaterialSlab(material, x);
0049   const auto qOverP = q / p;
0050   const auto absQ = std::abs(q);
0051   const auto absPdg = Acts::makeAbsolutePdgParticle(i);
0052
0053   auto dEBethe = computeEnergyLossBethe(slab, m, qOverP, absQ);
0054   auto dELandau = computeEnergyLossLandau(slab, m, qOverP, absQ);
0055   auto dELandauSigma = computeEnergyLossLandauSigma(slab, m, qOverP, absQ);
0056   auto dELandauSigmaQOverP =
0057       computeEnergyLossLandauSigmaQOverP(slab, m, qOverP, absQ);
0058   auto dERad = computeEnergyLossRadiative(slab, absPdg, m, qOverP, absQ);
0059   auto dEMean = computeEnergyLossMean(slab, absPdg, m, qOverP, absQ);
0060   auto dEMode = computeEnergyLossMode(slab, absPdg, m, qOverP, absQ);
0061
0062   BOOST_CHECK_LT(0, dEBethe);
0063   BOOST_CHECK_LT(0, dELandau);
0064   BOOST_CHECK_LT(0, dELandauSigma);
0065   BOOST_CHECK_LT(0, dELandauSigmaQOverP);
0066   BOOST_CHECK_LE(dELandauSigma, dEBethe);
0067   // radiative terms only kick above some threshold -> can be zero
0068   BOOST_CHECK_LE(0, dERad);
0069   BOOST_CHECK_LT(0, dEMean);
0070   BOOST_CHECK_LT(0, dEMode);
0071   BOOST_CHECK_LE((dEBethe + dERad), dEMean);
0072   // TODO verify mode/mean relation for full energy loss
0073   // BOOST_CHECK_LE(dEMode, dEMean);
0074 }
0075
0076 // consistency checks for multiple scattering
0077 BOOST_DATA_TEST_CASE(multiple_scattering_consistency,
0078                      thickness* particle* momentum, x, i, m, q, p) {
0079   const auto slab = Acts::MaterialSlab(material, x);
0080   const auto slabDoubled = Acts::MaterialSlab(material, 2 * x);
0081   const auto qOverP = q / p;
0082   const auto qOver2P = q / (2 * p);
0083   const auto absQ = std::abs(q);
0084   const auto absPdg = Acts::makeAbsolutePdgParticle(i);
0085
0086   auto t0 = computeMultipleScatteringTheta0(slab, absPdg, m, qOverP, absQ);
0087   BOOST_CHECK_LT(0, t0);
0088   // use the anti-particle -> same scattering
0089   auto tanti = computeMultipleScatteringTheta0(slab, absPdg, m, -qOverP, absQ);
0090   BOOST_CHECK_LT(0, tanti);
0091   BOOST_CHECK_EQUAL(t0, tanti);
0092   // double the material -> more scattering
0093   auto t2x =
0094       computeMultipleScatteringTheta0(slabDoubled, absPdg, m, qOverP, absQ);
0095   BOOST_CHECK_LT(0, t2x);
0096   BOOST_CHECK_LT(t0, t2x);
0097   // double the momentum -> less scattering
0098   auto t2p = computeMultipleScatteringTheta0(slab, absPdg, m, qOver2P, absQ);
0099   BOOST_CHECK_LT(0, t2p);
0100   BOOST_CHECK_LT(t2p, t0);
0101 }
0102
0103 // no material -> no interactions
0104 BOOST_DATA_TEST_CASE(vacuum, thickness* particle* momentum, x, i, m, q, p) {
0105   const auto vacuum = Acts::MaterialSlab(Acts::Material(), x);
0106   const auto qOverP = q / p;
0107   const auto absQ = std::abs(q);
0108   const auto absPdg = Acts::makeAbsolutePdgParticle(i);
0109
0110   BOOST_CHECK_EQUAL(computeEnergyLossBethe(vacuum, m, qOverP, absQ), 0);
0111   BOOST_CHECK_EQUAL(computeEnergyLossLandau(vacuum, m, qOverP, absQ), 0);
0112   BOOST_CHECK_EQUAL(computeEnergyLossLandauSigma(vacuum, m, qOverP, absQ), 0);
0113   BOOST_CHECK_EQUAL(computeEnergyLossLandauSigmaQOverP(vacuum, m, qOverP, absQ),
0114                     0);
0115   BOOST_CHECK_EQUAL(computeEnergyLossRadiative(vacuum, absPdg, m, qOverP, absQ),
0116                     0);
0117   BOOST_CHECK_EQUAL(computeEnergyLossMean(vacuum, absPdg, m, qOverP, absQ), 0);
0118   BOOST_CHECK_EQUAL(computeEnergyLossMode(vacuum, absPdg, m, qOverP, absQ), 0);
0119   BOOST_CHECK_EQUAL(
0120       computeMultipleScatteringTheta0(vacuum, absPdg, m, qOverP, absQ), 0);
0121 }
0122
0123 // Silicon Bethe Energy Loss Validation
0124 // PDG value from https://pdg.lbl.gov/2022/AtomicNuclearProperties
0125 static const double momentum[] = {0.1003_GeV, 1.101_GeV, 10.11_GeV, 100.1_GeV};
0126 static const double energy_loss[] = {2.608, 1.803, 2.177, 2.451};
0127
0128 BOOST_DATA_TEST_CASE(silicon_energy_loss,
0129                      data::make(momentum) ^ data::make(energy_loss), p, loss) {
0130   const Acts::Material silicon = Acts::Test::makeSilicon();
0131   const auto thickness = 1_cm;
0132   const auto slab = Acts::MaterialSlab(silicon, thickness);
0133   const auto m = 105.7_MeV;
0134   const auto q = -1_e;
0135   const auto qOverP = q / p;
0136   const auto absQ = std::abs(q);
0137
0138   // Difference is within 5% from PDG value
0139   BOOST_CHECK_CLOSE(computeEnergyLossBethe(slab, m, qOverP, absQ) / thickness /
0140                         slab.material().massDensity() / (1_MeV * 1_cm2 / 1_g),
0141                     loss, 5.);
0142 }
0143
0144 // Silicon Landau Energy Loss Validation
0145 BOOST_AUTO_TEST_CASE(silicon_landau) {
0146   const Acts::Material silicon = Acts::Test::makeSilicon();
0147   const auto thickness = 0.17_cm;
0148   const auto slab = Acts::MaterialSlab(silicon, thickness);
0149   const float m = 105.7_MeV;
0150   const float q = -1_e;
0151   const float qOverP = q / 10_GeV;
0152   const float absQ = std::abs(q);
0153
0154   // Difference is within 5% from PDG value
0155   const auto dE = computeEnergyLossLandau(slab, m, qOverP, absQ) / 1_MeV;
0156   BOOST_CHECK_CLOSE(dE, 0.525, 5.);
0157
0158   // Difference is within 10% from PDG value
0159   const auto fwhm =
0160       Acts::computeEnergyLossLandauFwhm(slab, m, qOverP, absQ) / 1_MeV;
0161   BOOST_CHECK_CLOSE(fwhm, 0.13, 10.);
0162 }
0163
0164 BOOST_AUTO_TEST_SUITE_END()