Acts/Propagator/EigenStepperDefaultExtension.hpp

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/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Utilities/VectorHelpers.hpp"
0014
0015 #include <array>
0016
0017 namespace Acts {
0018
0019 class IVolumeMaterial;
0020
0021 /// @brief Default evaluator of the k_i's and elements of the transport matrix
0022 /// D of the RKN4 stepping. This is a pure implementation by textbook.
0023 struct EigenStepperDefaultExtension {
0024   /// @brief Evaluator of the k_i's of the RKN4. For the case of i = 0 this
0025   /// step sets up qop, too.
0026   ///
0027   /// @tparam i Index of the k_i, i = [0, 3]
0028   /// @tparam stepper_t Type of the stepper
0029   ///
0030   /// @param [in] state State of the stepper
0031   /// @param [in] stepper Stepper of the propagation
0032   /// @param [in] volumeMaterial Material of the volume
0033   /// @param [out] knew Next k_i that is evaluated
0034   /// @param [in] bField B-Field at the evaluation position
0035   /// @param [out] kQoP k_i elements of the momenta
0036   /// @param [in] h Step size (= 0. ^ 0.5 * StepSize ^ StepSize)
0037   /// @param [in] kprev Evaluated k_{i - 1}
0038   ///
0039   /// @return Boolean flag if the calculation is valid
0040   template <int i, typename stepper_t>
0041   bool k(const typename stepper_t::State& state, const stepper_t& stepper,
0042          const IVolumeMaterial* volumeMaterial, Vector3& knew,
0043          const Vector3& bField, std::array<double, 4>& kQoP,
0044          const double h = 0., const Vector3& kprev = Vector3::Zero())
0045     requires(i >= 0 && i <= 3)
0046   {
0047     (void)volumeMaterial;
0048
0049     auto qop = stepper.qOverP(state);
0050     // First step does not rely on previous data
0051     if constexpr (i == 0) {
0052       knew = qop * stepper.direction(state).cross(bField);
0053       kQoP = {0., 0., 0., 0.};
0054     } else {
0055       knew = qop * (stepper.direction(state) + h * kprev).cross(bField);
0056     }
0057     return true;
0058   }
0059
0060   /// @brief Veto function after a RKN4 step was accepted by judging on the
0061   /// error of the step. Since the textbook does not deliver further vetos,
0062   /// this is a dummy function.
0063   ///
0064   /// @tparam stepper_t Type of the stepper
0065   ///
0066   /// @param [in] state State of the stepper
0067   /// @param [in] stepper Stepper of the propagation
0068   /// @param [in] volumeMaterial Material of the volume
0069   /// @param [in] h Step size
0070   ///
0071   /// @return Boolean flag if the calculation is valid
0072   template <typename stepper_t>
0073   bool finalize(typename stepper_t::State& state, const stepper_t& stepper,
0074                 const IVolumeMaterial* volumeMaterial, const double h) const {
0075     (void)volumeMaterial;
0076
0077     propagateTime(state, stepper, h);
0078     return true;
0079   }
0080
0081   /// @brief Veto function after a RKN4 step was accepted by judging on the
0082   /// error of the step. Since the textbook does not deliver further vetos,
0083   /// this is just for the evaluation of the transport matrix.
0084   ///
0085   /// @tparam stepper_t Type of the stepper
0086   ///
0087   /// @param [in] state State of the stepper
0088   /// @param [in] stepper Stepper of the propagation
0089   /// @param [in] volumeMaterial Material of the volume
0090   /// @param [in] h Step size
0091   /// @param [out] D Transport matrix
0092   ///
0093   /// @return Boolean flag if the calculation is valid
0094   template <typename stepper_t>
0095   bool finalize(typename stepper_t::State& state, const stepper_t& stepper,
0096                 const IVolumeMaterial* volumeMaterial, const double h,
0097                 FreeMatrix& D) const {
0098     (void)volumeMaterial;
0099
0100     propagateTime(state, stepper, h);
0101     return transportMatrix(state, stepper, h, D);
0102   }
0103
0104  private:
0105   /// @brief Propagation function for the time coordinate
0106   ///
0107   /// @tparam stepper_t Type of the stepper
0108   ///
0109   /// @param [in, out] state State of the stepper
0110   /// @param [in] stepper Stepper of the propagation
0111   /// @param [in] h Step size
0112   template <typename stepper_t>
0113   void propagateTime(typename stepper_t::State& state, const stepper_t& stepper,
0114                      const double h) const {
0115     /// This evaluation is based on dt/ds = 1/v = 1/(beta * c) with the velocity
0116     /// v, the speed of light c and beta = v/c. This can be re-written as dt/ds
0117     /// = sqrt(m^2/p^2 + c^{-2}) with the mass m and the momentum p.
0118     auto m = stepper.particleHypothesis(state).mass();
0119     auto p = stepper.absoluteMomentum(state);
0120     auto dtds = std::sqrt(1 + m * m / (p * p));
0121     state.pars[eFreeTime] += h * dtds;
0122     if (state.covTransport) {
0123       state.derivative(3) = dtds;
0124     }
0125   }
0126
0127   /// @brief Calculates the transport matrix D for the jacobian
0128   ///
0129   /// @tparam stepper_t Type of the stepper
0130   ///
0131   /// @param [in] state State of the stepper
0132   /// @param [in] stepper Stepper of the propagation
0133   /// @param [in] h Step size
0134   /// @param [out] D Transport matrix
0135   ///
0136   /// @return Boolean flag if evaluation is valid
0137   template <typename stepper_t>
0138   bool transportMatrix(typename stepper_t::State& state,
0139                        const stepper_t& stepper, const double h,
0140                        FreeMatrix& D) const {
0141     /// The calculations are based on ATL-SOFT-PUB-2009-002. The update of the
0142     /// Jacobian matrix is requires only the calculation of eq. 17 and 18.
0143     /// Since the terms of eq. 18 are currently 0, this matrix is not needed
0144     /// in the calculation. The matrix A from eq. 17 consists out of 3
0145     /// different parts. The first one is given by the upper left 3x3 matrix
0146     /// that are calculated by the derivatives dF/dT (called dFdT) and dG/dT
0147     /// (calls dGdT). The second is given by the top 3 lines of the rightmost
0148     /// column. This is calculated by dFdL and dGdL. The remaining non-zero term
0149     /// is calculated directly. The naming of the variables is explained in eq.
0150     /// 11 and are directly related to the initial problem in eq. 7.
0151     /// The evaluation is based by propagating the parameters T and lambda as
0152     /// given in eq. 16 and evaluating the derivations for matrix A.
0153     /// @note The translation for u_{n+1} in eq. 7 is in this case a
0154     /// 3-dimensional vector without a dependency of Lambda or lambda neither in
0155     /// u_n nor in u_n'. The second and fourth eq. in eq. 14 have the constant
0156     /// offset matrices h * Id and Id respectively. This involves that the
0157     /// constant offset does not exist for rectangular matrix dGdu' (due to the
0158     /// missing Lambda part) and only exists for dFdu' in dlambda/dlambda.
0159
0160     auto m = state.particleHypothesis.mass();
0161     auto& sd = state.stepData;
0162     auto dir = stepper.direction(state);
0163     auto qop = stepper.qOverP(state);
0164     auto p = stepper.absoluteMomentum(state);
0165     auto dtds = std::sqrt(1 + m * m / (p * p));
0166
0167     D = FreeMatrix::Identity();
0168
0169     double half_h = h * 0.5;
0170     // This sets the reference to the sub matrices
0171     // dFdx is already initialised as (3x3) idendity
0172     auto dFdT = D.block<3, 3>(0, 4);
0173     auto dFdL = D.block<3, 1>(0, 7);
0174     // dGdx is already initialised as (3x3) zero
0175     auto dGdT = D.block<3, 3>(4, 4);
0176     auto dGdL = D.block<3, 1>(4, 7);
0177
0178     SquareMatrix3 dk1dT = SquareMatrix3::Zero();
0179     SquareMatrix3 dk2dT = SquareMatrix3::Identity();
0180     SquareMatrix3 dk3dT = SquareMatrix3::Identity();
0181     SquareMatrix3 dk4dT = SquareMatrix3::Identity();
0182
0183     Vector3 dk1dL = Vector3::Zero();
0184     Vector3 dk2dL = Vector3::Zero();
0185     Vector3 dk3dL = Vector3::Zero();
0186     Vector3 dk4dL = Vector3::Zero();
0187
0188     // For the case without energy loss
0189     dk1dL = dir.cross(sd.B_first);
0190     dk2dL = (dir + half_h * sd.k1).cross(sd.B_middle) +
0191             qop * half_h * dk1dL.cross(sd.B_middle);
0192     dk3dL = (dir + half_h * sd.k2).cross(sd.B_middle) +
0193             qop * half_h * dk2dL.cross(sd.B_middle);
0194     dk4dL =
0195         (dir + h * sd.k3).cross(sd.B_last) + qop * h * dk3dL.cross(sd.B_last);
0196
0197     dk1dT(0, 1) = sd.B_first.z();
0198     dk1dT(0, 2) = -sd.B_first.y();
0199     dk1dT(1, 0) = -sd.B_first.z();
0200     dk1dT(1, 2) = sd.B_first.x();
0201     dk1dT(2, 0) = sd.B_first.y();
0202     dk1dT(2, 1) = -sd.B_first.x();
0203     dk1dT *= qop;
0204
0205     dk2dT += half_h * dk1dT;
0206     dk2dT = qop * VectorHelpers::cross(dk2dT, sd.B_middle);
0207
0208     dk3dT += half_h * dk2dT;
0209     dk3dT = qop * VectorHelpers::cross(dk3dT, sd.B_middle);
0210
0211     dk4dT += h * dk3dT;
0212     dk4dT = qop * VectorHelpers::cross(dk4dT, sd.B_last);
0213
0214     dFdT.setIdentity();
0215     dFdT += h / 6. * (dk1dT + dk2dT + dk3dT);
0216     dFdT *= h;
0217
0218     dFdL = (h * h) / 6. * (dk1dL + dk2dL + dk3dL);
0219
0220     dGdT += h / 6. * (dk1dT + 2. * (dk2dT + dk3dT) + dk4dT);
0221
0222     dGdL = h / 6. * (dk1dL + 2. * (dk2dL + dk3dL) + dk4dL);
0223
0224     D(3, 7) = h * m * m * qop / dtds;
0225     return true;
0226   }
0227 };
0228
0229 }  // namespace Acts