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