Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-02-22 09:55:27

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