fix: Implement DirectNavigator direction handling (acts-project#3702)

Our current implementation of the `DirectNavigator` does not handle the propagation direction which is essential if we for example reverse direction to smooth track in the Kalman filter.

Apart from that I added a search mechanism in the initialization for the starting surface.

---

This pull request to `Core/include/Acts/Propagator/DirectNavigator.hpp` includes changes to improve the navigation logic in the `DirectNavigator` class. The most important changes include adding new utility methods for surface navigation, updating the surface index handling, and enhancing logging for better debugging.

Enhancements to navigation logic:

* Added utility methods `navSurface`, `nextSurface`, `endOfSurfaces`, and `remainingSurfaces` to the `DirectNavigator` class to streamline surface navigation. (`Core/include/Acts/Propagator/DirectNavigator.hpp`) [[1]](diffhunk://#diff-56afaecc4b563cfe416f8615f35dd5ac21f706869dd2c164dd326087b3edd655R77-R100) [[2]](diffhunk://#diff-56afaecc4b563cfe416f8615f35dd5ac21f706869dd2c164dd326087b3edd655L138-R192) [[3]](diffhunk://#diff-56afaecc4b563cfe416f8615f35dd5ac21f706869dd2c164dd326087b3edd655L166-R215) [[4]](diffhunk://#diff-56afaecc4b563cfe416f8615f35dd5ac21f706869dd2c164dd326087b3edd655L186-R229) [[5]](diffhunk://#diff-56afaecc4b563cfe416f8615f35dd5ac21f706869dd2c164dd326087b3edd655L227-R279) [[6]](diffhunk://#diff-56afaecc4b563cfe416f8615f35dd5ac21f706869dd2c164dd326087b3edd655L255-R299)

Updates to surface index handling:

* Changed the type of `surfaceIndex` from `std::size_t` to `int` to support bidirectional navigation. (`Core/include/Acts/Propagator/DirectNavigator.hpp`)

Logging improvements:

* Enhanced logging in the `initialize` method to provide more detailed information about the start surface and surface sequence. (`Core/include/Acts/Propagator/DirectNavigator.hpp`)

Core/include/Acts/Propagator/DirectNavigator.hpp

@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include "Acts/Definitions/Direction.hpp"
 #include "Acts/Geometry/BoundarySurfaceT.hpp"
 #include "Acts/Geometry/Layer.hpp"
 #include "Acts/Geometry/TrackingGeometry.hpp"
@@ -17,6 +18,7 @@
 #include "Acts/Surfaces/BoundaryTolerance.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Intersection.hpp"
+#include "Acts/Utilities/Logger.hpp"
 
 #include <algorithm>
 #include <iterator>
@@ -63,7 +65,7 @@ class DirectNavigator {
     Options options;
 
     /// Index of the next surface to try
-    std::size_t surfaceIndex = 0;
+    int surfaceIndex = 0;
 
     /// Navigation state - external interface: the current surface
     const Surface* currentSurface = nullptr;
@@ -72,6 +74,35 @@ class DirectNavigator {
     bool targetReached = false;
     /// Navigation state - external interface: a break has been detected
     bool navigationBreak = false;
+
+    const Surface* navSurface() const {
+      return options.surfaces.at(surfaceIndex);
+    }
+
+    void nextSurface(Direction direction) {
+      if (direction == Direction::Forward) {
+        ++surfaceIndex;
+      } else {
+        --surfaceIndex;
+      }
+    }
+
+    bool endOfSurfaces() const {
+      return surfaceIndex < 0 ||
+             surfaceIndex >= static_cast<int>(options.surfaces.size());
+    }
+
+    int remainingSurfaces(Direction direction) const {
+      if (direction == Direction::Forward) {
+        return options.surfaces.size() - surfaceIndex;
+      }
+      return surfaceIndex + 1;
+    }
+
+    void resetSurfaceIndex(Direction direction) {
+      surfaceIndex =
+          direction == Direction::Forward ? 0 : options.surfaces.size() - 1;
+    }
   };
 
   explicit DirectNavigator(std::unique_ptr<const Logger> _logger =
@@ -135,17 +166,41 @@ class DirectNavigator {
   void initialize(propagator_state_t& state,
                   const stepper_t& /*stepper*/) const {
     ACTS_VERBOSE("Initialize. Surface sequence for navigation:");
-    for (auto surface : state.navigation.options.surfaces) {
+    for (const Surface* surface : state.navigation.options.surfaces) {
       ACTS_VERBOSE(surface->geometryId()
                    << " - " << surface->center(state.geoContext).transpose());
     }
 
     // We set the current surface to the start surface
     state.navigation.currentSurface = state.navigation.options.startSurface;
-    if (state.navigation.currentSurface) {
+    if (state.navigation.currentSurface != nullptr) {
       ACTS_VERBOSE("Current surface set to start surface "
                    << state.navigation.currentSurface->geometryId());
+    } else {
+      ACTS_VERBOSE("Current surface set to nullptr");
+    }
+
+    // Reset the surface index
+    state.navigation.resetSurfaceIndex(state.options.direction);
+    for (const Surface* surface : state.navigation.options.surfaces) {
+      // make sure we skip over the start surface
+      state.navigation.nextSurface(state.options.direction);
+      if (surface == state.navigation.currentSurface) {
+        break;
+      }
+    }
+    ACTS_VERBOSE("Start surface index set to "
+                 << state.navigation.surfaceIndex);
+    if (state.navigation.endOfSurfaces()) {
+      ACTS_DEBUG(
+          "Did not find the start surface in the sequence. Assuming it is not "
+          "part of the sequence. Trusting the correctness of the input "
+          "sequence. Resetting the surface index.");
+      state.navigation.resetSurfaceIndex(state.options.direction);
     }
+
+    state.navigation.navigationBreak = false;
+    state.navigation.targetReached = false;
   }
 
   /// @brief Navigator pre step call
@@ -157,21 +212,25 @@ class DirectNavigator {
   /// @param [in] stepper Stepper in use
   template <typename propagator_state_t, typename stepper_t>
   void preStep(propagator_state_t& state, const stepper_t& stepper) const {
+    if (state.navigation.navigationBreak) {
+      return;
+    }
+
     ACTS_VERBOSE("pre step");
 
     // Navigator target always resets the current surface
     state.navigation.currentSurface = nullptr;
 
     // Output the position in the sequence
-    ACTS_VERBOSE((state.navigation.options.surfaces.size() -
-                  state.navigation.surfaceIndex)
+    ACTS_VERBOSE(state.navigation.remainingSurfaces(state.options.direction)
                  << " out of " << state.navigation.options.surfaces.size()
                  << " surfaces remain to try.");
 
-    if (state.navigation.surfaceIndex >=
-        state.navigation.options.surfaces.size()) {
+    if (state.navigation.endOfSurfaces()) {
       // Set the navigation break
+      ACTS_VERBOSE("End of surfaces reached, navigation break.");
       state.navigation.navigationBreak = true;
+      stepper.releaseStepSize(state.stepping, ConstrainedStep::actor);
       // If no externally provided target is given, the target is reached
       if (state.navigation.options.targetSurface == nullptr) {
         state.navigation.targetReached = true;
@@ -183,8 +242,7 @@ class DirectNavigator {
 
     // Establish & update the surface status
     // TODO we do not know the intersection index - passing the closer one
-    const auto& surface =
-        *state.navigation.options.surfaces.at(state.navigation.surfaceIndex);
+    const auto& surface = *state.navigation.navSurface();
     const double farLimit = std::numeric_limits<double>::max();
     const auto index =
         chooseIntersection(
@@ -202,7 +260,7 @@ class DirectNavigator {
           "Surface not reachable anymore, switching to next one in "
           "sequence");
       // Move the sequence to the next surface
-      ++state.navigation.surfaceIndex;
+      state.navigation.nextSurface(state.options.direction);
     } else {
       ACTS_VERBOSE("Navigation stepSize set to "
                    << stepper.outputStepSize(state.stepping));
@@ -218,26 +276,27 @@ class DirectNavigator {
   /// @param [in] stepper Stepper in use
   template <typename propagator_state_t, typename stepper_t>
   void postStep(propagator_state_t& state, const stepper_t& stepper) const {
+    if (state.navigation.navigationBreak) {
+      return;
+    }
+
     ACTS_VERBOSE("post step");
 
     // Navigator post step always resets the current surface
     state.navigation.currentSurface = nullptr;
 
     // Output the position in the sequence
-    ACTS_VERBOSE((state.navigation.options.surfaces.size() -
-                  state.navigation.surfaceIndex)
+    ACTS_VERBOSE(state.navigation.remainingSurfaces(state.options.direction)
                  << " out of " << state.navigation.options.surfaces.size()
                  << " surfaces remain to try.");
 
-    if (state.navigation.surfaceIndex >=
-        state.navigation.options.surfaces.size()) {
+    if (state.navigation.endOfSurfaces()) {
       return;
     }
 
     // Establish the surface status
     // TODO we do not know the intersection index - passing the closer one
-    const auto& surface =
-        *state.navigation.options.surfaces.at(state.navigation.surfaceIndex);
+    const auto& surface = *state.navigation.navSurface();
     const double farLimit = std::numeric_limits<double>::max();
     const auto index =
         chooseIntersection(
@@ -252,14 +311,12 @@ class DirectNavigator {
         *m_logger);
     if (surfaceStatus == Intersection3D::Status::onSurface) {
       // Set the current surface
-      state.navigation.currentSurface =
-          state.navigation.options.surfaces.at(state.navigation.surfaceIndex);
+      state.navigation.currentSurface = state.navigation.navSurface();
       ACTS_VERBOSE("Current surface set to  "
                    << state.navigation.currentSurface->geometryId());
       // Move the sequence to the next surface
-      ++state.navigation.surfaceIndex;
-      if (state.navigation.surfaceIndex <
-          state.navigation.options.surfaces.size()) {
+      state.navigation.nextSurface(state.options.direction);
+      if (!state.navigation.endOfSurfaces()) {
         ACTS_VERBOSE("Next surface candidate is  "
                      << state.navigation.options.surfaces
                             .at(state.navigation.surfaceIndex)

Core/include/Acts/Propagator/MultiStepperAborters.hpp

@@ -71,8 +71,8 @@ struct MultiStepperSurfaceReached : public SurfaceReached {
       }
 
       ACTS_VERBOSE(
-          "MultiStepperSurfaceReached aborter | "
-          "Target intersection not found. Maybe next time?");
+          "MultiStepperSurfaceReached aborter | Average distance to target: "
+          << sIntersection.pathLength());
     }
 
     bool reached = true;
@@ -84,6 +84,7 @@ struct MultiStepperSurfaceReached : public SurfaceReached {
 
       if (!SurfaceReached::checkAbort(singleState, singleStepper, navigator,
                                       logger)) {
+        cmp.status() = Acts::Intersection3D::Status::reachable;
         reached = false;
       } else {
         cmp.status() = Acts::Intersection3D::Status::onSurface;

Core/include/Acts/Propagator/Navigator.hpp

@@ -205,8 +205,6 @@ class Navigator {
       : m_cfg{std::move(cfg)}, m_logger{std::move(_logger)} {}
 
   State makeState(const Options& options) const {
-    assert(options.startSurface != nullptr && "Start surface must be set");
-
     State state;
     state.options = options;
     state.startSurface = options.startSurface;

Core/include/Acts/Propagator/TryAllNavigator.hpp

@@ -96,8 +96,6 @@ class TryAllNavigatorBase {
       : m_cfg(std::move(cfg)), m_logger{std::move(_logger)} {}
 
   State makeState(const Options& options) const {
-    assert(options.startSurface != nullptr && "Start surface must be set");
-
     State state;
     state.options = options;
     state.startSurface = options.startSurface;
@@ -587,8 +585,6 @@ class TryAllOverstepNavigator : public TryAllNavigatorBase {
       : TryAllNavigatorBase(std::move(cfg), std::move(logger)) {}
 
   State makeState(const Options& options) const {
-    assert(options.startSurface != nullptr && "Start surface must be set");
-
     State state;
     state.options = options;
     state.startSurface = options.startSurface;

Core/include/Acts/TrackFitting/GaussianSumFitter.hpp

@@ -123,15 +123,11 @@ struct GaussianSumFitter {
       using Actors = ActorList<GsfActor>;
       using PropagatorOptions = typename propagator_t::template Options<Actors>;
 
-      std::vector<const Surface*> backwardSequence(
-          std::next(sSequence.rbegin()), sSequence.rend());
-      backwardSequence.push_back(opts.referenceSurface);
-
       PropagatorOptions propOptions(opts.geoContext, opts.magFieldContext);
 
       propOptions.setPlainOptions(opts.propagatorPlainOptions);
 
-      propOptions.navigation.surfaces = backwardSequence;
+      propOptions.navigation.surfaces = sSequence;
       propOptions.actorList.template get<GsfActor>()
           .m_cfg.bethe_heitler_approx = &m_betheHeitlerApproximation;