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