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