Acts/Propagator/Propagator.ipp

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/Propagator/Propagator.hpp"
0012
0013 #include "Acts/EventData/TrackParametersConcept.hpp"
0014 #include "Acts/Propagator/ActorList.hpp"
0015 #include "Acts/Propagator/ConstrainedStep.hpp"
0016 #include "Acts/Propagator/NavigationTarget.hpp"
0017 #include "Acts/Propagator/PropagatorError.hpp"
0018 #include "Acts/Propagator/StandardAborters.hpp"
0019 #include "Acts/Propagator/detail/LoopProtection.hpp"
0020 #include "Acts/Utilities/Intersection.hpp"
0021
0022 #include <concepts>
0023
0024 template <typename S, typename N>
0025 template <typename propagator_state_t>
0026 Acts::Result<void> Acts::Propagator<S, N>::propagate(
0027     propagator_state_t& state) const {
0028   ACTS_VERBOSE("Entering propagation.");
0029
0030   state.stage = PropagatorStage::prePropagation;
0031
0032   // Pre-Propagation: call to the actor list, abort condition check
0033   state.options.actorList.act(state, m_stepper, m_navigator, logger());
0034
0035   if (state.options.actorList.checkAbort(state, m_stepper, m_navigator,
0036                                          logger())) {
0037     ACTS_VERBOSE("Propagation terminated without going into stepping loop.");
0038
0039     state.stage = PropagatorStage::postPropagation;
0040
0041     state.options.actorList.act(state, m_stepper, m_navigator, logger());
0042
0043     return Result<void>::success();
0044   }
0045
0046   auto getNextTarget = [&]() -> Result<NavigationTarget> {
0047     for (unsigned int i = 0; i < state.options.maxTargetSkipping; ++i) {
0048       NavigationTarget nextTarget = m_navigator.nextTarget(
0049           state.navigation, state.position, state.direction);
0050       if (nextTarget.isNone()) {
0051         return NavigationTarget::None();
0052       }
0053       IntersectionStatus preStepSurfaceStatus = m_stepper.updateSurfaceStatus(
0054           state.stepping, nextTarget.surface(), nextTarget.intersectionIndex(),
0055           state.options.direction, nextTarget.boundaryTolerance(),
0056           state.options.surfaceTolerance, ConstrainedStep::Type::Navigator,
0057           logger());
0058       if (preStepSurfaceStatus == IntersectionStatus::onSurface) {
0059         // This indicates a geometry overlap which is not handled by the
0060         // navigator, so we skip this target.
0061         // This can also happen in a well-behaved geometry with external
0062         // surfaces.
0063         ACTS_VERBOSE("Pre-step surface status is onSurface, skipping target "
0064                      << nextTarget.surface().geometryId());
0065         continue;
0066       }
0067       if (preStepSurfaceStatus == IntersectionStatus::reachable) {
0068         return nextTarget;
0069       }
0070     }
0071
0072     ACTS_ERROR("getNextTarget failed to find a valid target surface after "
0073                << state.options.maxTargetSkipping << " attempts.");
0074     return Result<NavigationTarget>::failure(
0075         PropagatorError::NextTargetLimitReached);
0076   };
0077
0078   // priming error condition
0079   bool terminatedNormally = false;
0080
0081   // Pre-Stepping: target setting
0082   state.stage = PropagatorStage::preStep;
0083
0084   Result<NavigationTarget> nextTargetResult = getNextTarget();
0085   if (!nextTargetResult.ok()) {
0086     return nextTargetResult.error();
0087   }
0088   NavigationTarget nextTarget = *nextTargetResult;
0089
0090   ACTS_VERBOSE("Starting stepping loop.");
0091
0092   // Stepping loop
0093   for (; state.steps < state.options.maxSteps; ++state.steps) {
0094     // Perform a step
0095     Result<double> res =
0096         m_stepper.step(state.stepping, state.options.direction,
0097                        m_navigator.currentVolumeMaterial(state.navigation));
0098     if (!res.ok()) {
0099       ACTS_ERROR("Step failed with " << res.error() << ": "
0100                                      << res.error().message());
0101       // pass error to caller
0102       return res.error();
0103     }
0104     // Accumulate the path length
0105     state.pathLength += *res;
0106     // Update the position and direction
0107     state.position = m_stepper.position(state.stepping);
0108     state.direction =
0109         state.options.direction * m_stepper.direction(state.stepping);
0110
0111     ACTS_VERBOSE("Step with size " << *res << " performed. We are now at "
0112                                    << state.position.transpose()
0113                                    << " with direction "
0114                                    << state.direction.transpose());
0115
0116     // release actor and aborter constrains after step was performed
0117     m_stepper.releaseStepSize(state.stepping, ConstrainedStep::Type::Navigator);
0118     m_stepper.releaseStepSize(state.stepping, ConstrainedStep::Type::Actor);
0119
0120     // Post-stepping: check target status, call actors, check abort conditions
0121     state.stage = PropagatorStage::postStep;
0122
0123     if (!nextTarget.isNone()) {
0124       IntersectionStatus postStepSurfaceStatus = m_stepper.updateSurfaceStatus(
0125           state.stepping, nextTarget.surface(), nextTarget.intersectionIndex(),
0126           state.options.direction, nextTarget.boundaryTolerance(),
0127           state.options.surfaceTolerance, ConstrainedStep::Type::Navigator,
0128           logger());
0129       if (postStepSurfaceStatus == IntersectionStatus::onSurface) {
0130         m_navigator.handleSurfaceReached(state.navigation, state.position,
0131                                          state.direction, nextTarget.surface());
0132       }
0133       if (postStepSurfaceStatus != IntersectionStatus::reachable) {
0134         nextTarget = NavigationTarget::None();
0135       }
0136     }
0137
0138     state.options.actorList.act(state, m_stepper, m_navigator, logger());
0139
0140     if (state.options.actorList.checkAbort(state, m_stepper, m_navigator,
0141                                            logger())) {
0142       terminatedNormally = true;
0143       break;
0144     }
0145
0146     // Update the position and direction because actors might have changed it
0147     state.position = m_stepper.position(state.stepping);
0148     state.direction =
0149         state.options.direction * m_stepper.direction(state.stepping);
0150
0151     // Pre-Stepping: target setting
0152     state.stage = PropagatorStage::preStep;
0153
0154     if (!nextTarget.isNone() &&
0155         !m_navigator.checkTargetValid(state.navigation, state.position,
0156                                       state.direction)) {
0157       ACTS_VERBOSE("Target is not valid anymore.");
0158       nextTarget = NavigationTarget::None();
0159     }
0160
0161     if (nextTarget.isNone()) {
0162       // navigator step constraint is not valid anymore
0163       m_stepper.releaseStepSize(state.stepping,
0164                                 ConstrainedStep::Type::Navigator);
0165
0166       nextTargetResult = getNextTarget();
0167       if (!nextTargetResult.ok()) {
0168         return nextTargetResult.error();
0169       }
0170       nextTarget = *nextTargetResult;
0171     }
0172   }  // end of stepping loop
0173
0174   // check if we didn't terminate normally via aborters
0175   if (!terminatedNormally) {
0176     ACTS_ERROR("Propagation reached the step count limit of "
0177                << state.options.maxSteps << " (did " << state.steps
0178                << " steps)");
0179     return PropagatorError::StepCountLimitReached;
0180   }
0181
0182   ACTS_VERBOSE("Stepping loop done.");
0183
0184   state.stage = PropagatorStage::postPropagation;
0185
0186   // Post-stepping call to the actor list
0187   state.options.actorList.act(state, m_stepper, m_navigator, logger());
0188
0189   return Result<void>::success();
0190 }
0191
0192 template <typename S, typename N>
0193 template <typename parameters_t, typename propagator_options_t,
0194           typename path_aborter_t>
0195 auto Acts::Propagator<S, N>::propagate(const parameters_t& start,
0196                                        const propagator_options_t& options,
0197                                        bool createFinalParameters) const
0198     -> Result<
0199         actor_list_t_result_t<StepperBoundTrackParameters,
0200                               typename propagator_options_t::actor_list_type>> {
0201   static_assert(std::copy_constructible<StepperBoundTrackParameters>,
0202                 "return track parameter type must be copy-constructible");
0203
0204   auto state = makeState<propagator_options_t, path_aborter_t>(options);
0205
0206   auto initRes =
0207       initialize<decltype(state), parameters_t, path_aborter_t>(state, start);
0208   if (!initRes.ok()) {
0209     return initRes.error();
0210   }
0211
0212   // Perform the actual propagation
0213   auto propagationResult = propagate(state);
0214
0215   return makeResult(std::move(state), propagationResult, options,
0216                     createFinalParameters);
0217 }
0218
0219 template <typename S, typename N>
0220 template <typename parameters_t, typename propagator_options_t,
0221           typename target_aborter_t, typename path_aborter_t>
0222 auto Acts::Propagator<S, N>::propagate(
0223     const parameters_t& start, const Surface& target,
0224     const propagator_options_t& options) const
0225     -> Result<
0226         actor_list_t_result_t<StepperBoundTrackParameters,
0227                               typename propagator_options_t::actor_list_type>> {
0228   static_assert(BoundTrackParametersConcept<parameters_t>,
0229                 "Parameters do not fulfill bound parameters concept.");
0230
0231   auto state =
0232       makeState<propagator_options_t, target_aborter_t, path_aborter_t>(
0233           target, options);
0234
0235   auto initRes =
0236       initialize<decltype(state), parameters_t, path_aborter_t>(state, start);
0237   if (!initRes.ok()) {
0238     return initRes.error();
0239   }
0240
0241   // Perform the actual propagation
0242   auto propagationResult = propagate(state);
0243
0244   return makeResult(std::move(state), propagationResult, target, options);
0245 }
0246
0247 template <typename S, typename N>
0248 template <typename propagator_options_t, typename path_aborter_t>
0249 auto Acts::Propagator<S, N>::makeState(
0250     const propagator_options_t& options) const {
0251   // Type of track parameters produced by the propagation
0252   using ReturnParameterType = StepperBoundTrackParameters;
0253
0254   static_assert(std::copy_constructible<ReturnParameterType>,
0255                 "return track parameter type must be copy-constructible");
0256
0257   // Expand the actor list with a path aborter
0258   path_aborter_t pathAborter;
0259   pathAborter.internalLimit = options.pathLimit;
0260
0261   auto actorList = options.actorList.append(pathAborter);
0262
0263   // Create the extended options and declare their type
0264   auto eOptions = options.extend(actorList);
0265
0266   using OptionsType = decltype(eOptions);
0267   using StateType =
0268       actor_list_t_state_t<OptionsType,
0269                            typename propagator_options_t::actor_list_type>;
0270
0271   StateType state{eOptions, m_stepper.makeState(eOptions.stepping),
0272                   m_navigator.makeState(eOptions.navigation)};
0273
0274   return state;
0275 }
0276
0277 template <typename S, typename N>
0278 template <typename propagator_options_t, typename target_aborter_t,
0279           typename path_aborter_t>
0280 auto Acts::Propagator<S, N>::makeState(
0281     const Surface& target, const propagator_options_t& options) const {
0282   // Expand the actor list with a target and path aborter
0283   target_aborter_t targetAborter;
0284   targetAborter.surface = &target;
0285   path_aborter_t pathAborter;
0286   pathAborter.internalLimit = options.pathLimit;
0287
0288   auto actorList = options.actorList.append(targetAborter, pathAborter);
0289
0290   // Create the extended options and declare their type
0291   auto eOptions = options.extend(actorList);
0292   eOptions.navigation.targetSurface = &target;
0293
0294   using OptionsType = decltype(eOptions);
0295   using StateType =
0296       actor_list_t_state_t<OptionsType,
0297                            typename propagator_options_t::actor_list_type>;
0298
0299   StateType state{eOptions, m_stepper.makeState(eOptions.stepping),
0300                   m_navigator.makeState(eOptions.navigation)};
0301
0302   return state;
0303 }
0304
0305 template <typename S, typename N>
0306 template <typename propagator_state_t, typename parameters_t,
0307           typename path_aborter_t>
0308 Acts::Result<void> Acts::Propagator<S, N>::initialize(
0309     propagator_state_t& state, const parameters_t& start) const {
0310   static_assert(BoundTrackParametersConcept<parameters_t>,
0311                 "Parameters do not fulfill bound parameters concept.");
0312
0313   m_stepper.initialize(state.stepping, start.toBound());
0314
0315   state.position = m_stepper.position(state.stepping);
0316   state.direction =
0317       state.options.direction * m_stepper.direction(state.stepping);
0318
0319   state.navigation.options.startSurface = &start.referenceSurface();
0320
0321   // Navigator initialize state call
0322   auto navInitRes =
0323       m_navigator.initialize(state.navigation, state.position, state.direction,
0324                              state.options.direction);
0325   if (!navInitRes.ok()) {
0326     return navInitRes.error();
0327   }
0328
0329   // Apply the loop protection - it resets the internal path limit
0330   detail::setupLoopProtection(
0331       state, m_stepper, state.options.actorList.template get<path_aborter_t>(),
0332       false, logger());
0333
0334   return Result<void>::success();
0335 }
0336
0337 template <typename S, typename N>
0338 template <typename propagator_state_t, typename propagator_options_t>
0339 auto Acts::Propagator<S, N>::makeResult(propagator_state_t state,
0340                                         Result<void> propagationResult,
0341                                         const propagator_options_t& /*options*/,
0342                                         bool createFinalParameters) const
0343     -> Result<
0344         actor_list_t_result_t<StepperBoundTrackParameters,
0345                               typename propagator_options_t::actor_list_type>> {
0346   // Type of track parameters produced by the propagation
0347   using ReturnParameterType = StepperBoundTrackParameters;
0348
0349   static_assert(std::copy_constructible<ReturnParameterType>,
0350                 "return track parameter type must be copy-constructible");
0351
0352   // Type of the full propagation result, including output from actors
0353   using ResultType =
0354       actor_list_t_result_t<ReturnParameterType,
0355                             typename propagator_options_t::actor_list_type>;
0356
0357   if (!propagationResult.ok()) {
0358     return propagationResult.error();
0359   }
0360
0361   ResultType result{};
0362   moveStateToResult(state, result);
0363
0364   const Surface* currentSurface = m_navigator.currentSurface(state.navigation);
0365   if (createFinalParameters && currentSurface == nullptr) {
0366     if (!m_stepper.prepareCurvilinearState(state.stepping)) {
0367       // information to compute curvilinearState is incomplete.
0368       return propagationResult.error();
0369     }
0370     /// Convert into return type and fill the result object
0371     auto curvState = m_stepper.curvilinearState(state.stepping);
0372     // Fill the end parameters
0373     result.endParameters = std::get<StepperBoundTrackParameters>(curvState);
0374     // Only fill the transport jacobian when covariance transport was done
0375     if (state.stepping.covTransport) {
0376       result.transportJacobian = std::get<Jacobian>(curvState);
0377     }
0378   } else if (createFinalParameters && currentSurface != nullptr) {
0379     // We are at a surface, so we need to compute the bound state
0380     auto boundState =
0381         m_stepper.boundState(state.stepping, *currentSurface).value();
0382     // Fill the end parameters
0383     result.endParameters = std::get<StepperBoundTrackParameters>(boundState);
0384     // Only fill the transport jacobian when covariance transport was done
0385     if (state.stepping.covTransport) {
0386       result.transportJacobian = std::get<Jacobian>(boundState);
0387     }
0388   }
0389
0390   return Result<ResultType>::success(std::move(result));
0391 }
0392
0393 template <typename S, typename N>
0394 template <typename propagator_state_t, typename propagator_options_t>
0395 auto Acts::Propagator<S, N>::makeResult(
0396     propagator_state_t state, Result<void> propagationResult,
0397     const Surface& target, const propagator_options_t& /*options*/) const
0398     -> Result<
0399         actor_list_t_result_t<StepperBoundTrackParameters,
0400                               typename propagator_options_t::actor_list_type>> {
0401   // Type of track parameters produced at the end of the propagation
0402   using ReturnParameterType = StepperBoundTrackParameters;
0403
0404   static_assert(std::copy_constructible<ReturnParameterType>,
0405                 "return track parameter type must be copy-constructible");
0406
0407   // Type of the full propagation result, including output from actors
0408   using ResultType =
0409       actor_list_t_result_t<ReturnParameterType,
0410                             typename propagator_options_t::actor_list_type>;
0411
0412   if (!propagationResult.ok()) {
0413     return propagationResult.error();
0414   }
0415
0416   ResultType result{};
0417   moveStateToResult(state, result);
0418
0419   // Compute the final results and mark the propagation as successful
0420   auto bsRes = m_stepper.boundState(state.stepping, target);
0421   if (!bsRes.ok()) {
0422     return bsRes.error();
0423   }
0424   const auto& bs = *bsRes;
0425
0426   // Fill the end parameters
0427   result.endParameters = std::get<StepperBoundTrackParameters>(bs);
0428   // Only fill the transport jacobian when covariance transport was done
0429   if (state.stepping.covTransport) {
0430     result.transportJacobian = std::get<Jacobian>(bs);
0431   }
0432   return Result<ResultType>::success(std::move(result));
0433 }
0434
0435 template <typename S, typename N>
0436 template <typename propagator_state_t, typename result_t>
0437 void Acts::Propagator<S, N>::moveStateToResult(propagator_state_t& state,
0438                                                result_t& result) const {
0439   result.tuple() = std::move(state.tuple());
0440
0441   result.steps = state.steps;
0442   result.pathLength = state.pathLength;
0443
0444   result.statistics.stepping = state.stepping.statistics;
0445   result.statistics.navigation = state.navigation.statistics;
0446 }
0447
0448 template <typename derived_t>
0449 Acts::Result<Acts::BoundTrackParameters>
0450 Acts::detail::BasePropagatorHelper<derived_t>::propagateToSurface(
0451     const BoundTrackParameters& start, const Surface& target,
0452     const Options& options) const {
0453   using ResultType = Result<typename derived_t::template actor_list_t_result_t<
0454       BoundTrackParameters, ActorList<>>>;
0455   using DerivedOptions = typename derived_t::template Options<>;
0456
0457   DerivedOptions derivedOptions(options);
0458
0459   // dummy initialization
0460   ResultType res = ResultType::failure(PropagatorError::Failure);
0461
0462   // Due to the geometry of the perigee surface the overstepping tolerance
0463   // is sometimes not met.
0464   if (target.type() == Surface::SurfaceType::Perigee) {
0465     res = static_cast<const derived_t*>(this)
0466               ->template propagate<BoundTrackParameters, DerivedOptions,
0467                                    ForcedSurfaceReached, PathLimitReached>(
0468                   start, target, derivedOptions);
0469   } else {
0470     res = static_cast<const derived_t*>(this)
0471               ->template propagate<BoundTrackParameters, DerivedOptions,
0472                                    SurfaceReached, PathLimitReached>(
0473                   start, target, derivedOptions);
0474   }
0475
0476   if (!res.ok()) {
0477     return res.error();
0478   }
0479
0480   // Without errors we can expect a valid endParameters when propagating to a
0481   // target surface
0482   assert((*res).endParameters);
0483   return std::move((*res).endParameters.value());
0484 }