Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-10-07 08:48:03

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/EventData/BoundTrackParameters.hpp"
0014 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
0015 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0016 #include "Acts/Propagator/ConstrainedStep.hpp"
0017 #include "Acts/Propagator/NavigationTarget.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 
0023 #include <optional>
0024 #include <stdexcept>
0025 #include <string>
0026 
0027 namespace Acts {
0028 
0029 class IVolumeMaterial;
0030 
0031 /// @brief Stepper that moves the track along a helix in closed form
0032 ///
0033 /// Each step reads the magnetic field at the start position and moves the
0034 /// track along the exact helix in that field. The transport jacobian is also
0035 /// closed form. It is the exact derivative of the helix step, so it treats
0036 /// the field as constant over the step and has no field gradient term.
0037 ///
0038 /// The field at the end of a step is also the field at the start of the next
0039 /// step. The stepper uses it to estimate the position error from the field
0040 /// change along the step, and it adapts the step size to
0041 /// @ref StepperPlainOptions::stepTolerance. An accepted step needs no extra
0042 /// field lookup.
0043 ///
0044 /// The navigation estimates the distance to a surface with a straight line.
0045 /// A step of that length along the helix passes a surface perpendicular to the
0046 /// track by about |omega x T|^2 h^3 / 3, with omega = (q/p) B, and the next
0047 /// step goes back. A target aborter must accept an intersection that far behind
0048 /// the track. The default near limit of @ref SurfaceReached is 100 um.
0049 ///
0050 /// @note The stepper propagates in vacuum only. It ignores volume material.
0051 class HelixStepper final {
0052  public:
0053   /// Type alias for bound track parameters
0054   using BoundParameters = BoundTrackParameters;
0055   /// Type alias for jacobian matrix
0056   using Jacobian = BoundMatrix;
0057   /// Type alias for covariance matrix
0058   using Covariance = BoundMatrix;
0059 
0060   /// Configuration for the helix stepper.
0061   struct Config {
0062     /// Magnetic field provider
0063     std::shared_ptr<const MagneticFieldProvider> bField;
0064   };
0065 
0066   /// Runtime options for helix propagation.
0067   struct Options : public StepperPlainOptions {
0068     /// Constructor
0069     /// @param gctx Geometry context
0070     /// @param mctx Magnetic field context
0071     Options(const GeometryContext& gctx, const MagneticFieldContext& mctx)
0072         : StepperPlainOptions(gctx, mctx) {}
0073 
0074     /// Set plain stepper options
0075     /// @param options Plain stepper options
0076     void setPlainOptions(const StepperPlainOptions& options) {
0077       static_cast<StepperPlainOptions&>(*this) = options;
0078     }
0079   };
0080 
0081   /// @brief State for track parameter propagation
0082   ///
0083   /// It contains the stepping information and is provided thread local
0084   /// by the propagator
0085   struct State {
0086     /// Constructor from the options and the field cache
0087     ///
0088     /// @param [in] optionsIn is the configuration of the stepper
0089     /// @param [in] fieldCacheIn is the cache object for the magnetic field
0090     State(const Options& optionsIn, MagneticFieldProvider::Cache fieldCacheIn)
0091         : options(optionsIn), fieldCache(std::move(fieldCacheIn)) {}
0092 
0093     /// Configuration options for the stepper
0094     Options options;
0095 
0096     /// Internal free vector parameters
0097     FreeVector pars = FreeVector::Zero();
0098 
0099     /// The propagation derivative
0100     FreeVector derivative = FreeVector::Zero();
0101 
0102     /// Jacobian from local to the global frame
0103     BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
0104 
0105     /// Pure transport jacobian part from the helix steps
0106     FreeMatrix jacTransport = FreeMatrix::Identity();
0107 
0108     /// Covariance matrix for track parameter uncertainties, set if the
0109     /// covariance is transported
0110     std::optional<Covariance> cov;
0111 
0112     /// Particle hypothesis
0113     ParticleHypothesis particleHypothesis = ParticleHypothesis::pion();
0114 
0115     /// Adaptive step size of the helix steps
0116     ConstrainedStep stepSize;
0117 
0118     /// Last performed step (for overstep limit calculation)
0119     double previousStepSize = 0.;
0120 
0121     /// Magnetic field at the current position, reused as the next step's
0122     /// field. Reset whenever the position is set from outside.
0123     std::optional<Vector3> field;
0124 
0125     /// Accumulated path length state
0126     double pathAccumulated = 0.;
0127 
0128     /// Total number of performed steps
0129     std::size_t nSteps = 0;
0130 
0131     /// Statistics of the stepper
0132     StepperStatistics statistics;
0133 
0134     /// This caches the current magnetic field cell and stays
0135     /// (and interpolates) within it as long as this is valid.
0136     MagneticFieldProvider::Cache fieldCache;
0137   };
0138 
0139   /// Constructor requires knowledge of the detector's magnetic field
0140   /// @param bField The magnetic field provider
0141   explicit HelixStepper(std::shared_ptr<const MagneticFieldProvider> bField);
0142 
0143   /// @brief Constructor with configuration
0144   /// @param config The configuration of the stepper
0145   explicit HelixStepper(const Config& config);
0146 
0147   /// Create a state object
0148   /// @param options Stepper options
0149   /// @return State object
0150   State makeState(const Options& options) const;
0151 
0152   /// Initialize the state from bound track parameters
0153   /// @param state The state to initialize
0154   /// @param par The bound track parameters
0155   void initialize(State& state, const BoundParameters& par) const;
0156 
0157   /// Initialize the state from bound parameters
0158   /// @param state The state to initialize
0159   /// @param boundParams Bound track parameters vector
0160   /// @param cov Covariance matrix
0161   /// @param particleHypothesis Particle hypothesis
0162   /// @param surface Reference surface
0163   void initialize(State& state, const BoundVector& boundParams,
0164                   const std::optional<BoundMatrix>& cov,
0165                   ParticleHypothesis particleHypothesis,
0166                   const Surface& surface) const;
0167 
0168   /// Get the field for the stepping, it checks first if the access is still
0169   /// within the Cell, and updates the cell if necessary.
0170   ///
0171   /// @param [in,out] state is the propagation state associated with the track
0172   ///                 the magnetic field cell is used (and potentially updated)
0173   /// @param [in] pos is the field position
0174   /// @return Magnetic field vector
0175   Result<Vector3> getField(State& state, const Vector3& pos) const {
0176     return m_bField->getField(pos, state.fieldCache);
0177   }
0178 
0179   /// Global particle position accessor
0180   ///
0181   /// @param state [in] The stepping state (thread-local cache)
0182   /// @return Position vector
0183   Vector3 position(const State& state) const {
0184     return state.pars.template segment<3>(eFreePos0);
0185   }
0186 
0187   /// Momentum direction accessor
0188   ///
0189   /// @param state [in] The stepping state (thread-local cache)
0190   /// @return Direction vector
0191   Vector3 direction(const State& state) const {
0192     return state.pars.template segment<3>(eFreeDir0);
0193   }
0194 
0195   /// QoP direction accessor
0196   ///
0197   /// @param state [in] The stepping state (thread-local cache)
0198   /// @return Charge over momentum
0199   double qOverP(const State& state) const { return state.pars[eFreeQOverP]; }
0200 
0201   /// Absolute momentum accessor
0202   ///
0203   /// @param state [in] The stepping state (thread-local cache)
0204   /// @return Absolute momentum
0205   double absoluteMomentum(const State& state) const {
0206     return particleHypothesis(state).extractMomentum(qOverP(state));
0207   }
0208 
0209   /// Momentum accessor
0210   ///
0211   /// @param state [in] The stepping state (thread-local cache)
0212   /// @return Momentum vector
0213   Vector3 momentum(const State& state) const {
0214     return absoluteMomentum(state) * direction(state);
0215   }
0216 
0217   /// Charge access
0218   ///
0219   /// @param state [in] The stepping state (thread-local cache)
0220   /// @return Particle charge
0221   double charge(const State& state) const {
0222     return particleHypothesis(state).extractCharge(qOverP(state));
0223   }
0224 
0225   /// Particle hypothesis
0226   ///
0227   /// @param state [in] The stepping state (thread-local cache)
0228   /// @return Particle hypothesis
0229   const ParticleHypothesis& particleHypothesis(const State& state) const {
0230     return state.particleHypothesis;
0231   }
0232 
0233   /// Time access
0234   ///
0235   /// @param state [in] The stepping state (thread-local cache)
0236   /// @return Time
0237   double time(const State& state) const { return state.pars[eFreeTime]; }
0238 
0239   /// Update surface status
0240   ///
0241   /// It checks the status to the reference surface & updates
0242   /// the step size accordingly
0243   ///
0244   /// @param [in,out] state The stepping state (thread-local cache)
0245   /// @param [in] surface The surface provided
0246   /// @param [in] index The surface intersection index
0247   /// @param [in] navDir The navigation direction
0248   /// @param [in] boundaryTolerance The boundary check for this status update
0249   /// @param [in] surfaceTolerance Surface tolerance used for intersection
0250   /// @param [in] stype The step size type to be set
0251   /// @param [in] logger A @c Logger instance
0252   /// @return Intersection status
0253   IntersectionStatus updateSurfaceStatus(
0254       State& state, const Surface& surface, std::uint8_t index,
0255       Direction navDir, const BoundaryTolerance& boundaryTolerance,
0256       double surfaceTolerance, ConstrainedStep::Type stype,
0257       const Logger& logger = getDummyLogger()) const {
0258     return detail::updateSingleSurfaceStatus<HelixStepper>(
0259         *this, state, surface, index, navDir, boundaryTolerance,
0260         surfaceTolerance, stype, logger);
0261   }
0262 
0263   /// Update step size
0264   ///
0265   /// This method intersects the provided surface and update the navigation
0266   /// step estimation accordingly (hence it changes the state). It also
0267   /// returns the status of the intersection to trigger onSurface in case
0268   /// the surface is reached.
0269   ///
0270   /// @param state [in,out] The stepping state (thread-local cache)
0271   /// @param target [in] The NavigationTarget
0272   /// @param direction [in] The propagation direction
0273   /// @param stype [in] The step size type to be set
0274   void updateStepSize(State& state, const NavigationTarget& target,
0275                       Direction direction, ConstrainedStep::Type stype) const {
0276     static_cast<void>(direction);
0277     double stepSize = target.pathLength();
0278     updateStepSize(state, stepSize, stype);
0279   }
0280 
0281   /// Update step size - explicitly with a double
0282   ///
0283   /// @param state [in,out] The stepping state (thread-local cache)
0284   /// @param stepSize [in] The step size value
0285   /// @param stype [in] The step size type to be set
0286   void updateStepSize(State& state, double stepSize,
0287                       ConstrainedStep::Type stype) const {
0288     state.previousStepSize = state.stepSize.value();
0289     state.stepSize.update(stepSize, stype);
0290   }
0291 
0292   /// Get the step size
0293   ///
0294   /// @param state [in] The stepping state (thread-local cache)
0295   /// @param stype [in] The step size type to be returned
0296   /// @return Step size
0297   double getStepSize(const State& state, ConstrainedStep::Type stype) const {
0298     return state.stepSize.value(stype);
0299   }
0300 
0301   /// Release the Step size
0302   ///
0303   /// @param state [in,out] The stepping state (thread-local cache)
0304   /// @param [in] stype The step size type to be released
0305   void releaseStepSize(State& state, ConstrainedStep::Type stype) const {
0306     state.stepSize.release(stype);
0307   }
0308 
0309   /// Output the Step Size - single component
0310   ///
0311   /// @param state [in,out] The stepping state (thread-local cache)
0312   /// @return String representation of step size
0313   std::string outputStepSize(const State& state) const {
0314     return state.stepSize.toString();
0315   }
0316 
0317   /// Get the step size constraints
0318   ///
0319   /// @param state [in] The stepping state (thread-local cache)
0320   /// @return The step size constraints
0321   const ConstrainedStep& stepSize(const State& state) const {
0322     return state.stepSize;
0323   }
0324 
0325   /// Get the stepper statistics
0326   ///
0327   /// @param state [in] The stepping state (thread-local cache)
0328   /// @return The statistics since the last initialization
0329   const StepperStatistics& statistics(const State& state) const {
0330     return state.statistics;
0331   }
0332 
0333   /// Get the path length
0334   ///
0335   /// @param state [in] The stepping state (thread-local cache)
0336   /// @return The path length since the last initialization
0337   double pathLength(const State& state) const { return state.pathAccumulated; }
0338 
0339   /// Check if the state carries a covariance
0340   ///
0341   /// @param state [in] The stepping state (thread-local cache)
0342   /// @return True if the covariance is transported
0343   bool hasCovariance(const State& state) const { return state.cov.has_value(); }
0344 
0345   /// Get the covariance at the anchor
0346   ///
0347   /// The anchor is the frame of the last initialization, transport or update.
0348   ///
0349   /// @param state [in] The stepping state (thread-local cache)
0350   /// @return The covariance at the anchor, or no value if the state does not
0351   ///         carry a covariance
0352   const std::optional<Covariance>& covariance(const State& state) const {
0353     return state.cov;
0354   }
0355 
0356   /// Set the covariance at the anchor
0357   ///
0358   /// The state must already carry a covariance, because the stepper only
0359   /// transports the Jacobian while it has one.
0360   ///
0361   /// @param [in,out] state The stepping state (thread-local cache)
0362   /// @param [in] covariance The new covariance at the anchor
0363   /// @throws std::logic_error if the state does not carry a covariance
0364   void setCovariance(State& state, const Covariance& covariance) const {
0365     if (!state.cov.has_value()) {
0366       throw std::logic_error(
0367           "Cannot set the covariance of a state without a covariance");
0368     }
0369     state.cov = covariance;
0370   }
0371 
0372   /// Get the bound parameters at the current position
0373   ///
0374   /// The parameters carry the covariance at the anchor if the state has one.
0375   ///
0376   /// @note It does not check if the state is on @p surface or anchored on it
0377   ///
0378   /// @param [in] state The stepping state (thread-local cache)
0379   /// @param [in] surface The surface of the parameters
0380   /// @return The bound parameters, or a failure if the position cannot be
0381   ///         expressed on @p surface
0382   Result<BoundParameters> boundParameters(const State& state,
0383                                           const Surface& surface) const;
0384 
0385   /// @brief If necessary fill additional members needed for
0386   /// transportToCurvilinear
0387   ///
0388   /// Compute path length derivatives in case they have not been computed
0389   /// yet, which is the case if no step has been executed yet.
0390   ///
0391   /// @param [in, out] state State of the stepper
0392   /// @return true if nothing is missing after this call, false otherwise.
0393   bool prepareCurvilinearState(State& state) const;
0394 
0395   /// Get the curvilinear parameters at the current position
0396   ///
0397   /// The parameters carry the covariance at the anchor if the state has one.
0398   ///
0399   /// @param [in] state The stepping state (thread-local cache)
0400   /// @return The curvilinear parameters
0401   BoundParameters curvilinearParameters(const State& state) const;
0402 
0403   /// Method to update a stepper state to the some parameters
0404   ///
0405   /// This anchors the state on @p surface.
0406   ///
0407   /// @param [in,out] state State object that will be updated
0408   /// @param [in] freeParams Free parameters that will be written into @p state
0409   /// @param [in] boundParams Corresponding bound parameters used to update jacToGlobal in @p state
0410   /// @param [in] covariance The covariance that will be written into @p state
0411   ///                        if the state carries a covariance
0412   /// @param [in] surface The surface used to update the jacToGlobal
0413   void update(State& state, const FreeVector& freeParams,
0414               const BoundVector& boundParams, const Covariance& covariance,
0415               const Surface& surface) const;
0416 
0417   /// Method to update the stepper state
0418   ///
0419   /// @param [in,out] state State object that will be updated
0420   /// @param [in] uposition the updated position
0421   /// @param [in] udirection the updated direction
0422   /// @param [in] qop the updated qop value
0423   /// @param [in] time the updated time value
0424   void update(State& state, const Vector3& uposition, const Vector3& udirection,
0425               double qop, double time) const;
0426 
0427   /// Transport the covariance to the curvilinear frame at the current position
0428   ///
0429   /// This anchors the state on the curvilinear frame. Without a covariance
0430   /// the state does not change.
0431   ///
0432   /// @param [in,out] state State of the stepper
0433   /// @return The jacobian from the previous anchor to the curvilinear frame
0434   Jacobian transportToCurvilinear(State& state) const;
0435 
0436   /// Transport the covariance to a surface at the current position
0437   ///
0438   /// This anchors the state on @p surface. Without a covariance the state
0439   /// does not change. The state must be on @p surface, and it stays unchanged
0440   /// if it is not.
0441   ///
0442   /// @param [in,out] state State of the stepper
0443   /// @param [in] surface The surface to transport the covariance to
0444   /// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound
0445   /// @return The jacobian from the previous anchor to @p surface, or a failure
0446   ///         if the state is not on @p surface
0447   Result<Jacobian> transportToBound(
0448       State& state, const Surface& surface,
0449       const FreeToBoundCorrection& freeToBoundCorrection =
0450           FreeToBoundCorrection(false)) const;
0451 
0452   /// Perform a helix track parameter propagation step
0453   ///
0454   /// @param [in,out] state State of the stepper
0455   /// @param propDir is the direction of propagation
0456   /// @param material is ignored, the stepper propagates in vacuum only
0457   /// @return the result of the step
0458   ///
0459   /// @note The state contains the desired step size. It can be negative during
0460   ///       backwards track propagation, and since the step size is adaptive,
0461   ///       it can be modified by the stepper class during propagation.
0462   Result<double> step(State& state, Direction propDir,
0463                       const IVolumeMaterial* material) const;
0464 
0465  private:
0466   /// Magnetic field inside of the detector
0467   std::shared_ptr<const MagneticFieldProvider> m_bField;
0468 };
0469 
0470 template <>
0471 struct SupportsBoundParameters<HelixStepper> : public std::true_type {};
0472 
0473 }  // namespace Acts