Fatras/src/FatrasSimulation.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017-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 "ActsExamples/Fatras/FatrasSimulation.hpp"
0010
0011 #include "Acts/Definitions/Direction.hpp"
0012 #include "Acts/EventData/TrackParameters.hpp"
0013 #include "Acts/Geometry/GeometryContext.hpp"
0014 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0015 #include "Acts/Propagator/EigenStepper.hpp"
0016 #include "Acts/Propagator/Navigator.hpp"
0017 #include "Acts/Propagator/Propagator.hpp"
0018 #include "Acts/Propagator/StraightLineStepper.hpp"
0019 #include "Acts/Surfaces/Surface.hpp"
0020 #include "Acts/Utilities/Logger.hpp"
0021 #include "Acts/Utilities/Result.hpp"
0022 #include "ActsExamples/EventData/SimHit.hpp"
0023 #include "ActsExamples/EventData/SimParticle.hpp"
0024 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0025 #include "ActsExamples/Framework/IAlgorithm.hpp"
0026 #include "ActsExamples/Framework/RandomNumbers.hpp"
0027 #include "ActsFatras/EventData/Barcode.hpp"
0028 #include "ActsFatras/EventData/Particle.hpp"
0029 #include "ActsFatras/Kernel/InteractionList.hpp"
0030 #include "ActsFatras/Kernel/Simulation.hpp"
0031 #include "ActsFatras/Physics/Decay/NoDecay.hpp"
0032 #include "ActsFatras/Physics/ElectroMagnetic/PhotonConversion.hpp"
0033 #include "ActsFatras/Physics/NuclearInteraction/NuclearInteractionParameters.hpp"
0034 #include "ActsFatras/Physics/StandardInteractions.hpp"
0035 #include "ActsFatras/Selectors/SelectorHelpers.hpp"
0036 #include "ActsFatras/Selectors/SurfaceSelectors.hpp"
0037
0038 #include <algorithm>
0039 #include <array>
0040 #include <map>
0041 #include <ostream>
0042 #include <stdexcept>
0043 #include <system_error>
0044 #include <utility>
0045 #include <vector>
0046
0047 #include <boost/version.hpp>
0048
0049 namespace {
0050
0051 /// Simple struct to select surfaces where hits should be generated.
0052 struct HitSurfaceSelector {
0053   bool sensitive = false;
0054   bool material = false;
0055   bool passive = false;
0056
0057   /// Check if the surface should be used.
0058   bool operator()(const Acts::Surface &surface) const {
0059     // sensitive/material are not mutually exclusive
0060     bool isSensitive = surface.associatedDetectorElement() != nullptr;
0061     bool isMaterial = surface.surfaceMaterial() != nullptr;
0062     // passive should be an orthogonal category
0063     bool isPassive = !(isSensitive || isMaterial);
0064     return (isSensitive && sensitive) || (isMaterial && material) ||
0065            (isPassive && passive);
0066   }
0067 };
0068
0069 }  // namespace
0070
0071 // Same interface as `ActsFatras::Simulation` but with concrete types.
0072 struct ActsExamples::detail::FatrasSimulation {
0073   virtual ~FatrasSimulation() = default;
0074   virtual Acts::Result<std::vector<ActsFatras::FailedParticle>> simulate(
0075       const Acts::GeometryContext &, const Acts::MagneticFieldContext &,
0076       ActsExamples::RandomEngine &, const ActsExamples::SimParticleContainer &,
0077       ActsExamples::SimParticleContainer::sequence_type &,
0078       ActsExamples::SimParticleContainer::sequence_type &,
0079       ActsExamples::SimHitContainer::sequence_type &) const = 0;
0080 };
0081
0082 namespace {
0083
0084 // Magnetic-field specific PIMPL implementation.
0085 //
0086 // This always uses the EigenStepper with default extensions for charged
0087 // particle propagation and is thus limited to propagation in vacuum at the
0088 // moment.
0089 // @TODO: Remove this, unneeded after #675
0090 struct FatrasSimulationT final : ActsExamples::detail::FatrasSimulation {
0091   using CutPMin = ActsFatras::Min<ActsFatras::Casts::P>;
0092
0093   // typedefs for charge particle simulation
0094   // propagate charged particles numerically in the given magnetic field
0095   using ChargedStepper = Acts::EigenStepper<>;
0096   using ChargedPropagator = Acts::Propagator<ChargedStepper, Acts::Navigator>;
0097   // charged particles w/ standard em physics list and selectable hits
0098   using ChargedSelector = CutPMin;
0099   using ChargedSimulation = ActsFatras::SingleParticleSimulation<
0100       ChargedPropagator, ActsFatras::StandardChargedElectroMagneticInteractions,
0101       HitSurfaceSelector, ActsFatras::NoDecay>;
0102
0103   // typedefs for neutral particle simulation
0104   // propagate neutral particles with just straight lines
0105   using NeutralStepper = Acts::StraightLineStepper;
0106   using NeutralPropagator = Acts::Propagator<NeutralStepper, Acts::Navigator>;
0107   // neutral particles w/ photon conversion and no hits
0108   using NeutralSelector = CutPMin;
0109   using NeutralInteractions =
0110       ActsFatras::InteractionList<ActsFatras::PhotonConversion>;
0111   using NeutralSimulation = ActsFatras::SingleParticleSimulation<
0112       NeutralPropagator, NeutralInteractions, ActsFatras::NoSurface,
0113       ActsFatras::NoDecay>;
0114
0115   // combined simulation type
0116   using Simulation = ActsFatras::Simulation<ChargedSelector, ChargedSimulation,
0117                                             NeutralSelector, NeutralSimulation>;
0118
0119   Simulation simulation;
0120
0121   FatrasSimulationT(const ActsExamples::FatrasSimulation::Config &cfg,
0122                     Acts::Logging::Level lvl)
0123       : simulation(
0124             ChargedSimulation(
0125                 ChargedPropagator(ChargedStepper(cfg.magneticField),
0126                                   Acts::Navigator{{cfg.trackingGeometry}}),
0127                 Acts::getDefaultLogger("Simulation", lvl)),
0128             NeutralSimulation(
0129                 NeutralPropagator(NeutralStepper(),
0130                                   Acts::Navigator{{cfg.trackingGeometry}}),
0131                 Acts::getDefaultLogger("Simulation", lvl))) {
0132     using namespace ActsFatras;
0133     using namespace ActsFatras::detail;
0134     // apply the configuration
0135
0136     // minimal p cut on input particles and as is-alive check for interactions
0137     simulation.selectCharged.valMin = cfg.pMin;
0138     simulation.selectNeutral.valMin = cfg.pMin;
0139     simulation.charged.interactions =
0140         makeStandardChargedElectroMagneticInteractions(cfg.pMin);
0141
0142     // processes are enabled by default
0143     if (!cfg.emScattering) {
0144       simulation.charged.interactions.template disable<StandardScattering>();
0145     }
0146     if (!cfg.emEnergyLossIonisation) {
0147       simulation.charged.interactions.template disable<StandardBetheBloch>();
0148     }
0149     if (!cfg.emEnergyLossRadiation) {
0150       simulation.charged.interactions.template disable<StandardBetheHeitler>();
0151     }
0152     if (!cfg.emPhotonConversion) {
0153       simulation.neutral.interactions.template disable<PhotonConversion>();
0154     }
0155
0156     // configure hit surfaces for charged particles
0157     simulation.charged.selectHitSurface.sensitive = cfg.generateHitsOnSensitive;
0158     simulation.charged.selectHitSurface.material = cfg.generateHitsOnMaterial;
0159     simulation.charged.selectHitSurface.passive = cfg.generateHitsOnPassive;
0160
0161     simulation.charged.maxStepSize = cfg.maxStepSize;
0162     simulation.charged.pathLimit = cfg.pathLimit;
0163     simulation.neutral.maxStepSize = cfg.maxStepSize;
0164     simulation.neutral.pathLimit = cfg.pathLimit;
0165   }
0166   ~FatrasSimulationT() final = default;
0167
0168   Acts::Result<std::vector<ActsFatras::FailedParticle>> simulate(
0169       const Acts::GeometryContext &geoCtx,
0170       const Acts::MagneticFieldContext &magCtx, ActsExamples::RandomEngine &rng,
0171       const ActsExamples::SimParticleContainer &inputParticles,
0172       ActsExamples::SimParticleContainer::sequence_type
0173           &simulatedParticlesInitial,
0174       ActsExamples::SimParticleContainer::sequence_type
0175           &simulatedParticlesFinal,
0176       ActsExamples::SimHitContainer::sequence_type &simHits) const final {
0177     return simulation.simulate(geoCtx, magCtx, rng, inputParticles,
0178                                simulatedParticlesInitial,
0179                                simulatedParticlesFinal, simHits);
0180   }
0181 };
0182
0183 }  // namespace
0184
0185 ActsExamples::FatrasSimulation::FatrasSimulation(Config cfg,
0186                                                  Acts::Logging::Level lvl)
0187     : ActsExamples::IAlgorithm("FatrasSimulation", lvl), m_cfg(std::move(cfg)) {
0188   ACTS_DEBUG("hits on sensitive surfaces: " << m_cfg.generateHitsOnSensitive);
0189   ACTS_DEBUG("hits on material surfaces: " << m_cfg.generateHitsOnMaterial);
0190   ACTS_DEBUG("hits on passive surfaces: " << m_cfg.generateHitsOnPassive);
0191
0192   if (!m_cfg.generateHitsOnSensitive && !m_cfg.generateHitsOnMaterial &&
0193       !m_cfg.generateHitsOnPassive) {
0194     ACTS_WARNING("FatrasSimulation not configured to generate any hits!");
0195   }
0196
0197   if (!m_cfg.trackingGeometry) {
0198     throw std::invalid_argument{"Missing tracking geometry"};
0199   }
0200   if (!m_cfg.magneticField) {
0201     throw std::invalid_argument{"Missing magnetic field"};
0202   }
0203   if (!m_cfg.randomNumbers) {
0204     throw std::invalid_argument("Missing random numbers tool");
0205   }
0206
0207   // construct the simulation for the specific magnetic field
0208   m_sim = std::make_unique<FatrasSimulationT>(m_cfg, lvl);
0209
0210   m_inputParticles.initialize(m_cfg.inputParticles);
0211   m_outputParticlesInitial.initialize(m_cfg.outputParticlesInitial);
0212   m_outputParticlesFinal.initialize(m_cfg.outputParticlesFinal);
0213   m_outputSimHits.initialize(m_cfg.outputSimHits);
0214 }
0215
0216 // explicit destructor needed for the PIMPL implementation to work
0217 ActsExamples::FatrasSimulation::~FatrasSimulation() = default;
0218
0219 ActsExamples::ProcessCode ActsExamples::FatrasSimulation::execute(
0220     const AlgorithmContext &ctx) const {
0221   // read input containers
0222   const auto &inputParticles = m_inputParticles(ctx);
0223
0224   ACTS_DEBUG(inputParticles.size() << " input particles");
0225
0226   // prepare output containers
0227   SimParticleContainer::sequence_type particlesInitialUnordered;
0228   SimParticleContainer::sequence_type particlesFinalUnordered;
0229   SimHitContainer::sequence_type 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                              inputParticles, 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());
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   // order output containers
0263 #if BOOST_VERSION >= 107800
0264   SimParticleContainer particlesInitial(particlesInitialUnordered.begin(),
0265                                         particlesInitialUnordered.end());
0266   SimParticleContainer particlesFinal(particlesFinalUnordered.begin(),
0267                                       particlesFinalUnordered.end());
0268   SimHitContainer simHits(simHitsUnordered.begin(), simHitsUnordered.end());
0269 #else
0270   // working around a nasty boost bug
0271   // https://github.com/boostorg/container/issues/244
0272
0273   SimParticleContainer particlesInitial;
0274   SimParticleContainer particlesFinal;
0275   SimHitContainer simHits;
0276
0277   particlesInitial.reserve(particlesInitialUnordered.size());
0278   particlesFinal.reserve(particlesFinalUnordered.size());
0279   simHits.reserve(simHitsUnordered.size());
0280
0281   for (const auto &p : particlesInitialUnordered) {
0282     particlesInitial.insert(p);
0283   }
0284   for (const auto &p : particlesFinalUnordered) {
0285     particlesFinal.insert(p);
0286   }
0287   for (const auto &h : simHitsUnordered) {
0288     simHits.insert(h);
0289   }
0290 #endif
0291
0292   // store ordered output containers
0293   m_outputParticlesInitial(ctx, std::move(particlesInitial));
0294   m_outputParticlesFinal(ctx, std::move(particlesFinal));
0295   m_outputSimHits(ctx, std::move(simHits));
0296
0297   return ActsExamples::ProcessCode::SUCCESS;
0298 }