feat: Add ability for grouped hit requirements to track selector (#2871)

With this PR, `TrackSelector` allows selecting tracks based on minimum hit requirements where hits are required to match a set of geometry ID masks. With this pattern, you can implement **minimum hit requirements per subdetector**, i.e. at least 3 hits in the Pixel detector, etc.

It uses `GeometryIdHierarchyMap` to implement the matching, and then combines hit counts from all matches in the map.

From python, you would set this up like:

```python
measurementCounter = acts.TrackSelector.MeasurementCounter()
# At least 2 hits in volume 3
measurementCounter.addCounter([acts.GeometryIdentifier().setVolume(3)], 2)
```

and then set it to `TrackSelectorConfig`:

```python
TrackSelectorConfig(
  pt=(1.0 * u.GeV if ttbar else 0.0, None),
  absEta=(None, 3.0),
  loc0=(-4.0 * u.mm, 4.0 * u.mm),
  nMeasurementsMin=7,
  nMeasurementsGroupMin=measurementCounter
)
```

You can also select e.g. layer IDs, where due to the OR behavior within a single map, you can model cuts like "at least 2 hits on the two innermost pixel layers".

Core/include/Acts/TrackFinding/TrackSelector.hpp

@@ -8,12 +8,19 @@
 
 #pragma once
 
+#include "Acts/EventData/TrackStateType.hpp"
+#include "Acts/Geometry/GeometryHierarchyMap.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+
 #include <cmath>
 #include <functional>
 #include <limits>
+#include <numeric>
 #include <ostream>
 #include <vector>
 
+#include <boost/container/small_vector.hpp>
+
 namespace Acts {
 
 /// Class which performs filtering of tracks. It accepts an input and an output
@@ -23,6 +30,26 @@ class TrackSelector {
   static constexpr double inf = std::numeric_limits<double>::infinity();
 
  public:
+  struct MeasurementCounter {
+    // Combination of a geometry hierarchy map and a minimum hit count
+    using CounterElement =
+        std::pair<GeometryHierarchyMap<unsigned int>, unsigned int>;
+
+    boost::container::small_vector<CounterElement, 4> counters;
+
+    template <typename track_proxy_t>
+    bool isValidTrack(const track_proxy_t& track) const;
+
+    void addCounter(const std::vector<GeometryIdentifier>& identifiers,
+                    unsigned int threshold) {
+      std::vector<GeometryHierarchyMap<unsigned int>::InputElement> elements;
+      for (const auto& id : identifiers) {
+        elements.emplace_back(id, 0);
+      }
+      counters.emplace_back(std::move(elements), threshold);
+    }
+  };
+
   /// Configuration of a set of cuts for a single eta bin
   /// Default construction yields a set of cuts that accepts everything.
   struct Config {
@@ -51,6 +78,9 @@ class TrackSelector {
     std::size_t maxSharedHits = std::numeric_limits<std::size_t>::max();
     double maxChi2 = inf;
 
+    // Defaults to: no cut
+    MeasurementCounter measurementCounter;
+
     // Helper factory functions to produce a populated config object more
     // conveniently
 
@@ -136,7 +166,8 @@ class TrackSelector {
 
     /// Auto-converting constructor from a single cut configuration.
     /// Results in a single absolute eta bin from 0 to infinity.
-    EtaBinnedConfig(Config cutSet) : cutSets{cutSet}, absEtaEdges{{0, inf}} {}
+    EtaBinnedConfig(Config cutSet)
+        : cutSets{std::move(cutSet)}, absEtaEdges{{0, inf}} {}
 
     /// Add a new eta bin with the given upper bound.
     /// @param etaMax Upper bound of the new eta bin
@@ -395,7 +426,8 @@ bool TrackSelector::isValidTrack(const track_proxy_t& track) const {
          checkMax(track.nHoles(), cuts.maxHoles) &&
          checkMax(track.nOutliers(), cuts.maxOutliers) &&
          checkMax(track.nSharedHits(), cuts.maxSharedHits) &&
-         checkMax(track.chi2(), cuts.maxChi2);
+         checkMax(track.chi2(), cuts.maxChi2) &&
+         cuts.measurementCounter.isValidTrack(track);
 }
 
 inline TrackSelector::TrackSelector(
@@ -433,4 +465,43 @@ inline TrackSelector::TrackSelector(
 inline TrackSelector::TrackSelector(const Config& config)
     : TrackSelector{EtaBinnedConfig{config}} {}
 
+template <typename track_proxy_t>
+bool TrackSelector::MeasurementCounter::isValidTrack(
+    const track_proxy_t& track) const {
+  // No hit cuts, accept everything
+  if (counters.empty()) {
+    return true;
+  }
+
+  boost::container::small_vector<unsigned int, 4> counterValues;
+  counterValues.resize(counters.size(), 0);
+
+  for (const auto& ts : track.trackStatesReversed()) {
+    if (!ts.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) {
+      continue;
+    }
+
+    const auto geoId = ts.referenceSurface().geometryId();
+
+    for (std::size_t i = 0; i < counters.size(); i++) {
+      const auto& [counterMap, threshold] = counters[i];
+      const auto it = counterMap.find(geoId);
+      if (it == counterMap.end()) {
+        continue;
+      }
+
+      counterValues[i]++;
+    }
+  }
+
+  for (std::size_t i = 0; i < counters.size(); i++) {
+    const auto& [counterMap, threshold] = counters[i];
+    const unsigned int value = counterValues[i];
+    if (value < threshold) {
+      return false;
+    }
+  }
+
+  return true;
+}
 }  // namespace Acts

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

@@ -113,10 +113,21 @@
     defaults=[(None, None)],
 )
 
+
 TrackSelectorConfig = namedtuple(
     "TrackSelectorConfig",
-    ["loc0", "loc1", "time", "eta", "absEta", "pt", "phi", "nMeasurementsMin"],
-    defaults=[(None, None)] * 7 + [None],
+    [
+        "loc0",
+        "loc1",
+        "time",
+        "eta",
+        "absEta",
+        "pt",
+        "phi",
+        "nMeasurementsMin",
+        "nMeasurementsGroupMin",
+    ],
+    defaults=[(None, None)] * 7 + [None] * 2,
 )
 
 CkfConfig = namedtuple(
@@ -1177,6 +1188,7 @@ def addCKFTracks(
                 ptMin=trackSelectorConfig.pt[0],
                 ptMax=trackSelectorConfig.pt[1],
                 minMeasurements=trackSelectorConfig.nMeasurementsMin,
+                measurementCounter=trackSelectorConfig.nMeasurementsGroupMin,
             )
         )
         if trackSelectorConfig is not None

Examples/Python/src/ExampleAlgorithms.cpp

@@ -80,6 +80,14 @@ void addExampleAlgorithms(Context& ctx) {
                     .def(py::init<const Config&>(), py::arg("config"))
                     .def(py::init<const EtaBinnedConfig&>(), py::arg("config"));
 
+    {
+      auto mc = py::class_<Acts::TrackSelector::MeasurementCounter>(
+                    tool, "MeasurementCounter")
+                    .def(py::init<>())
+                    .def("addCounter",
+                         &Acts::TrackSelector::MeasurementCounter::addCounter);
+    }
+
     {
       auto c = py::class_<Config>(tool, "Config").def(py::init<>());
 
@@ -101,6 +109,7 @@ void addExampleAlgorithms(Context& ctx) {
       ACTS_PYTHON_MEMBER(ptMin);
       ACTS_PYTHON_MEMBER(ptMax);
       ACTS_PYTHON_MEMBER(minMeasurements);
+      ACTS_PYTHON_MEMBER(measurementCounter);
       ACTS_PYTHON_STRUCT_END();
 
       pythonRangeProperty(c, "loc0", &Config::loc0Min, &Config::loc0Max);

Tests/UnitTests/Core/TrackFinding/TrackSelectorTests.cpp

@@ -9,6 +9,13 @@
 #include <boost/test/data/test_case.hpp>
 #include <boost/test/unit_test.hpp>
 
+#include "Acts/EventData/MultiTrajectory.hpp"
+#include "Acts/EventData/TrackContainer.hpp"
+#include "Acts/EventData/TrackStateType.hpp"
+#include "Acts/EventData/VectorMultiTrajectory.hpp"
+#include "Acts/EventData/VectorTrackContainer.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/Surfaces/PerigeeSurface.hpp"
 #include "Acts/TrackFinding/TrackSelector.hpp"
 
 #include <limits>
@@ -41,6 +48,32 @@ struct MockTrack {
   std::size_t nOutliers() const { return m_nOutliers; }
   std::size_t nSharedHits() const { return m_nSharedHits; }
   float chi2() const { return m_chi2; }
+
+  // To comply with concept, not actually used
+ private:
+  struct MockTrackState {
+    const Surface& referenceSurface() const {
+      static const auto srf =
+          Surface::makeShared<PlaneSurface>(Vector3::Zero(), Vector3::UnitZ());
+      return *srf;
+    }
+
+    ConstTrackStateType typeFlags() const {
+      static const ConstTrackStateType::raw_type raw{0};
+      return {raw};
+    }
+  };
+
+  struct TrackStateRange {
+    auto begin() const { return m_trackStates.begin(); }
+    auto end() const { return m_trackStates.end(); }
+
+   private:
+    std::vector<MockTrackState> m_trackStates;
+  };
+
+ public:
+  TrackStateRange trackStatesReversed() const { return {}; }
 };
 
 double thetaFromEta(double eta) {
@@ -571,4 +604,141 @@ BOOST_AUTO_TEST_CASE(TestConstructor) {
   }
 }
 
+BOOST_AUTO_TEST_CASE(SubsetHitCountCut) {
+  auto makeSurface = [](GeometryIdentifier id) {
+    auto srf =
+        Surface::makeShared<PlaneSurface>(Vector3::Zero(), Vector3::UnitZ());
+
+    srf->assignGeometryId(id);
+    return srf;
+  };
+
+  auto addTrackState = [](auto& track, const auto& surface,
+                          TrackStateFlag flag) {
+    auto ts = track.appendTrackState();
+    ts.setReferenceSurface(surface);
+    ts.typeFlags().set(flag);
+    return ts;
+  };
+
+  auto addMeasurement = [&](auto& track, const auto& surface) {
+    return addTrackState(track, surface, TrackStateFlag::MeasurementFlag);
+  };
+
+  auto addMaterial = [&](auto& track, const auto& surface) {
+    return addTrackState(track, surface, TrackStateFlag::MaterialFlag);
+  };
+
+  TrackContainer tc{VectorTrackContainer{}, VectorMultiTrajectory{}};
+
+  auto makeTrack = [&]() {
+    auto track = tc.makeTrack();
+
+    using namespace Acts::UnitLiterals;
+    track.parameters() << 0, 0, M_PI / 2, M_PI / 2, 1 / 1_GeV, 0;
+    auto perigee = Surface::makeShared<PerigeeSurface>(Vector3::Zero());
+    track.setReferenceSurface(perigee);
+    return track;
+  };
+
+  auto vol7_lay3_sen2 = makeSurface(
+      GeometryIdentifier{}.setVolume(7).setLayer(3).setSensitive(2));
+  auto vol7_lay4 = makeSurface(GeometryIdentifier{}.setVolume(7).setLayer(4));
+  auto vol7_lay3_sen8 = makeSurface(
+      GeometryIdentifier{}.setVolume(7).setLayer(3).setSensitive(8));
+  auto vol7_lay5_sen11 = makeSurface(
+      GeometryIdentifier{}.setVolume(7).setLayer(5).setSensitive(11));
+  auto vol7_lay5_sen12 = makeSurface(
+      GeometryIdentifier{}.setVolume(7).setLayer(5).setSensitive(12));
+  auto vol7_lay6_sen3 = makeSurface(
+      GeometryIdentifier{}.setVolume(7).setLayer(6).setSensitive(3));
+
+  auto vol8_lay8_sen1 = makeSurface(
+      GeometryIdentifier{}.setVolume(8).setLayer(8).setSensitive(1));
+  auto vol8_lay8_sen2 = makeSurface(
+      GeometryIdentifier{}.setVolume(8).setLayer(8).setSensitive(2));
+  auto vol8_lay9_sen1 = makeSurface(
+      GeometryIdentifier{}.setVolume(8).setLayer(9).setSensitive(1));
+
+  TrackSelector::Config cfgVol7;
+  cfgVol7.measurementCounter.addCounter({GeometryIdentifier{}.setVolume(7)}, 3);
+  TrackSelector selectorVol7{cfgVol7};
+
+  auto trackVol7 = makeTrack();
+
+  BOOST_CHECK(!selectorVol7.isValidTrack(trackVol7));
+
+  // 1 hit in vol7
+  addMeasurement(trackVol7, vol7_lay3_sen2);
+  addMaterial(trackVol7, vol7_lay4);
+
+  BOOST_CHECK(!selectorVol7.isValidTrack(trackVol7));
+  addMeasurement(trackVol7, vol7_lay5_sen11);
+  BOOST_CHECK(!selectorVol7.isValidTrack(trackVol7));
+
+  // Now we should have enough hits
+  addMeasurement(trackVol7, vol7_lay6_sen3);
+  BOOST_CHECK(selectorVol7.isValidTrack(trackVol7));
+
+  TrackSelector::Config cfgVol8;
+  cfgVol8.measurementCounter.addCounter({GeometryIdentifier{}.setVolume(8)}, 2);
+  TrackSelector selectorVol8{cfgVol8};
+
+  // Previous trackVol7 has no measurements in volume 8
+  BOOST_CHECK(!selectorVol8.isValidTrack(trackVol7));
+
+  auto trackVol8 = makeTrack();
+  BOOST_CHECK(!selectorVol8.isValidTrack(trackVol8));
+
+  addMeasurement(trackVol8, vol8_lay8_sen1);
+  BOOST_CHECK(!selectorVol8.isValidTrack(trackVol8));
+  addMeasurement(trackVol8, vol8_lay8_sen2);
+  BOOST_CHECK(selectorVol8.isValidTrack(trackVol8));
+  addMeasurement(trackVol8, vol8_lay9_sen1);
+  BOOST_CHECK(selectorVol8.isValidTrack(trackVol8));
+
+  TrackSelector::Config cfgVol7Lay5;
+  cfgVol7Lay5.measurementCounter.addCounter(
+      {GeometryIdentifier{}.setVolume(7).setLayer(5)}, 2);
+  TrackSelector selectorVol7Lay5{cfgVol7Lay5};
+
+  // Only one hit on volume 7 layer 5
+  BOOST_CHECK(!selectorVol7Lay5.isValidTrack(trackVol7));
+  addMeasurement(trackVol7, vol7_lay5_sen12);
+  BOOST_CHECK(selectorVol7Lay5.isValidTrack(trackVol7));
+
+  // Check requirement on volume 7 OR 8
+  TrackSelector::Config cfgVol7Or8;
+  cfgVol7Or8.measurementCounter.addCounter(
+      {GeometryIdentifier{}.setVolume(7), GeometryIdentifier{}.setVolume(8)},
+      4);
+  TrackSelector selectorVol7Or8{cfgVol7Or8};
+
+  // threshold is 4
+  // this track has enough hits in volume 7 only
+  BOOST_CHECK(selectorVol7Or8.isValidTrack(trackVol7));
+  // this track does not have enough hits in volume 8 only
+  BOOST_CHECK(!selectorVol7Or8.isValidTrack(trackVol8));
+
+  // add 1 hit in volume 7 to push it over the threshold
+  addMeasurement(trackVol8, vol7_lay3_sen8);
+  // now it passes
+  BOOST_CHECK(selectorVol7Or8.isValidTrack(trackVol8));
+
+  TrackSelector::Config cfgVol7And8;
+  cfgVol7And8.measurementCounter.addCounter({GeometryIdentifier{}.setVolume(7)},
+                                            4);
+  cfgVol7And8.measurementCounter.addCounter({GeometryIdentifier{}.setVolume(8)},
+                                            2);
+  TrackSelector selectorVol7And8{cfgVol7And8};
+
+  // this track has enough hits in vol 7 but not enough in vol 8
+  BOOST_CHECK(!selectorVol7And8.isValidTrack(trackVol7));
+
+  addMeasurement(trackVol7, vol8_lay8_sen1);
+  addMeasurement(trackVol7, vol8_lay8_sen2);
+
+  BOOST_CHECK(selectorVol7And8.isValidTrack(trackVol7));
+}
+
 BOOST_AUTO_TEST_SUITE_END()