IntegrationTests/Fatras/FatrasSimulationTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020-2021 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 http://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/EventData/TrackParameters.hpp"
0014 #include "Acts/MagneticField/ConstantBField.hpp"
0015 #include "Acts/Propagator/EigenStepper.hpp"
0016 #include "Acts/Propagator/Navigator.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::sort(container.begin(), container.end(),
0127             [](const auto& lhs, const auto& rhs) {
0128               return lhs.particleId() < rhs.particleId();
0129             });
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(FatrasSimulation, dataset, pdg, phi, eta, p,
0154                      numParticles) {
0155   using namespace Acts::UnitLiterals;
0156
0157   Acts::GeometryContext geoCtx;
0158   Acts::MagneticFieldContext magCtx;
0159   Acts::Logging::Level logLevel = Acts::Logging::Level::DEBUG;
0160
0161   // construct the example detector
0162   Acts::Test::CylindricalTrackingGeometry geoBuilder(geoCtx);
0163   auto trackingGeometry = geoBuilder();
0164
0165   // construct the propagators
0166   Navigator navigator({trackingGeometry});
0167   ChargedStepper chargedStepper(
0168       std::make_shared<Acts::ConstantBField>(Acts::Vector3{0, 0, 1_T}));
0169   ChargedPropagator chargedPropagator(std::move(chargedStepper), navigator);
0170   NeutralPropagator neutralPropagator(NeutralStepper(), navigator);
0171
0172   // construct the simulator
0173   ChargedSimulation simulatorCharged(
0174       std::move(chargedPropagator),
0175       Acts::getDefaultLogger("ChargedSimulation", logLevel));
0176   NeutralSimulation simulatorNeutral(
0177       std::move(neutralPropagator),
0178       Acts::getDefaultLogger("NeutralSimulation", logLevel));
0179   Simulation simulator(std::move(simulatorCharged),
0180                        std::move(simulatorNeutral));
0181
0182   // prepare simulation call parameters
0183   // random number generator
0184   Generator generator;
0185   // input/ output particle and hits containers
0186   std::vector<ActsFatras::Particle> input;
0187   std::vector<ActsFatras::Particle> simulatedInitial;
0188   std::vector<ActsFatras::Particle> simulatedFinal;
0189   std::vector<ActsFatras::Hit> hits;
0190
0191   // create input particles. particle number should ne non-zero.
0192   for (auto i = numParticles; 0 < i; --i) {
0193     const auto pid = ActsFatras::Barcode().setVertexPrimary(42).setParticle(i);
0194     const auto particle =
0195         ActsFatras::Particle(pid, pdg)
0196             .setDirection(Acts::makeDirectionFromPhiEta(phi, eta))
0197             .setAbsoluteMomentum(p);
0198     input.push_back(std::move(particle));
0199   }
0200   BOOST_TEST_INFO(input.front());
0201   BOOST_CHECK_EQUAL(input.size(), numParticles);
0202
0203   // run the simulation
0204   auto result = simulator.simulate(geoCtx, magCtx, generator, input,
0205                                    simulatedInitial, simulatedFinal, hits);
0206
0207   // should always succeed
0208   BOOST_CHECK(result.ok());
0209
0210   // ensure simulated particle containers have consistent content
0211   BOOST_CHECK_EQUAL(simulatedInitial.size(), simulatedFinal.size());
0212   for (std::size_t i = 0; i < simulatedInitial.size(); ++i) {
0213     const auto& initialParticle = simulatedInitial[i];
0214     const auto& finalParticle = simulatedFinal[i];
0215     // particle identify should not change during simulation
0216     BOOST_CHECK_EQUAL(initialParticle.particleId(), finalParticle.particleId());
0217     BOOST_CHECK_EQUAL(initialParticle.process(), finalParticle.process());
0218     BOOST_CHECK_EQUAL(initialParticle.pdg(), finalParticle.pdg());
0219     BOOST_CHECK_EQUAL(initialParticle.charge(), finalParticle.charge());
0220     BOOST_CHECK_EQUAL(initialParticle.mass(), finalParticle.mass());
0221   }
0222
0223   // we have no particle cuts and should not lose any particles.
0224   // might end up with more due to secondaries
0225   BOOST_CHECK_LE(input.size(), simulatedInitial.size());
0226   BOOST_CHECK_LE(input.size(), simulatedFinal.size());
0227   // there should be some hits if we started with a charged particle
0228   if (Acts::findCharge(pdg) != 0) {
0229     BOOST_CHECK_LT(0u, hits.size());
0230   }
0231
0232   // sort all outputs by particle id to simply further tests
0233   sortByParticleId(input);
0234   sortByParticleId(simulatedInitial);
0235   sortByParticleId(simulatedFinal);
0236   sortByParticleId(hits);
0237
0238   // check that all particle ids are unique
0239   BOOST_CHECK(areParticleIdsUnique(input));
0240   BOOST_CHECK(areParticleIdsUnique(simulatedInitial));
0241   BOOST_CHECK(areParticleIdsUnique(simulatedFinal));
0242   // hits must necessarily contain particle id duplicates
0243   // check that every input particles is simulated
0244   for (const auto& particle : input) {
0245     BOOST_CHECK(containsParticleId(simulatedInitial, particle));
0246     BOOST_CHECK(containsParticleId(simulatedFinal, particle));
0247   }
0248   // check that all hits can be associated to a particle
0249   for (const auto& hit : hits) {
0250     BOOST_CHECK(containsParticleId(simulatedInitial, hit));
0251     BOOST_CHECK(containsParticleId(simulatedFinal, hit));
0252   }
0253 }