refactor: Refactor SurfaceReached aborter

Core/include/Acts/Propagator/MaterialInteractor.hpp

@@ -12,6 +12,7 @@
 #include "Acts/Geometry/TrackingVolume.hpp"
 #include "Acts/Material/MaterialInteraction.hpp"
 #include "Acts/Material/MaterialSlab.hpp"
+#include "Acts/Propagator/Propagator.hpp"
 #include "Acts/Propagator/detail/PointwiseMaterialInteraction.hpp"
 #include "Acts/Propagator/detail/VolumeMaterialInteraction.hpp"
 #include "Acts/Surfaces/Surface.hpp"
@@ -55,74 +56,85 @@ struct MaterialInteractor {
   void operator()(propagator_state_t& state, const stepper_t& stepper,
                   const navigator_t& navigator, result_type& result,
                   const Logger& logger) const {
-    // In case of Volume material update the result of the previous step
-    if (recordInteractions && !result.materialInteractions.empty() &&
-        !result.materialInteractions.back().volume.empty() &&
-        result.materialInteractions.back().updatedVolumeStep == false) {
-      updateResult(state, stepper, result);
-    }
-
-    // If we are on target, everything should have been done
-    if (navigator.targetReached(state.navigation)) {
+    if (state.stage == PropagatorStage::postPropagation) {
       return;
     }
+
     // Do nothing if nothing is what is requested.
     if (!(multipleScattering || energyLoss || recordInteractions)) {
       return;
     }
-    // We only have material interactions if there is potential material
-    const Surface* surface = navigator.currentSurface(state.navigation);
-    const TrackingVolume* volume = navigator.currentVolume(state.navigation);
 
-    if (!(surface && surface->surfaceMaterial()) &&
-        !(volume && volume->volumeMaterial())) {
-      return;
-    }
+    // Handle surface material
 
+    // Note that start and target surface conditions are handled in the
+    // interaction code
+    const Surface* surface = navigator.currentSurface(state.navigation);
+
+    // We only have material interactions if there is potential material
     if (surface && surface->surfaceMaterial()) {
+      ACTS_VERBOSE("MaterialInteractor | "
+                   << "Found material on surface " << surface->geometryId());
+
       // Prepare relevant input particle properties
       detail::PointwiseMaterialInteraction d(surface, state, stepper);
 
       // Determine the effective traversed material and its properties
       // Material exists but it's not real, i.e. vacuum; there is nothing to do
-      if (!d.evaluateMaterialSlab(state, navigator)) {
-        return;
+      if (d.evaluateMaterialSlab(state, navigator)) {
+        // Evaluate the material effects
+        d.evaluatePointwiseMaterialInteraction(multipleScattering, energyLoss);
+
+        if (energyLoss) {
+          using namespace UnitLiterals;
+          ACTS_VERBOSE("MaterialInteractor | "
+                       << d.slab << " absPdg=" << d.absPdg
+                       << " mass=" << d.mass / 1_MeV << "MeV"
+                       << " momentum=" << d.momentum / 1_GeV << "GeV"
+                       << " energyloss=" << d.Eloss / 1_MeV << "MeV");
+        }
+
+        // To integrate process noise, we need to transport
+        // the covariance to the current position in space
+        if (d.performCovarianceTransport) {
+          stepper.transportCovarianceToCurvilinear(state.stepping);
+        }
+        // Change the noise updater depending on the navigation direction
+        NoiseUpdateMode mode = (state.options.direction == Direction::Forward)
+                                   ? addNoise
+                                   : removeNoise;
+        // Apply the material interactions
+        d.updateState(state, stepper, mode);
+
+        // Record the result
+        recordResult(d, result);
       }
+    }
 
-      // Evaluate the material effects
-      d.evaluatePointwiseMaterialInteraction(multipleScattering, energyLoss);
+    // Handle volume material
 
-      if (energyLoss) {
-        using namespace UnitLiterals;
-        ACTS_VERBOSE(d.slab << " absPdg=" << d.absPdg
-                            << " mass=" << d.mass / 1_MeV << "MeV"
-                            << " momentum=" << d.momentum / 1_GeV << "GeV"
-                            << " energyloss=" << d.Eloss / 1_MeV << "MeV");
-      }
+    // In case of Volume material update the result of the previous step
+    if (!result.materialInteractions.empty() &&
+        !result.materialInteractions.back().volume.empty() &&
+        result.materialInteractions.back().updatedVolumeStep == false) {
+      updateResult(state, stepper, result);
+    }
+
+    const TrackingVolume* volume = navigator.currentVolume(state.navigation);
+
+    // We only have material interactions if there is potential material
+    if (volume && volume->volumeMaterial()) {
+      ACTS_VERBOSE("MaterialInteractor | "
+                   << "Found material in volume " << volume->geometryId());
 
-      // To integrate process noise, we need to transport
-      // the covariance to the current position in space
-      if (d.performCovarianceTransport) {
-        stepper.transportCovarianceToCurvilinear(state.stepping);
-      }
-      // Change the noise updater depending on the navigation direction
-      NoiseUpdateMode mode = (state.options.direction == Direction::Forward)
-                                 ? addNoise
-                                 : removeNoise;
-      // Apply the material interactions
-      d.updateState(state, stepper, mode);
-      // Record the result
-      recordResult(d, result);
-    } else if (recordInteractions && volume && volume->volumeMaterial()) {
       // Prepare relevant input particle properties
       detail::VolumeMaterialInteraction d(volume, state, stepper);
       // Determine the effective traversed material and its properties
       // Material exists but it's not real, i.e. vacuum; there is nothing to do
-      if (!d.evaluateMaterialSlab(state, navigator)) {
-        return;
+      if (d.evaluateMaterialSlab(state, navigator)) {
+        // Record the result
+        recordResult(d, result);
       }
-      // Record the result
-      recordResult(d, result);
     }
   }
 
@@ -159,16 +171,18 @@ struct MaterialInteractor {
   /// @param [in, out] result Result storage
   void recordResult(const detail::VolumeMaterialInteraction& d,
                     result_type& result) const {
-    // Record the interaction
-    MaterialInteraction mi;
-    mi.position = d.pos;
-    mi.time = d.time;
-    mi.direction = d.dir;
-    mi.surface = nullptr;
-    mi.volume = d.volume;
-    mi.pathCorrection = d.pathCorrection;
-    mi.materialSlab = d.slab;
-    result.materialInteractions.push_back(std::move(mi));
+    // Record the interaction if requested
+    if (recordInteractions) {
+      MaterialInteraction mi;
+      mi.position = d.pos;
+      mi.time = d.time;
+      mi.direction = d.dir;
+      mi.surface = nullptr;
+      mi.volume = d.volume;
+      mi.pathCorrection = d.pathCorrection;
+      mi.materialSlab = d.slab;
+      result.materialInteractions.push_back(std::move(mi));
+    }
   }
 
   /// @brief This function update the previous material step
@@ -179,19 +193,21 @@ struct MaterialInteractor {
   template <typename propagator_state_t, typename stepper_t>
   void updateResult(propagator_state_t& state, const stepper_t& stepper,
                     result_type& result) const {
-    // Update the previous interaction
-    Vector3 shift = stepper.position(state.stepping) -
-                    result.materialInteractions.back().position;
-    double momentum = stepper.direction(state.stepping).norm();
-    result.materialInteractions.back().deltaP =
-        momentum - result.materialInteractions.back().direction.norm();
-    result.materialInteractions.back().materialSlab.scaleThickness(
-        shift.norm());
-    result.materialInteractions.back().updatedVolumeStep = true;
-    result.materialInX0 +=
-        result.materialInteractions.back().materialSlab.thicknessInX0();
-    result.materialInL0 +=
-        result.materialInteractions.back().materialSlab.thicknessInL0();
+    // Update the previous interaction if requested
+    if (recordInteractions) {
+      Vector3 shift = stepper.position(state.stepping) -
+                      result.materialInteractions.back().position;
+      double momentum = stepper.direction(state.stepping).norm();
+      result.materialInteractions.back().deltaP =
+          momentum - result.materialInteractions.back().direction.norm();
+      result.materialInteractions.back().materialSlab.scaleThickness(
+          shift.norm());
+      result.materialInteractions.back().updatedVolumeStep = true;
+      result.materialInX0 +=
+          result.materialInteractions.back().materialSlab.thicknessInX0();
+      result.materialInL0 +=
+          result.materialInteractions.back().materialSlab.thicknessInL0();
+    }
   }
 };
 

Core/include/Acts/Propagator/MultiStepperAborters.hpp

@@ -13,10 +13,7 @@
 
 namespace Acts {
 
-/// This
-struct MultiStepperSurfaceReached {
-  MultiStepperSurfaceReached() = default;
-
+struct MultiStepperSurfaceReached : public SurfaceReached {
   /// If this is set, we are also happy if the mean of the components is on the
   /// surface. How the averaging is performed depends on the stepper
   /// implementation
@@ -27,6 +24,9 @@ struct MultiStepperSurfaceReached {
   /// false
   double averageOnSurfaceTolerance = 0.2;
 
+  MultiStepperSurfaceReached() = default;
+  MultiStepperSurfaceReached(double oLimit) : SurfaceReached(oLimit) {}
+
   /// boolean operator for abort condition without using the result
   ///
   /// @tparam propagator_state_t Type of the propagator state
@@ -41,56 +41,27 @@ struct MultiStepperSurfaceReached {
             typename navigator_t>
   bool operator()(propagator_state_t& state, const stepper_t& stepper,
                   const navigator_t& navigator, const Logger& logger) const {
-    return (*this)(state, stepper, navigator,
-                   *navigator.targetSurface(state.navigation), logger);
-  }
-
-  /// boolean operator for abort condition without using the result
-  ///
-  /// @tparam propagator_state_t Type of the propagator state
-  /// @tparam stepper_t Type of the stepper
-  /// @tparam navigator_t Type of the navigator
-  ///
-  /// @param [in,out] state The propagation state object
-  /// @param [in] stepper Stepper used for the propagation
-  /// @param [in] navigator Navigator used for the propagation
-  /// @param [in] targetSurface The target surface
-  /// @param logger a logger instance
-  template <typename propagator_state_t, typename stepper_t,
-            typename navigator_t>
-  bool operator()(propagator_state_t& state, const stepper_t& stepper,
-                  const navigator_t& navigator, const Surface& targetSurface,
-                  const Logger& logger) const {
-    bool reached = true;
-    const auto oldCurrentSurface = navigator.currentSurface(state.navigation);
-
-    for (auto cmp : stepper.componentIterable(state.stepping)) {
-      auto singleState = cmp.singleState(state);
-      const auto& singleStepper = cmp.singleStepper(stepper);
-
-      if (!SurfaceReached{}(singleState, singleStepper, navigator,
-                            targetSurface, logger)) {
-        reached = false;
-      } else {
-        cmp.status() = Acts::Intersection3D::Status::onSurface;
-      }
+    if (surface == nullptr) {
+      ACTS_VERBOSE(
+          "MultiStepperSurfaceReached aborter | "
+          "No target surface set.");
+      return false;
     }
 
     // However, if mean of all is on surface, we are happy as well
     if (averageOnSurface) {
       const auto sIntersection =
-          targetSurface
-              .intersect(
+          surface
+              ->intersect(
                   state.geoContext, stepper.position(state.stepping),
                   state.options.direction * stepper.direction(state.stepping),
-                  BoundaryCheck(true), averageOnSurfaceTolerance)
+                  BoundaryCheck(boundaryCheck), 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);
-
+        ACTS_VERBOSE(
+            "MultiStepperSurfaceReached aborter | "
+            "Reached target in average mode");
         for (auto cmp : stepper.componentIterable(state.stepping)) {
           cmp.status() = Intersection3D::Status::onSurface;
         }
@@ -99,14 +70,23 @@ struct MultiStepperSurfaceReached {
       }
     }
 
-    // These values are changed by the single component aborters but must be
-    // reset if not all components are on the target
-    if (!reached) {
-      navigator.currentSurface(state.navigation, oldCurrentSurface);
-      navigator.targetReached(state.navigation, false);
+    bool reached = true;
+
+    for (auto cmp : stepper.componentIterable(state.stepping)) {
+      // note that this is not copying anything heavy
+      auto singleState = cmp.singleState(state);
+      const auto& singleStepper = cmp.singleStepper(stepper);
+
+      if (!SurfaceReached::operator()(singleState, singleStepper, navigator,
+                                      logger)) {
+        reached = false;
+      } else {
+        cmp.status() = Acts::Intersection3D::Status::onSurface;
+      }
     }
 
     return reached;
   }
 };
+
 }  // namespace Acts

Core/include/Acts/Propagator/Propagator.ipp

@@ -230,6 +230,7 @@ auto Acts::Propagator<S, N>::propagate(
 
   // Type of provided options
   target_aborter_t targetAborter;
+  targetAborter.surface = &target;
   path_aborter_t pathAborter;
   pathAborter.internalLimit = options.pathLimit;
   auto abortList = options.abortList.append(targetAborter, pathAborter);