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/ConstrainedStep.hpp"
0015 #include "Acts/Propagator/NavigationTarget.hpp"
0016 #include "Acts/Propagator/PropagatorError.hpp"
0017 #include "Acts/Propagator/StandardAborters.hpp"
0018 #include "Acts/Propagator/detail/LoopProtection.hpp"
0019 #include "Acts/Utilities/Intersection.hpp"
0020
0021 #include <concepts>
0022
0023 namespace Acts {
0024
0025 template <typename S, typename N>
0026 template <typename propagator_state_t>
0027 Result<void> Propagator<S, N>::propagate(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     return state.options.actorList.act(state, m_stepper, m_navigator, logger());
0042   }
0043
0044   auto getNextTarget = [&]() -> Result<NavigationTarget> {
0045     for (unsigned int i = 0; i < state.options.maxTargetSkipping; ++i) {
0046       NavigationTarget nextTarget = m_navigator.nextTarget(
0047           state.navigation, state.position, state.direction);
0048       if (nextTarget.isNone()) {
0049         return NavigationTarget::None();
0050       }
0051       IntersectionStatus preStepSurfaceStatus = m_stepper.updateSurfaceStatus(
0052           state.stepping, nextTarget.surface(), nextTarget.intersectionIndex(),
0053           state.options.direction, nextTarget.boundaryTolerance(),
0054           state.options.surfaceTolerance, ConstrainedStep::Type::Navigator,
0055           logger());
0056       if (preStepSurfaceStatus == IntersectionStatus::onSurface) {
0057         // This indicates a geometry overlap which is not handled by the
0058         // navigator, so we skip this target.
0059         // This can also happen in a well-behaved geometry with external
0060         // surfaces.
0061         ACTS_VERBOSE("Pre-step surface status is onSurface, skipping target "
0062                      << nextTarget.surface().geometryId());
0063         continue;
0064       }
0065       if (preStepSurfaceStatus == IntersectionStatus::reachable) {
0066         return nextTarget;
0067       }
0068     }
0069
0070     ACTS_DEBUG("getNextTarget failed to find a valid target surface after "
0071                << state.options.maxTargetSkipping << " attempts.");
0072     return Result<NavigationTarget>::failure(
0073         PropagatorError::NextTargetLimitReached);
0074   };
0075
0076   // priming error condition
0077   bool terminatedNormally = false;
0078
0079   // Pre-Stepping: target setting
0080   state.stage = PropagatorStage::preStep;
0081
0082   Result<NavigationTarget> nextTargetResult = getNextTarget();
0083   if (!nextTargetResult.ok()) {
0084     ACTS_DEBUG("Failed to get next target: "
0085                << nextTargetResult.error() << ": "
0086                << nextTargetResult.error().message());
0087     return nextTargetResult.error();
0088   }
0089   NavigationTarget nextTarget = *nextTargetResult;
0090
0091   ACTS_VERBOSE("Starting stepping loop.");
0092
0093   // Stepping loop
0094   for (; state.steps < state.options.maxSteps; ++state.steps) {
0095     // Perform a step
0096     Result<double> res =
0097         m_stepper.step(state.stepping, state.options.direction,
0098                        m_navigator.currentVolumeMaterial(state.navigation));
0099     if (!res.ok()) {
0100       ACTS_DEBUG("Step failed with " << res.error() << ": "
0101                                      << res.error().message());
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     Result<void> actResult =
0139         state.options.actorList.act(state, m_stepper, m_navigator, logger());
0140     if (!actResult.ok()) {
0141       ACTS_DEBUG("Actor call failed: " << actResult.error() << ": "
0142                                        << actResult.error().message());
0143       return actResult.error();
0144     }
0145
0146     if (state.options.actorList.checkAbort(state, m_stepper, m_navigator,
0147                                            logger())) {
0148       terminatedNormally = true;
0149       break;
0150     }
0151
0152     // Update the position and direction because actors might have changed it
0153     state.position = m_stepper.position(state.stepping);
0154     state.direction =
0155         state.options.direction * m_stepper.direction(state.stepping);
0156
0157     // Pre-Stepping: target setting
0158     state.stage = PropagatorStage::preStep;
0159
0160     if (!nextTarget.isNone() &&
0161         !m_navigator.checkTargetValid(state.navigation, state.position,
0162                                       state.direction)) {
0163       ACTS_VERBOSE("Target is not valid anymore.");
0164       nextTarget = NavigationTarget::None();
0165     }
0166
0167     if (nextTarget.isNone()) {
0168       // navigator step constraint is not valid anymore
0169       m_stepper.releaseStepSize(state.stepping,
0170                                 ConstrainedStep::Type::Navigator);
0171
0172       nextTargetResult = getNextTarget();
0173       if (!nextTargetResult.ok()) {
0174         ACTS_DEBUG("Failed to get next target: "
0175                    << nextTargetResult.error() << ": "
0176                    << nextTargetResult.error().message());
0177         return nextTargetResult.error();
0178       }
0179       nextTarget = *nextTargetResult;
0180     }
0181   }  // end of stepping loop
0182
0183   // check if we didn't terminate normally via aborters
0184   if (!terminatedNormally) {
0185     ACTS_DEBUG("Propagation reached the step count limit of "
0186                << state.options.maxSteps << " (did " << state.steps
0187                << " steps)");
0188     return PropagatorError::StepCountLimitReached;
0189   }
0190
0191   ACTS_VERBOSE("Stepping loop done.");
0192
0193   state.stage = PropagatorStage::postPropagation;
0194
0195   // Post-stepping call to the actor list
0196   if (auto postPropagationResult =
0197           state.options.actorList.act(state, m_stepper, m_navigator, logger());
0198       !postPropagationResult.ok()) {
0199     ACTS_DEBUG("Post-propagation actor call failed: "
0200                << postPropagationResult.error() << ": "
0201                << postPropagationResult.error().message());
0202     return postPropagationResult.error();
0203   }
0204   return Result<void>::success();
0205 }
0206
0207 template <typename S, typename N>
0208 template <typename parameters_t, typename propagator_options_t,
0209           typename path_aborter_t>
0210 auto Propagator<S, N>::propagate(const parameters_t& start,
0211                                  const propagator_options_t& options,
0212                                  bool createFinalParameters) const
0213     -> Result<ResultType<propagator_options_t>> {
0214   static_assert(std::copy_constructible<StepperBoundTrackParameters>,
0215                 "return track parameter type must be copy-constructible");
0216
0217   auto state = makeState<propagator_options_t, path_aborter_t>(options);
0218
0219   auto initRes =
0220       initialize<decltype(state), parameters_t, path_aborter_t>(state, start);
0221   if (!initRes.ok()) {
0222     ACTS_DEBUG("Initialization failed: " << initRes.error() << ": "
0223                                          << initRes.error().message());
0224     return initRes.error();
0225   }
0226
0227   // Perform the actual propagation
0228   auto propagationResult = propagate(state);
0229
0230   return makeResult(std::move(state), propagationResult, options,
0231                     createFinalParameters);
0232 }
0233
0234 template <typename S, typename N>
0235 template <typename parameters_t, typename propagator_options_t,
0236           typename target_aborter_t, typename path_aborter_t>
0237 auto Propagator<S, N>::propagate(const parameters_t& start,
0238                                  const Surface& target,
0239                                  const propagator_options_t& options) const
0240     -> Result<ResultType<propagator_options_t>> {
0241   static_assert(BoundTrackParametersConcept<parameters_t>,
0242                 "Parameters do not fulfill bound parameters concept.");
0243
0244   auto state =
0245       makeState<propagator_options_t, target_aborter_t, path_aborter_t>(
0246           target, options);
0247
0248   auto initRes =
0249       initialize<decltype(state), parameters_t, path_aborter_t>(state, start);
0250   if (!initRes.ok()) {
0251     ACTS_DEBUG("Initialization failed: " << initRes.error() << ": "
0252                                          << initRes.error().message());
0253     return initRes.error();
0254   }
0255
0256   // Perform the actual propagation
0257   auto propagationResult = propagate(state);
0258
0259   return makeResult(std::move(state), propagationResult, target, options);
0260 }
0261
0262 template <typename S, typename N>
0263 template <typename propagator_options_t, typename path_aborter_t>
0264 auto Propagator<S, N>::makeState(const propagator_options_t& options) const {
0265   // Expand the actor list with a path aborter
0266   path_aborter_t pathAborter;
0267   pathAborter.internalLimit = options.pathLimit;
0268
0269   auto actorList = options.actorList.append(pathAborter);
0270
0271   // Create the extended options and declare their type
0272   auto eOptions = options.extend(actorList);
0273
0274   using OptionsType = decltype(eOptions);
0275   using StateType = State<OptionsType>;
0276
0277   StateType state{eOptions, m_stepper.makeState(eOptions.stepping),
0278                   m_navigator.makeState(eOptions.navigation)};
0279
0280   return state;
0281 }
0282
0283 template <typename S, typename N>
0284 template <typename propagator_options_t, typename target_aborter_t,
0285           typename path_aborter_t>
0286 auto Propagator<S, N>::makeState(const Surface& target,
0287                                  const propagator_options_t& options) const {
0288   // Expand the actor list with a target and path aborter
0289   target_aborter_t targetAborter;
0290   targetAborter.surface = &target;
0291   path_aborter_t pathAborter;
0292   pathAborter.internalLimit = options.pathLimit;
0293
0294   auto actorList = options.actorList.append(targetAborter, pathAborter);
0295
0296   // Create the extended options and declare their type
0297   auto eOptions = options.extend(actorList);
0298   eOptions.navigation.targetSurface = &target;
0299
0300   using OptionsType = decltype(eOptions);
0301   using StateType = State<OptionsType>;
0302
0303   StateType state{eOptions, m_stepper.makeState(eOptions.stepping),
0304                   m_navigator.makeState(eOptions.navigation)};
0305
0306   return state;
0307 }
0308
0309 template <typename S, typename N>
0310 template <typename propagator_state_t, typename parameters_t,
0311           typename path_aborter_t>
0312 Result<void> Propagator<S, N>::initialize(propagator_state_t& state,
0313                                           const parameters_t& start) const {
0314   static_assert(BoundTrackParametersConcept<parameters_t>,
0315                 "Parameters do not fulfill bound parameters concept.");
0316
0317   m_stepper.initialize(state.stepping, start.toBound());
0318
0319   state.position = m_stepper.position(state.stepping);
0320   state.direction =
0321       state.options.direction * m_stepper.direction(state.stepping);
0322
0323   state.navigation.options.startSurface = &start.referenceSurface();
0324
0325   // Navigator initialize state call
0326   auto navInitRes =
0327       m_navigator.initialize(state.navigation, state.position, state.direction,
0328                              state.options.direction);
0329   if (!navInitRes.ok()) {
0330     ACTS_DEBUG("Navigator initialization failed: "
0331                << navInitRes.error() << ": " << navInitRes.error().message());
0332     return navInitRes.error();
0333   }
0334
0335   // Apply the loop protection - it resets the internal path limit
0336   detail::setupLoopProtection(
0337       state, m_stepper, state.options.actorList.template get<path_aborter_t>(),
0338       false, logger());
0339
0340   return Result<void>::success();
0341 }
0342
0343 template <typename S, typename N>
0344 template <typename propagator_state_t, typename propagator_options_t>
0345 auto Propagator<S, N>::makeResult(propagator_state_t state,
0346                                   Result<void> propagationResult,
0347                                   const propagator_options_t& /*options*/,
0348                                   bool createFinalParameters) const
0349     -> Result<ResultType<propagator_options_t>> {
0350   // Type of the full propagation result, including output from actors
0351   using ThisResultType = ResultType<propagator_options_t>;
0352
0353   if (!propagationResult.ok()) {
0354     ACTS_DEBUG("Propagation failed: " << propagationResult.error() << ": "
0355                                       << propagationResult.error().message());
0356     return propagationResult.error();
0357   }
0358
0359   ThisResultType result{};
0360   moveStateToResult(state, result);
0361
0362   if (const Surface* currentSurface =
0363           m_navigator.currentSurface(state.navigation);
0364       createFinalParameters && currentSurface == nullptr) {
0365     if (!m_stepper.prepareCurvilinearState(state.stepping)) {
0366       // information to compute curvilinearState is incomplete.
0367       ACTS_DEBUG("Failed to prepare curvilinear state.");
0368       return PropagatorError::Failure;
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 = m_stepper.boundState(state.stepping, *currentSurface);
0381     if (!boundState.ok()) {
0382       ACTS_DEBUG("Failed to get bound state at current surface: "
0383                  << boundState.error() << ": " << boundState.error().message());
0384       return boundState.error();
0385     }
0386     // Fill the end parameters
0387     result.endParameters = std::get<StepperBoundTrackParameters>(*boundState);
0388     // Only fill the transport jacobian when covariance transport was done
0389     if (state.stepping.covTransport) {
0390       result.transportJacobian = std::get<Jacobian>(*boundState);
0391     }
0392   }
0393
0394   return Result<ThisResultType>::success(std::move(result));
0395 }
0396
0397 template <typename S, typename N>
0398 template <typename propagator_state_t, typename propagator_options_t>
0399 auto Propagator<S, N>::makeResult(propagator_state_t state,
0400                                   Result<void> propagationResult,
0401                                   const Surface& target,
0402                                   const propagator_options_t& /*options*/) const
0403     -> Result<ResultType<propagator_options_t>> {
0404   // Type of the full propagation result, including output from actors
0405   using ThisResultType = ResultType<propagator_options_t>;
0406
0407   if (!propagationResult.ok()) {
0408     ACTS_DEBUG("Propagation failed: " << propagationResult.error() << ": "
0409                                       << propagationResult.error().message());
0410     return propagationResult.error();
0411   }
0412
0413   // Get the bound state at the target surface
0414   auto bsRes = m_stepper.boundState(state.stepping, target);
0415   if (!bsRes.ok()) {
0416     // This likely indicates that the position is outside the surface bounds
0417     ACTS_DEBUG("Failed to get bound state at target surface: "
0418                << bsRes.error() << ": " << bsRes.error().message());
0419     return bsRes.error();
0420   }
0421   const auto& bs = *bsRes;
0422
0423   ThisResultType result{};
0424   moveStateToResult(state, result);
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<ThisResultType>::success(std::move(result));
0433 }
0434
0435 template <typename S, typename N>
0436 template <typename propagator_state_t, typename propagator_result_t>
0437 void Propagator<S, N>::moveStateToResult(propagator_state_t& state,
0438                                          propagator_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 Result<BoundTrackParameters>
0450 detail::BasePropagatorHelper<derived_t>::propagateToSurface(
0451     const BoundTrackParameters& start, const Surface& target,
0452     const Options& options) const {
0453   using DerivedOptions = typename derived_t::template Options<>;
0454   using DerivedResult = typename derived_t::template ResultType<DerivedOptions>;
0455
0456   DerivedOptions derivedOptions(options);
0457
0458   // dummy initialization
0459   Result<DerivedResult> res =
0460       Result<DerivedResult>::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 }
0485
0486 }  // namespace Acts