Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-11 08:00:35

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