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