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 "ActsFatras/EventData/Particle.hpp"
0015 #include "ActsFatras/Physics/ElectroMagnetic/Scattering.hpp"
0016 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0017 #include "ActsTests/CommonHelpers/PredefinedMaterials.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 ActsTests {
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, 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
0057 BOOST_AUTO_TEST_SUITE(PhysicsSuite)
0058
0059 BOOST_DATA_TEST_CASE(GeneralMixture, Dataset::parameters, pdg, phi, theta, p,
0060                      seed) {
0061   test(ActsFatras::GeneralMixtureScattering(), seed,
0062        Dataset::makeParticle(pdg, phi, theta, p));
0063 }
0064
0065 BOOST_DATA_TEST_CASE(GaussianMixture, Dataset::parameters, pdg, phi, theta, p,
0066                      seed) {
0067   test(ActsFatras::GaussianMixtureScattering(), seed,
0068        Dataset::makeParticle(pdg, phi, theta, p));
0069 }
0070
0071 BOOST_DATA_TEST_CASE(Highland, Dataset::parameters, pdg, phi, theta, p, seed) {
0072   test(ActsFatras::HighlandScattering(), seed,
0073        Dataset::makeParticle(pdg, phi, theta, p));
0074 }
0075
0076 BOOST_AUTO_TEST_CASE(HighlandRms) {
0077   auto scattering = ActsFatras::HighlandScattering();
0078   auto particle = Dataset::makeParticle(Acts::PdgParticle::eMuon, 0, 0, 1);
0079   auto materialSlab = makePercentSlab();
0080
0081   auto theta0 = Acts::computeMultipleScatteringTheta0(
0082       materialSlab, particle.absolutePdg(), particle.mass(), particle.qOverP(),
0083       particle.absoluteCharge());
0084
0085   std::ranlux48 gen(0);
0086
0087   std::vector<double> thetaYZs;
0088   std::vector<double> theta3Ds;
0089
0090   for (std::size_t i = 0; i < 10000; i++) {
0091     auto newParticle = particle;
0092     scattering(gen, materialSlab, newParticle);
0093
0094     double thetaYZ =
0095         std::atan2(newParticle.direction().y(), newParticle.direction().z());
0096     double theta3d =
0097         std::acos(newParticle.direction().dot(particle.direction()));
0098
0099     thetaYZs.push_back(thetaYZ);
0100     theta3Ds.push_back(theta3d);
0101   }
0102
0103   double rmsThetaYZ = rms(thetaYZs, 0);
0104   double rmsTheta3D = rms(theta3Ds, 0);
0105
0106   CHECK_CLOSE_REL(rmsThetaYZ, theta0, 0.02);
0107   CHECK_CLOSE_REL(rmsTheta3D, std::numbers::sqrt2 * theta0, 0.02);
0108 }
0109
0110 BOOST_AUTO_TEST_SUITE_END()
0111
0112 }  // namespace ActsTests