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