Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-12-15 09:24:04

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/TrackParameters.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/MaterialEffectsAccumulator.hpp"
0022 #include "Acts/Propagator/detail/SteppingHelper.hpp"
0023 
0024 namespace Acts {
0025 
0026 class IVolumeMaterial;
0027 
0028 class SympyStepper {
0029  public:
0030   /// Jacobian, Covariance and State definitions
0031   using Jacobian = BoundMatrix;
0032   using Covariance = BoundSquareMatrix;
0033   using BoundState = std::tuple<BoundTrackParameters, Jacobian, double>;
0034 
0035   struct Config {
0036     std::shared_ptr<const MagneticFieldProvider> bField;
0037   };
0038 
0039   struct Options : public StepperPlainOptions {
0040     bool doDense = true;
0041     double maxXOverX0Step = 1;
0042 
0043     Options(const GeometryContext& gctx, const MagneticFieldContext& mctx)
0044         : StepperPlainOptions(gctx, mctx) {}
0045 
0046     void setPlainOptions(const StepperPlainOptions& options) {
0047       static_cast<StepperPlainOptions&>(*this) = options;
0048     }
0049   };
0050 
0051   /// @brief State for track parameter propagation
0052   ///
0053   /// It contains the stepping information and is provided thread local
0054   /// by the propagator
0055   struct State {
0056     /// Constructor from the initial bound track parameters
0057     ///
0058     /// @param [in] optionsIn is the configuration of the stepper
0059     /// @param [in] fieldCacheIn is the cache object for the magnetic field
0060     ///
0061     /// @note the covariance matrix is copied when needed
0062     State(const Options& optionsIn, MagneticFieldProvider::Cache fieldCacheIn)
0063         : options(optionsIn), fieldCache(std::move(fieldCacheIn)) {}
0064 
0065     /// Configuration options for the stepper
0066     Options options;
0067 
0068     /// Internal free vector parameters
0069     FreeVector pars = FreeVector::Zero();
0070 
0071     /// Particle hypothesis
0072     ParticleHypothesis particleHypothesis = ParticleHypothesis::pion();
0073 
0074     /// Covariance matrix (and indicator)
0075     /// associated with the initial error on track parameters
0076     bool covTransport = false;
0077     /// Covariance matrix for error propagation
0078     Covariance cov = Covariance::Zero();
0079 
0080     /// The full jacobian of the transport entire transport
0081     Jacobian jacobian = Jacobian::Identity();
0082 
0083     /// Jacobian from local to the global frame
0084     BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
0085 
0086     /// Pure transport jacobian part from runge kutta integration
0087     FreeMatrix jacTransport = FreeMatrix::Identity();
0088 
0089     /// The propagation derivative
0090     FreeVector derivative = FreeVector::Zero();
0091 
0092     /// Accummulated path length state
0093     double pathAccumulated = 0.;
0094 
0095     /// Total number of performed steps
0096     std::size_t nSteps = 0;
0097 
0098     /// Totoal number of attempted steps
0099     std::size_t nStepTrials = 0;
0100 
0101     /// Adaptive step size of the runge-kutta integration
0102     ConstrainedStep stepSize;
0103 
0104     /// Last performed step (for overstep limit calculation)
0105     double previousStepSize = 0.;
0106 
0107     /// This caches the current magnetic field cell and stays
0108     /// (and interpolates) within it as long as this is valid.
0109     /// See step() code for details.
0110     MagneticFieldProvider::Cache fieldCache;
0111 
0112     /// Statistics of the stepper
0113     StepperStatistics statistics;
0114 
0115     /// Accumulator for material effects along the trajectory
0116     detail::MaterialEffectsAccumulator materialEffectsAccumulator;
0117   };
0118 
0119   /// Constructor requires knowledge of the detector's magnetic field
0120   /// @param bField The magnetic field provider
0121   explicit SympyStepper(std::shared_ptr<const MagneticFieldProvider> bField);
0122 
0123   /// @brief Constructor with configuration
0124   /// @param config The configuration of the stepper
0125   explicit SympyStepper(const Config& config);
0126 
0127   State makeState(const Options& options) const;
0128 
0129   void initialize(State& state, const BoundTrackParameters& par) const;
0130 
0131   void initialize(State& state, const BoundVector& boundParams,
0132                   const std::optional<BoundMatrix>& cov,
0133                   ParticleHypothesis particleHypothesis,
0134                   const Surface& surface) const;
0135 
0136   /// Get the field for the stepping, it checks first if the access is still
0137   /// within the Cell, and updates the cell if necessary.
0138   ///
0139   /// @param [in,out] state is the propagation state associated with the track
0140   ///                 the magnetic field cell is used (and potentially updated)
0141   /// @param [in] pos is the field position
0142   Result<Vector3> getField(State& state, const Vector3& pos) const {
0143     // get the field from the cell
0144     return m_bField->getField(pos, state.fieldCache);
0145   }
0146 
0147   /// Global particle position accessor
0148   ///
0149   /// @param state [in] The stepping state (thread-local cache)
0150   Vector3 position(const State& state) const {
0151     return state.pars.template segment<3>(eFreePos0);
0152   }
0153 
0154   /// Momentum direction accessor
0155   ///
0156   /// @param state [in] The stepping state (thread-local cache)
0157   Vector3 direction(const State& state) const {
0158     return state.pars.template segment<3>(eFreeDir0);
0159   }
0160 
0161   /// QoP direction accessor
0162   ///
0163   /// @param state [in] The stepping state (thread-local cache)
0164   double qOverP(const State& state) const { return state.pars[eFreeQOverP]; }
0165 
0166   /// Absolute momentum accessor
0167   ///
0168   /// @param state [in] The stepping state (thread-local cache)
0169   double absoluteMomentum(const State& state) const {
0170     return particleHypothesis(state).extractMomentum(qOverP(state));
0171   }
0172 
0173   /// Momentum accessor
0174   ///
0175   /// @param state [in] The stepping state (thread-local cache)
0176   Vector3 momentum(const State& state) const {
0177     return absoluteMomentum(state) * direction(state);
0178   }
0179 
0180   /// Charge access
0181   ///
0182   /// @param state [in] The stepping state (thread-local cache)
0183   double charge(const State& state) const {
0184     return particleHypothesis(state).extractCharge(qOverP(state));
0185   }
0186 
0187   /// Particle hypothesis
0188   ///
0189   /// @param state [in] The stepping state (thread-local cache)
0190   const ParticleHypothesis& particleHypothesis(const State& state) const {
0191     return state.particleHypothesis;
0192   }
0193 
0194   /// Time access
0195   ///
0196   /// @param state [in] The stepping state (thread-local cache)
0197   double time(const State& state) const { return state.pars[eFreeTime]; }
0198 
0199   /// Update surface status
0200   ///
0201   /// It checks the status to the reference surface & updates
0202   /// the step size accordingly
0203   ///
0204   /// @param [in,out] state The stepping state (thread-local cache)
0205   /// @param [in] surface The surface provided
0206   /// @param [in] index The surface intersection index
0207   /// @param [in] navDir The navigation direction
0208   /// @param [in] boundaryTolerance The boundary check for this status update
0209   /// @param [in] surfaceTolerance Surface tolerance used for intersection
0210   /// @param [in] stype The step size type to be set
0211   /// @param [in] logger A @c Logger instance
0212   IntersectionStatus updateSurfaceStatus(
0213       State& state, const Surface& surface, std::uint8_t index,
0214       Direction navDir, const BoundaryTolerance& boundaryTolerance,
0215       double surfaceTolerance, ConstrainedStep::Type stype,
0216       const Logger& logger = getDummyLogger()) const {
0217     return detail::updateSingleSurfaceStatus<SympyStepper>(
0218         *this, state, surface, index, navDir, boundaryTolerance,
0219         surfaceTolerance, stype, logger);
0220   }
0221 
0222   /// Update step size
0223   ///
0224   /// This method intersects the provided surface and update the navigation
0225   /// step estimation accordingly (hence it changes the state). It also
0226   /// returns the status of the intersection to trigger onSurface in case
0227   /// the surface is reached.
0228   ///
0229   /// @param state [in,out] The stepping state (thread-local cache)
0230   /// @param target [in] The NavigationTarget
0231   /// @param direction [in] The propagation direction
0232   /// @param stype [in] The step size type to be set
0233   void updateStepSize(State& state, const NavigationTarget& target,
0234                       Direction direction, ConstrainedStep::Type stype) const {
0235     (void)direction;
0236     double stepSize = target.pathLength();
0237     updateStepSize(state, stepSize, stype);
0238   }
0239 
0240   /// Update step size - explicitly with a double
0241   ///
0242   /// @param state [in,out] The stepping state (thread-local cache)
0243   /// @param stepSize [in] The step size value
0244   /// @param stype [in] The step size type to be set
0245   void updateStepSize(State& state, double stepSize,
0246                       ConstrainedStep::Type stype) const {
0247     state.previousStepSize = state.stepSize.value();
0248     state.stepSize.update(stepSize, stype);
0249   }
0250 
0251   /// Get the step size
0252   ///
0253   /// @param state [in] The stepping state (thread-local cache)
0254   /// @param stype [in] The step size type to be returned
0255   double getStepSize(const State& state, ConstrainedStep::Type stype) const {
0256     return state.stepSize.value(stype);
0257   }
0258 
0259   /// Release the Step size
0260   ///
0261   /// @param state [in,out] The stepping state (thread-local cache)
0262   /// @param [in] stype The step size type to be released
0263   void releaseStepSize(State& state, ConstrainedStep::Type stype) const {
0264     state.stepSize.release(stype);
0265   }
0266 
0267   /// Output the Step Size - single component
0268   ///
0269   /// @param state [in,out] The stepping state (thread-local cache)
0270   std::string outputStepSize(const State& state) const {
0271     return state.stepSize.toString();
0272   }
0273 
0274   /// Create and return the bound state at the current position
0275   ///
0276   /// @brief This transports (if necessary) the covariance
0277   /// to the surface and creates a bound state. It does not check
0278   /// if the transported state is at the surface, this needs to
0279   /// be guaranteed by the propagator
0280   ///
0281   /// @param [in] state State that will be presented as @c BoundState
0282   /// @param [in] surface The surface to which we bind the state
0283   /// @param [in] transportCov Flag steering covariance transport
0284   /// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound
0285   ///
0286   /// @return A bound state:
0287   ///   - the parameters at the surface
0288   ///   - the stepwise jacobian towards it (from last bound)
0289   ///   - and the path length (from start - for ordering)
0290   Result<BoundState> boundState(
0291       State& state, const Surface& surface, bool transportCov = true,
0292       const FreeToBoundCorrection& freeToBoundCorrection =
0293           FreeToBoundCorrection(false)) const;
0294 
0295   /// @brief If necessary fill additional members needed for curvilinearState
0296   ///
0297   /// Compute path length derivatives in case they have not been computed
0298   /// yet, which is the case if no step has been executed yet.
0299   ///
0300   /// @param [in, out] state State of the stepper
0301   /// @return true if nothing is missing after this call, false otherwise.
0302   bool prepareCurvilinearState(State& state) const;
0303 
0304   /// Create and return a curvilinear state at the current position
0305   ///
0306   /// @brief This transports (if necessary) the covariance
0307   /// to the current position and creates a curvilinear state.
0308   ///
0309   /// @param [in] state State that will be presented as @c CurvilinearState
0310   /// @param [in] transportCov Flag steering covariance transport
0311   ///
0312   /// @return A curvilinear state:
0313   ///   - the curvilinear parameters at given position
0314   ///   - the stepweise jacobian towards it (from last bound)
0315   ///   - and the path length (from start - for ordering)
0316   BoundState curvilinearState(State& state, bool transportCov = true) const;
0317 
0318   /// Method to update a stepper state to the some parameters
0319   ///
0320   /// @param [in,out] state State object that will be updated
0321   /// @param [in] freeParams Free parameters that will be written into @p state
0322   /// @param [in] boundParams Corresponding bound parameters used to update jacToGlobal in @p state
0323   /// @param [in] covariance The covariance that will be written into @p state
0324   /// @param [in] surface The surface used to update the jacToGlobal
0325   void update(State& state, const FreeVector& freeParams,
0326               const BoundVector& boundParams, const Covariance& covariance,
0327               const Surface& surface) const;
0328 
0329   /// Method to update the stepper state
0330   ///
0331   /// @param [in,out] state State object that will be updated
0332   /// @param [in] uposition the updated position
0333   /// @param [in] udirection the updated direction
0334   /// @param [in] qOverP the updated qOverP value
0335   /// @param [in] time the updated time value
0336   void update(State& state, const Vector3& uposition, const Vector3& udirection,
0337               double qOverP, double time) const;
0338 
0339   /// Method for on-demand transport of the covariance
0340   /// to a new curvilinear frame at current  position,
0341   /// or direction of the state
0342   ///
0343   /// @param [in,out] state State of the stepper
0344   void transportCovarianceToCurvilinear(State& state) const;
0345 
0346   /// Method for on-demand transport of the covariance
0347   /// to a new curvilinear frame at current position,
0348   /// or direction of the state
0349   ///
0350   /// @tparam surface_t the Surface type
0351   ///
0352   /// @param [in,out] state State of the stepper
0353   /// @param [in] surface is the surface to which the covariance is forwarded to
0354   /// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound
0355   /// @note no check is done if the position is actually on the surface
0356   void transportCovarianceToBound(
0357       State& state, const Surface& surface,
0358       const FreeToBoundCorrection& freeToBoundCorrection =
0359           FreeToBoundCorrection(false)) const;
0360 
0361   /// Perform a Runge-Kutta track parameter propagation step
0362   ///
0363   /// @param [in,out] state State of the stepper
0364   /// @param propDir is the direction of propagation
0365   /// @param material is the optional volume material we are stepping through.
0366   //         This is simply ignored if `nullptr`.
0367   /// @return the result of the step
0368   ///
0369   /// @note The state contains the desired step size. It can be negative during
0370   ///       backwards track propagation, and since we're using an adaptive
0371   ///       algorithm, it can be modified by the stepper class during
0372   ///       propagation.
0373   Result<double> step(State& state, Direction propDir,
0374                       const IVolumeMaterial* material) const;
0375 
0376   /// Method that reset the Jacobian to the Identity for when no bound state are
0377   /// available
0378   ///
0379   /// @param [in,out] state State of the stepper
0380   void setIdentityJacobian(State& state) const;
0381 
0382  protected:
0383   /// Magnetic field inside of the detector
0384   std::shared_ptr<const MagneticFieldProvider> m_bField;
0385 };
0386 
0387 template <>
0388 struct SupportsBoundParameters<SympyStepper> : public std::true_type {};
0389 
0390 }  // namespace Acts