Acts/Propagator/StepperExtensionList.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 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/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Propagator/detail/Auctioneer.hpp"
0014 #include "Acts/Utilities/detail/Extendable.hpp"
0015 #include "Acts/Utilities/detail/MPL/all_of.hpp"
0016 #include "Acts/Utilities/detail/MPL/has_duplicates.hpp"
0017
0018 #include <array>
0019
0020 namespace Acts {
0021
0022 /// @brief Container of extensions used in the stepper of the propagation. This
0023 /// struct allows a broadcast of function calls for each element in the list.
0024 /// The broadcasts occur for a certain function at each step in a specific
0025 /// order.
0026 /// The first function is an evaluater if an extension is or how many extensions
0027 /// are applicable for an upcoming step.
0028 /// The next functions called are the evaluations of the k_1 - k_4 or the RKN4
0029 /// integration.
0030 /// The last function call in a step is the finalize() method. This method is an
0031 /// overloaded function (optionally propagates the covariance).
0032 /// Each method has the possibility to break the evaluation of a given step if
0033 /// an extension reports that something went wrong (e.g. a particle lost too
0034 /// much momentum during the step)
0035 /// @tparam extensions Types of the extensions
0036 template <typename... extensions>
0037 struct StepperExtensionList : private detail::Extendable<extensions...> {
0038  private:
0039   // Checkout for duplicates in the extensions
0040   static_assert(!detail::has_duplicates_v<extensions...>,
0041                 "same extension type specified several times");
0042
0043   static constexpr unsigned int nExtensions = sizeof...(extensions);
0044
0045   static_assert(nExtensions != 0, "no extension type specified");
0046
0047   // Access to all extensions
0048   using detail::Extendable<extensions...>::tuple;
0049
0050   // Vector of valid extensions for a step
0051   std::array<bool, nExtensions> validExtensions{};
0052
0053  public:
0054   // Access to an extension
0055   using detail::Extendable<extensions...>::get;
0056
0057   /// @brief Evaluation function to set valid extensions for an upcoming
0058   /// integration step
0059   ///
0060   /// @tparam propagator_state_t Type of the state of the propagator
0061   /// @tparam stepper_t Type of the stepper
0062   /// @tparam navigtor_t Type of the navigator
0063   ///
0064   /// @param [in] state State of the propagator
0065   /// @param [in] stepper Stepper of the propagation
0066   /// @param [in] navigator Navigator of the propagation
0067   template <typename propagator_state_t, typename stepper_t,
0068             typename navigtor_t>
0069   bool validExtensionForStep(const propagator_state_t& state,
0070                              const stepper_t& stepper,
0071                              const navigtor_t& navigator) {
0072     const auto bids = std::apply(
0073         [&](const auto&... ext) {
0074           return std::array<int, nExtensions>{
0075               ext.bid(state, stepper, navigator)...};
0076         },
0077         tuple());
0078
0079     validExtensions = state.stepping.auctioneer(std::move(bids));
0080
0081     return (std::find(validExtensions.begin(), validExtensions.end(), true) !=
0082             validExtensions.end());
0083   }
0084
0085   /// @brief This functions broadcasts the call for evaluating a generic k. It
0086   /// collects all arguments and extensions, test their validity for the
0087   /// evaluation and passes them forward for evaluation and returns a boolean as
0088   /// indicator if the evaluation is valid.
0089   template <typename propagator_state_t, typename stepper_t,
0090             typename navigator_t>
0091   bool k(const propagator_state_t& state, const stepper_t& stepper,
0092          const navigator_t& navigator, Vector3& knew, const Vector3& bField,
0093          std::array<double, 4>& kQoP, const int i, const double h = 0.,
0094          const Vector3& kprev = Vector3::Zero()) {
0095     // TODO replace with integer-templated lambda with C++20
0096     auto impl = [&, i, h](auto intType, auto& implRef) {
0097       constexpr int N = decltype(intType)::value;
0098
0099       if constexpr (N == 0) {
0100         return true;
0101       } else {
0102         // If element is invalid: continue
0103         if (!std::get<N - 1>(validExtensions)) {
0104           return implRef(std::integral_constant<int, N - 1>{}, implRef);
0105         }
0106         // Continue as long as evaluations are 'true'
0107         if (std::get<N - 1>(this->tuple())
0108                 .template k(state, stepper, navigator, knew, bField, kQoP, i, h,
0109                             kprev)) {
0110           return implRef(std::integral_constant<int, N - 1>{}, implRef);
0111         } else {
0112           // Break at false
0113           return false;
0114         }
0115       }
0116     };
0117
0118     return impl(std::integral_constant<int, nExtensions>{}, impl);
0119   }
0120
0121   /// @brief This functions broadcasts the call for evaluating k1. It collects
0122   /// all arguments and extensions, test their validity for the evaluation and
0123   /// passes them forward for evaluation and returns a boolean as indicator if
0124   /// the evaluation is valid.
0125   template <typename propagator_state_t, typename stepper_t,
0126             typename navigator_t>
0127   bool k1(const propagator_state_t& state, const stepper_t& stepper,
0128           const navigator_t& navigator, Vector3& knew, const Vector3& bField,
0129           std::array<double, 4>& kQoP) {
0130     return k(state, stepper, navigator, knew, bField, kQoP, 0);
0131   }
0132
0133   /// @brief This functions broadcasts the call for evaluating k2. It collects
0134   /// all arguments and extensions and passes them forward for evaluation and
0135   /// returns a boolean as indicator if the evaluation is valid.
0136   template <typename propagator_state_t, typename stepper_t,
0137             typename navigator_t>
0138   bool k2(const propagator_state_t& state, const stepper_t& stepper,
0139           const navigator_t& navigator, Vector3& knew, const Vector3& bField,
0140           std::array<double, 4>& kQoP, const double h, const Vector3& kprev) {
0141     return k(state, stepper, navigator, knew, bField, kQoP, 1, h, kprev);
0142   }
0143
0144   /// @brief This functions broadcasts the call for evaluating k3. It collects
0145   /// all arguments and extensions and passes them forward for evaluation and
0146   /// returns a boolean as indicator if the evaluation is valid.
0147   template <typename propagator_state_t, typename stepper_t,
0148             typename navigator_t>
0149   bool k3(const propagator_state_t& state, const stepper_t& stepper,
0150           const navigator_t& navigator, Vector3& knew, const Vector3& bField,
0151           std::array<double, 4>& kQoP, const double h, const Vector3& kprev) {
0152     return k(state, stepper, navigator, knew, bField, kQoP, 2, h, kprev);
0153   }
0154
0155   /// @brief This functions broadcasts the call for evaluating k4. It collects
0156   /// all arguments and extensions and passes them forward for evaluation and
0157   /// returns a boolean as indicator if the evaluation is valid.
0158   template <typename propagator_state_t, typename stepper_t,
0159             typename navigator_t>
0160   bool k4(const propagator_state_t& state, const stepper_t& stepper,
0161           const navigator_t& navigator, Vector3& knew, const Vector3& bField,
0162           std::array<double, 4>& kQoP, const double h, const Vector3& kprev) {
0163     return k(state, stepper, navigator, knew, bField, kQoP, 3, h, kprev);
0164   }
0165
0166   /// @brief This functions broadcasts the call of the method finalize(). It
0167   /// collects all extensions and arguments and passes them forward for
0168   /// evaluation and returns a boolean.
0169   template <typename propagator_state_t, typename stepper_t,
0170             typename navigator_t>
0171   bool finalize(propagator_state_t& state, const stepper_t& stepper,
0172                 const navigator_t& navigator, const double h, FreeMatrix& D) {
0173     // TODO replace with integer-templated lambda with C++20
0174     auto impl = [&, h](auto intType, auto& implRef) {
0175       constexpr int N = decltype(intType)::value;
0176
0177       if constexpr (N == 0) {
0178         return true;
0179       } else {
0180         // If element is invalid: continue
0181         if (!std::get<N - 1>(validExtensions)) {
0182           return implRef(std::integral_constant<int, N - 1>{}, implRef);
0183         }
0184         // Continue as long as evaluations are 'true'
0185         if (std::get<N - 1>(this->tuple())
0186                 .finalize(state, stepper, navigator, h, D)) {
0187           return implRef(std::integral_constant<int, N - 1>{}, implRef);
0188         } else {
0189           // Break at false
0190           return false;
0191         }
0192       }
0193     };
0194
0195     return impl(std::integral_constant<int, nExtensions>{}, impl);
0196   }
0197
0198   /// @brief This functions broadcasts the call of the method finalize(). It
0199   /// collects all extensions and arguments and passes them forward for
0200   /// evaluation and returns a boolean.
0201   template <typename propagator_state_t, typename stepper_t,
0202             typename navigator_t>
0203   bool finalize(propagator_state_t& state, const stepper_t& stepper,
0204                 const navigator_t& navigator, const double h) {
0205     // TODO replace with integer-templated lambda with C++20
0206     auto impl = [&, h](auto intType, auto& implRef) {
0207       constexpr int N = decltype(intType)::value;
0208
0209       if constexpr (N == 0) {
0210         return true;
0211       } else {
0212         // If element is invalid: continue
0213         if (!std::get<N - 1>(validExtensions)) {
0214           return implRef(std::integral_constant<int, N - 1>{}, implRef);
0215         }
0216
0217         // Continue as long as evaluations are 'true'
0218         if (std::get<N - 1>(this->tuple())
0219                 .finalize(state, stepper, navigator, h)) {
0220           return implRef(std::integral_constant<int, N - 1>{}, implRef);
0221         } else {
0222           // Break at false
0223           return false;
0224         }
0225       }
0226     };
0227
0228     return impl(std::integral_constant<int, nExtensions>{}, impl);
0229   }
0230 };
0231
0232 }  // namespace Acts