Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-10-30 07:53:45

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