Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:12:13

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 #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   { p.generatePathLimits(rng, prt) } -> std::same_as<std::pair<double, double>>;
0087 };
0088 
0089 template <typename process_t>
0090 concept ContinuousProcessConcept = !PointLikeProcessConcept<process_t>;
0091 
0092 template <typename process_t>
0093 struct PointLikeProcessTrait {
0094   static constexpr bool value = PointLikeProcessConcept<process_t>;
0095 };
0096 
0097 template <typename process_t>
0098 struct ContinuousProcessTrait {
0099   static constexpr bool value = ContinuousProcessConcept<process_t>;
0100 };
0101 
0102 template <typename processes_t>
0103 using ContinuousIndices = TupleFilter<ContinuousProcessTrait, processes_t>;
0104 template <typename processes_t>
0105 using PointLikeIndices = TupleFilter<PointLikeProcessTrait, processes_t>;
0106 
0107 }  // namespace detail
0108 
0109 /// Compile-time set of interaction processes for the simulation.
0110 ///
0111 /// Two different type of interaction processes are supported: continuous and
0112 /// point-like interactions.
0113 ///
0114 /// Continuous processes scale with the passed material. They tpyically
0115 /// describe effective results of a large number of small interactions such as
0116 /// multiple scattering or ionisation. Continuous process types **must** provide
0117 /// a call operator with the following signature:
0118 ///
0119 ///     template <typename generator_t>
0120 ///     bool
0121 ///     operator()(
0122 ///         generator_t& rng,
0123 ///         const Acts::MaterialSlab& slab,
0124 ///         Particle& particle,
0125 ///         std::vector<Particle>& generatedParticles) const
0126 ///
0127 /// If multiple continuous processes are defined, they are executed serially in
0128 /// the order in which they are given.
0129 ///
0130 /// For point-like processes, the passed material only affects the probability
0131 /// with which they occur but not the interaction itself, e.g. photon conversion
0132 /// into electron pairs. They are simulated by first drawing a limit on the
0133 /// material paths and then executing the interaction with the shortest limit
0134 /// when the drawn amount of material has been passed. Point-like process
0135 /// types **must** provide the following two member functions:
0136 ///
0137 ///     // generate X0/L0 limits
0138 ///     template <typename generator_t>
0139 ///     std::pair<Scalar, Scalar>
0140 ///     generatePathLimits(
0141 ///         generator& rng,
0142 ///         const Particle& particle) const
0143 ///
0144 ///     // run the process simulation
0145 ///     template <typename generator_t>
0146 ///     bool
0147 ///     run(
0148 ///         generator_t& rng,
0149 ///         Particle& particle,
0150 ///         std::vector<Particle>& generatedParticles) const
0151 ///
0152 /// For both continuous and point-like interactions, the output particle is
0153 /// modified in-place (if needed) and the return value indicates a break
0154 /// condition in the simulation, i.e. the particle is dead (true) or alive
0155 /// (false) after the interaction.
0156 ///
0157 /// @note If an interaction destroys the incoming particle, the process
0158 ///   simulation should indicate this via the break condition only and not
0159 ///   by reducing the particle momentum to zero. The incoming particle should
0160 ///   retain its initial kinematic state; the final kinematic state before
0161 ///   destruction is typically of more interest to the user and this simplifies
0162 ///   validation.
0163 ///
0164 /// The physics processes are extendable by the user to accommodate their
0165 /// specific requirements. While the set of available physics processes must be
0166 /// configured at compile-time, within that set, processes can again be
0167 /// selectively disabled at run-time. By default all processes are applied.
0168 template <typename... processes_t>
0169 class InteractionList {
0170   using Mask = std::bitset<sizeof...(processes_t)>;
0171   using Processes = std::tuple<processes_t...>;
0172   using ContinuousIndices = detail::ContinuousIndices<Processes>;
0173   using PointLikeIndices = detail::PointLikeIndices<Processes>;
0174 
0175  public:
0176   /// Point-like interaction selection.
0177   struct Selection {
0178     double x0Limit = std::numeric_limits<double>::infinity();
0179     double l0Limit = std::numeric_limits<double>::infinity();
0180     std::size_t x0Process = std::numeric_limits<std::size_t>::max();
0181     std::size_t l0Process = std::numeric_limits<std::size_t>::max();
0182   };
0183 
0184   /// Disable a specific process identified by index.
0185   void disable(std::size_t process) { m_mask.set(process); }
0186   /// Disable a specific process identified by type.
0187   ///
0188   /// @note Disables only the first element, if multiple elements of the same
0189   ///   type exist.
0190   template <typename process_t>
0191   void disable() {
0192     m_mask.set(detail::TupleIndexOf<process_t, Processes>::value);
0193   }
0194 
0195   /// Access a specific process identified by index.
0196   template <std::size_t kProcess>
0197   std::tuple_element_t<kProcess, Processes>& get() {
0198     return std::get<kProcess>(m_processes);
0199   }
0200   /// Access a specific process identified by type.
0201   ///
0202   /// @warning This function only works if all configured processes have
0203   ///   different types.
0204   template <typename process_t>
0205   process_t& get() {
0206     return std::get<process_t>(m_processes);
0207   }
0208 
0209   /// Simulate the combined effects from all continuous interactions.
0210   ///
0211   /// @tparam generator_t must be a RandomNumberEngine
0212   /// @param[in]     rng       is the random number generator
0213   /// @param[in]     slab      is the passed material
0214   /// @param[in,out] particle  is the particle being updated
0215   /// @param[out]    generated is the container of generated particles
0216   /// @return Break condition, i.e. whether a process stopped the propagation
0217   template <typename generator_t>
0218   bool runContinuous(generator_t& rng, const Acts::MaterialSlab& slab,
0219                      Particle& particle,
0220                      std::vector<Particle>& generated) const {
0221     return runContinuousImpl(rng, slab, particle, generated,
0222                              ContinuousIndices());
0223   }
0224 
0225   /// Arm the point-like interactions by generating limits and select processes.
0226   ///
0227   /// @tparam generator_t must be a RandomNumberEngine
0228   /// @param[in] rng      is the random number generator
0229   /// @param[in] particle is the initial particle state
0230   /// @return X0/L0 limits for the particle and the process index that should be
0231   ///   executed once the limit has been reached.
0232   template <typename generator_t>
0233   Selection armPointLike(generator_t& rng, const Particle& particle) const {
0234     Selection selection;
0235     armPointLikeImpl(rng, particle, selection, PointLikeIndices());
0236     return selection;
0237   }
0238 
0239   /// Simulate the effects from a single point-like interaction.
0240   ///
0241   /// @tparam generator_t must be a RandomNumberEngine
0242   /// @param[in]     rng          is the random number generator
0243   /// @param[in]     processIndex is the index of the process to be executed
0244   /// @param[in,out] particle     is the particle being updated
0245   /// @param[out]    generated    is the container of generated particles
0246   /// @return Break condition, i.e. whether a process killed the particle
0247   ///
0248   /// The process index is expected to originate from a previous
0249   /// `armPointLike(...)` call, but this is not enforced. How to select the
0250   /// correct process requires more information that is not available here.
0251   template <typename generator_t>
0252   bool runPointLike(generator_t& rng, std::size_t processIndex,
0253                     Particle& particle,
0254                     std::vector<Particle>& generated) const {
0255     return runPointLikeImpl(rng, processIndex, particle, generated,
0256                             PointLikeIndices());
0257   }
0258 
0259  private:
0260   // allow processes to be masked. defaults to zeros -> no masked processes
0261   Mask m_mask;
0262   Processes m_processes;
0263 
0264   // for the `runContinuous` call, we need to iterate over all available
0265   // processes and apply the ones that implement the continuous process
0266   // interface. this is done using an index-based compile-time recursive call.
0267   template <typename generator_t, std::size_t kI0, std::size_t... kIs>
0268   bool runContinuousImpl(generator_t& rng, const Acts::MaterialSlab& slab,
0269                          Particle& particle, std::vector<Particle>& generated,
0270                          std::index_sequence<kI0, kIs...> /*indices*/) const {
0271     const auto& process = std::get<kI0>(m_processes);
0272     // only call process if it is not masked
0273     if (!m_mask[kI0] && process(rng, slab, particle, generated)) {
0274       // exit early in case the process signals an abort
0275       return true;
0276     }
0277     return runContinuousImpl(rng, slab, particle, generated,
0278                              std::index_sequence<kIs...>());
0279   }
0280   template <typename generator_t>
0281   bool runContinuousImpl(generator_t& /*rng*/,
0282                          const Acts::MaterialSlab& /*slab*/,
0283                          Particle& /*particle*/,
0284                          std::vector<Particle>& /*generated*/,
0285                          std::index_sequence<> /*indices*/) const {
0286     return false;
0287   }
0288 
0289   // for the `armPointLike` call, we need to iterate over all available
0290   // processes and select the ones that generate the smallest limits. this is
0291   // done using an index-based compile-time recursive call.
0292   template <typename generator_t, std::size_t kI0, std::size_t... kIs>
0293   void armPointLikeImpl(generator_t& rng, const Particle& particle,
0294                         Selection& selection,
0295                         std::index_sequence<kI0, kIs...> /*indices*/) const {
0296     // only arm the process if it is not masked
0297     if (!m_mask[kI0]) {
0298       auto [x0Limit, l0Limit] =
0299           std::get<kI0>(m_processes).generatePathLimits(rng, particle);
0300       if (x0Limit < selection.x0Limit) {
0301         selection.x0Limit = x0Limit;
0302         selection.x0Process = kI0;
0303       }
0304       if (l0Limit < selection.l0Limit) {
0305         selection.l0Limit = l0Limit;
0306         selection.l0Process = kI0;
0307       }
0308     }
0309     // continue with the remaining processes
0310     armPointLikeImpl(rng, particle, selection, std::index_sequence<kIs...>());
0311   }
0312   template <typename generator_t>
0313   void armPointLikeImpl(generator_t& /*rng*/, const Particle& /*particle*/,
0314                         Selection& /*selection*/,
0315                         std::index_sequence<> /*indices*/) const {}
0316 
0317   // for the `runPointLike` call we need to call just one process. since we can
0318   // not select a tuple element with a runtime index, we need to iterate over
0319   // all processes with a compile-time recursive function until we reach the
0320   // requested one.
0321   template <typename generator_t, std::size_t kI0, std::size_t... kIs>
0322   bool runPointLikeImpl(generator_t& rng, std::size_t processIndex,
0323                         Particle& particle, std::vector<Particle>& generated,
0324                         std::index_sequence<kI0, kIs...> /*indices*/) const {
0325     if (kI0 == processIndex) {
0326       if (m_mask[kI0]) {
0327         // the selected process is masked. since nothing is executed the
0328         // particle continues to be alive; not a break condition.
0329         return false;
0330       }
0331       return std::get<kI0>(m_processes).run(rng, particle, generated);
0332     }
0333     // continue the iteration with the remaining processes
0334     return runPointLikeImpl(rng, processIndex, particle, generated,
0335                             std::index_sequence<kIs...>());
0336   }
0337   template <typename generator_t>
0338   bool runPointLikeImpl(generator_t& /*rng*/, std::size_t /*processIndex*/,
0339                         Particle& /*particle*/,
0340                         std::vector<Particle>& /*generated*/,
0341                         std::index_sequence<> /*indices*/) const {
0342     // the requested process index is outside the possible range. **do not**
0343     // treat this as an error to simplify the case of an empty physics lists or
0344     // a default process index.
0345     return false;
0346   }
0347 };
0348 
0349 }  // namespace ActsFatras