Simulator.hpp

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

// This file is part of the Acts project.
//
// Copyright (C) 2018-2020 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 <algorithm>
#include <cassert>
#include <iterator>
#include <memory>
#include <vector>

#include "Acts/EventData/ChargePolicy.hpp"
#include "Acts/EventData/SingleCurvilinearTrackParameters.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/Propagator/AbortList.hpp"
#include "Acts/Propagator/ActionList.hpp"
#include "Acts/Propagator/DebugOutputActor.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/Interactor.hpp"
#include "ActsFatras/Kernel/detail/SimulatorError.hpp"

namespace ActsFatras {

template <typename propagator_t, typename physics_list_t,
          typename hit_surface_selector_t>
struct ParticleSimulator {
  propagator_t propagator;
  physics_list_t physics;
  hit_surface_selector_t selectHitSurface;
  std::shared_ptr<const Acts::Logger> localLogger = nullptr;

  ParticleSimulator(propagator_t &&propagator_, Acts::Logging::Level lvl)
      : propagator(propagator_),
        localLogger(Acts::getDefaultLogger("Simulator", lvl)) {}

  const Acts::Logger &logger() const { return *localLogger; }

  template <typename generator_t>
  Acts::Result<InteractorResult> simulate(
      const Acts::GeometryContext &geoCtx,
      const Acts::MagneticFieldContext &magCtx, generator_t &generator,
      const Particle &particle) const {
    assert(localLogger and "Missing local logger");

    // propagator-related additional types
    using Interact =
        Interactor<generator_t, physics_list_t, hit_surface_selector_t>;
    using Actions = Acts::ActionList<Interact, Acts::DebugOutputActor>;
    using Abort = Acts::AbortList<typename Interact::ParticleNotAlive,
                                  Acts::EndOfWorldReached>;
    using PropagatorOptions = Acts::PropagatorOptions<Actions, Abort>;

    // Construct per-call options.
    PropagatorOptions options(geoCtx, magCtx);
    options.absPdgCode = particle.pdg();
    options.mass = particle.mass();
    options.debug = localLogger->doPrint(Acts::Logging::Level::DEBUG);
    // setup the interactor as part of the propagator options
    auto &interactor = options.actionList.template get<Interact>();
    interactor.generator = &generator;
    interactor.physics = physics;
    interactor.selectHitSurface = selectHitSurface;
    interactor.particle = particle;

    // run with a start parameter type depending on the particle charge.
    // TODO make track parameters consistently constructible regardless
    //      of the charge policy and template the class on the parameter.
    if (particle.charge() != 0) {
      Acts::SingleCurvilinearTrackParameters<Acts::ChargedPolicy> start(
          std::nullopt, particle.position(),
          particle.absMomentum() * particle.unitDirection(), particle.charge(),
          particle.time());
      auto result = propagator.propagate(start, options);
      if (result.ok()) {
        return result.value().template get<InteractorResult>();
      } else {
        return result.error();
      }
    } else {
      Acts::SingleCurvilinearTrackParameters<Acts::NeutralPolicy> start(
          std::nullopt, particle.position(),
          particle.absMomentum() * particle.unitDirection(), particle.time());
      auto result = propagator.propagate(start, options);
      if (result.ok()) {
        return result.value().template get<InteractorResult>();
      } else {
        return result.error();
      }
    }
  }
};

template <typename charged_selector_t, typename charged_simulator_t,
          typename neutral_selector_t, typename neutral_simulator_t>
struct Simulator {
  struct FailedParticle {
    Particle particle;
    std::error_code error;
  };

  charged_selector_t selectCharged;
  neutral_selector_t selectNeutral;
  charged_simulator_t charged;
  neutral_simulator_t neutral;

  Simulator(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 ParticleSimulatorResult = Acts::Result<InteractorResult>;

    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::SimulatorError::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.
        // they must therefore be either charged or neutral.
        ParticleSimulatorResult result = ParticleSimulatorResult::success({});
        if (selectCharged(initialParticle)) {
          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 {
    const bool isValidCharged = selectCharged(particle);
    const bool isValidNeutral = selectNeutral(particle);
    assert(not(isValidCharged and isValidNeutral) and
           "Charge selectors are not mutually exclusive");
    return isValidCharged xor isValidNeutral;
  }

  template <typename particles_t, typename hits_t>
  void copyOutputs(const InteractorResult &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