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