Simulation.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Simulation.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2018-2021 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include "Acts/EventData/Charge.hpp"
#include "Acts/EventData/GenericCurvilinearTrackParameters.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/Propagator/AbortList.hpp"
#include "Acts/Propagator/ActionList.hpp"
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/Propagator/StandardAborters.hpp"
#include "Acts/Utilities/Logger.hpp"
#include "Acts/Utilities/Result.hpp"
#include "ActsFatras/EventData/Hit.hpp"
#include "ActsFatras/EventData/Particle.hpp"
#include "ActsFatras/Kernel/SimulationResult.hpp"
#include "ActsFatras/Kernel/detail/SimulationActor.hpp"
#include "ActsFatras/Kernel/detail/SimulationError.hpp"

#include <algorithm>
#include <cassert>
#include <iterator>
#include <memory>
#include <vector>

namespace ActsFatras {

template <typename propagator_t, typename interactions_t,
          typename hit_surface_selector_t, typename decay_t>
struct SingleParticleSimulation {
  propagator_t propagator;
  decay_t decay;
  interactions_t interactions;
  hit_surface_selector_t selectHitSurface;
  std::unique_ptr<const Acts::Logger> logger;

  SingleParticleSimulation(propagator_t &&propagator_,
                           std::unique_ptr<const Acts::Logger> _logger)
      : propagator(propagator_), logger(std::move(_logger)) {}

  template <typename generator_t>
  Acts::Result<SimulationResult> simulate(
      const Acts::GeometryContext &geoCtx,
      const Acts::MagneticFieldContext &magCtx, generator_t &generator,
      const Particle &particle) const {
    // propagator-related additional types
    using Actor = detail::SimulationActor<generator_t, decay_t, interactions_t,
                                          hit_surface_selector_t>;
    using Aborter = typename Actor::ParticleNotAlive;
    using Result = typename Actor::result_type;
    using Actions = Acts::ActionList<Actor>;
    using Abort = Acts::AbortList<Aborter, Acts::EndOfWorldReached>;
    using PropagatorOptions = Acts::PropagatorOptions<Actions, Abort>;

    // Construct per-call options.
    PropagatorOptions options(geoCtx, magCtx);
    // setup the interactor as part of the propagator options
    auto &actor = options.actionList.template get<Actor>();
    actor.generator = &generator;
    actor.decay = decay;
    actor.interactions = interactions;
    actor.selectHitSurface = selectHitSurface;
    actor.initialParticle = particle;

    if (particle.hasReferenceSurface()) {
      auto result = propagator.propagate(
          particle.boundParameters(geoCtx).value(), options);
      if (not result.ok()) {
        return result.error();
      }
      auto &value = result.value().template get<Result>();
      return std::move(value);
    }

    auto result =
        propagator.propagate(particle.curvilinearParameters(), options);
    if (not result.ok()) {
      return result.error();
    }
    auto &value = result.value().template get<Result>();
    return std::move(value);
  }
};

struct FailedParticle {
  Particle particle;
  std::error_code error;
};

template <typename charged_selector_t, typename charged_simulator_t,
          typename neutral_selector_t, typename neutral_simulator_t>
struct Simulation {
  charged_selector_t selectCharged;
  neutral_selector_t selectNeutral;
  charged_simulator_t charged;
  neutral_simulator_t neutral;

  Simulation(charged_simulator_t &&charged_, neutral_simulator_t &&neutral_)
      : charged(std::move(charged_)), neutral(std::move(neutral_)) {}

