refactor: Keep unbound surfaces in the navigation stream

Core/include/Acts/Propagator/NavigationTarget.hpp

@@ -189,7 +189,11 @@ class NavigationTarget {
   constexpr const BoundaryTolerance& boundaryTolerance() const noexcept {
     return m_boundaryTolerance;
   }
-
+  /// Returns whether the target has been reached by the propagator.
+  /// @return the target passed flag
+  constexpr bool isReached() const noexcept { return m_reached; }
+  /// Marks the target as reached
+  constexpr void targetReached() noexcept { m_reached = true; }
   /// Returns whether the intersection was successful or not
   /// @return true if the intersection is valid
   constexpr bool isValid() const noexcept { return m_intersection.isValid(); }
@@ -279,11 +283,13 @@ class NavigationTarget {
   const Surface* m_surfaceRepresentation = nullptr;
   /// The boundary tolerance used for this intersection
   BoundaryTolerance m_boundaryTolerance = BoundaryTolerance::None();
+  /// State toggling whether the target has been reached by the propgator
+  bool m_reached{false};
 
   /// Default constructor creating a none target
   constexpr NavigationTarget() = default;
 
-  /// @brief print method
+  /// Print method
   /// @param ostr: Stream to which the object is printed
   void print(std::ostream& ostr) const;
 };

Core/include/Acts/Propagator/Navigator.hpp

@@ -241,7 +241,6 @@ class Navigator final {
       navLayerIndex.reset();
       navBoundaries.clear();
       navBoundaryIndex.reset();
-      navCandidates.clear();
       navCandidateIndex.reset();
 
       currentLayer = nullptr;

Core/include/Acts/Propagator/NavigatorOptions.hpp

@@ -44,6 +44,10 @@ struct NavigatorPlainOptions {
   /// The far limit to resolve surfaces
   double farLimit = std::numeric_limits<double>::max();
 
+  /// Switch deciding whether surfaces without boundary
+  /// are cleared after a volume switch although they've not been reached yet
+  bool eraseUnboundVolChange{true};
+
   /// Delegate to decide whether free surfaces are appended to the navigation
   /// stream given the current volume and the track coordinates. If the
   /// delegate is set, it is called in each candidate resolution step

Core/include/Acts/Surfaces/BoundaryTolerance.hpp

@@ -224,7 +224,20 @@ class BoundaryTolerance {
   bool isTolerated(const Vector2& boundDelta,
                    const SquareMatrix2& boundToCartesian) const;
 
+  /// Define the ostream operator to print the object
+  /// @param ostr: Reference to the ostream
+  /// @param tolerance: Reference to the tolerance to print
+  friend std::ostream& operator<<(std::ostream& ostr,
+                                  const BoundaryTolerance& tolerance) {
+    tolerance.print(ostr);
+    return ostr;
+  }
+
  private:
+  /// Print method
+  /// @param ostr: Stream to which the object is printed
+  void print(std::ostream& ostr) const;
+
   Variant m_variant;
 
   /// Check if the boundary check is of a specific type.

Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp

@@ -30,13 +30,15 @@
 #include "Acts/TrackFitting/detail/VoidFitterComponents.hpp"
 #include "Acts/Utilities/CalibrationContext.hpp"
 #include "Acts/Utilities/Delegate.hpp"
+#include "Acts/Utilities/Helpers.hpp"
 #include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/Result.hpp"
 #include "Acts/Utilities/TrackHelpers.hpp"
 
 #include <functional>
 #include <limits>
 #include <memory>
+#include <span>
 #include <type_traits>
 #include <unordered_map>
 
@@ -185,6 +187,10 @@ struct Gx2FitterOptions {
   /// Max number of iterations during the fit (abort condition)
   std::size_t nUpdateMax = 5;
 
+  /// Number of iterations before missed measurements are tried to
+  /// be added back to the track
+  std::size_t nCoolDownMissed = 10;
+
   /// Check for convergence (abort condition). Set to 0 to skip.
   double relChi2changeCutOff = 1e-7;
 };
@@ -231,9 +237,8 @@ struct Gx2FitterResult {
   /// Measurement surfaces without hits
   std::vector<const Surface*> missedActiveSurfaces;
 
-  /// Measurement surfaces handled in both forward and
-  /// backward filtering
-  std::vector<const Surface*> passedAgainSurfaces;
+  /// Hit measurement surfaces
+  std::vector<const Surface*> hitActiveSurfaces;
 
   /// Count how many surfaces have been hit
   std::size_t surfaceCount = 0;
@@ -756,6 +761,8 @@ class Gx2Fitter {
 
     /// Allows retrieving measurements for a surface
     const std::unordered_map<const Surface*, SourceLink>* inputMeasurements{};
+    /// Allows to skip measurements that cannot be considered as active
+    std::span<const Surface*> vetoedActiveSurfaces{};
 
     /// Whether to consider multiple scattering.
     bool multipleScattering = false;
@@ -900,8 +907,10 @@ class Gx2Fitter {
 
       // Here we handle all measurements
       if (auto sourceLinkIt = inputMeasurements->find(surface);
-          sourceLinkIt != inputMeasurements->end()) {
+          sourceLinkIt != inputMeasurements->end() &&
+          !Acts::rangeContainsValue(vetoedActiveSurfaces, surface)) {
         ACTS_DEBUG("    The surface contains a measurement.");
+        result.hitActiveSurfaces.push_back(surface);
 
         // Transport the covariance to the surface
         stepper.transportCovarianceToBound(state.stepping, *surface,
@@ -1210,6 +1219,9 @@ class Gx2Fitter {
     ACTS_VERBOSE("Preparing " << std::distance(it, end)
                               << " input measurements");
     std::unordered_map<const Surface*, SourceLink> inputMeasurements{};
+    // Keep track of the hits that were not hit during the propagation
+    std::vector<const Surface*> missedMeasurements{};
+    std::size_t lastMissedIter{0u};
 
     for (; it != end; ++it) {
       inputMeasurements.try_emplace(gx2fOptions.extensions.surfaceAccessor(*it),
@@ -1266,6 +1278,12 @@ class Gx2Fitter {
     for (nUpdate = 0; nUpdate < gx2fOptions.nUpdateMax; nUpdate++) {
       ACTS_DEBUG("nUpdate = " << nUpdate + 1 << "/" << gx2fOptions.nUpdateMax);
 
+      if (!missedMeasurements.empty() &&
+          lastMissedIter + gx2fOptions.nCoolDownMissed <= nUpdate) {
+        ACTS_DEBUG("Try to add back " << missedMeasurements.size()
+                                      << " measurements");
+        missedMeasurements.clear();
+      }
       // set up propagator and co
       PropagatorOptions propagatorOptions{gx2fOptions.propagatorPlainOptions};
 
@@ -1277,6 +1295,7 @@ class Gx2Fitter {
 
       auto& gx2fActor = propagatorOptions.actorList.template get<GX2FActor>();
       gx2fActor.inputMeasurements = &inputMeasurements;
+      gx2fActor.vetoedActiveSurfaces = missedMeasurements;
       gx2fActor.multipleScattering = false;
       gx2fActor.extensions = gx2fOptions.extensions;
       gx2fActor.calibrationContext = &gx2fOptions.calibrationContext.get();
@@ -1296,6 +1315,7 @@ class Gx2Fitter {
 
       auto& r = propagatorState.template get<Gx2FitterResult<traj_t>>();
       r.fittedStates = &trajectoryTempBackend;
+      r.hitActiveSurfaces.reserve(inputMeasurements.size());
 
       // Clear the track container. It could be more performant to update the
       // existing states, but this needs some more thinking.
@@ -1316,11 +1336,31 @@ class Gx2Fitter {
       // TODO Improve Propagator + Actor [allocate before loop], rewrite
       // makeMeasurements
       auto& propRes = *result;
+
       GX2FResult gx2fResult = std::move(propRes.template get<GX2FResult>());
 
       auto track = trackContainerTemp.makeTrack();
       tipIndex = gx2fResult.lastMeasurementIndex;
 
+      if (gx2fResult.hitActiveSurfaces.size() + missedMeasurements.size() !=
+              inputMeasurements.size()) {
+        ACTS_DEBUG("Propagation only hit "
+                   << gx2fResult.hitActiveSurfaces.size() << "/"
+                   << (inputMeasurements.size() - missedMeasurements.size())
+                   << " surfaces. Remove non hit surfaces from the list");
+        for (const auto& [surface, _] : inputMeasurements) {
+          if (!Acts::rangeContainsValue(gx2fResult.hitActiveSurfaces,
+                                        surface)) {
+            missedMeasurements.push_back(surface);
+            ACTS_VERBOSE(" -- Remove hit associated to "
+                         << surface->geometryId() << ".");
+          }
+        }
+      }
+      ACTS_VERBOSE("Found " << gx2fResult.hitActiveSurfaces.size()
+                            << " measurements, "
+                            << ", " << gx2fResult.measurementHoles << " holes");
+
       // It could happen, that no measurements were found. Then the track would
       // be empty and the following operations would be invalid. Usually, this
       // only happens during the first iteration, due to bad initial parameters.
@@ -1372,6 +1412,9 @@ class Gx2Fitter {
         ACTS_INFO("Not enough measurements. Require "
                   << extendedSystem.findRequiredNdf() + 1 << ", but only "
                   << extendedSystem.ndf() << " could be used.");
+        for (const Surface* surface : gx2fResult.hitActiveSurfaces) {
+          ACTS_INFO(" --- " << surface->geometryId());
+        }
         return Experimental::GlobalChiSquareFitterError::NotEnoughMeasurements;
       }
 
@@ -1543,7 +1586,8 @@ class Gx2Fitter {
 
       ACTS_VERBOSE("aMatrix:\n"
                    << extendedSystem.aMatrix() << "\n"
-                   << "bVector:\n"
+                   << "determinant: " << extendedSystem.aMatrix().determinant()
+                   << ", bVector:\n"
                    << extendedSystem.bVector() << "\n"
                    << "deltaParamsExtended:\n"
                    << deltaParamsExtended << "\n"

Core/src/Navigation/NavigationStream.cpp

@@ -84,8 +84,9 @@ bool NavigationStream::initialize(const GeometryContext& gctx,
   }
 
   // Append the multi intersection candidates
-  m_candidates.insert(m_candidates.end(), additionalCandidates.begin(),
-                      additionalCandidates.end());
+  m_candidates.insert(m_candidates.end(),
+                      std::make_move_iterator(additionalCandidates.begin()),
+                      std::make_move_iterator(additionalCandidates.end()));
 
   // Sort the candidates by path length
   std::ranges::sort(m_candidates, NavigationTarget::pathLengthOrder);

Core/src/Navigation/Navigator.cpp

@@ -313,6 +313,8 @@ void Navigator::handleSurfaceReached(State& state, const Vector3& position,
 
   // handling portals in gen3 configuration
   if (m_geometryVersion == GeometryVersion::Gen3) {
+    state.navCandidate().targetReached();
+
     if (state.navCandidate().isPortalTarget() &&
         &state.navCandidate().surface() == &surface) {
       ACTS_VERBOSE(volInfo(state) << "Handling portal status.");
@@ -581,12 +583,6 @@ void Navigator::resolveCandidates(State& state, const Vector3& position,
         "Navigator: No navigation policy found for current volume.");
   }
 
-  auto policyState = state.policyStateManager.currentState();
-  state.currentVolume->initializeNavigationCandidates(
-      state.options.geoContext, args, policyState, appendOnly, logger());
-
-  ACTS_VERBOSE(volInfo(state) << "Found " << state.stream.candidates().size()
-                              << " navigation candidates.");
   for (const Surface* surface : state.options.externalSurfaces) {
     const GeometryIdentifier geoId = surface->geometryId();
     // Don't add any surface which is not in the same volume (volume bits)
@@ -598,6 +594,28 @@ void Navigator::resolveCandidates(State& state, const Vector3& position,
                                 << " surface " << geoId);
     appendOnly.addSurfaceCandidate(*surface, BoundaryTolerance::Infinite());
   }
+
+  if (!state.options.eraseUnboundVolChange) {
+    // unbound targets are copied over to the next surface if they've not been
+    // yet reached by the propagator
+    for (const NavigationTarget& target : state.navCandidates) {
+      if (target.isSurfaceTarget() && target.boundaryTolerance().isInfinite() &&
+          !target.isReached()) {
+        appendOnly.addSurfaceCandidate(target.surface(),
+                                       BoundaryTolerance::Infinite());
+        ACTS_VERBOSE(volInfo(state) << "Target " << target
+                                    << " not been reached yet. Try now");
+      }
+    }
+  }
+
+  auto policyState = state.policyStateManager.currentState();
+  state.currentVolume->initializeNavigationCandidates(
+      state.options.geoContext, args, policyState, appendOnly, logger());
+
+  ACTS_VERBOSE(volInfo(state) << "Found " << state.stream.candidates().size()
+                              << " navigation candidates.");
+
   bool pruneFreeCand{false};
   if (!state.freeCandidates.empty()) {
     for (const auto& [surface, wasReached] : state.freeCandidates) {
@@ -615,8 +633,11 @@ void Navigator::resolveCandidates(State& state, const Vector3& position,
         appendOnly.addSurfaceCandidate(*surface, BoundaryTolerance::Infinite());
         pruneFreeCand = !state.options.freeSurfaceSelector.connected();
       }
-    };
+    }
   }
+
+  state.navCandidates.clear();
+
   state.stream.initialize(state.options.geoContext, {position, direction},
                           BoundaryTolerance::None(),
                           state.options.surfaceTolerance);
@@ -626,8 +647,6 @@ void Navigator::resolveCandidates(State& state, const Vector3& position,
                << " navigation candidates after initialization.\n"
                << state.stream.candidates());
 
-  state.navCandidates.clear();
-
   double farLimit = state.options.farLimit;
   // If the user has not provided the selection delegate, then
   // just apply a simple candidate pruning. Constrain the maximum

Core/src/Propagator/NavigationTarget.cpp

@@ -32,6 +32,8 @@ void NavigationTarget::print(std::ostream& ostr) const {
         }
       },
       m_target);
-  ostr << ", path length: " << pathLength();
+  ostr << ", path length: " << pathLength()
+       << ", reached: " << (isReached() ? "yes" : "no") << ", "
+       << boundaryTolerance();
 }
 }  // namespace Acts

Core/src/Surfaces/BoundaryTolerance.cpp

@@ -97,4 +97,20 @@ bool BoundaryTolerance::isTolerated(
   throw std::logic_error("Unsupported tolerance type");
 }
 
+void BoundaryTolerance::print(std::ostream& ostr) const {
+  if (isInfinite()) {
+    ostr << "no boundary limit applied";
+  } else if (isNone()) {
+    ostr << "strict boundary limit";
+  } else if (hasAbsoluteEuclidean()) {
+    ostr << "boundary limit  dX < " << asAbsoluteEuclidean().tolerance;
+  } else if (hasChi2Bound()) {
+    ostr << "local chi2 " << asChi2Bound().maxChi2 << " with weights:\n"
+         << asChi2Bound().weightMatrix() << "\n";
+  } else if (hasChi2Cartesian()) {
+    ostr << "local chi2 " << asChi2Cartesian().maxChi2 << " with weights:\n"
+         << asChi2Cartesian().weightMatrix() << "\n";
+  }
+}
+
 }  // namespace Acts