Merge branch 'main' into feat-detector-python-infrastructure

Core/include/Acts/Navigation/SurfaceCandidatesUpdators.hpp

@@ -45,8 +45,9 @@ inline static void updateCandidates(const GeometryContext& gctx,
         (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);
+    // TODO surface tolerance
+    auto sIntersection = sRep.intersect(gctx, position, direction,
+                                        c.boundaryCheck, s_onSurfaceTolerance);
     // Re-order and swap if necessary
     if (sIntersection.intersection.pathLength + s_onSurfaceTolerance <
             nState.overstepTolerance and

Core/include/Acts/Propagator/AtlasStepper.hpp

@@ -74,10 +74,10 @@ class AtlasStepper {
       const auto pos = pars.position(gctx);
       const auto Vp = pars.parameters();
 
-      double Sf = sin(Vp[eBoundPhi]);
-      double Cf = cos(Vp[eBoundPhi]);
-      double Se = sin(Vp[eBoundTheta]);
-      double Ce = cos(Vp[eBoundTheta]);
+      double Sf = std::sin(Vp[eBoundPhi]);
+      double Cf = std::cos(Vp[eBoundPhi]);
+      double Se = std::sin(Vp[eBoundTheta]);
+      double Ce = std::cos(Vp[eBoundTheta]);
 
       pVector[0] = pos[ePos0];
       pVector[1] = pos[ePos1];
@@ -168,8 +168,8 @@ class AtlasStepper {
         const auto& surface = pars.referenceSurface();
         // the disc needs polar coordinate adaptations
         if (surface.type() == Surface::Disc) {
-          double lCf = cos(Vp[1]);
-          double lSf = sin(Vp[1]);
+          double lCf = std::cos(Vp[1]);
+          double lSf = std::sin(Vp[1]);
           double Ax[3] = {transform(0, 0), transform(1, 0), transform(2, 0)};
           double Ay[3] = {transform(0, 1), transform(1, 1), transform(2, 1)};
           double d0 = lCf * Ax[0] + lSf * Ay[0];
@@ -385,11 +385,13 @@ class AtlasStepper {
   /// @param [in] surface The surface provided
   /// @param [in] bcheck The boundary check for this status update
   /// @param [in] logger Logger instance to use
+  /// @param [in] surfaceTolerance Surface tolerance used for intersection
   Intersection3D::Status updateSurfaceStatus(
       State& state, const Surface& surface, const BoundaryCheck& bcheck,
-      const Logger& logger = getDummyLogger()) const {
+      const Logger& logger = getDummyLogger(),
+      ActsScalar surfaceTolerance = s_onSurfaceTolerance) const {
     return detail::updateSingleSurfaceStatus<AtlasStepper>(
-        *this, state, surface, bcheck, logger);
+        *this, state, surface, bcheck, logger, surfaceTolerance);
   }
 
   /// Update step size
@@ -568,10 +570,10 @@ class AtlasStepper {
     Vector3 pos(state.pVector[0], state.pVector[1], state.pVector[2]);
     Vector3 mom(state.pVector[4], state.pVector[5], state.pVector[6]);
 
-    double Sf = sin(boundParams[eBoundPhi]);
-    double Cf = cos(boundParams[eBoundPhi]);
-    double Se = sin(boundParams[eBoundTheta]);
-    double Ce = cos(boundParams[eBoundTheta]);
+    double Sf = std::sin(boundParams[eBoundPhi]);
+    double Cf = std::cos(boundParams[eBoundPhi]);
+    double Se = std::sin(boundParams[eBoundTheta]);
+    double Ce = std::cos(boundParams[eBoundTheta]);
 
     const auto transform = surface.referenceFrame(state.geoContext, pos, mom);
 
@@ -639,8 +641,8 @@ class AtlasStepper {
     // special treatment for surface types
     // the disc needs polar coordinate adaptations
     if (surface.type() == Surface::Disc) {
-      double lCf = cos(boundParams[eBoundLoc1]);
-      double lSf = sin(boundParams[eBoundLoc1]);
+      double lCf = std::cos(boundParams[eBoundLoc1]);
+      double lSf = std::sin(boundParams[eBoundLoc1]);
       double Ax[3] = {transform(0, 0), transform(1, 0), transform(2, 0)};
       double Ay[3] = {transform(0, 1), transform(1, 1), transform(2, 1)};
       double d0 = lCf * Ax[0] + lSf * Ay[0];
@@ -1220,7 +1222,7 @@ class AtlasStepper {
           2. * h *
           (std::abs((A1 + A6) - (A3 + A4)) + std::abs((B1 + B6) - (B3 + B4)) +
            std::abs((C1 + C6) - (C3 + C4)));
-      if (EST > state.options.tolerance) {
+      if (std::abs(EST) > std::abs(state.options.tolerance)) {
         h = h * .5;
         state.stepping.stepSize.setValue(h);
         //        dltm = 0.;

Core/include/Acts/Propagator/DirectNavigator.hpp

@@ -222,7 +222,8 @@ class DirectNavigator {
     if (state.navigation.navSurfaceIter != state.navigation.navSurfaces.end()) {
       // Establish & update the surface status
       auto surfaceStatus = stepper.updateSurfaceStatus(
-          state.stepping, **state.navigation.navSurfaceIter, false);
+          state.stepping, **state.navigation.navSurfaceIter, false, *m_logger,
+          state.options.targetTolerance);
       if (surfaceStatus == Intersection3D::Status::unreachable) {
         ACTS_VERBOSE(
             "Surface not reachable anymore, switching to next one in "
@@ -269,7 +270,8 @@ class DirectNavigator {
     if (state.navigation.navSurfaceIter != state.navigation.navSurfaces.end()) {
       // Establish the surface status
       auto surfaceStatus = stepper.updateSurfaceStatus(
-          state.stepping, **state.navigation.navSurfaceIter, false);
+          state.stepping, **state.navigation.navSurfaceIter, false, *m_logger,
+          state.options.targetTolerance);
       if (surfaceStatus == Intersection3D::Status::onSurface) {
         // Set the current surface
         state.navigation.currentSurface = *state.navigation.navSurfaceIter;

Core/include/Acts/Propagator/EigenStepper.hpp

@@ -11,6 +11,8 @@
 // Workaround for building on clang+libstdc++
 #include "Acts/Utilities/detail/ReferenceWrapperAnyCompat.hpp"
 
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/Tolerance.hpp"
 #include "Acts/Definitions/Units.hpp"
 #include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
@@ -68,20 +70,17 @@ class EigenStepper {
     /// @param [in] par The track parameters at start
     /// @param [in] ndir The navigation direction w.r.t momentum
     /// @param [in] ssize is the maximum step size
-    /// @param [in] stolerance is the stepping tolerance
     ///
     /// @note the covariance matrix is copied when needed
     template <typename charge_t>
     explicit State(const GeometryContext& gctx,
                    MagneticFieldProvider::Cache fieldCacheIn,
                    const SingleBoundTrackParameters<charge_t>& par,
                    Direction ndir = Direction::Forward,
-                   double ssize = std::numeric_limits<double>::max(),
-                   double stolerance = s_onSurfaceTolerance)
+                   double ssize = std::numeric_limits<double>::max())
         : absCharge(std::abs(par.charge())),
           navDir(ndir),
           stepSize(ndir * std::abs(ssize)),
-          tolerance(stolerance),
           fieldCache(std::move(fieldCacheIn)),
           geoContext(gctx) {
       pars.template segment<3>(eFreePos0) = par.position(gctx);
@@ -135,9 +134,6 @@ class EigenStepper {
     /// Last performed step (for overstep limit calculation)
     double previousStepSize = 0.;
 
-    /// The tolerance for the stepping
-    double tolerance = s_onSurfaceTolerance;
-
     /// This caches the current magnetic field cell and stays
     /// (and interpolates) within it as long as this is valid.
     /// See step() code for details.
@@ -172,8 +168,7 @@ class EigenStepper {
                   std::reference_wrapper<const MagneticFieldContext> mctx,
                   const SingleBoundTrackParameters<charge_t>& par,
                   Direction navDir = Direction::Forward,
-                  double ssize = std::numeric_limits<double>::max(),
-                  double stolerance = s_onSurfaceTolerance) const;
+                  double ssize = std::numeric_limits<double>::max()) const;
 
   /// @brief Resets the state
   ///
@@ -254,11 +249,13 @@ class EigenStepper {
   /// @param [in] surface The surface provided
   /// @param [in] bcheck The boundary check for this status update
   /// @param [in] logger A @c Logger instance
+  /// @param [in] surfaceTolerance Surface tolerance used for intersection
   Intersection3D::Status updateSurfaceStatus(
       State& state, const Surface& surface, const BoundaryCheck& bcheck,
-      const Logger& logger = getDummyLogger()) const {
+      const Logger& logger = getDummyLogger(),
+      ActsScalar surfaceTolerance = s_onSurfaceTolerance) const {
     return detail::updateSingleSurfaceStatus<EigenStepper>(
-        *this, state, surface, bcheck, logger);
+        *this, state, surface, bcheck, logger, surfaceTolerance);
   }
 
   /// Update step size

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -20,8 +20,8 @@ auto Acts::EigenStepper<E, A>::makeState(
     std::reference_wrapper<const GeometryContext> gctx,
     std::reference_wrapper<const MagneticFieldContext> mctx,
     const SingleBoundTrackParameters<charge_t>& par, Direction navDir,
-    double ssize, double stolerance) const -> State {
-  return State{gctx, m_bField->makeCache(mctx), par, navDir, ssize, stolerance};
+    double ssize) const -> State {
+  return State{gctx, m_bField->makeCache(mctx), par, navDir, ssize};
 }
 
 template <typename E, typename A>

Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp

@@ -299,7 +299,6 @@ class MultiEigenStepperLoop
     /// @param [in] multipars The track multi-component track-parameters at start
     /// @param [in] ndir The navigation direction w.r.t momentum
     /// @param [in] ssize is the maximum step size
-    /// @param [in] stolerance is the stepping tolerance
     ///
     /// @note the covariance matrix is copied when needed
     template <typename charge_t>
@@ -308,8 +307,7 @@ class MultiEigenStepperLoop
         const std::shared_ptr<const MagneticFieldProvider>& bfield,
         const MultiComponentBoundTrackParameters<charge_t>& multipars,
         Direction ndir = Direction::Forward,
-        double ssize = std::numeric_limits<double>::max(),
-        double stolerance = s_onSurfaceTolerance)
+        double ssize = std::numeric_limits<double>::max())
         : navDir(ndir), geoContext(gctx), magContext(mctx) {
       if (multipars.components().empty()) {
         throw std::invalid_argument(
@@ -321,10 +319,9 @@ class MultiEigenStepperLoop
 
       for (auto i = 0ul; i < multipars.components().size(); ++i) {
         const auto [weight, singlePars] = multipars[i];
-        components.push_back(
-            {SingleState(gctx, bfield->makeCache(mctx), std::move(singlePars),
-                         ndir, ssize, stolerance),
-             weight, Intersection3D::Status::onSurface});
+        components.push_back({SingleState(gctx, bfield->makeCache(mctx),
+                                          std::move(singlePars), ndir, ssize),
+                              weight, Intersection3D::Status::onSurface});
       }
 
       if (std::get<2>(multipars.components().front())) {
@@ -339,10 +336,8 @@ class MultiEigenStepperLoop
                   std::reference_wrapper<const MagneticFieldContext> mctx,
                   const MultiComponentBoundTrackParameters<charge_t>& par,
                   Direction navDir = Direction::Forward,
-                  double ssize = std::numeric_limits<double>::max(),
-                  double stolerance = s_onSurfaceTolerance) const {
-    return State(gctx, mctx, SingleStepper::m_bField, par, navDir, ssize,
-                 stolerance);
+                  double ssize = std::numeric_limits<double>::max()) const {
+    return State(gctx, mctx, SingleStepper::m_bField, par, navDir, ssize);
   }
 
   /// @brief Resets the state
@@ -690,17 +685,19 @@ class MultiEigenStepperLoop
   /// @param state [in,out] The stepping state (thread-local cache)
   /// @param surface [in] The surface provided
   /// @param bcheck [in] The boundary check for this status update
-  /// @param logger [in] A @c Loggerinstance
+  /// @param logger [in] A @c Logger instance
+  /// @param [in] surfaceTolerance Surface tolerance used for intersection
   Intersection3D::Status updateSurfaceStatus(
       State& state, const Surface& surface, const BoundaryCheck& bcheck,
-      const Logger& logger = getDummyLogger()) const {
+      const Logger& logger = getDummyLogger(),
+      ActsScalar surfaceTolerance = s_onSurfaceTolerance) const {
     using Status = Intersection3D::Status;
 
     std::array<int, 4> counts = {0, 0, 0, 0};
 
     for (auto& component : state.components) {
       component.status = detail::updateSingleSurfaceStatus<SingleStepper>(
-          *this, component.state, surface, bcheck, logger);
+          *this, component.state, surface, bcheck, logger, surfaceTolerance);
       ++counts[static_cast<std::size_t>(component.status)];
     }
 

Core/include/Acts/Propagator/MultiStepperAborters.hpp

@@ -80,11 +80,11 @@ struct MultiStepperSurfaceReached {
     if (averageOnSurface) {
       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,
+          averageOnSurfaceTolerance);
 
       if (sIntersection.intersection.status ==
-              Intersection3D::Status::onSurface or
-          sIntersection.intersection.pathLength < averageOnSurfaceTolerance) {
+          Intersection3D::Status::onSurface) {
         ACTS_VERBOSE("Reached target in average mode");
         navigator.currentSurface(state.navigation, &targetSurface);
         navigator.targetReached(state.navigation, true);

Core/include/Acts/Propagator/Navigator.hpp

@@ -623,7 +623,8 @@ class Navigator {
     // If we are on the surface pointed at by the index, we can make
     // it the current one to pass it to the other actors
     auto surfaceStatus =
-        stepper.updateSurfaceStatus(state.stepping, *surface, true, logger());
+        stepper.updateSurfaceStatus(state.stepping, *surface, true, logger(),
+                                    state.options.targetTolerance);
     if (surfaceStatus == Intersection3D::Status::onSurface) {
       ACTS_VERBOSE(volInfo(state)
                    << "Status Surface successfully hit, storing it.");
@@ -712,8 +713,9 @@ class Navigator {
           break;
         }
       }
-      auto surfaceStatus = stepper.updateSurfaceStatus(state.stepping, *surface,
-                                                       boundaryCheck, logger());
+      auto surfaceStatus =
+          stepper.updateSurfaceStatus(state.stepping, *surface, boundaryCheck,
+                                      logger(), state.options.targetTolerance);
       if (surfaceStatus == Intersection3D::Status::reachable) {
         ACTS_VERBOSE(volInfo(state)
                      << "Surface reachable, step size updated to "
@@ -876,8 +878,9 @@ class Navigator {
         }
       }
       // Try to step towards it
-      auto layerStatus = stepper.updateSurfaceStatus(
-          state.stepping, *layerSurface, true, logger());
+      auto layerStatus =
+          stepper.updateSurfaceStatus(state.stepping, *layerSurface, true,
+                                      logger(), state.options.targetTolerance);
       if (layerStatus == Intersection3D::Status::reachable) {
         ACTS_VERBOSE(volInfo(state) << "Layer reachable, step size updated to "
                                     << stepper.outputStepSize(state.stepping));
@@ -1017,8 +1020,9 @@ class Navigator {
       // That is the current boundary surface
       auto boundarySurface = state.navigation.navBoundary().representation;
       // Step towards the boundary surfrace
-      auto boundaryStatus = stepper.updateSurfaceStatus(
-          state.stepping, *boundarySurface, true, logger());
+      auto boundaryStatus =
+          stepper.updateSurfaceStatus(state.stepping, *boundarySurface, true,
+                                      logger(), state.options.targetTolerance);
       if (boundaryStatus == Intersection3D::Status::reachable) {
         ACTS_VERBOSE(volInfo(state)
                      << "Boundary reachable, step size updated to "
@@ -1169,7 +1173,7 @@ class Navigator {
         stepper.getStepSize(state.stepping, ConstrainedStep::aborter);
     // No overstepping on start layer, otherwise ask the stepper
     navOpts.overstepLimit = (cLayer != nullptr)
-                                ? s_onSurfaceTolerance
+                                ? state.options.targetTolerance
                                 : stepper.overstepLimit(state.stepping);
 
     // get the surfaces
@@ -1315,7 +1319,8 @@ class Navigator {
         return true;
       }
       auto targetStatus = stepper.updateSurfaceStatus(
-          state.stepping, *state.navigation.targetSurface, true, logger());
+          state.stepping, *state.navigation.targetSurface, true, logger(),
+          state.options.targetTolerance);
       // the only advance could have been to the target
       if (targetStatus == Intersection3D::Status::onSurface) {
         // set the target surface

Core/include/Acts/Propagator/Propagator.ipp

@@ -146,8 +146,7 @@ auto Acts::Propagator<S, N>::propagate(
   StateType state{
       eOptions,
       m_stepper.makeState(eOptions.geoContext, eOptions.magFieldContext, start,
-                          eOptions.direction, eOptions.maxStepSize,
-                          eOptions.tolerance),
+                          eOptions.direction, eOptions.maxStepSize),
       m_navigator.makeState(&start.referenceSurface(), nullptr)};
 
   static_assert(
@@ -233,8 +232,7 @@ auto Acts::Propagator<S, N>::propagate(
   StateType state{
       eOptions,
       m_stepper.makeState(eOptions.geoContext, eOptions.magFieldContext, start,
-                          eOptions.direction, eOptions.maxStepSize,
-                          eOptions.tolerance),
+                          eOptions.direction, eOptions.maxStepSize),
       m_navigator.makeState(&start.referenceSurface(), &target)};
 
   static_assert(

Core/include/Acts/Propagator/StepperConcept.hpp

@@ -128,7 +128,7 @@ constexpr bool MultiStepperStateConcept= require<
         constexpr static bool covariance_transport_exists = require<has_method<const S, void, covariance_transport_curvilinear_t, state&>,
                                                                     has_method<const S, void, covariance_transport_bound_t, state&, const Surface&, const FreeToBoundCorrection&>>;
         static_assert(covariance_transport_exists, "covarianceTransport method not found");
-        constexpr static bool update_surface_exists = has_method<const S, Intersection3D::Status, update_surface_status_t, state&, const Surface&, const BoundaryCheck&, const Logger&>;
+        constexpr static bool update_surface_exists = has_method<const S, Intersection3D::Status, update_surface_status_t, state&, const Surface&, const BoundaryCheck&, const Logger&, ActsScalar>;
         static_assert(update_surface_exists, "updateSurfaceStatus method not found");
         constexpr static bool set_step_size_exists = has_method<const S, void, set_step_size_t, state&, double, ConstrainedStep::Type, bool>;
         static_assert(set_step_size_exists, "setStepSize method not found");

Core/include/Acts/Propagator/StraightLineStepper.hpp

@@ -244,11 +244,13 @@ class StraightLineStepper {
   /// @param [in] surface The surface provided
   /// @param [in] bcheck The boundary check for this status update
   /// @param [in] logger A logger instance
+  /// @param [in] surfaceTolerance Surface tolerance used for intersection
   Intersection3D::Status updateSurfaceStatus(
       State& state, const Surface& surface, const BoundaryCheck& bcheck,
-      const Logger& logger = getDummyLogger()) const {
+      const Logger& logger = getDummyLogger(),
+      ActsScalar surfaceTolerance = s_onSurfaceTolerance) const {
     return detail::updateSingleSurfaceStatus<StraightLineStepper>(
-        *this, state, surface, bcheck, logger);
+        *this, state, surface, bcheck, logger, surfaceTolerance);
   }
 
   /// Update step size