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