feat: Allow volume constrain for propagation (#3470)

Add optional volume constrains to the propagation which are executed by the `VolumeConstraintAborter`. This is directly exercised in the track finding and wired to python.

This is generally useful if track finding / fitting should be done in a sub detector region. Same goes for simulation and extrapolation.

Core/include/Acts/Propagator/PropagatorOptions.hpp

@@ -33,13 +33,22 @@ struct PurePropagatorPlainOptions {
   /// Absolute maximum path length
   double pathLimit = std::numeric_limits<double>::max();
 
-  /// Required tolerance to reach surface
-  double surfaceTolerance = s_onSurfaceTolerance;
-
   /// Loop protection step, it adapts the pathLimit
   bool loopProtection = true;
   /// Allowed loop fraction, 1 is a full loop
   double loopFraction = 0.5;
+
+  /// Required tolerance to reach surface
+  double surfaceTolerance = s_onSurfaceTolerance;
+
+  /// Constrain the propagation to selected volumes
+  /// @note ignored if empty
+  /// @note requires `VolumeConstraintAborter` aborter
+  std::vector<std::uint32_t> constrainToVolumeIds;
+  /// Additional volumes to be considered as end of world
+  /// @note ignored if empty
+  /// @note requires `VolumeConstraintAborter` aborter
+  std::vector<std::uint32_t> endOfWorldVolumeIds;
 };
 
 }  // namespace detail

Core/include/Acts/Propagator/StandardAborters.hpp

@@ -161,7 +161,7 @@ struct ForcedSurfaceReached : SurfaceReached {
       : SurfaceReached(std::numeric_limits<double>::lowest()) {}
 };
 
-/// This is the condition that the end of World has been reached
+/// This is the condition that the end of world has been reached
 /// it then triggers an propagation abort
 struct EndOfWorldReached {
   /// boolean operator for abort condition without using the result
@@ -181,6 +181,60 @@ struct EndOfWorldReached {
   }
 };
 
