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