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