IntegrationTests/Fatras/FatrasSimulationTests.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/ParticleData.hpp"
0013 #include "Acts/MagneticField/ConstantBField.hpp"
0014 #include "Acts/Propagator/EigenStepper.hpp"
0015 #include "Acts/Propagator/Navigator.hpp"
0016 #include "Acts/Propagator/Propagator.hpp"
0017 #include "Acts/Propagator/StraightLineStepper.hpp"
0018 #include "Acts/Tests/CommonHelpers/CylindricalTrackingGeometry.hpp"
0019 #include "Acts/Utilities/Logger.hpp"
0020 #include "Acts/Utilities/UnitVectors.hpp"
0021 #include "ActsFatras/Kernel/InteractionList.hpp"
0022 #include "ActsFatras/Kernel/Simulation.hpp"
0023 #include "ActsFatras/Physics/Decay/NoDecay.hpp"
0024 #include "ActsFatras/Physics/StandardInteractions.hpp"
0025 #include "ActsFatras/Selectors/ParticleSelectors.hpp"
0026 #include "ActsFatras/Selectors/SurfaceSelectors.hpp"
0027
0028 #include <algorithm>
0029 #include <random>
0030
0031 using namespace Acts::UnitLiterals;
0032
0033 namespace {
0034
0035 /// Mock-up process that splits a particle into two above a momentum threshold.
0036 struct SplitEnergyLoss {
0037   double splitMomentumMin = 5_GeV;
0038
0039   template <typename generator_t>
0040   bool operator()(generator_t& /*generator*/,
0041                   const Acts::MaterialSlab& /*slab*/,
0042                   ActsFatras::Particle& particle,
0043                   std::vector<ActsFatras::Particle>& generated) const {
0044     const auto p = particle.absoluteMomentum();
0045     if (splitMomentumMin < p) {
0046       particle.setAbsoluteMomentum(0.5 * p);
0047       const auto pid = particle.particleId().makeDescendant();
0048       generated.push_back(
0049           particle.withParticleId(pid).setAbsoluteMomentum(0.5 * p));
0050     }
0051     // never break
0052     return false;
0053   }
0054 };
0055
0056 // setup propagator-related types
0057 // use the default navigation
0058 using Navigator = Acts::Navigator;
0059 using MagneticField = Acts::ConstantBField;
0060 // propagate charged particles numerically in a constant magnetic field
0061 using ChargedStepper = Acts::EigenStepper<>;
0062 using ChargedPropagator = Acts::Propagator<ChargedStepper, Navigator>;
0063 // propagate neutral particles with just straight lines
0064 using NeutralStepper = Acts::StraightLineStepper;
0065 using NeutralPropagator = Acts::Propagator<NeutralStepper, Navigator>;
0066
0067 // setup simulator-related types
0068 // the random number generator type
0069 using Generator = std::ranlux48;
0070 // all charged particles w/ a mock-up physics list and hits everywhere
0071 using ChargedSelector = ActsFatras::EveryParticle;
0072 using ChargedInteractions =
0073     ActsFatras::InteractionList<ActsFatras::detail::StandardScattering,
0074                                 SplitEnergyLoss>;
0075 using ChargedSimulation =
0076     ActsFatras::SingleParticleSimulation<ChargedPropagator, ChargedInteractions,
0077                                          ActsFatras::EverySurface,
0078                                          ActsFatras::NoDecay>;
0079 // all neutral particles w/o physics and no hits
0080 using NeutralSelector = ActsFatras::EveryParticle;
0081 using NeutralInteractions = ActsFatras::InteractionList<>;
0082 using NeutralSimulation =
0083     ActsFatras::SingleParticleSimulation<NeutralPropagator, NeutralInteractions,
0084                                          ActsFatras::NoSurface,
0085                                          ActsFatras::NoDecay>;
0086 // full simulator type for charged and neutrals
0087 using Simulation = ActsFatras::Simulation<ChargedSelector, ChargedSimulation,
0088                                           NeutralSelector, NeutralSimulation>;
0089
0090 // parameters for data-driven test cases
0091
0092 const auto rangePdg =
0093     boost::unit_test::data::make(std::vector<Acts::PdgParticle>{
0094         Acts::PdgParticle::eElectron,
0095         Acts::PdgParticle::eMuon,
0096         Acts::PdgParticle::ePionPlus,
0097         Acts::PdgParticle::ePionZero,
0098     });
0099 const auto rangePhi = boost::unit_test::data::make(std::vector<double>{
0100     -135_degree,
0101     -45_degree,
0102     45_degree,
0103     135_degree,
0104 });
0105 const auto rangeEta = boost::unit_test::data::make(std::vector<double>{
0106     -1.0,
0107     0.0,
0108     1.0,
0109     3.0,
0110 });
0111 const auto rangeP = boost::unit_test::data::make(std::vector<double>{
0112     1_GeV,
0113     10_GeV,
0114 });
0115 const auto rangeNumParticles = boost::unit_test::data::make(std::vector<int>{
0116     1,
0117     2,
0118 });
0119
0120 const auto dataset =
0121     rangePdg * rangePhi * rangeEta * rangeP * rangeNumParticles;
0122
0123 // helper functions for tests
0124 template <typename Container>
0125 void sortByParticleId(Container& container) {
0126   std::ranges::sort(container, std::less{},
0127                     [](const auto& c) { return c.particleId(); });
0128 }
0129 template <typename Container>
0130 bool areParticleIdsUnique(const Container& sortedByParticleId) {
0131   // assumes the container is sorted by particle id
0132   auto ret =
0133       std::adjacent_find(sortedByParticleId.begin(), sortedByParticleId.end(),
0134                          [](const auto& lhs, const auto& rhs) {
0135                            return lhs.particleId() == rhs.particleId();
0136                          });
0137   return ret == sortedByParticleId.end();
0138 }
0139 template <typename Container, typename Value>
0140 bool containsParticleId(const Container& sortedByParticleId,
0141                         const Value& value) {
0142   return std::binary_search(sortedByParticleId.begin(),
0143                             sortedByParticleId.end(), value,
0144                             [](const auto& lhs, const auto& rhs) {
0145                               return lhs.particleId() < rhs.particleId();
0146                             });
0147 }
0148
0149 }  // namespace
0150
0151 BOOST_DATA_TEST_CASE(FatrasSimulation, dataset, pdg, phi, eta, p,
0152                      numParticles) {
0153   Acts::GeometryContext geoCtx;
0154   Acts::MagneticFieldContext magCtx;
0155   Acts::Logging::Level logLevel = Acts::Logging::Level::DEBUG;
0156
0157   // construct the example detector
0158   Acts::Test::CylindricalTrackingGeometry geoBuilder(geoCtx);
0159   auto trackingGeometry = geoBuilder();
0160
0161   // construct the propagators
0162   Navigator navigator({trackingGeometry});
0163   ChargedStepper chargedStepper(
0164       std::make_shared<Acts::ConstantBField>(Acts::Vector3{0, 0, 1_T}));
0165   ChargedPropagator chargedPropagator(std::move(chargedStepper), navigator);
0166   NeutralPropagator neutralPropagator(NeutralStepper(), navigator);
0167
0168   // construct the simulator
0169   ChargedSimulation simulatorCharged(
0170       std::move(chargedPropagator),
0171       Acts::getDefaultLogger("ChargedSimulation", logLevel));
0172   NeutralSimulation simulatorNeutral(
0173       std::move(neutralPropagator),
0174       Acts::getDefaultLogger("NeutralSimulation", logLevel));
0175   Simulation simulator(std::move(simulatorCharged),
0176                        std::move(simulatorNeutral));
0177
0178   // prepare simulation call parameters
0179   // random number generator
0180   Generator generator;
0181   // input/ output particle and hits containers
0182   std::vector<ActsFatras::Particle> input;
0183   std::vector<ActsFatras::Particle> simulatedInitial;
0184   std::vector<ActsFatras::Particle> simulatedFinal;
0185   std::vector<ActsFatras::Hit> hits;
0186
0187   // create input particles. particle number should ne non-zero.
0188   for (auto i = numParticles; 0 < i; --i) {
0189     const auto pid = ActsFatras::Barcode().setVertexPrimary(42).setParticle(i);
0190     const auto particle =
0191         ActsFatras::Particle(pid, pdg)
0192             .setDirection(Acts::makeDirectionFromPhiEta(phi, eta))
0193             .setAbsoluteMomentum(p);
0194     input.push_back(std::move(particle));
0195   }
0196   BOOST_TEST_INFO(input.front());
0197   BOOST_CHECK_EQUAL(input.size(), numParticles);
0198
0199   // run the simulation
0200   auto result = simulator.simulate(geoCtx, magCtx, generator, input,
0201                                    simulatedInitial, simulatedFinal, hits);
0202
0203   // should always succeed
0204   BOOST_CHECK(result.ok());
0205
0206   // ensure simulated particle containers have consistent content
0207   BOOST_CHECK_EQUAL(simulatedInitial.size(), simulatedFinal.size());
0208   for (std::size_t i = 0; i < simulatedInitial.size(); ++i) {
0209     const auto& initialParticle = simulatedInitial[i];
0210     const auto& finalParticle = simulatedFinal[i];
0211     // particle identify should not change during simulation
0212     BOOST_CHECK_EQUAL(initialParticle.particleId(), finalParticle.particleId());
0213     BOOST_CHECK_EQUAL(initialParticle.process(), finalParticle.process());
0214     BOOST_CHECK_EQUAL(initialParticle.pdg(), finalParticle.pdg());
0215     BOOST_CHECK_EQUAL(initialParticle.charge(), finalParticle.charge());
0216     BOOST_CHECK_EQUAL(initialParticle.mass(), finalParticle.mass());
0217   }
0218
0219   // we have no particle cuts and should not lose any particles.
0220   // might end up with more due to secondaries
0221   BOOST_CHECK_LE(input.size(), simulatedInitial.size());
0222   BOOST_CHECK_LE(input.size(), simulatedFinal.size());
0223   // there should be some hits if we started with a charged particle
0224   if (Acts::findCharge(pdg) != 0) {
0225     BOOST_CHECK_LT(0u, hits.size());
0226   }
0227
0228   // sort all outputs by particle id to simply further tests
0229   sortByParticleId(input);
0230   sortByParticleId(simulatedInitial);
0231   sortByParticleId(simulatedFinal);
0232   sortByParticleId(hits);
0233
0234   // check that all particle ids are unique
0235   BOOST_CHECK(areParticleIdsUnique(input));
0236   BOOST_CHECK(areParticleIdsUnique(simulatedInitial));
0237   BOOST_CHECK(areParticleIdsUnique(simulatedFinal));
0238   // hits must necessarily contain particle id duplicates
0239   // check that every input particles is simulated
0240   for (const auto& particle : input) {
0241     BOOST_CHECK(containsParticleId(simulatedInitial, particle));
0242     BOOST_CHECK(containsParticleId(simulatedFinal, particle));
0243   }
0244   // check that all hits can be associated to a particle
0245   for (const auto& hit : hits) {
0246     BOOST_CHECK(containsParticleId(simulatedInitial, hit));
0247     BOOST_CHECK(containsParticleId(simulatedFinal, hit));
0248   }
0249 }