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