Fatras/Physics/ScatteringTests.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/Material/Interactions.hpp"
0014 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0015 #include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp"
0016 #include "ActsFatras/EventData/Particle.hpp"
0017 #include "ActsFatras/Physics/ElectroMagnetic/Scattering.hpp"
0018
0019 #include <cmath>
0020 #include <cstdint>
0021 #include <limits>
0022 #include <numbers>
0023 #include <random>
0024
0025 #include "Dataset.hpp"
0026
0027 namespace {
0028
0029 constexpr auto eps = std::numeric_limits<double>::epsilon();
0030
0031 double rms(const std::vector<double>& values, double mean) {
0032   return std::sqrt(std::accumulate(values.begin(), values.end(), 0.0,
0033                                    [mean](double sum, double value) {
0034                                      return sum +
0035                                             (value - mean) * (value - mean);
0036                                    }) /
0037                    values.size());
0038 }
0039
0040 // Common test method that will be instantiated for each scattering model.
0041 template <typename Scattering>
0042 void test(const Scattering& scattering, std::uint32_t seed,
0043           const ActsFatras::Particle& before) {
0044   std::ranlux48 gen(seed);
0045   ActsFatras::Particle after = before;
0046
0047   const auto outgoing = scattering(gen, Acts::Test::makePercentSlab(), after);
0048   // scattering leaves absolute energy/momentum unchanged
0049   CHECK_CLOSE_REL(after.absoluteMomentum(), before.absoluteMomentum(), eps);
0050   CHECK_CLOSE_REL(after.energy(), before.energy(), eps);
0051   // scattering has changed the direction
0052   BOOST_CHECK_LT(before.direction().dot(after.direction()), 1);
0053   // scattering creates no new particles
0054   BOOST_CHECK(outgoing.empty());
0055 }
0056 }  // namespace
0057
0058 BOOST_AUTO_TEST_SUITE(FatrasScattering)
0059
0060 BOOST_DATA_TEST_CASE(GeneralMixture, Dataset::parameters, pdg, phi, theta, p,
0061                      seed) {
0062   test(ActsFatras::GeneralMixtureScattering(), seed,
0063        Dataset::makeParticle(pdg, phi, theta, p));
0064 }
0065
0066 BOOST_DATA_TEST_CASE(GaussianMixture, Dataset::parameters, pdg, phi, theta, p,
0067                      seed) {
0068   test(ActsFatras::GaussianMixtureScattering(), seed,
0069        Dataset::makeParticle(pdg, phi, theta, p));
0070 }
0071
0072 BOOST_DATA_TEST_CASE(Highland, Dataset::parameters, pdg, phi, theta, p, seed) {
0073   test(ActsFatras::HighlandScattering(), seed,
0074        Dataset::makeParticle(pdg, phi, theta, p));
0075 }
0076
0077 BOOST_AUTO_TEST_CASE(HighlandRms) {
0078   auto scattering = ActsFatras::HighlandScattering();
0079   auto particle = Dataset::makeParticle(Acts::PdgParticle::eMuon, 0, 0, 1);
0080   auto materialSlab = Acts::Test::makePercentSlab();
0081
0082   auto theta0 = Acts::computeMultipleScatteringTheta0(
0083       materialSlab, particle.absolutePdg(), particle.mass(), particle.qOverP(),
0084       particle.absoluteCharge());
0085
0086   std::ranlux48 gen(0);
0087
0088   std::vector<double> thetaYZs;
0089   std::vector<double> theta3Ds;
0090
0091   for (std::size_t i = 0; i < 10000; i++) {
0092     auto newParticle = particle;
0093     scattering(gen, materialSlab, newParticle);
0094
0095     double thetaYZ =
0096         std::atan2(newParticle.direction().y(), newParticle.direction().z());
0097     double theta3d =
0098         std::acos(newParticle.direction().dot(particle.direction()));
0099
0100     thetaYZs.push_back(thetaYZ);
0101     theta3Ds.push_back(theta3d);
0102   }
0103
0104   double rmsThetaYZ = rms(thetaYZs, 0);
0105   double rmsTheta3D = rms(theta3Ds, 0);
0106
0107   CHECK_CLOSE_REL(rmsThetaYZ, theta0, 0.02);
0108   CHECK_CLOSE_REL(rmsTheta3D, std::numbers::sqrt2 * theta0, 0.02);
0109 }
0110
0111 BOOST_AUTO_TEST_SUITE_END()