|
|
|||
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
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|