refactor: Make surface intersection tolerance configurable (#2106)

`s_onSurfaceTolerance` is used in a bunch of places where the user might actually expect a configurable tolerance.

in this PR I add a tolerance parameter to `Surface::intersect` which can then be used by the navigator or target reached aborter

Core/include/Acts/Navigation/NextNavigator.hpp

@@ -45,9 +45,6 @@ class NextNavigator {
     bool resolveMaterial = true;
     /// stop at every surface regardless what it is
     bool resolvePassive = false;
-
-    /// The tolerance used to defined "reached"
-    double tolerance = s_onSurfaceTolerance;
   };
 
   /// Nested State struct

Core/include/Acts/Propagator/DirectNavigator.hpp

@@ -40,9 +40,6 @@ class DirectNavigator {
                       getDefaultLogger("DirectNavigator", Logging::INFO))
       : m_logger{std::move(_logger)} {}
 
-  /// The tolerance used to define "surface reached"
-  double tolerance = s_onSurfaceTolerance;
-
   /// Nested Actor struct, called Initializer
   ///
   /// This is needed for the initialization of the

Core/include/Acts/Propagator/Navigator.hpp

@@ -150,9 +150,6 @@ class Navigator {
     /// stop at every surface regardless what it is
     bool resolvePassive = false;
 
-    /// The tolerance used to defined "reached"
-    double tolerance = s_onSurfaceTolerance;
-
     /// Wether to perform boundary checks for layer resolving (improves
     /// navigation for bended tracks)
     BoundaryCheck boundaryCheckLayerResolving = true;
@@ -1118,7 +1115,8 @@ class Navigator {
       // target volume and layer search through global search
       auto targetIntersection = state.navigation.targetSurface->intersect(
           state.geoContext, stepper.position(state.stepping),
-          state.stepping.navDir * stepper.direction(state.stepping), false);
+          state.stepping.navDir * stepper.direction(state.stepping), false,
+          state.options.targetTolerance);
       if (targetIntersection) {
         ACTS_VERBOSE(volInfo(state)
                      << "Target estimate position ("

Core/include/Acts/Propagator/StandardAborters.hpp

@@ -140,7 +140,8 @@ struct SurfaceReached {
     const double tolerance = state.options.targetTolerance;
     const auto sIntersection = targetSurface.intersect(
         state.geoContext, stepper.position(state.stepping),
-        state.stepping.navDir * stepper.direction(state.stepping), true);
+        state.stepping.navDir * stepper.direction(state.stepping), true,
+        tolerance);
 
     // The target is reached
     bool targetReached = (sIntersection.intersection.status ==

Core/include/Acts/Surfaces/ConeSurface.hpp

@@ -167,14 +167,15 @@ class ConeSurface : public Surface {
   /// @param position The position to start from
   /// @param direction The direction at start
   /// @param bcheck the Boundary Check
+  /// @param tolerance the tolerance used for the intersection
   ///
   /// If possible returns both solutions for the cylinder
   ///
   /// @return SurfaceIntersection object (contains intersection & surface)
-  SurfaceIntersection intersect(const GeometryContext& gctx,
-                                const Vector3& position,
-                                const Vector3& direction,
-                                const BoundaryCheck& bcheck) const final;
+  SurfaceIntersection intersect(
+      const GeometryContext& gctx, const Vector3& position,
+      const Vector3& direction, const BoundaryCheck& bcheck = false,
+      double tolerance = s_onSurfaceTolerance) const final;
 
   /// The pathCorrection for derived classes with thickness
   ///

Core/include/Acts/Surfaces/CylinderSurface.hpp

@@ -177,14 +177,15 @@ class CylinderSurface : public Surface {
   /// @param position The position to start from
   /// @param direction The direction at start
   /// @param bcheck the Boundary Check
+  /// @param tolerance the tolerance used for the intersection
   ///
   /// If possible returns both solutions for the cylinder
   ///
   /// @return SurfaceIntersection object (contains intersection & surface)
-  SurfaceIntersection intersect(const GeometryContext& gctx,
-                                const Vector3& position,
-                                const Vector3& direction,
-                                const BoundaryCheck& bcheck) const final;
+  SurfaceIntersection intersect(
+      const GeometryContext& gctx, const Vector3& position,
+      const Vector3& direction, const BoundaryCheck& bcheck = false,
+      ActsScalar tolerance = s_onSurfaceTolerance) const final;
 
   /// Path correction due to incident of the track
   ///

Core/include/Acts/Surfaces/DiscSurface.hpp

@@ -243,6 +243,7 @@ class DiscSurface : public Surface {
   /// @param direction The global direction at the starting point
   ///        @note expected to be normalized (no checking)
   /// @param bcheck The boundary check prescription
+  /// @param tolerance the tolerance used for the intersection
   ///
   ///  <b>mathematical motivation:</b>
   ///
@@ -263,8 +264,8 @@ class DiscSurface : public Surface {
   /// @return The SurfaceIntersection object
   SurfaceIntersection intersect(
       const GeometryContext& gctx, const Vector3& position,
-      const Vector3& direction,
-      const BoundaryCheck& bcheck = false) const final;
+      const Vector3& direction, const BoundaryCheck& bcheck = false,
+      ActsScalar tolerance = s_onSurfaceTolerance) const final;
 
   /// Implement the binningValue
   ///

Core/include/Acts/Surfaces/LineSurface.hpp

@@ -197,6 +197,7 @@ class LineSurface : public Surface {
   /// @param direction The global direction at the starting point
   ///        @note exptected to be normalized
   /// @param bcheck The boundary check directive for the estimate
+  /// @param tolerance the tolerance used for the intersection
   ///
   ///   <b>mathematical motivation:</b>
   ///   Given two lines in parameteric form:<br>
@@ -228,8 +229,8 @@ class LineSurface : public Surface {
   /// @return is the intersection object
   SurfaceIntersection intersect(
       const GeometryContext& gctx, const Vector3& position,
-      const Vector3& direction,
-      const BoundaryCheck& bcheck = false) const final;
+      const Vector3& direction, const BoundaryCheck& bcheck = false,
+      ActsScalar tolerance = s_onSurfaceTolerance) const final;
 
   /// the pathCorrection for derived classes with thickness
   /// is by definition 1 for LineSurfaces

Core/include/Acts/Surfaces/PlaneSurface.hpp

@@ -156,6 +156,7 @@ class PlaneSurface : public Surface {
   /// @param direction The direction of the interesection attempt,
   /// (@note expected to be normalized)
   /// @param bcheck The boundary check directive
+  /// @param tolerance the tolerance used for the intersection
   ///
   /// <b>mathematical motivation:</b>
   ///
@@ -177,8 +178,8 @@ class PlaneSurface : public Surface {
   /// @return the SurfaceIntersection object
   SurfaceIntersection intersect(
       const GeometryContext& gctx, const Vector3& position,
-      const Vector3& direction,
-      const BoundaryCheck& bcheck = false) const final;
+      const Vector3& direction, const BoundaryCheck& bcheck = false,
+      ActsScalar tolerance = s_onSurfaceTolerance) const final;
 
   /// Return a Polyhedron for the surfaces
   ///

Core/include/Acts/Surfaces/Surface.hpp

@@ -398,12 +398,13 @@ class Surface : public virtual GeometryObject,
   /// @param position The position to start from
   /// @param direction The direction at start
   /// @param bcheck the Boundary Check
+  /// @param tolerance the tolerance used for the intersection
   ///
   /// @return SurfaceIntersection object (contains intersection & surface)
-  virtual SurfaceIntersection intersect(const GeometryContext& gctx,
-                                        const Vector3& position,
-                                        const Vector3& direction,
-                                        const BoundaryCheck& bcheck) const = 0;
+  virtual SurfaceIntersection intersect(
+      const GeometryContext& gctx, const Vector3& position,
+      const Vector3& direction, const BoundaryCheck& bcheck = false,
+      ActsScalar tolerance = s_onSurfaceTolerance) const = 0;
 
   /// Output Method for std::ostream, to be overloaded by child classes
   ///

Core/include/Acts/Surfaces/detail/PlanarHelper.hpp

@@ -25,7 +25,8 @@ namespace PlanarHelper {
 /// @return The intersection
 inline Intersection3D intersect(const Transform3& transform,
                                 const Vector3& position,
-                                const Vector3& direction) {
+                                const Vector3& direction,
+                                ActsScalar tolerance) {
   // Get the matrix from the transform (faster access)
   const auto& tMatrix = transform.matrix();
   const Vector3 pnormal = tMatrix.block<3, 1>(0, 2).transpose();
@@ -36,10 +37,9 @@ inline Intersection3D intersect(const Transform3& transform,
     // Translate that into a path
     ActsScalar path = (pnormal.dot((pcenter - position))) / (denom);
     // Is valid hence either on surface or reachable
-    Intersection3D::Status status =
-        std::abs(path) < std::abs(s_onSurfaceTolerance)
-            ? Intersection3D::Status::onSurface
-            : Intersection3D::Status::reachable;
+    Intersection3D::Status status = std::abs(path) < std::abs(tolerance)
+                                        ? Intersection3D::Status::onSurface
+                                        : Intersection3D::Status::reachable;
     // Return the intersection
     return Intersection3D{(position + path * direction), path, status};
   }

Core/src/Surfaces/ConeSurface.cpp

@@ -280,7 +280,8 @@ Acts::detail::RealQuadraticEquation Acts::ConeSurface::intersectionSolver(
 
 Acts::SurfaceIntersection Acts::ConeSurface::intersect(
     const GeometryContext& gctx, const Vector3& position,
-    const Vector3& direction, const BoundaryCheck& bcheck) const {
+    const Vector3& direction, const BoundaryCheck& bcheck,
+    ActsScalar tolerance) const {
   // Solve the quadratic equation
   auto qe = intersectionSolver(gctx, position, direction);
 
@@ -291,21 +292,19 @@ Acts::SurfaceIntersection Acts::ConeSurface::intersect(
 
   // Check the validity of the first solution
   Vector3 solution1 = position + qe.first * direction;
-  Intersection3D::Status status1 =
-      std::abs(qe.first) < std::abs(s_onSurfaceTolerance)
-          ? Intersection3D::Status::onSurface
-          : Intersection3D::Status::reachable;
+  Intersection3D::Status status1 = std::abs(qe.first) < std::abs(tolerance)
+                                       ? Intersection3D::Status::onSurface
+                                       : Intersection3D::Status::reachable;
 
   if (bcheck and not isOnSurface(gctx, solution1, direction, bcheck)) {
     status1 = Intersection3D::Status::missed;
   }
 
   // Check the validity of the second solution
   Vector3 solution2 = position + qe.first * direction;
-  Intersection3D::Status status2 =
-      std::abs(qe.second) < std::abs(s_onSurfaceTolerance)
-          ? Intersection3D::Status::onSurface
-          : Intersection3D::Status::reachable;
+  Intersection3D::Status status2 = std::abs(qe.second) < std::abs(tolerance)
+                                       ? Intersection3D::Status::onSurface
+                                       : Intersection3D::Status::reachable;
   if (bcheck and not isOnSurface(gctx, solution2, direction, bcheck)) {
     status2 = Intersection3D::Status::missed;
   }

Core/src/Surfaces/CylinderSurface.cpp

@@ -210,7 +210,8 @@ Acts::detail::RealQuadraticEquation Acts::CylinderSurface::intersectionSolver(
 
 Acts::SurfaceIntersection Acts::CylinderSurface::intersect(
     const GeometryContext& gctx, const Vector3& position,
-    const Vector3& direction, const BoundaryCheck& bcheck) const {
+    const Vector3& direction, const BoundaryCheck& bcheck,
+    ActsScalar tolerance) const {
   const auto& gctxTransform = transform(gctx);
 
   // Solve the quadratic equation
@@ -223,10 +224,9 @@ Acts::SurfaceIntersection Acts::CylinderSurface::intersect(
 
   // Check the validity of the first solution
   Vector3 solution1 = position + qe.first * direction;
-  Intersection3D::Status status1 =
-      std::abs(qe.first) < std::abs(s_onSurfaceTolerance)
-          ? Intersection3D::Status::onSurface
-          : Intersection3D::Status::reachable;
+  Intersection3D::Status status1 = std::abs(qe.first) < std::abs(tolerance)
+                                       ? Intersection3D::Status::onSurface
+                                       : Intersection3D::Status::reachable;
 
   // Helper method for boundary check
   auto boundaryCheck =
@@ -245,8 +245,8 @@ Acts::SurfaceIntersection Acts::CylinderSurface::intersect(
       // Create the reference vector in local
       const Vector3 vecLocal(solution - tMatrix.block<3, 1>(0, 3));
       double cZ = vecLocal.dot(tMatrix.block<3, 1>(0, 2));
-      double tolerance = s_onSurfaceTolerance + bcheck.tolerance()[eBoundLoc1];
-      double hZ = cBounds.get(CylinderBounds::eHalfLengthZ) + tolerance;
+      double modifiedTolerance = tolerance + bcheck.tolerance()[eBoundLoc1];
+      double hZ = cBounds.get(CylinderBounds::eHalfLengthZ) + modifiedTolerance;
       return std::abs(cZ) < std::abs(hZ) ? status
                                          : Intersection3D::Status::missed;
     }
@@ -264,10 +264,9 @@ Acts::SurfaceIntersection Acts::CylinderSurface::intersect(
   }
   // Check the validity of the second solution
   Vector3 solution2 = position + qe.second * direction;
-  Intersection3D::Status status2 =
-      std::abs(qe.second) < std::abs(s_onSurfaceTolerance)
-          ? Intersection3D::Status::onSurface
-          : Intersection3D::Status::reachable;
+  Intersection3D::Status status2 = std::abs(qe.second) < std::abs(tolerance)
+                                       ? Intersection3D::Status::onSurface
+                                       : Intersection3D::Status::reachable;
   // Check first solution for boundary compatiblity
   status2 = boundaryCheck(solution2, status2);
   Intersection3D second(solution2, qe.second, status2);

Core/src/Surfaces/DiscSurface.cpp

@@ -287,12 +287,13 @@ Acts::FreeToBoundMatrix Acts::DiscSurface::freeToBoundJacobian(
 
 Acts::SurfaceIntersection Acts::DiscSurface::intersect(
     const GeometryContext& gctx, const Vector3& position,
-    const Vector3& direction, const BoundaryCheck& bcheck) const {
+    const Vector3& direction, const BoundaryCheck& bcheck,
+    ActsScalar tolerance) const {
   // Get the contextual transform
   auto gctxTransform = transform(gctx);
   // Use the intersection helper for planar surfaces
   auto intersection =
-      PlanarHelper::intersect(gctxTransform, position, direction);
+      PlanarHelper::intersect(gctxTransform, position, direction, tolerance);
   // Evaluate boundary check if requested (and reachable)
   if (intersection.status != Intersection3D::Status::unreachable and bcheck and
       m_bounds != nullptr) {
@@ -302,9 +303,9 @@ Acts::SurfaceIntersection Acts::DiscSurface::intersect(
     const Vector2 lcartesian = tMatrix.block<3, 2>(0, 0).transpose() * vecLocal;
     if (bcheck.type() == BoundaryCheck::Type::eAbsolute and
         m_bounds->coversFullAzimuth()) {
-      double tolerance = s_onSurfaceTolerance + bcheck.tolerance()[eBoundLoc0];
+      double modifiedTolerance = tolerance + bcheck.tolerance()[eBoundLoc0];
       if (not m_bounds->insideRadialBounds(VectorHelpers::perp(lcartesian),
-                                           tolerance)) {
+                                           modifiedTolerance)) {
         intersection.status = Intersection3D::Status::missed;
       }
     } else if (not insideBounds(localCartesianToPolar(lcartesian), bcheck)) {

Core/src/Surfaces/LineSurface.cpp

@@ -133,7 +133,8 @@ const Acts::SurfaceBounds& Acts::LineSurface::bounds() const {
 
 Acts::SurfaceIntersection Acts::LineSurface::intersect(
     const GeometryContext& gctx, const Vector3& position,
-    const Vector3& direction, const BoundaryCheck& bcheck) const {
+    const Vector3& direction, const BoundaryCheck& bcheck,
+    ActsScalar tolerance) const {
   // following nominclature found in header file and doxygen documentation
   // line one is the straight track
   const Vector3& ma = position;
@@ -148,10 +149,10 @@ Acts::SurfaceIntersection Acts::LineSurface::intersect(
   double denom = 1 - eaTeb * eaTeb;
   // validity parameter
   Intersection3D::Status status = Intersection3D::Status::unreachable;
-  if (std::abs(denom) > std::abs(s_onSurfaceTolerance)) {
+  if (std::abs(denom) > std::abs(tolerance)) {
     double u = (mab.dot(ea) - mab.dot(eb) * eaTeb) / denom;
     // Check if we are on the surface already
-    status = std::abs(u) < std::abs(s_onSurfaceTolerance)
+    status = std::abs(u) < std::abs(tolerance)
                  ? Intersection3D::Status::onSurface
                  : Intersection3D::Status::reachable;
     Vector3 result = (ma + u * ea);
@@ -161,11 +162,10 @@ Acts::SurfaceIntersection Acts::LineSurface::intersect(
       // At closest approach: check inside R or and inside Z
       const Vector3 vecLocal(result - mb);
       double cZ = vecLocal.dot(eb);
-      double hZ =
-          m_bounds->get(LineBounds::eHalfLengthZ) + s_onSurfaceTolerance;
+      double hZ = m_bounds->get(LineBounds::eHalfLengthZ) + tolerance;
       if ((std::abs(cZ) > std::abs(hZ)) or
           ((vecLocal - cZ * eb).norm() >
-           m_bounds->get(LineBounds::eR) + s_onSurfaceTolerance)) {
+           m_bounds->get(LineBounds::eR) + tolerance)) {
         status = Intersection3D::Status::missed;
       }
     }

Core/src/Surfaces/PlaneSurface.cpp

@@ -172,12 +172,13 @@ double Acts::PlaneSurface::pathCorrection(const GeometryContext& gctx,
 
 Acts::SurfaceIntersection Acts::PlaneSurface::intersect(
     const GeometryContext& gctx, const Vector3& position,
-    const Vector3& direction, const BoundaryCheck& bcheck) const {
+    const Vector3& direction, const BoundaryCheck& bcheck,
+    ActsScalar tolerance) const {
   // Get the contextual transform
   const auto& gctxTransform = transform(gctx);
   // Use the intersection helper for planar surfaces
   auto intersection =
-      PlanarHelper::intersect(gctxTransform, position, direction);
+      PlanarHelper::intersect(gctxTransform, position, direction, tolerance);
   // Evaluate boundary check if requested (and reachable)
   if (intersection.status != Intersection3D::Status::unreachable and bcheck) {
     // Built-in local to global for speed reasons

Tests/UnitTests/Core/Propagator/NavigatorTests.cpp

@@ -208,6 +208,8 @@ struct PropagatorState {
     size_t debugMsgWidth = 50;
 
     const Acts::Logger& logger = Acts::getDummyLogger();
+
+    ActsScalar targetTolerance = s_onSurfaceTolerance;
   };
 
   /// Navigation cache: the start surface

Tests/UnitTests/Core/Surfaces/SurfaceStub.hpp

@@ -83,7 +83,8 @@ class SurfaceStub : public Surface {
   SurfaceIntersection intersect(const GeometryContext& /*gctx*/,
                                 const Vector3& /*position*/,
                                 const Vector3& /*direction*/,
-                                const BoundaryCheck& /*bcheck*/) const final {
+                                const BoundaryCheck& /*bcheck*/,
+                                const ActsScalar /*tolerance*/) const final {
     Intersection3D stubIntersection(Vector3(20., 0., 0.), 20.,
                                     Intersection3D::Status::reachable);
     return SurfaceIntersection(stubIntersection, this);