Fatras/src/FatrasSimulation.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 "ActsExamples/Fatras/FatrasSimulation.hpp"
0010
0011 #include "Acts/Geometry/GeometryContext.hpp"
0012 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0013 #include "Acts/Propagator/Navigator.hpp"
0014 #include "Acts/Propagator/Propagator.hpp"
0015 #include "Acts/Propagator/StraightLineStepper.hpp"
0016 #include "Acts/Propagator/SympyStepper.hpp"
0017 #include "Acts/Surfaces/Surface.hpp"
0018 #include "Acts/Utilities/Logger.hpp"
0019 #include "Acts/Utilities/Result.hpp"
0020 #include "ActsExamples/EventData/SimHit.hpp"
0021 #include "ActsExamples/EventData/SimParticle.hpp"
0022 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0023 #include "ActsExamples/Framework/IAlgorithm.hpp"
0024 #include "ActsExamples/Framework/RandomNumbers.hpp"
0025 #include "ActsFatras/EventData/Hit.hpp"
0026 #include "ActsFatras/EventData/Particle.hpp"
0027 #include "ActsFatras/Kernel/InteractionList.hpp"
0028 #include "ActsFatras/Kernel/Simulation.hpp"
0029 #include "ActsFatras/Physics/Decay/NoDecay.hpp"
0030 #include "ActsFatras/Physics/ElectroMagnetic/PhotonConversion.hpp"
0031 #include "ActsFatras/Physics/StandardInteractions.hpp"
0032 #include "ActsFatras/Selectors/SelectorHelpers.hpp"
0033 #include "ActsFatras/Selectors/SurfaceSelectors.hpp"
0034
0035 #include <algorithm>
0036 #include <ostream>
0037 #include <stdexcept>
0038 #include <system_error>
0039 #include <utility>
0040 #include <vector>
0041
0042 namespace {
0043
0044 /// Simple struct to select surfaces where hits should be generated.
0045 struct HitSurfaceSelector {
0046   bool sensitive = false;
0047   bool material = false;
0048   bool passive = false;
0049
0050   /// Check if the surface should be used.
0051   bool operator()(const Acts::Surface &surface) const {
0052     // sensitive/material are not mutually exclusive
0053     bool isSensitive = surface.associatedDetectorElement() != nullptr;
0054     bool isMaterial = surface.surfaceMaterial() != nullptr;
0055     // passive should be an orthogonal category
0056     bool isPassive = !(isSensitive || isMaterial);
0057     return (isSensitive && sensitive) || (isMaterial && material) ||
0058            (isPassive && passive);
0059   }
0060 };
0061
0062 }  // namespace
0063
0064 // Same interface as `ActsFatras::Simulation` but with concrete types.
0065 struct ActsExamples::detail::FatrasSimulation {
0066   virtual ~FatrasSimulation() = default;
0067   virtual Acts::Result<std::vector<ActsFatras::FailedParticle>> simulate(
0068       const Acts::GeometryContext &, const Acts::MagneticFieldContext &,
0069       ActsExamples::RandomEngine &, const std::vector<ActsFatras::Particle> &,
0070       std::vector<ActsFatras::Particle> &, std::vector<ActsFatras::Particle> &,
0071       std::vector<ActsFatras::Hit> &) const = 0;
0072 };
0073
0074 namespace {
0075
0076 // Magnetic-field specific PIMPL implementation.
0077 //
0078 // This always uses the SympyStepper for charged particle propagation and is
0079 // thus limited to propagation in vacuum at the moment.
0080 struct FatrasSimulationT final : ActsExamples::detail::FatrasSimulation {
0081   using CutPMin = ActsFatras::Min<ActsFatras::Casts::P>;
0082
0083   // typedefs for charge particle simulation
0084   // propagate charged particles numerically in the given magnetic field
0085   using ChargedStepper = Acts::SympyStepper;
0086   using ChargedPropagator = Acts::Propagator<ChargedStepper, Acts::Navigator>;
0087   // charged particles w/ standard em physics list and selectable hits
0088   using ChargedSelector = CutPMin;
0089   using ChargedSimulation = ActsFatras::SingleParticleSimulation<
0090       ChargedPropagator, ActsFatras::StandardChargedElectroMagneticInteractions,
0091       HitSurfaceSelector, ActsFatras::NoDecay>;
0092
0093   // typedefs for neutral particle simulation
0094   // propagate neutral particles with just straight lines
0095   using NeutralStepper = Acts::StraightLineStepper;
0096   using NeutralPropagator = Acts::Propagator<NeutralStepper, Acts::Navigator>;
0097   // neutral particles w/ photon conversion and no hits
0098   using NeutralSelector = CutPMin;
0099   using NeutralInteractions =
0100       ActsFatras::InteractionList<ActsFatras::PhotonConversion>;
0101   using NeutralSimulation = ActsFatras::SingleParticleSimulation<
0102       NeutralPropagator, NeutralInteractions, ActsFatras::NoSurface,
0103       ActsFatras::NoDecay>;
0104
0105   // combined simulation type
0106   using Simulation = ActsFatras::Simulation<ChargedSelector, ChargedSimulation,
0107                                             NeutralSelector, NeutralSimulation>;
0108
0109   Simulation simulation;
0110
0111   FatrasSimulationT(const ActsExamples::FatrasSimulation::Config &cfg,
0112                     Acts::Logging::Level lvl)
0113       : simulation(
0114             ChargedSimulation(
0115                 ChargedPropagator(
0116                     ChargedStepper(cfg.magneticField),
0117                     Acts::Navigator({cfg.trackingGeometry},
0118                                     Acts::getDefaultLogger("SimNav", lvl)),
0119                     Acts::getDefaultLogger("SimProp", lvl)),
0120                 Acts::getDefaultLogger("Simulation", lvl)),
0121             NeutralSimulation(
0122                 NeutralPropagator(
0123                     NeutralStepper(),
0124                     Acts::Navigator({cfg.trackingGeometry},
0125                                     Acts::getDefaultLogger("SimNav", lvl)),
0126                     Acts::getDefaultLogger("SimProp", lvl)),
0127                 Acts::getDefaultLogger("Simulation", lvl))) {
0128     using namespace ActsFatras;
0129     using namespace ActsFatras::detail;
0130     // apply the configuration
0131
0132     // minimal p cut on input particles and as is-alive check for interactions
0133     simulation.selectCharged.valMin = cfg.pMin;
0134     simulation.selectNeutral.valMin = cfg.pMin;
0135     simulation.charged.interactions =
0136         makeStandardChargedElectroMagneticInteractions(cfg.pMin);
0137
0138     // processes are enabled by default
0139     if (!cfg.emScattering) {
0140       simulation.charged.interactions.template disable<StandardScattering>();
0141     }
0142     if (!cfg.emEnergyLossIonisation) {
0143       simulation.charged.interactions.template disable<StandardBetheBloch>();
0144     }
0145     if (!cfg.emEnergyLossRadiation) {
0146       simulation.charged.interactions.template disable<StandardBetheHeitler>();
0147     }
0148     if (!cfg.emPhotonConversion) {
0149       simulation.neutral.interactions.template disable<PhotonConversion>();
0150     }
0151
0152     // configure hit surfaces for charged particles
0153     simulation.charged.selectHitSurface.sensitive = cfg.generateHitsOnSensitive;
0154     simulation.charged.selectHitSurface.material = cfg.generateHitsOnMaterial;
0155     simulation.charged.selectHitSurface.passive = cfg.generateHitsOnPassive;
0156
0157     simulation.charged.maxStepSize = cfg.maxStepSize;
0158     simulation.charged.pathLimit = cfg.pathLimit;
0159     simulation.neutral.maxStepSize = cfg.maxStepSize;
0160     simulation.neutral.pathLimit = cfg.pathLimit;
0161   }
0162   ~FatrasSimulationT() final = default;
0163
0164   Acts::Result<std::vector<ActsFatras::FailedParticle>> simulate(
0165       const Acts::GeometryContext &geoCtx,
0166       const Acts::MagneticFieldContext &magCtx, ActsExamples::RandomEngine &rng,
0167       const std::vector<ActsFatras::Particle> &inputParticles,
0168       std::vector<ActsFatras::Particle> &simulatedParticlesInitial,
0169       std::vector<ActsFatras::Particle> &simulatedParticlesFinal,
0170       std::vector<ActsFatras::Hit> &simHits) const final {
0171     return simulation.simulate(geoCtx, magCtx, rng, inputParticles,
0172                                simulatedParticlesInitial,
0173                                simulatedParticlesFinal, simHits);
0174   }
0175 };
0176
0177 }  // namespace
0178
0179 ActsExamples::FatrasSimulation::FatrasSimulation(Config cfg,
0180                                                  Acts::Logging::Level lvl)
0181     : ActsExamples::IAlgorithm("FatrasSimulation", lvl), m_cfg(std::move(cfg)) {
0182   ACTS_DEBUG("hits on sensitive surfaces: " << m_cfg.generateHitsOnSensitive);
0183   ACTS_DEBUG("hits on material surfaces: " << m_cfg.generateHitsOnMaterial);
0184   ACTS_DEBUG("hits on passive surfaces: " << m_cfg.generateHitsOnPassive);
0185
0186   if (!m_cfg.generateHitsOnSensitive && !m_cfg.generateHitsOnMaterial &&
0187       !m_cfg.generateHitsOnPassive) {
0188     ACTS_WARNING("FatrasSimulation not configured to generate any hits!");
0189   }
0190
0191   if (!m_cfg.trackingGeometry) {
0192     throw std::invalid_argument{"Missing tracking geometry"};
0193   }
0194   if (!m_cfg.magneticField) {
0195     throw std::invalid_argument{"Missing magnetic field"};
0196   }
0197   if (!m_cfg.randomNumbers) {
0198     throw std::invalid_argument("Missing random numbers tool");
0199   }
0200
0201   // construct the simulation for the specific magnetic field
0202   m_sim = std::make_unique<FatrasSimulationT>(m_cfg, lvl);
0203
0204   m_inputParticles.initialize(m_cfg.inputParticles);
0205   m_outputParticles.initialize(m_cfg.outputParticles);
0206   m_outputSimHits.initialize(m_cfg.outputSimHits);
0207 }
0208
0209 // explicit destructor needed for the PIMPL implementation to work
0210 ActsExamples::FatrasSimulation::~FatrasSimulation() = default;
0211
0212 ActsExamples::ProcessCode ActsExamples::FatrasSimulation::execute(
0213     const AlgorithmContext &ctx) const {
0214   // read input containers
0215   const auto &inputParticles = m_inputParticles(ctx);
0216
0217   ACTS_DEBUG(inputParticles.size() << " input particles");
0218
0219   // prepare input container
0220   std::vector<ActsFatras::Particle> particlesInput;
0221   particlesInput.reserve(inputParticles.size());
0222   for (const auto &p : inputParticles) {
0223     particlesInput.push_back(p.initialState());
0224   }
0225
0226   // prepare output containers
0227   std::vector<ActsFatras::Particle> particlesInitialUnordered;
0228   std::vector<ActsFatras::Particle> particlesFinalUnordered;
0229   std::vector<ActsFatras::Hit> simHitsUnordered;
0230   // reserve appropriate resources
0231   particlesInitialUnordered.reserve(inputParticles.size());
0232   particlesFinalUnordered.reserve(inputParticles.size());
0233   simHitsUnordered.reserve(inputParticles.size() *
0234                            m_cfg.averageHitsPerParticle);
0235
0236   // run the simulation w/ a local random generator
0237   auto rng = m_cfg.randomNumbers->spawnGenerator(ctx);
0238   auto ret = m_sim->simulate(ctx.geoContext, ctx.magFieldContext, rng,
0239                              particlesInput, particlesInitialUnordered,
0240                              particlesFinalUnordered, simHitsUnordered);
0241   // fatal error leads to panic
0242   if (!ret.ok()) {
0243     ACTS_FATAL("event " << ctx.eventNumber << " simulation failed with error "
0244                         << ret.error().message());
0245     return ProcessCode::ABORT;
0246   }
0247   // failed particles are just logged. assumes that failed particles are due
0248   // to edge-cases representing a tiny fraction of the event; not due to a
0249   // fundamental issue.
0250   for (const auto &failed : ret.value()) {
0251     ACTS_ERROR("event " << ctx.eventNumber << " particle " << failed.particle
0252                         << " failed to simulate with error " << failed.error
0253                         << ": " << failed.error.message());
0254   }
0255
0256   ACTS_DEBUG(particlesInitialUnordered.size()
0257              << " simulated particles (initial state)");
0258   ACTS_DEBUG(particlesFinalUnordered.size()
0259              << " simulated particles (final state)");
0260   ACTS_DEBUG(simHitsUnordered.size() << " simulated hits");
0261
0262   if (particlesInitialUnordered.size() != particlesFinalUnordered.size()) {
0263     ACTS_ERROR("number of initial and final state particles differ");
0264   }
0265
0266   // order output containers
0267   SimParticleStateContainer particlesInitial(particlesInitialUnordered.begin(),
0268                                              particlesInitialUnordered.end());
0269   SimParticleStateContainer particlesFinal(particlesFinalUnordered.begin(),
0270                                            particlesFinalUnordered.end());
0271   SimHitContainer simHits(simHitsUnordered.begin(), simHitsUnordered.end());
0272
0273   SimParticleContainer particlesSimulated;
0274   particlesSimulated.reserve(particlesInitial.size());
0275   for (const auto &particleInitial : particlesInitial) {
0276     SimParticle particleSimulated(particleInitial, particleInitial);
0277
0278     if (auto it = particlesFinal.find(particleInitial.particleId());
0279         it != particlesFinal.end()) {
0280       particleSimulated.finalState() = *it;
0281     } else {
0282       ACTS_ERROR("particle " << particleInitial.particleId()
0283                              << " has no final state");
0284     }
0285
0286     particlesSimulated.insert(particleSimulated);
0287   }
0288
0289   // store ordered output containers
0290   m_outputParticles(ctx, std::move(particlesSimulated));
0291   m_outputSimHits(ctx, std::move(simHits));
0292
0293   return ActsExamples::ProcessCode::SUCCESS;
0294 }