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