+/// This is the condition that the end of world has been reached
+/// it then triggers a propagation abort
+struct VolumeConstraintAborter {
+  /// boolean operator for abort condition without using the result
+  ///
+  /// @tparam propagator_state_t Type of the propagator state
+  /// @tparam navigator_t Type of the navigator
+  ///
+  /// @param [in,out] state The propagation state object
+  /// @param [in] navigator The navigator object
+  /// @param logger a logger instance
+  template <typename propagator_state_t, typename stepper_t,
+            typename navigator_t>
+  bool checkAbort(propagator_state_t& state, const stepper_t& /*stepper*/,
+                  const navigator_t& navigator, const Logger& logger) const {
+    const auto& constrainToVolumeIds = state.options.constrainToVolumeIds;
+    const auto& endOfWorldVolumeIds = state.options.endOfWorldVolumeIds;
+
+    if (constrainToVolumeIds.empty() && endOfWorldVolumeIds.empty()) {
+      return false;
+    }
+    const auto* currentVolume = navigator.currentVolume(state.navigation);
+
+    // We need a volume to check its ID
+    if (currentVolume == nullptr) {
+      return false;
+    }
+
+    const auto currentVolumeId =
+        static_cast<std::uint32_t>(currentVolume->geometryId().volume());
+
+    if (!constrainToVolumeIds.empty() &&
+        std::find(constrainToVolumeIds.begin(), constrainToVolumeIds.end(),
+                  currentVolumeId) == constrainToVolumeIds.end()) {
+      ACTS_VERBOSE(
+          "VolumeConstraintAborter aborter | Abort with volume constrain "
+          << currentVolumeId);
+      return true;
+    }
+
+    if (!endOfWorldVolumeIds.empty() &&
+        std::find(endOfWorldVolumeIds.begin(), endOfWorldVolumeIds.end(),
+                  currentVolumeId) != endOfWorldVolumeIds.end()) {
+      ACTS_VERBOSE(
+          "VolumeConstraintAborter aborter | Abort with additional end of "
+          "world volume "
+          << currentVolumeId);
+      return true;
+    }
+
+    return false;
+  }
+};
+
 /// Aborter that checks if the propagation has reached any surface
 struct AnySurfaceReached {
   template <typename propagator_state_t, typename stepper_t,

Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp

@@ -596,15 +596,19 @@ class CombinatorialKalmanFilter {
 
       const bool isEndOfWorldReached =
           endOfWorldReached.checkAbort(state, stepper, navigator, logger());
+      const bool isVolumeConstraintReached = volumeConstraintAborter.checkAbort(
+          state, stepper, navigator, logger());
       const bool isPathLimitReached = result.pathLimitReached.checkAbort(
           state, stepper, navigator, logger());
       const bool isTargetReached =
           targetReached.checkAbort(state, stepper, navigator, logger());
       const bool allBranchesStopped = result.activeBranches.empty();
       if (isEndOfWorldReached || isPathLimitReached || isTargetReached ||
-          allBranchesStopped) {
+          allBranchesStopped || isVolumeConstraintReached) {
         if (isEndOfWorldReached) {
           ACTS_VERBOSE("End of world reached");
+        } else if (isVolumeConstraintReached) {
+          ACTS_VERBOSE("Volume constraint reached");
         } else if (isPathLimitReached) {
           ACTS_VERBOSE("Path limit reached");
         } else if (isTargetReached) {
@@ -1173,6 +1177,9 @@ class CombinatorialKalmanFilter {
     /// End of world aborter
     EndOfWorldReached endOfWorldReached;
 
+    /// Volume constraint aborter
+    VolumeConstraintAborter volumeConstraintAborter;
+
     /// Actor logger instance
     const Logger* actorLogger{nullptr};
     /// Updater logger instance

Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp

@@ -130,12 +130,17 @@ class TrackFindingAlgorithm final : public IAlgorithm {
     bool trimTracks = true;
 
     // Pixel and strip volume ids to be used for maxPixel/StripHoles cuts
-    std::set<Acts::GeometryIdentifier::Value> pixelVolumes;
-    std::set<Acts::GeometryIdentifier::Value> stripVolumes;
+    std::vector<std::uint32_t> pixelVolumeIds;
+    std::vector<std::uint32_t> stripVolumeIds;
 
-    /// additional track selector settings
+    // additional track selector settings
     std::size_t maxPixelHoles = std::numeric_limits<std::size_t>::max();
     std::size_t maxStripHoles = std::numeric_limits<std::size_t>::max();
+
+    /// The volume ids to constrain the track finding to
+    std::vector<std::uint32_t> constrainToVolumeIds;
+    /// The volume ids to stop the track finding at
+    std::vector<std::uint32_t> endOfWorldVolumeIds;
   };
 
   /// Constructor of the track finding algorithm

Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp

@@ -217,16 +217,18 @@ class BranchStopper {
     }
 
     bool tooManyHolesPS = false;
-    if (!(m_cfg.pixelVolumes.empty() && m_cfg.stripVolumes.empty())) {
+    if (!(m_cfg.pixelVolumeIds.empty() && m_cfg.stripVolumeIds.empty())) {
       auto& branchState = branchStateAccessor(track);
       // count both holes and outliers as holes for pixel/strip counts
       if (trackState.typeFlags().test(Acts::TrackStateFlag::HoleFlag) ||
           trackState.typeFlags().test(Acts::TrackStateFlag::OutlierFlag)) {
-        if (m_cfg.pixelVolumes.contains(
-                trackState.referenceSurface().geometryId().volume())) {
+        auto volumeId = trackState.referenceSurface().geometryId().volume();
+        if (std::find(m_cfg.pixelVolumeIds.begin(), m_cfg.pixelVolumeIds.end(),
+                      volumeId) != m_cfg.pixelVolumeIds.end()) {
           ++branchState.nPixelHoles;
-        } else if (m_cfg.stripVolumes.contains(
-                       trackState.referenceSurface().geometryId().volume())) {
+        } else if (std::find(m_cfg.stripVolumeIds.begin(),
+                             m_cfg.stripVolumeIds.end(),
+                             volumeId) != m_cfg.stripVolumeIds.end()) {
           ++branchState.nStripHoles;
         }
       }
@@ -350,11 +352,15 @@ ProcessCode TrackFindingAlgorithm::execute(const AlgorithmContext& ctx) const {
   firstPropOptions.maxSteps = m_cfg.maxSteps;
   firstPropOptions.direction = m_cfg.reverseSearch ? Acts::Direction::Backward
                                                    : Acts::Direction::Forward;
+  firstPropOptions.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
+  firstPropOptions.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
 
   Acts::PropagatorPlainOptions secondPropOptions(ctx.geoContext,
                                                  ctx.magFieldContext);
   secondPropOptions.maxSteps = m_cfg.maxSteps;
   secondPropOptions.direction = firstPropOptions.direction.invert();
+  secondPropOptions.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
+  secondPropOptions.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
 
   // Set the CombinatorialKalmanFilter options
   TrackFinderOptions firstOptions(ctx.geoContext, ctx.magFieldContext,
@@ -379,6 +385,8 @@ ProcessCode TrackFindingAlgorithm::execute(const AlgorithmContext& ctx) const {
       logger().cloneWithSuffix("Propagator"));
 
   ExtrapolatorOptions extrapolationOptions(ctx.geoContext, ctx.magFieldContext);
+  extrapolationOptions.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
+  extrapolationOptions.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
 
   // Perform the track finding for all initial parameters
   ACTS_DEBUG("Invoke track finding with " << initialParameters.size()

Examples/Python/python/acts/examples/reconstruction.py

@@ -147,8 +147,24 @@
         "maxPixelHoles",
         "maxStripHoles",
         "trimTracks",
+        "constrainToVolumes",
+        "endOfWorldVolumes",
+    ],
+    defaults=[
+        15.0,
+        25.0,
+        10,
+        None,
+        None,
+        None,
+        None,
+        None,
+        None,
+        None,
+        None,
+        None,
+        None,
     ],
-    defaults=[15.0, 25.0, 10, None, None, None, None, None, None, None, None],
 )
 
 AmbiguityResolutionConfig = namedtuple(
@@ -1531,11 +1547,13 @@ def addCKFTracks(
             reverseSearch=reverseSearch,
             seedDeduplication=ckfConfig.seedDeduplication,
             stayOnSeed=ckfConfig.stayOnSeed,
-            pixelVolumes=ckfConfig.pixelVolumes,
-            stripVolumes=ckfConfig.stripVolumes,
+            pixelVolumeIds=ckfConfig.pixelVolumes,
+            stripVolumeIds=ckfConfig.stripVolumes,
             maxPixelHoles=ckfConfig.maxPixelHoles,
             maxStripHoles=ckfConfig.maxStripHoles,
             trimTracks=ckfConfig.trimTracks,
+            constrainToVolumeIds=ckfConfig.constrainToVolumes,
+            endOfWorldVolumeIds=ckfConfig.endOfWorldVolumes,
         ),
     )
     s.addAlgorithm(trackFinder)

Examples/Python/src/TrackFinding.cpp

@@ -337,11 +337,13 @@ void addTrackFinding(Context& ctx) {
     ACTS_PYTHON_MEMBER(reverseSearch);
     ACTS_PYTHON_MEMBER(seedDeduplication);
     ACTS_PYTHON_MEMBER(stayOnSeed);
-    ACTS_PYTHON_MEMBER(pixelVolumes);
-    ACTS_PYTHON_MEMBER(stripVolumes);
+    ACTS_PYTHON_MEMBER(pixelVolumeIds);
+    ACTS_PYTHON_MEMBER(stripVolumeIds);
     ACTS_PYTHON_MEMBER(maxPixelHoles);
     ACTS_PYTHON_MEMBER(maxStripHoles);
     ACTS_PYTHON_MEMBER(trimTracks);
+    ACTS_PYTHON_MEMBER(constrainToVolumeIds);
+    ACTS_PYTHON_MEMBER(endOfWorldVolumeIds);
     ACTS_PYTHON_STRUCT_END();
   }
 

Examples/Scripts/Python/full_chain_itk.py

@@ -126,8 +126,8 @@
         seedDeduplication=True,
         stayOnSeed=True,
         # ITk volumes from Noemi's plot
-        pixelVolumes={8, 9, 10, 13, 14, 15, 16, 18, 19, 20},
-        stripVolumes={22, 23, 24},
+        pixelVolumes=[8, 9, 10, 13, 14, 15, 16, 18, 19, 20],
+        stripVolumes=[22, 23, 24],
         maxPixelHoles=1,
         maxStripHoles=2,
     ),

Examples/Scripts/Python/full_chain_odd.py

@@ -363,10 +363,27 @@
             numMeasurementsCutOff=10,
             seedDeduplication=True,
             stayOnSeed=True,
-            pixelVolumes={16, 17, 18},
-            stripVolumes={23, 24, 25},
+            pixelVolumes=[16, 17, 18],
+            stripVolumes=[23, 24, 25],
             maxPixelHoles=1,
             maxStripHoles=2,
+            constrainToVolumes=[
+                2,  # beam pipe
+                32,
+                4,  # beam pip gap
+                16,
+                17,
+                18,  # pixel
+                20,  # PST
+                23,
+                24,
+                25,  # short strip
+                26,
+                8,  # long strip gap
+                28,
+                29,
+                30,  # long strip
+            ],
         ),
         outputDirRoot=outputDir if args.output_root else None,
         outputDirCsv=outputDir if args.output_csv else None,

Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp

@@ -12,6 +12,7 @@
 #include "Acts/Definitions/PdgParticle.hpp"
 #include "Acts/Definitions/Units.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Geometry/TrackingVolume.hpp"
 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
 #include "Acts/Material/MaterialSlab.hpp"
 #include "Acts/Propagator/ConstrainedStep.hpp"
@@ -148,6 +149,11 @@ struct MockNavigator {
     return state.currentSurface;
   }
 
+  const Acts::TrackingVolume *currentVolume(
+      const MockNavigatorState & /*state*/) const {
+    return nullptr;
+  }
+
   bool endOfWorldReached(const MockNavigatorState & /*state*/) const {
     return false;
   }
@@ -158,6 +164,10 @@ struct MockPropagatorState {
   MockStepperState stepping;
   Acts::GeometryContext geoContext;
   Acts::PropagatorStage stage = Acts::PropagatorStage::invalid;
+
+  struct {
+    std::vector<std::uint32_t> constrainToVolumeIds;
+  } options;
 };
 
 template <typename SurfaceSelector>