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