Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
InteractionsTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file InteractionsTests.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <boost/test/data/test_case.hpp>
10 #include <boost/test/unit_test.hpp>
11 
12 #include "Acts/Definitions/PdgParticle.hpp"
13 #include "Acts/Definitions/Units.hpp"
14 #include "Acts/Material/Interactions.hpp"
15 #include "Acts/Material/Material.hpp"
16 #include "Acts/Material/MaterialSlab.hpp"
17 #include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp"
18 
19 #include <utility>
20 
21 namespace data = boost::unit_test::data;
22 using namespace Acts::UnitLiterals;
23 
24 // fixed material
25 static const Acts::Material material = Acts::Test::makeSilicon();
26 // variable values for other parameters
27 // thickness
28 static const double valuesThickness[] = {200_um, 1_mm};
29 static auto thickness = data::make(valuesThickness);
30 // particle type, mass, and charge
31 static const Acts::PdgParticle pdg[] = {Acts::eElectron, Acts::eMuon,
32  Acts::ePionPlus, Acts::eProton};
33 static const double mass[] = {511_keV, 105.7_MeV, 139.6_MeV, 938.3_MeV};
34 static const double charge[] = {-1_e, -1_e, 1_e, 1_e};
35 static const auto particle =
36  data::make(pdg) ^ data::make(mass) ^ data::make(charge);
37 // momentum range
38 static const auto momentum_low = data::xrange(100_MeV, 10_GeV, 100_MeV);
39 static const auto momentum_med = data::xrange(10_GeV, 100_GeV, 10_GeV);
40 static const auto momentum_high = data::xrange(100_GeV, 10_TeV, 100_GeV);
41 static const auto momentum = momentum_low + momentum_med + momentum_high;
42 
43 BOOST_AUTO_TEST_SUITE(interactions)
44 
45 // consistency checks for the energy loss values
46 BOOST_DATA_TEST_CASE(energy_loss_consistency, thickness* particle* momentum, x,
47  i, m, q, p) {
48  const auto slab = Acts::MaterialSlab(material, x);
49  const auto qOverP = q / p;
50  const auto absQ = std::abs(q);
51  const auto absPdg = Acts::makeAbsolutePdgParticle(i);
52 
53  auto dEBethe = computeEnergyLossBethe(slab, m, qOverP, absQ);
54  auto dELandau = computeEnergyLossLandau(slab, m, qOverP, absQ);
55  auto dELandauSigma = computeEnergyLossLandauSigma(slab, m, qOverP, absQ);
56  auto dELandauSigmaQOverP =
57  computeEnergyLossLandauSigmaQOverP(slab, m, qOverP, absQ);
58  auto dERad = computeEnergyLossRadiative(slab, absPdg, m, qOverP, absQ);
59  auto dEMean = computeEnergyLossMean(slab, absPdg, m, qOverP, absQ);
60  auto dEMode = computeEnergyLossMode(slab, absPdg, m, qOverP, absQ);
61 
62  BOOST_CHECK_LT(0, dEBethe);
63  BOOST_CHECK_LT(0, dELandau);
64  BOOST_CHECK_LT(0, dELandauSigma);
65  BOOST_CHECK_LT(0, dELandauSigmaQOverP);
66  BOOST_CHECK_LE(dELandauSigma, dEBethe);
67  // radiative terms only kick above some threshold -> can be zero
68  BOOST_CHECK_LE(0, dERad);
69  BOOST_CHECK_LT(0, dEMean);
70  BOOST_CHECK_LT(0, dEMode);
71  BOOST_CHECK_LE((dEBethe + dERad), dEMean);
72  // TODO verify mode/mean relation for full energy loss
73  // BOOST_CHECK_LE(dEMode, dEMean);
74 }
75 
76 // consistency checks for multiple scattering
77 BOOST_DATA_TEST_CASE(multiple_scattering_consistency,
78  thickness* particle* momentum, x, i, m, q, p) {
79  const auto slab = Acts::MaterialSlab(material, x);
80  const auto slabDoubled = Acts::MaterialSlab(material, 2 * x);
81  const auto qOverP = q / p;
82  const auto qOver2P = q / (2 * p);
83  const auto absQ = std::abs(q);
84  const auto absPdg = Acts::makeAbsolutePdgParticle(i);
85 
86  auto t0 = computeMultipleScatteringTheta0(slab, absPdg, m, qOverP, absQ);
87  BOOST_CHECK_LT(0, t0);
88  // use the anti-particle -> same scattering
89  auto tanti = computeMultipleScatteringTheta0(slab, absPdg, m, -qOverP, absQ);
90  BOOST_CHECK_LT(0, tanti);
91  BOOST_CHECK_EQUAL(t0, tanti);
92  // double the material -> more scattering
93  auto t2x =
94  computeMultipleScatteringTheta0(slabDoubled, absPdg, m, qOverP, absQ);
95  BOOST_CHECK_LT(0, t2x);
96  BOOST_CHECK_LT(t0, t2x);
97  // double the momentum -> less scattering
98  auto t2p = computeMultipleScatteringTheta0(slab, absPdg, m, qOver2P, absQ);
99  BOOST_CHECK_LT(0, t2p);
100  BOOST_CHECK_LT(t2p, t0);
101 }
102 
103 // no material -> no interactions
104 BOOST_DATA_TEST_CASE(vacuum, thickness* particle* momentum, x, i, m, q, p) {
105  const auto vacuum = Acts::MaterialSlab(Acts::Material(), x);
106  const auto qOverP = q / p;
107  const auto absQ = std::abs(q);
108  const auto absPdg = Acts::makeAbsolutePdgParticle(i);
109 
110  BOOST_CHECK_EQUAL(computeEnergyLossBethe(vacuum, m, qOverP, absQ), 0);
111  BOOST_CHECK_EQUAL(computeEnergyLossLandau(vacuum, m, qOverP, absQ), 0);
112  BOOST_CHECK_EQUAL(computeEnergyLossLandauSigma(vacuum, m, qOverP, absQ), 0);
113  BOOST_CHECK_EQUAL(computeEnergyLossLandauSigmaQOverP(vacuum, m, qOverP, absQ),
114  0);
115  BOOST_CHECK_EQUAL(computeEnergyLossRadiative(vacuum, absPdg, m, qOverP, absQ),
116  0);
117  BOOST_CHECK_EQUAL(computeEnergyLossMean(vacuum, absPdg, m, qOverP, absQ), 0);
118  BOOST_CHECK_EQUAL(computeEnergyLossMode(vacuum, absPdg, m, qOverP, absQ), 0);
119  BOOST_CHECK_EQUAL(
120  computeMultipleScatteringTheta0(vacuum, absPdg, m, qOverP, absQ), 0);
121 }
122 
123 // Silicon Bethe Energy Loss Validation
124 // PDG value from https://pdg.lbl.gov/2022/AtomicNuclearProperties
125 static const double momentum[] = {0.1003_GeV, 1.101_GeV, 10.11_GeV, 100.1_GeV};
126 static const double energy_loss[] = {2.608, 1.803, 2.177, 2.451};
127 
128 BOOST_DATA_TEST_CASE(silicon_energy_loss,
129  data::make(momentum) ^ data::make(energy_loss), p, loss) {
130  const Acts::Material silicon = Acts::Test::makeSilicon();
131  const auto thickness = 1_cm;
132  const auto slab = Acts::MaterialSlab(silicon, thickness);
133  const auto m = 105.7_MeV;
134  const auto q = -1_e;
135  const auto qOverP = q / p;
136  const auto absQ = std::abs(q);
137 
138  // Difference is within 5% from PDG value
139  BOOST_CHECK_CLOSE(computeEnergyLossBethe(slab, m, qOverP, absQ) / thickness /
140  slab.material().massDensity() / (1_MeV * 1_cm2 / 1_g),
141  loss, 5.);
142 }
143 
144 // Silicon Landau Energy Loss Validation
145 BOOST_AUTO_TEST_CASE(silicon_landau) {
146  const Acts::Material silicon = Acts::Test::makeSilicon();
147  const auto thickness = 0.17_cm;
148  const auto slab = Acts::MaterialSlab(silicon, thickness);
149  const float m = 105.7_MeV;
150  const float q = -1_e;
151  const float qOverP = q / 10_GeV;
152  const float absQ = std::abs(q);
153 
154  // Difference is within 5% from PDG value
155  const auto dE = computeEnergyLossLandau(slab, m, qOverP, absQ) / 1_MeV;
156  BOOST_CHECK_CLOSE(dE, 0.525, 5.);
157 
158  // Difference is within 10% from PDG value
159  const auto fwhm =
160  Acts::computeEnergyLossLandauFwhm(slab, m, qOverP, absQ) / 1_MeV;
161  BOOST_CHECK_CLOSE(fwhm, 0.13, 10.);
162 }
163 
164 BOOST_AUTO_TEST_SUITE_END()