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