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