Acts/Propagator/StraightLineStepper.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2016-2020 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 // Workaround for building on clang+libstdc++
0012 #include "Acts/Utilities/detail/ReferenceWrapperAnyCompat.hpp"
0013
0014 #include "Acts/Definitions/Algebra.hpp"
0015 #include "Acts/Definitions/Direction.hpp"
0016 #include "Acts/Definitions/Tolerance.hpp"
0017 #include "Acts/Definitions/TrackParametrization.hpp"
0018 #include "Acts/EventData/TrackParameters.hpp"
0019 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
0020 #include "Acts/Geometry/GeometryContext.hpp"
0021 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0022 #include "Acts/MagneticField/NullBField.hpp"
0023 #include "Acts/Propagator/ConstrainedStep.hpp"
0024 #include "Acts/Propagator/PropagatorTraits.hpp"
0025 #include "Acts/Propagator/detail/SteppingHelper.hpp"
0026 #include "Acts/Surfaces/BoundaryCheck.hpp"
0027 #include "Acts/Surfaces/Surface.hpp"
0028 #include "Acts/Utilities/Intersection.hpp"
0029 #include "Acts/Utilities/Logger.hpp"
0030 #include "Acts/Utilities/Result.hpp"
0031
0032 #include <algorithm>
0033 #include <cmath>
0034 #include <functional>
0035 #include <limits>
0036 #include <string>
0037 #include <tuple>
0038
0039 namespace Acts {
0040
0041 /// @brief straight line stepper based on Surface intersection
0042 ///
0043 /// The straight line stepper is a simple navigation stepper
0044 /// to be used to navigate through the tracking geometry. It can be
0045 /// used for simple material mapping, navigation validation
0046 class StraightLineStepper {
0047  public:
0048   using Jacobian = BoundMatrix;
0049   using Covariance = BoundSquareMatrix;
0050   using BoundState = std::tuple<BoundTrackParameters, Jacobian, double>;
0051   using CurvilinearState =
0052       std::tuple<CurvilinearTrackParameters, Jacobian, double>;
0053   using BField = NullBField;
0054
0055   /// State for track parameter propagation
0056   ///
0057   struct State {
0058     State() = delete;
0059
0060     /// Constructor from the initial bound track parameters
0061     ///
0062     /// @param [in] gctx is the context object for the geometry
0063     /// @param [in] par The track parameters at start
0064     /// @param [in] ssize is the maximum step size
0065     /// @param [in] stolerance is the stepping tolerance
0066     ///
0067     /// @note the covariance matrix is copied when needed
0068     explicit State(const GeometryContext& gctx,
0069                    const MagneticFieldContext& /*mctx*/,
0070                    const BoundTrackParameters& par,
0071                    double ssize = std::numeric_limits<double>::max(),
0072                    double stolerance = s_onSurfaceTolerance)
0073         : particleHypothesis(par.particleHypothesis()),
0074           stepSize(ssize),
0075           tolerance(stolerance),
0076           geoContext(gctx) {
0077       Vector3 position = par.position(gctx);
0078       Vector3 direction = par.direction();
0079       pars.template segment<3>(eFreePos0) = position;
0080       pars.template segment<3>(eFreeDir0) = direction;
0081       pars[eFreeTime] = par.time();
0082       pars[eFreeQOverP] = par.parameters()[eBoundQOverP];
0083       if (par.covariance()) {
0084         // Get the reference surface for navigation
0085         const auto& surface = par.referenceSurface();
0086         // set the covariance transport flag to true and copy
0087         covTransport = true;
0088         cov = BoundSquareMatrix(*par.covariance());
0089         jacToGlobal = surface.boundToFreeJacobian(gctx, position, direction);
0090       }
0091     }
0092
0093     /// Jacobian from local to the global frame
0094     BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
0095
0096     /// Pure transport jacobian part from runge kutta integration
0097     FreeMatrix jacTransport = FreeMatrix::Identity();
0098
0099     /// The full jacobian of the transport entire transport
0100     Jacobian jacobian = Jacobian::Identity();
0101
0102     /// The propagation derivative
0103     FreeVector derivative = FreeVector::Zero();
0104
0105     /// Internal free vector parameters
0106     FreeVector pars = FreeVector::Zero();
0107
0108     /// Particle hypothesis
0109     ParticleHypothesis particleHypothesis = ParticleHypothesis::pion();
0110
0111     /// Boolean to indicate if you need covariance transport
0112     bool covTransport = false;
0113     Covariance cov = Covariance::Zero();
0114
0115     /// accummulated path length state
0116     double pathAccumulated = 0.;
0117
0118     /// Total number of performed steps
0119     std::size_t nSteps = 0;
0120
0121     /// Totoal number of attempted steps
0122     std::size_t nStepTrials = 0;
0123
0124     /// adaptive step size of the runge-kutta integration
0125     ConstrainedStep stepSize;
0126
0127     // Previous step size for overstep estimation (ignored for SL stepper)
0128     double previousStepSize = 0.;
0129
0130     /// The tolerance for the stepping
0131     double tolerance = s_onSurfaceTolerance;
0132
0133     // Cache the geometry context of this propagation
0134     std::reference_wrapper<const GeometryContext> geoContext;
0135   };
0136
0137   /// Always use the same propagation state type, independently of the initial
0138   /// track parameter type and of the target surface
0139   using state_type = State;
0140
0141   StraightLineStepper() = default;
0142
0143   State makeState(std::reference_wrapper<const GeometryContext> gctx,
0144                   std::reference_wrapper<const MagneticFieldContext> mctx,
0145                   const BoundTrackParameters& par,
0146                   double ssize = std::numeric_limits<double>::max(),
0147                   double stolerance = s_onSurfaceTolerance) const {
0148     return State{gctx, mctx, par, ssize, stolerance};
0149   }
0150
0151   /// @brief Resets the state
0152   ///
0153   /// @param [in, out] state State of the stepper
0154   /// @param [in] boundParams Parameters in bound parametrisation
0155   /// @param [in] cov Covariance matrix
0156   /// @param [in] surface The reset @c State will be on this surface
0157   /// @param [in] stepSize Step size
0158   void resetState(
0159       State& state, const BoundVector& boundParams,
0160       const BoundSquareMatrix& cov, const Surface& surface,
0161       const double stepSize = std::numeric_limits<double>::max()) const;
0162
0163   /// Get the field for the stepping, this gives back a zero field
0164   Result<Vector3> getField(State& /*state*/, const Vector3& /*pos*/) const {
0165     // get the field from the cell
0166     return Result<Vector3>::success({0., 0., 0.});
0167   }
0168
0169   /// Global particle position accessor
0170   ///
0171   /// @param state [in] The stepping state (thread-local cache)
0172   Vector3 position(const State& state) const {
0173     return state.pars.template segment<3>(eFreePos0);
0174   }
0175
0176   /// Momentum direction accessor
0177   ///
0178   /// @param state [in] The stepping state (thread-local cache)
0179   Vector3 direction(const State& state) const {
0180     return state.pars.template segment<3>(eFreeDir0);
0181   }
0182
0183   /// QoP direction accessor
0184   ///
0185   /// @param state [in] The stepping state (thread-local cache)
0186   double qOverP(const State& state) const { return state.pars[eFreeQOverP]; }
0187
0188   /// Absolute momentum accessor
0189   ///
0190   /// @param state [in] The stepping state (thread-local cache)
0191   double absoluteMomentum(const State& state) const {
0192     return particleHypothesis(state).extractMomentum(qOverP(state));
0193   }
0194
0195   /// Momentum accessor
0196   ///
0197   /// @param state [in] The stepping state (thread-local cache)
0198   Vector3 momentum(const State& state) const {
0199     return absoluteMomentum(state) * direction(state);
0200   }
0201
0202   /// Charge access
0203   ///
0204   /// @param state [in] The stepping state (thread-local cache)
0205   double charge(const State& state) const {
0206     return particleHypothesis(state).extractCharge(qOverP(state));
0207   }
0208
0209   /// Particle hypothesis
0210   ///
0211   /// @param state [in] The stepping state (thread-local cache)
0212   const ParticleHypothesis& particleHypothesis(const State& state) const {
0213     return state.particleHypothesis;
0214   }
0215
0216   /// Time access
0217   ///
0218   /// @param state [in] The stepping state (thread-local cache)
0219   double time(const State& state) const { return state.pars[eFreeTime]; }
0220
0221   /// Update surface status
0222   ///
0223   /// This method intersects the provided surface and update the navigation
0224   /// step estimation accordingly (hence it changes the state). It also
0225   /// returns the status of the intersection to trigger onSurface in case
0226   /// the surface is reached.
0227   ///
0228   /// @param [in,out] state The stepping state (thread-local cache)
0229   /// @param [in] surface The surface provided
0230   /// @param [in] index The surface intersection index
0231   /// @param [in] navDir The navigation direction
0232   /// @param [in] bcheck The boundary check for this status update
0233   /// @param [in] surfaceTolerance Surface tolerance used for intersection
0234   /// @param [in] logger A logger instance
0235   Intersection3D::Status updateSurfaceStatus(
0236       State& state, const Surface& surface, std::uint8_t index,
0237       Direction navDir, const BoundaryCheck& bcheck,
0238       ActsScalar surfaceTolerance = s_onSurfaceTolerance,
0239       const Logger& logger = getDummyLogger()) const {
0240     return detail::updateSingleSurfaceStatus<StraightLineStepper>(
0241         *this, state, surface, index, navDir, bcheck, surfaceTolerance, logger);
0242   }
0243
0244   /// Update step size
0245   ///
0246   /// It checks the status to the reference surface & updates
0247   /// the step size accordingly
0248   ///
0249   /// @param state [in,out] The stepping state (thread-local cache)
0250   /// @param oIntersection [in] The ObjectIntersection to layer, boundary, etc
0251   /// @param release [in] boolean to trigger step size release
0252   template <typename object_intersection_t>
0253   void updateStepSize(State& state, const object_intersection_t& oIntersection,
0254                       Direction /*direction*/, bool release = true) const {
0255     detail::updateSingleStepSize<StraightLineStepper>(state, oIntersection,
0256                                                       release);
0257   }
0258
0259   /// Update step size - explicitly with a double
0260   ///
0261   /// @param state [in,out] The stepping state (thread-local cache)
0262   /// @param stepSize [in] The step size value
0263   /// @param stype [in] The step size type to be set
0264   /// @param release [in] Do we release the step size?
0265   void updateStepSize(State& state, double stepSize,
0266                       ConstrainedStep::Type stype = ConstrainedStep::actor,
0267                       bool release = true) const {
0268     state.previousStepSize = state.stepSize.value();
0269     state.stepSize.update(stepSize, stype, release);
0270   }
0271
0272   /// Get the step size
0273   ///
0274   /// @param state [in] The stepping state (thread-local cache)
0275   /// @param stype [in] The step size type to be returned
0276   double getStepSize(const State& state, ConstrainedStep::Type stype) const {
0277     return state.stepSize.value(stype);
0278   }
0279
0280   /// Release the Step size
0281   ///
0282   /// @param [in,out] state The stepping state (thread-local cache)
0283   /// @param [in] stype The step size type to be released
0284   void releaseStepSize(State& state, ConstrainedStep::Type stype) const {
0285     state.stepSize.release(stype);
0286   }
0287
0288   /// Output the Step Size - single component
0289   ///
0290   /// @param state [in,out] The stepping state (thread-local cache)
0291   std::string outputStepSize(const State& state) const {
0292     return state.stepSize.toString();
0293   }
0294
0295   /// Create and return the bound state at the current position
0296   ///
0297   /// @brief It does not check if the transported state is at the surface, this
0298   /// needs to be guaranteed by the propagator
0299   ///
0300   /// @param [in] state State that will be presented as @c BoundState
0301   /// @param [in] surface The surface to which we bind the state
0302   /// @param [in] transportCov Flag steering covariance transport
0303   /// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound
0304   ///
0305   /// @return A bound state:
0306   ///   - the parameters at the surface
0307   ///   - the stepwise jacobian towards it (from last bound)
0308   ///   - and the path length (from start - for ordering)
0309   Result<BoundState> boundState(
0310       State& state, const Surface& surface, bool transportCov = true,
0311       const FreeToBoundCorrection& freeToBoundCorrection =
0312           FreeToBoundCorrection(false)) const;
0313
0314   /// @brief If necessary fill additional members needed for curvilinearState
0315   ///
0316   /// Compute path length derivatives in case they have not been computed
0317   /// yet, which is the case if no step has been executed yet.
0318   ///
0319   /// @param [in, out] prop_state State that will be presented as @c BoundState
0320   /// @param [in] navigator the navigator of the propagation
0321   /// @return true if nothing is missing after this call, false otherwise.
0322   template <typename propagator_state_t, typename navigator_t>
0323   bool prepareCurvilinearState(
0324       [[maybe_unused]] propagator_state_t& prop_state,
0325       [[maybe_unused]] const navigator_t& navigator) const {
0326     // test whether the accumulated path has still its initial value.
0327     if (prop_state.stepping.pathAccumulated == 0.) {
0328       // dr/ds :
0329       prop_state.stepping.derivative.template head<3>() =
0330           direction(prop_state.stepping);
0331       // dt / ds
0332       prop_state.stepping.derivative(eFreeTime) =
0333           std::hypot(1., prop_state.stepping.particleHypothesis.mass() /
0334                              absoluteMomentum(prop_state.stepping));
0335       // d (dr/ds) / ds : == 0
0336       prop_state.stepping.derivative.template segment<3>(4) =
0337           Acts::Vector3::Zero().transpose();
0338       // d qop / ds  == 0
0339       prop_state.stepping.derivative(eFreeQOverP) = 0.;
0340     }
0341     return true;
0342   }
0343
0344   /// Create and return a curvilinear state at the current position
0345   ///
0346   /// @brief This creates a curvilinear state.
0347   ///
0348   /// @param [in] state State that will be presented as @c CurvilinearState
0349   /// @param [in] transportCov Flag steering covariance transport
0350   ///
0351   /// @return A curvilinear state:
0352   ///   - the curvilinear parameters at given position
0353   ///   - the stepweise jacobian towards it (from last bound)
0354   ///   - and the path length (from start - for ordering)
0355   CurvilinearState curvilinearState(State& state,
0356                                     bool transportCov = true) const;
0357
0358   /// Method to update a stepper state to the some parameters
0359   ///
0360   /// @param [in,out] state State object that will be updated
0361   /// @param [in] freeParams Free parameters that will be written into @p state
0362   /// @param [in] boundParams Corresponding bound parameters used to update jacToGlobal in @p state
0363   /// @param [in] covariance Covariance that will be written into @p state
0364   /// @param [in] surface The surface used to update the jacToGlobal
0365   void update(State& state, const FreeVector& freeParams,
0366               const BoundVector& boundParams, const Covariance& covariance,
0367               const Surface& surface) const;
0368
0369   /// Method to update the stepper state
0370   ///
0371   /// @param [in,out] state State object that will be updated
0372   /// @param [in] uposition the updated position
0373   /// @param [in] udirection the updated direction
0374   /// @param [in] qop the updated qop value
0375   /// @param [in] time the updated time value
0376   void update(State& state, const Vector3& uposition, const Vector3& udirection,
0377               double qop, double time) const;
0378
0379   /// Method for on-demand transport of the covariance
0380   /// to a new curvilinear frame at current  position,
0381   /// or direction of the state - for the moment a dummy method
0382   ///
0383   /// @param [in,out] state State of the stepper
0384   void transportCovarianceToCurvilinear(State& state) const;
0385
0386   /// Method for on-demand transport of the covariance
0387   /// to a new curvilinear frame at current  position,
0388   /// or direction of the state - for the moment a dummy method
0389   ///
0390   /// @tparam surface_t the surface type - ignored here
0391   ///
0392   /// @param [in,out] state The stepper state
0393   /// @param [in] surface is the surface to which the covariance is
0394   ///        forwarded to
0395   /// @note no check is done if the position is actually on the surface
0396   /// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound
0397   ///
0398   void transportCovarianceToBound(
0399       State& state, const Surface& surface,
0400       const FreeToBoundCorrection& freeToBoundCorrection =
0401           FreeToBoundCorrection(false)) const;
0402
0403   /// Perform a straight line propagation step
0404   ///
0405   /// @param [in,out] state is the propagation state associated with the track
0406   /// parameters that are being propagated.
0407   ///                The state contains the desired step size,
0408   ///                it can be negative during backwards track propagation,
0409   ///                and since we're using an adaptive algorithm, it can
0410   ///                be modified by the stepper class during propagation.
0411   ///
0412   /// @return the step size taken
0413   template <typename propagator_state_t, typename navigator_t>
0414   Result<double> step(propagator_state_t& state,
0415                       const navigator_t& /*navigator*/) const {
0416     // use the adjusted step size
0417     const auto h = state.stepping.stepSize.value() * state.options.direction;
0418     const auto m = state.stepping.particleHypothesis.mass();
0419     const auto p = absoluteMomentum(state.stepping);
0420     // time propagates along distance as 1/b = sqrt(1 + m²/p²)
0421     const auto dtds = std::hypot(1., m / p);
0422     // Update the track parameters according to the equations of motion
0423     Vector3 dir = direction(state.stepping);
0424     state.stepping.pars.template segment<3>(eFreePos0) += h * dir;
0425     state.stepping.pars[eFreeTime] += h * dtds;
0426     // Propagate the jacobian
0427     if (state.stepping.covTransport) {
0428       // The step transport matrix in global coordinates
0429       FreeMatrix D = FreeMatrix::Identity();
0430       D.block<3, 3>(0, 4) = ActsSquareMatrix<3>::Identity() * h;
0431       // Extend the calculation by the time propagation
0432       // Evaluate dt/dlambda
0433       D(3, 7) = h * m * m * state.stepping.pars[eFreeQOverP] / dtds;
0434       // Set the derivative factor the time
0435       state.stepping.derivative(3) = dtds;
0436       // Update jacobian and derivative
0437       state.stepping.jacTransport = D * state.stepping.jacTransport;
0438       state.stepping.derivative.template head<3>() = dir;
0439     }
0440
0441     // state the path length
0442     state.stepping.pathAccumulated += h;
0443     ++state.stepping.nSteps;
0444     ++state.stepping.nStepTrials;
0445
0446     // return h
0447     return h;
0448   }
0449 };
0450
0451 template <typename navigator_t>
0452 struct SupportsBoundParameters<StraightLineStepper, navigator_t>
0453     : public std::true_type {};
0454
0455 }  // namespace Acts