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