  template <typename generator_t, typename input_particles_t,
            typename output_particles_t, typename hits_t>
  Acts::Result<std::vector<FailedParticle>> simulate(
      const Acts::GeometryContext &geoCtx,
      const Acts::MagneticFieldContext &magCtx, generator_t &generator,
      const input_particles_t &inputParticles,
      output_particles_t &simulatedParticlesInitial,
      output_particles_t &simulatedParticlesFinal, hits_t &hits) const {
    assert(
        (simulatedParticlesInitial.size() == simulatedParticlesFinal.size()) and
        "Inconsistent initial sizes of the simulated particle containers");

    using SingleParticleSimulationResult = Acts::Result<SimulationResult>;

    std::vector<FailedParticle> failedParticles;

    for (const Particle &inputParticle : inputParticles) {
      // only consider simulatable particles
      if (not selectParticle(inputParticle)) {
        continue;
      }
      // required to allow correct particle id numbering for secondaries later
      if ((inputParticle.particleId().generation() != 0u) or
          (inputParticle.particleId().subParticle() != 0u)) {
        return detail::SimulationError::eInvalidInputParticleId;
      }

      // Do a *depth-first* simulation of the particle and its secondaries,
      // i.e. we simulate all secondaries, tertiaries, ... before simulating
      // the next primary particle. Use the end of the output container as
      // a queue to store particles that should be simulated.
      //
      // WARNING the initial particle state output container will be modified
      //         during iteration. New secondaries are added to and failed
      //         particles might be removed. To avoid issues, access must always
      //         occur via indices.
      auto iinitial = simulatedParticlesInitial.size();
      simulatedParticlesInitial.push_back(inputParticle);
      for (; iinitial < simulatedParticlesInitial.size(); ++iinitial) {
        const auto &initialParticle = simulatedParticlesInitial[iinitial];

        // only simulatable particles are pushed to the container and here we
        // only need to switch between charged/neutral.
        SingleParticleSimulationResult result =
            SingleParticleSimulationResult::success({});
        if (initialParticle.charge() != Particle::Scalar(0)) {
          result = charged.simulate(geoCtx, magCtx, generator, initialParticle);
        } else {
          result = neutral.simulate(geoCtx, magCtx, generator, initialParticle);
        }

        if (not result.ok()) {
          // remove particle from output container since it was not simulated.
          simulatedParticlesInitial.erase(
              std::next(simulatedParticlesInitial.begin(), iinitial));
          // record the particle as failed
          failedParticles.push_back({initialParticle, result.error()});
          continue;
        }

        copyOutputs(result.value(), simulatedParticlesInitial,
                    simulatedParticlesFinal, hits);
        // since physics processes are independent, there can be particle id
        // collisions within the generated secondaries. they can be resolved by
        // renumbering within each sub-particle generation. this must happen
        // before the particle is simulated since the particle id is used to
        // associate generated hits back to the particle.
        renumberTailParticleIds(simulatedParticlesInitial, iinitial);
      }
    }

    // the overall function call succeeded, i.e. no fatal errors occured.
    // yet, there might have been some particle for which the propagation
    // failed. thus, the successful result contains a list of failed particles.
    // sounds a bit weird, but that is the way it is.
    return failedParticles;
  }

 private:
  bool selectParticle(const Particle &particle) const {
    if (particle.charge() != Particle::Scalar(0)) {
      return selectCharged(particle);
    } else {
      return selectNeutral(particle);
    }
  }

  template <typename particles_t, typename hits_t>
  void copyOutputs(const SimulationResult &result,
                   particles_t &particlesInitial, particles_t &particlesFinal,
                   hits_t &hits) const {
    // initial particle state was already pushed to the container before
    // store final particle state at the end of the simulation
    particlesFinal.push_back(result.particle);
    // move generated secondaries that should be simulated to the output
    std::copy_if(
        result.generatedParticles.begin(), result.generatedParticles.end(),
        std::back_inserter(particlesInitial),
        [this](const Particle &particle) { return selectParticle(particle); });
    std::copy(result.hits.begin(), result.hits.end(), std::back_inserter(hits));
  }

  template <typename particles_t>
  static void renumberTailParticleIds(particles_t &particles,
                                      std::size_t lastValid) {
    // iterate over adjacent pairs; potentially modify the second element.
    // assume e.g. a primary particle 2 with generation=subparticle=0 that
    // generates two secondaries during simulation. we have the following
    // ids (decoded as vertex|particle|generation|subparticle)
    //
    //     v|2|0|0, v|2|1|0, v|2|1|0
    //
    // where v represents the vertex numbers. this will be renumbered to
    //
    //     v|2|0|0, v|2|1|0, v|2|1|1
    //
    // if each secondary then generates two tertiaries we could have e.g.
    //
    //     v|2|0|0, v|2|1|0, v|2|1|1, v|2|2|0, v|2|2|1, v|2|2|0, v|2|2|1
    //
    // which would then be renumbered to
    //
    //     v|2|0|0, v|2|1|0, v|2|1|1, v|2|2|0, v|2|2|1, v|2|2|2, v|2|2|3
    //
    for (auto j = lastValid; (j + 1u) < particles.size(); ++j) {
      const auto prevId = particles[j].particleId();
      auto currId = particles[j + 1u].particleId();
      // NOTE primary/secondary vertex and particle are assumed to be equal
      // only renumber within one generation
      if (prevId.generation() != currId.generation()) {
        continue;
      }
      // ensure the sub-particle is strictly monotonic within a generation
      if (prevId.subParticle() < currId.subParticle()) {
        continue;
      }
      // sub-particle numbering must be non-zero
      currId.setSubParticle(prevId.subParticle() + 1u);
      particles[j + 1u] = particles[j + 1u].withParticleId(currId);
    }
  }
};

}  // namespace ActsFatras