Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-05 08:26:41

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2018-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 #pragma once
0010 
0011 #include "Acts/Material/MaterialSlab.hpp"
0012 #include "ActsFatras/EventData/Particle.hpp"
0013 
0014 #include <bitset>
0015 #include <tuple>
0016 #include <type_traits>
0017 #include <utility>
0018 
0019 namespace ActsFatras {
0020 namespace detail {
0021 
0022 /// Retrieve index of the first matching type in the tuple.
0023 ///
0024 /// Taken from https://stackoverflow.com/a/18063608.
0025 template <class T, class Tuple>
0026 struct TupleIndexOf;
0027 template <class T, class... Types>
0028 struct TupleIndexOf<T, std::tuple<T, Types...>> {
0029   static constexpr std::size_t value = 0u;
0030 };
0031 template <class T, class U, class... Types>
0032 struct TupleIndexOf<T, std::tuple<U, Types...>> {
0033   static constexpr std::size_t value =
0034       1u + TupleIndexOf<T, std::tuple<Types...>>::value;
0035 };
0036 
0037 // Construct an index sequence for a subset of the tuple elements.
0038 //
0039 // Whether an element is part of the subset is defined by the predicate
0040 // template type. It must take the element type as its only template parameter
0041 // and must provide a static `value` member value. If the value evaluates to
0042 // `true`, then the corresponding index will be part of the index sequence.
0043 //
0044 // Example: The tuple contains four elements, where all but the third one (i=2)
0045 // should be selected. This leads to the following recursive expansion
0046 // where the index sequence of the subset is filled from the front.
0047 //
0048 //        TupleFilterImpl<..., kCounter=4>          // select index=3
0049 //     -> TupleFilterImpl<..., kCounter=3, 3>       // skip   index=2
0050 //     -> TupleFilterImpl<..., kCounter=2, 3>       // select index=1
0051 //     -> TupleFilterImpl<..., kCounter=1, 1, 3>    // select index=0
0052 //     -> TupleFilterImpl<..., kCounter=0, 0, 1, 3> // terminate
0053 //
0054 template <template <typename> typename predicate_t, typename tuple_t,
0055           std::size_t kCounter, std::size_t... kIndices>
0056 struct TupleFilterImpl {
0057   static constexpr auto kIndex = kCounter - 1u;
0058   static constexpr bool kElementSelection =
0059       predicate_t<std::tuple_element_t<kIndex, tuple_t>>::value;
0060   // recursive type if the element would be selected
0061   using SelectElement = typename TupleFilterImpl<predicate_t, tuple_t, kIndex,
0062                                                  kIndex, kIndices...>::Type;
0063   // recursive type if the element would be skipped
0064   using SkipElement =
0065       typename TupleFilterImpl<predicate_t, tuple_t, kIndex, kIndices...>::Type;
0066   // select recursive type based on the selector decision
0067   using Type =
0068       std::conditional_t<kElementSelection, SelectElement, SkipElement>;
0069 };
0070 template <template <typename> typename predicate_t, typename tuple_t,
0071           std::size_t... kIndices>
0072 struct TupleFilterImpl<predicate_t, tuple_t, 0u, kIndices...> {
0073   using Type = std::index_sequence<kIndices...>;
0074 };
0075 template <template <typename> typename predicate_t, typename tuple_t>
0076 using TupleFilter = typename TupleFilterImpl<predicate_t, tuple_t,
0077                                              std::tuple_size_v<tuple_t>>::Type;
0078 
0079 /// Check if the given type is a point-like process.
0080 ///
0081 /// Only checks for the existence of the templated `generatePathLimits` method
0082 template <typename process_t>
0083 concept PointLikeProcessConcept = requires(
0084     const process_t& p, std::uniform_int_distribution<unsigned int>& rng,
0085     const Particle& prt) {
0086   {
0087     p.generatePathLimits(rng, prt)
0088   } -> std::same_as<std::pair<Particle::Scalar, Particle::Scalar>>;
0089 };
0090 
0091 template <typename process_t>
0092 concept ContinuousProcessConcept = !PointLikeProcessConcept<process_t>;
0093 
0094 template <typename process_t>
0095 struct PointLikeProcessTrait {
0096   static constexpr bool value = PointLikeProcessConcept<process_t>;
0097 };
0098 
0099 template <typename process_t>
0100 struct ContinuousProcessTrait {
0101   static constexpr bool value = ContinuousProcessConcept<process_t>;
0102 };
0103 
0104 template <typename processes_t>
0105 using ContinuousIndices = TupleFilter<ContinuousProcessTrait, processes_t>;
0106 template <typename processes_t>
0107 using PointLikeIndices = TupleFilter<PointLikeProcessTrait, processes_t>;
0108 
0109 }  // namespace detail
0110 
0111 /// Compile-time set of interaction processes for the simulation.
0112 ///
0113 /// Two different type of interaction processes are supported: continuous and
0114 /// point-like interactions.
0115 ///
0116 /// Continuous processes scale with the passed material. They tpyically
0117 /// describe effective results of a large number of small interactions such as
0118 /// multiple scattering or ionisation. Continuous process types **must** provide
0119 /// a call operator with the following signature:
0120 ///
0121 ///     template <typename generator_t>
0122 ///     bool
0123 ///     operator()(
0124 ///         generator_t& rng,
0125 ///         const Acts::MaterialSlab& slab,
0126 ///         Particle& particle,
0127 ///         std::vector<Particle>& generatedParticles) const
0128 ///
0129 /// If multiple continuous processes are defined, they are executed serially in
0130 /// the order in which they are given.
0131 ///
0132 /// For point-like processes, the passed material only affects the probability
0133 /// with which they occur but not the interaction itself, e.g. photon conversion
0134 /// into electron pairs. They are simulated by first drawing a limit on the
0135 /// material paths and then executing the interaction with the shortest limit
0136 /// when the drawn amount of material has been passed. Point-like process
0137 /// types **must** provide the following two member functions:
0138 ///
0139 ///     // generate X0/L0 limits
0140 ///     template <typename generator_t>
0141 ///     std::pair<Scalar, Scalar>
0142 ///     generatePathLimits(
0143 ///         generator& rng,
0144 ///         const Particle& particle) const
0145 ///
0146 ///     // run the process simulation
0147 ///     template <typename generator_t>
0148 ///     bool
0149 ///     run(
0150 ///         generator_t& rng,
0151 ///         Particle& particle,
0152 ///         std::vector<Particle>& generatedParticles) const
0153 ///
0154 /// For both continuous and point-like interactions, the output particle is
0155 /// modified in-place (if needed) and the return value indicates a break
0156 /// condition in the simulation, i.e. the particle is dead (true) or alive
0157 /// (false) after the interaction.
0158 ///
0159 /// @note If an interaction destroys the incoming particle, the process
0160 ///   simulation should indicate this via the break condition only and not
0161 ///   by reducing the particle momentum to zero. The incoming particle should
0162 ///   retain its initial kinematic state; the final kinematic state before
0163 ///   destruction is typically of more interest to the user and this simplifies
0164 ///   validation.
0165 ///
0166 /// The physics processes are extendable by the user to accommodate their
0167 /// specific requirements. While the set of available physics processes must be
0168 /// configured at compile-time, within that set, processes can again be
0169 /// selectively disabled at run-time. By default all processes are applied.
0170 template <typename... processes_t>
0171 class InteractionList {
0172   using Mask = std::bitset<sizeof...(processes_t)>;
0173   using Processes = std::tuple<processes_t...>;
0174   using ContinuousIndices = detail::ContinuousIndices<Processes>;
0175   using PointLikeIndices = detail::PointLikeIndices<Processes>;
0176 
0177  public:
0178   /// Point-like interaction selection.
0179   struct Selection {
0180     Particle::Scalar x0Limit =
0181         std::numeric_limits<Particle::Scalar>::infinity();
0182     Particle::Scalar l0Limit =
0183         std::numeric_limits<Particle::Scalar>::infinity();
0184     std::size_t x0Process = std::numeric_limits<std::size_t>::max();
0185     std::size_t l0Process = std::numeric_limits<std::size_t>::max();
0186   };
0187 
0188   /// Disable a specific process identified by index.
0189   void disable(std::size_t process) { m_mask.set(process); }
0190   /// Disable a specific process identified by type.
0191   ///
0192   /// @note Disables only the first element, if multiple elements of the same
0193   ///   type exist.
0194   template <typename process_t>
0195   void disable() {
0196     m_mask.set(detail::TupleIndexOf<process_t, Processes>::value);
0197   }
0198 
0199   /// Access a specific process identified by index.
0200   template <std::size_t kProcess>
0201   std::tuple_element_t<kProcess, Processes>& get() {
0202     return std::get<kProcess>(m_processes);
0203   }
0204   /// Access a specific process identified by type.
0205   ///
0206   /// @warning This function only works if all configured processes have
0207   ///   different types.
0208   template <typename process_t>
0209   process_t& get() {
0210     return std::get<process_t>(m_processes);
0211   }
0212 
0213   /// Simulate the combined effects from all continuous interactions.
0214   ///
0215   /// @tparam generator_t must be a RandomNumberEngine
0216   /// @param[in]     rng       is the random number generator
0217   /// @param[in]     slab      is the passed material
0218   /// @param[in,out] particle  is the particle being updated
0219   /// @param[out]    generated is the container of generated particles
0220   /// @return Break condition, i.e. whether a process stopped the propagation
0221   template <typename generator_t>
0222   bool runContinuous(generator_t& rng, const Acts::MaterialSlab& slab,
0223                      Particle& particle,
0224                      std::vector<Particle>& generated) const {
0225     return runContinuousImpl(rng, slab, particle, generated,
0226                              ContinuousIndices());
0227   }
0228 
0229   /// Arm the point-like interactions by generating limits and select processes.
0230   ///
0231   /// @tparam generator_t must be a RandomNumberEngine
0232   /// @param[in] rng      is the random number generator
0233   /// @param[in] particle is the initial particle state
0234   /// @return X0/L0 limits for the particle and the process index that should be
0235   ///   executed once the limit has been reached.
0236   template <typename generator_t>
0237   Selection armPointLike(generator_t& rng, const Particle& particle) const {
0238     Selection selection;
0239     armPointLikeImpl(rng, particle, selection, PointLikeIndices());
0240     return selection;
0241   }
0242 
0243   /// Simulate the effects from a single point-like interaction.
0244   ///
0245   /// @tparam generator_t must be a RandomNumberEngine
0246   /// @param[in]     rng          is the random number generator
0247   /// @param[in]     processIndex is the index of the process to be executed
0248   /// @param[in,out] particle     is the particle being updated
0249   /// @param[out]    generated    is the container of generated particles
0250   /// @return Break condition, i.e. whether a process killed the particle
0251   ///
0252   /// The process index is expected to originate from a previous
0253   /// `armPointLike(...)` call, but this is not enforced. How to select the
0254   /// correct process requires more information that is not available here.
0255   template <typename generator_t>
0256   bool runPointLike(generator_t& rng, std::size_t processIndex,
0257                     Particle& particle,
0258                     std::vector<Particle>& generated) const {
0259     return runPointLikeImpl(rng, processIndex, particle, generated,
0260                             PointLikeIndices());
0261   }
0262 
0263  private:
0264   // allow processes to be masked. defaults to zeros -> no masked processes
0265   Mask m_mask;
0266   Processes m_processes;
0267 
0268   // for the `runContinuous` call, we need to iterate over all available
0269   // processes and apply the ones that implement the continuous process
0270   // interface. this is done using an index-based compile-time recursive call.
0271   template <typename generator_t, std::size_t kI0, std::size_t... kIs>
0272   bool runContinuousImpl(generator_t& rng, const Acts::MaterialSlab& slab,
0273                          Particle& particle, std::vector<Particle>& generated,
0274                          std::index_sequence<kI0, kIs...> /*indices*/) const {
0275     const auto& process = std::get<kI0>(m_processes);
0276     // only call process if it is not masked
0277     if (!m_mask[kI0] && process(rng, slab, particle, generated)) {
0278       // exit early in case the process signals an abort
0279       return true;
0280     }
0281     return runContinuousImpl(rng, slab, particle, generated,
0282                              std::index_sequence<kIs...>());
0283   }
0284   template <typename generator_t>
0285   bool runContinuousImpl(generator_t& /*rng*/,
0286                          const Acts::MaterialSlab& /*slab*/,
0287                          Particle& /*particle*/,
0288                          std::vector<Particle>& /*generated*/,
0289                          std::index_sequence<> /*indices*/) const {
0290     return false;
0291   }
0292 
0293   // for the `armPointLike` call, we need to iterate over all available
0294   // processes and select the ones that generate the smallest limits. this is
0295   // done using an index-based compile-time recursive call.
0296   template <typename generator_t, std::size_t kI0, std::size_t... kIs>
0297   void armPointLikeImpl(generator_t& rng, const Particle& particle,
0298                         Selection& selection,
0299                         std::index_sequence<kI0, kIs...> /*indices*/) const {
0300     // only arm the process if it is not masked
0301     if (!m_mask[kI0]) {
0302       auto [x0Limit, l0Limit] =
0303           std::get<kI0>(m_processes).generatePathLimits(rng, particle);
0304       if (x0Limit < selection.x0Limit) {
0305         selection.x0Limit = x0Limit;
0306         selection.x0Process = kI0;
0307       }
0308       if (l0Limit < selection.l0Limit) {
0309         selection.l0Limit = l0Limit;
0310         selection.l0Process = kI0;
0311       }
0312     }
0313     // continue with the remaining processes
0314     armPointLikeImpl(rng, particle, selection, std::index_sequence<kIs...>());
0315   }
0316   template <typename generator_t>
0317   void armPointLikeImpl(generator_t& /*rng*/, const Particle& /*particle*/,
0318                         Selection& /*selection*/,
0319                         std::index_sequence<> /*indices*/) const {}
0320 
0321   // for the `runPointLike` call we need to call just one process. since we can
0322   // not select a tuple element with a runtime index, we need to iterate over
0323   // all processes with a compile-time recursive function until we reach the
0324   // requested one.
0325   template <typename generator_t, std::size_t kI0, std::size_t... kIs>
0326   bool runPointLikeImpl(generator_t& rng, std::size_t processIndex,
0327                         Particle& particle, std::vector<Particle>& generated,
0328                         std::index_sequence<kI0, kIs...> /*indices*/) const {
0329     if (kI0 == processIndex) {
0330       if (m_mask[kI0]) {
0331         // the selected process is masked. since nothing is executed the
0332         // particle continues to be alive; not a break condition.
0333         return false;
0334       }
0335       return std::get<kI0>(m_processes).run(rng, particle, generated);
0336     }
0337     // continue the iteration with the remaining processes
0338     return runPointLikeImpl(rng, processIndex, particle, generated,
0339                             std::index_sequence<kIs...>());
0340   }
0341   template <typename generator_t>
0342   bool runPointLikeImpl(generator_t& /*rng*/, std::size_t /*processIndex*/,
0343                         Particle& /*particle*/,
0344                         std::vector<Particle>& /*generated*/,
0345                         std::index_sequence<> /*indices*/) const {
0346     // the requested process index is outside the possible range. **do not**
0347     // treat this as an error to simplify the case of an empty physics lists or
0348     // a default process index.
0349     return false;
0350   }
0351 };
0352 
0353 }  // namespace ActsFatras