refactor!: Disambiguate surface intersection solutions (#2336)

Currently our intersection call can give an alternative solution but it
is not clear which one is to be preferred in which situation and the
results are unstable meaning that it could be swapped depending on where
you are on the ray.

Here I try to make the interface cleaner and the intersection order
stable.

blocked by
- #2366
- #2368

---------

Co-authored-by: kodiakhq[bot] <49736102+kodiakhq[bot]@users.noreply.github.com>
Co-authored-by: Benjamin Huth <[email protected]>
Co-authored-by: Benjamin Huth <[email protected]>
Co-authored-by: Paul Gessinger <[email protected]>

Core/include/Acts/Geometry/ApproachDescriptor.hpp

@@ -10,6 +10,7 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Intersection.hpp"
 
 #include <vector>
@@ -42,11 +43,17 @@ class ApproachDescriptor {
   /// @param position is the position from start of the search
   /// @param direction is the direction at the start of the search
   /// @param bcheck is the boundary check directive
+  /// @param pLimit The path limit
+  /// @param oLimit The overstep limit
+  /// @param tolerance The surface tolerance
   ///
   /// @return is a surface intersection
-  virtual ObjectIntersection<Surface> approachSurface(
-      const GeometryContext& gctx, const Vector3& position,
-      const Vector3& direction, const BoundaryCheck& bcheck) const = 0;
+  virtual SurfaceIntersection approachSurface(const GeometryContext& gctx,
+                                              const Vector3& position,
+                                              const Vector3& direction,
+                                              const BoundaryCheck& bcheck,
+                                              double pLimit, double oLimit,
+                                              double tolerance) const = 0;
 
   /// Get all the contained surfaces
   /// @return all contained surfaces of this approach descriptor

Core/include/Acts/Geometry/GenericApproachDescriptor.hpp

@@ -12,6 +12,7 @@
 #include "Acts/Geometry/ApproachDescriptor.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Surfaces/BoundaryCheck.hpp"
+#include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Helpers.hpp"
 #include "Acts/Utilities/Intersection.hpp"
 
@@ -28,10 +29,6 @@ class Surface;
 ///
 /// Class to decide and return which approaching surface to be taken,
 /// it's a generic descriptor for n surfaces
-///
-/// It is templated in order to allow for BoundarySurfaces from
-/// representing volumes of layers to be re-used
-
 class GenericApproachDescriptor : public ApproachDescriptor {
  public:
   /// A generic approach descriptor for new Acts::Surface objects
@@ -60,11 +57,17 @@ class GenericApproachDescriptor : public ApproachDescriptor {
   /// @param position The global position to start the approach from
   /// @param direction The momentum vector
   /// @param bcheck The boundary check prescription
+  /// @param pLimit The path limit
+  /// @param oLimit The overstep limit
+  /// @param tolerance The surface tolerance
   ///
-  /// @return : a SurfaceIntersection
-  ObjectIntersection<Surface> approachSurface(
-      const GeometryContext& gctx, const Vector3& position,
-      const Vector3& direction, const BoundaryCheck& bcheck) const override;
+  /// @return : a @c SurfaceIntersection
+  SurfaceIntersection approachSurface(const GeometryContext& gctx,
+                                      const Vector3& position,
+                                      const Vector3& direction,
+                                      const BoundaryCheck& bcheck,
+                                      double pLimit, double oLimit,
+                                      double tolerance) const override;
 
   /// return all contained surfaces of this approach descriptor
   const std::vector<const Surface*>& containedSurfaces() const override;

Core/include/Acts/Geometry/TrackingVolume.hpp

@@ -68,13 +68,18 @@ using MutableTrackingVolumeVector = std::vector<MutableTrackingVolumePtr>;
 using LayerArray = BinnedArray<LayerPtr>;
 using LayerVector = std::vector<LayerPtr>;
 
-// Intersection with Layer
+/// Intersection with @c Layer
 using LayerIntersection = ObjectIntersection<Layer, Surface>;
+/// Multi-intersection with @c Layer
+using LayerMultiIntersection = ObjectMultiIntersection<Layer, Surface>;
 
 /// BoundarySurface of a volume
 using BoundarySurface = BoundarySurfaceT<TrackingVolume>;
 /// Intersection with a @c BoundarySurface
 using BoundaryIntersection = ObjectIntersection<BoundarySurface, Surface>;
+/// Multi-intersection with a @c BoundarySurface
+using BoundaryMultiIntersection =
+    ObjectMultiIntersection<BoundarySurface, Surface>;
 
 /// @class TrackingVolume
 ///

Core/include/Acts/Navigation/NavigationStateFillers.hpp

@@ -51,7 +51,7 @@ struct SurfacesFiller {
                           const std::vector<const Surface*>& surfaces) {
     std::for_each(surfaces.begin(), surfaces.end(), [&](const auto& s) {
       nState.surfaceCandidates.push_back(NavigationState::SurfaceCandidate{
-          ObjectIntersection<Surface>{}, s, nullptr,
+          ObjectIntersection<Surface>::invalid(), s, nullptr,
           nState.surfaceBoundaryCheck});
     });
   }
@@ -68,7 +68,7 @@ struct PortalsFiller {
                           const std::vector<const Portal*>& portals) {
     std::for_each(portals.begin(), portals.end(), [&](const auto& p) {
       nState.surfaceCandidates.push_back(NavigationState::SurfaceCandidate{
-          ObjectIntersection<Surface>{}, nullptr, p, true});
+          ObjectIntersection<Surface>::invalid(), nullptr, p, true});
     });
   }
 };

Core/include/Acts/Navigation/SurfaceCandidatesUpdators.hpp

@@ -43,26 +43,21 @@ inline static void updateCandidates(const GeometryContext& gctx,
   for (auto& c : nCandidates) {
     // Get the surface representation: either native surfcae of portal
     const Surface& sRep =
-        (c.surface != nullptr) ? (*c.surface) : (c.portal->surface());
+        c.surface != nullptr ? *c.surface : c.portal->surface();
 
     // Get the intersection @todo make a templated intersector
-    auto sIntersection =
-        sRep.intersect(gctx, position, direction, c.boundaryCheck);
-    // Re-order and swap if necessary
-    if (sIntersection.intersection.pathLength + s_onSurfaceTolerance <
-            nState.overstepTolerance and
-        sIntersection.alternative.status >= Intersection3D::Status::reachable) {
-      sIntersection.swapSolutions();
-    }
-    c.objectIntersection = sIntersection;
+    // TODO surface tolerance
+    auto sIntersection = sRep.intersect(gctx, position, direction,
+                                        c.boundaryCheck, s_onSurfaceTolerance);
+    c.objectIntersection = sIntersection[c.objectIntersection.index()];
   }
   // Sort and stuff non-allowed solutions to the end
   std::sort(
       nCandidates.begin(), nCandidates.end(),
       [&](const auto& a, const auto& b) {
         // The two path lengths
-        ActsScalar pathToA = a.objectIntersection.intersection.pathLength;
-        ActsScalar pathToB = b.objectIntersection.intersection.pathLength;
+        ActsScalar pathToA = a.objectIntersection.pathLength();
+        ActsScalar pathToB = b.objectIntersection.pathLength();
         if (pathToA + s_onSurfaceTolerance < nState.overstepTolerance or
             std::abs(pathToA) < s_onSurfaceTolerance) {
           return false;

Core/include/Acts/Propagator/AtlasStepper.hpp

@@ -375,7 +375,10 @@ class AtlasStepper {
   }
 
   /// Overstep limit
-  double overstepLimit(const State& /*state*/) const {
+  ///
+  /// @param state The stepping state (thread-local cache)
+  double overstepLimit(const State& state) const {
+    (void)state;
     return -m_overstepLimit;
   }
 
@@ -414,7 +417,7 @@ class AtlasStepper {
   /// @param release [in] boolean to trigger step size release
   template <typename object_intersection_t>
   void updateStepSize(State& state, const object_intersection_t& oIntersection,
-                      bool release = true) const {
+                      Direction /*direction*/, bool release = true) const {
     detail::updateSingleStepSize<AtlasStepper>(state, oIntersection, release);
   }
 

Core/include/Acts/Propagator/EigenStepper.hpp

@@ -267,7 +267,7 @@ class EigenStepper {
   /// @param release [in] boolean to trigger step size release
   template <typename object_intersection_t>
   void updateStepSize(State& state, const object_intersection_t& oIntersection,
-                      bool release = true) const {
+                      Direction /*direction*/, bool release = true) const {
     detail::updateSingleStepSize<EigenStepper>(state, oIntersection, release);
   }
 
@@ -307,7 +307,10 @@ class EigenStepper {
   }
 
   /// Overstep limit
-  double overstepLimit(const State& /*state*/) const {
+  ///
+  /// @param state The stepping state (thread-local cache)
+  double overstepLimit(const State& state) const {
+    (void)state;
     // A dynamic overstep limit could sit here
     return -m_overstepLimit;
   }

Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp

@@ -23,12 +23,14 @@
 #include "Acts/Propagator/EigenStepper.hpp"
 #include "Acts/Propagator/EigenStepperError.hpp"
 #include "Acts/Propagator/Propagator.hpp"
+#include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/GaussianMixtureReduction.hpp"
 #include "Acts/Utilities/Intersection.hpp"
 #include "Acts/Utilities/Result.hpp"
 
 #include <algorithm>
 #include <cmath>
+#include <cstddef>
 #include <functional>
 #include <limits>
 #include <numeric>
@@ -770,30 +772,21 @@ class MultiEigenStepperLoop
   ///
   /// @param state [in,out] The stepping state (thread-local cache)
   /// @param oIntersection [in] The ObjectIntersection to layer, boundary, etc
+  /// @param direction [in] The propagation direction
   /// @param release [in] boolean to trigger step size release
   template <typename object_intersection_t>
   void updateStepSize(State& state, const object_intersection_t& oIntersection,
-                      bool release = true) const {
-    const Surface& surface = *oIntersection.representation;
+                      Direction direction, bool release = true) const {
+    const Surface& surface = *oIntersection.representation();
 
     for (auto& component : state.components) {
       auto intersection = surface.intersect(
           component.state.geoContext, SingleStepper::position(component.state),
-          SingleStepper::direction(component.state), true);
+          direction * SingleStepper::direction(component.state),
+          true)[oIntersection.index()];
 
-      // We don't know whatever was done to manipulate the intersection before
-      // (e.g. in Layer.ipp:240), so we trust and just adjust the sign
-      if (std::signbit(oIntersection.intersection.pathLength) !=
-          std::signbit(intersection.intersection.pathLength)) {
-        intersection.intersection.pathLength *= -1;
-      }
-
-      if (std::signbit(oIntersection.alternative.pathLength) !=
-          std::signbit(intersection.alternative.pathLength)) {
-        intersection.alternative.pathLength *= -1;
-      }
-
-      SingleStepper::updateStepSize(component.state, intersection, release);
+      SingleStepper::updateStepSize(component.state, intersection, direction,
+                                    release);
     }
   }
 

Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp

@@ -39,7 +39,8 @@ auto MultiEigenStepperLoop<E, R, A>::boundState(
             .intersect(state.geoContext,
                        cmpState.pars.template segment<3>(eFreePos0),
                        cmpState.pars.template segment<3>(eFreeDir0), false)
-            .intersection.position;
+            .closest()
+            .position();
 
     auto bs = SingleStepper::boundState(cmpState, surface, transportCov,
                                         freeToBoundCorrection);

Core/include/Acts/Propagator/MultiStepperAborters.hpp

@@ -78,13 +78,15 @@ struct MultiStepperSurfaceReached {
 
     // However, if mean of all is on surface, we are happy as well
     if (averageOnSurface) {
-      const auto sIntersection = targetSurface.intersect(
-          state.geoContext, stepper.position(state.stepping),
-          state.options.direction * stepper.direction(state.stepping), true,
-          averageOnSurfaceTolerance);
-
-      if (sIntersection.intersection.status ==
-          Intersection3D::Status::onSurface) {
+      const auto sIntersection =
+          targetSurface
+              .intersect(
+                  state.geoContext, stepper.position(state.stepping),
+                  state.options.direction * stepper.direction(state.stepping),
+                  true, averageOnSurfaceTolerance)
+              .closest();
+
+      if (sIntersection.status() == Intersection3D::Status::onSurface) {
         ACTS_VERBOSE("Reached target in average mode");
         navigator.currentSurface(state.navigation, &targetSurface);
         navigator.targetReached(state.navigation, true);