feat: Max chi2 for outliers in Core MeasurementSelector (#3475)

- Add `chi2CutOffOutlier` which allows to cut for holes
- Refactor cut selection
- Refactor return logic

This should allow testing #3474 properly in ACTS standalone

blocked by
- #3474

CI/physmon/workflows/physmon_trackfinding_1muon.py

@@ -141,6 +141,9 @@ def run_ckf_tracking(label, seeding):
                 maxOutliers=2,
             ),
             CkfConfig(
+                chi2CutOffMeasurement=15.0,
+                chi2CutOffOutlier=25.0,
+                numMeasurementsCutOff=10,
                 seedDeduplication=(
                     True if seeding != SeedingAlgorithm.TruthSmeared else False
                 ),

CI/physmon/workflows/physmon_trackfinding_4muon_50vertices.py

@@ -123,6 +123,9 @@
             maxOutliers=2,
         ),
         CkfConfig(
+            chi2CutOffMeasurement=15.0,
+            chi2CutOffOutlier=25.0,
+            numMeasurementsCutOff=10,
             seedDeduplication=True,
             stayOnSeed=True,
         ),

CI/physmon/workflows/physmon_trackfinding_ttbar_pu200.py

@@ -117,6 +117,9 @@
             maxOutliers=2,
         ),
         CkfConfig(
+            chi2CutOffMeasurement=15.0,
+            chi2CutOffOutlier=25.0,
+            numMeasurementsCutOff=10,
             seedDeduplication=True,
             stayOnSeed=True,
         ),

Core/include/Acts/TrackFinding/MeasurementSelector.hpp

@@ -36,10 +36,12 @@ namespace Acts {
 struct MeasurementSelectorCuts {
   /// bins in |eta| to specify variable selections
   std::vector<double> etaBins{};
-  /// Maximum local chi2 contribution.
+  /// Maximum local chi2 contribution to classify as measurement.
   std::vector<double> chi2CutOff{15};
   /// Maximum number of associated measurements on a single surface.
   std::vector<std::size_t> numMeasurementsCutOff{1};
+  /// Maximum local chi2 contribution to classify as outlier.
+  std::vector<double> chi2CutOffOutlier{};
 };
 
 /// @brief Measurement selection struct selecting those measurements compatible
@@ -73,7 +75,7 @@ class MeasurementSelector {
   /// @brief Constructor with config
   ///
   /// @param config a config instance
-  explicit MeasurementSelector(Config config);
+  explicit MeasurementSelector(const Config& config);
 
   /// @brief Function that select the measurements compatible with
   /// the given track parameter on a surface
@@ -93,11 +95,25 @@ class MeasurementSelector {
          bool& isOutlier, const Logger& logger) const;
 
  private:
-  template <typename traj_t, typename cut_value_t>
-  static cut_value_t getCut(
-      const typename traj_t::TrackStateProxy& trackState,
-      const Acts::MeasurementSelector::Config::Iterator selector,
-      const std::vector<cut_value_t>& cuts, const Logger& logger);
+  struct InternalCutBin {
+    double maxTheta{};
+    std::size_t maxNumMeasurements{};
+    double maxChi2Measurement{};
+    double maxChi2Outlier{};
+  };
+  using InternalCutBins = std::vector<InternalCutBin>;
+  using InternalConfig = Acts::GeometryHierarchyMap<InternalCutBins>;
+
+  struct Cuts {
+    std::size_t numMeasurements{};
+    double chi2Measurement{};
+    double chi2Outlier{};
+  };
+
+  static InternalCutBins convertCutBins(const MeasurementSelectorCuts& config);
+
+  static Cuts getCutsByTheta(const InternalCutBins& config, double theta);
+  Result<Cuts> getCuts(const GeometryIdentifier& geoID, double theta) const;
 
   double calculateChi2(
       const double* fullCalibrated, const double* fullCalibratedCovariance,
@@ -109,7 +125,7 @@ class MeasurementSelector {
                        false>::Projector projector,
       unsigned int calibratedSize) const;
 
-  Config m_config;
+  InternalConfig m_config;
 };
 
 }  // namespace Acts

Core/include/Acts/TrackFinding/MeasurementSelector.ipp

@@ -27,25 +27,17 @@ MeasurementSelector::select(
   }
 
   // Get geoID of this surface
-  auto geoID = candidates.front().referenceSurface().geometryId();
+  GeometryIdentifier geoID = candidates.front().referenceSurface().geometryId();
+  // Get the theta of the first track state
+  const double theta = candidates.front().predicted()[eBoundTheta];
   // Find the appropriate cuts
-  auto cuts = m_config.find(geoID);
-  if (cuts == m_config.end()) {
-    // for now we consider missing cuts an unrecoverable error
-    // TODO consider other options e.g. do not add measurements at all (not
-    // even as outliers)
-    return CombinatorialKalmanFilterError::MeasurementSelectionFailed;
+  const auto cutsResult = getCuts(geoID, theta);
+  if (!cutsResult.ok()) {
+    return cutsResult.error();
   }
+  const Cuts& cuts = *cutsResult;
 
-  assert(!cuts->chi2CutOff.empty());
-  const std::vector<double>& chi2CutOff = cuts->chi2CutOff;
-  const double maxChi2Cut =
-      std::min(*std::max_element(chi2CutOff.begin(), chi2CutOff.end()),
-               getCut<traj_t>(candidates.front(), cuts, chi2CutOff, logger));
-  const std::size_t numMeasurementsCut = getCut<traj_t>(
-      candidates.front(), cuts, cuts->numMeasurementsCutOff, logger);
-
-  if (numMeasurementsCut == 0ul) {
+  if (cuts.numMeasurements == 0ul) {
     return CombinatorialKalmanFilterError::MeasurementSelectionFailed;
   }
 
@@ -80,7 +72,7 @@ MeasurementSelector::select(
     }
 
     // only consider track states which pass the chi2 cut
-    if (chi2 >= maxChi2Cut) {
+    if (chi2 >= cuts.chi2Measurement) {
       continue;
     }
 
@@ -93,17 +85,25 @@ MeasurementSelector::select(
     ++passedCandidates;
   }
 
-  // If there are no measurements below the chi2 cut off, return the
-  // measurement with the best chi2 and tag it as an outlier
+  // Handle if there are no measurements below the chi2 cut off
   if (passedCandidates == 0ul) {
-    ACTS_VERBOSE("No measurement candidate. Return an outlier measurement chi2="
-                 << minChi2);
-    isOutlier = true;
+    if (minChi2 < cuts.chi2Outlier) {
+      ACTS_VERBOSE(
+          "No measurement candidate. Return an outlier measurement chi2="
+          << minChi2);
+      isOutlier = true;
+      return Result::success(std::make_pair(candidates.begin() + minIndex,
+                                            candidates.begin() + minIndex + 1));
+    } else {
+      ACTS_VERBOSE("No measurement candidate. Return empty chi2=" << minChi2);
+      return Result::success(
+          std::make_pair(candidates.begin(), candidates.begin()));
+    }
   }
 
-  if (passedCandidates <= 1ul || numMeasurementsCut == 1ul) {
+  if (passedCandidates <= 1ul) {
     // return single item range, no sorting necessary
-    ACTS_VERBOSE("Returning only 1 element");
+    ACTS_VERBOSE("Returning only 1 element chi2=" << minChi2);
     return Result::success(std::make_pair(candidates.begin() + minIndex,
                                           candidates.begin() + minIndex + 1));
   }
@@ -112,43 +112,12 @@ MeasurementSelector::select(
       candidates.begin(), candidates.begin() + passedCandidates,
       [](const auto& tsa, const auto& tsb) { return tsa.chi2() < tsb.chi2(); });
 
-  if (passedCandidates <= numMeasurementsCut) {
-    ACTS_VERBOSE("Number of selected measurements: "
-                 << passedCandidates << ", max: " << numMeasurementsCut);
-    return Result::success(std::make_pair(
-        candidates.begin(), candidates.begin() + passedCandidates));
-  }
-
   ACTS_VERBOSE("Number of selected measurements: "
-               << numMeasurementsCut << ", max: " << numMeasurementsCut);
+               << passedCandidates << ", max: " << cuts.numMeasurements);
 
   return Result::success(std::make_pair(
-      candidates.begin(), candidates.begin() + numMeasurementsCut));
-}
-
-template <typename traj_t, typename cut_value_t>
-cut_value_t MeasurementSelector::getCut(
-    const typename traj_t::TrackStateProxy& trackState,
-    const Acts::MeasurementSelector::Config::Iterator selector,
-    const std::vector<cut_value_t>& cuts, const Logger& logger) {
-  const auto& etaBins = selector->etaBins;
-  if (etaBins.empty()) {
-    return cuts[0];  // shortcut if no etaBins
-  }
-  const auto eta = std::atanh(std::cos(trackState.predicted()[eBoundTheta]));
-  const auto abseta = std::abs(eta);
-  std::size_t bin = 0;
-  for (auto etaBin : etaBins) {
-    if (etaBin >= abseta) {
-      break;
-    }
-    bin++;
-  }
-  if (bin >= cuts.size()) {
-    bin = cuts.size() - 1;
-  }
-  ACTS_VERBOSE("Get cut for eta=" << eta << ": " << cuts[bin]);
-  return cuts[bin];
+      candidates.begin(),
+      candidates.begin() + std::min(cuts.numMeasurements, passedCandidates)));
 }
 
 }  // namespace Acts

Core/include/Acts/TrackFinding/TrackSelector.hpp

@@ -75,6 +75,7 @@ class TrackSelector {
     std::size_t minMeasurements = 0;
     std::size_t maxHoles = std::numeric_limits<std::size_t>::max();
     std::size_t maxOutliers = std::numeric_limits<std::size_t>::max();
+    std::size_t maxHolesAndOutliers = std::numeric_limits<std::size_t>::max();
     std::size_t maxSharedHits = std::numeric_limits<std::size_t>::max();
     double maxChi2 = inf;
 
@@ -311,6 +312,7 @@ inline std::ostream& operator<<(std::ostream& os,
   print("pt", cuts.ptMin, cuts.ptMax);
   print("nHoles", 0, cuts.maxHoles);
   print("nOutliers", 0, cuts.maxOutliers);
+  print("nHoles + nOutliers", 0, cuts.maxHolesAndOutliers);
   print("nSharedHits", 0, cuts.maxSharedHits);
   print("chi2", 0.0, cuts.maxChi2);
   os << " - " << cuts.minMeasurements << " <= nMeasurements\n";
@@ -434,6 +436,8 @@ bool TrackSelector::isValidTrack(const track_proxy_t& track) const {
          checkMin(track.nMeasurements(), cuts.minMeasurements) &&
          checkMax(track.nHoles(), cuts.maxHoles) &&
          checkMax(track.nOutliers(), cuts.maxOutliers) &&
+         checkMax(track.nHoles() + track.nOutliers(),
+                  cuts.maxHolesAndOutliers) &&
          checkMax(track.nSharedHits(), cuts.maxSharedHits) &&
          checkMax(track.chi2(), cuts.maxChi2) &&
          cuts.measurementCounter.isValidTrack(track);

Core/src/TrackFinding/MeasurementSelector.cpp

@@ -10,19 +10,64 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/EventData/MeasurementHelpers.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
 
 #include <algorithm>
+#include <cstddef>
+#include <limits>
 
 namespace Acts {
 
 MeasurementSelector::MeasurementSelector()
-    : m_config{{GeometryIdentifier(), MeasurementSelectorCuts{}}} {}
+    : MeasurementSelector({{GeometryIdentifier(), MeasurementSelectorCuts{}}}) {
+}
 
 MeasurementSelector::MeasurementSelector(const MeasurementSelectorCuts& cuts)
-    : m_config{{GeometryIdentifier(), cuts}} {}
+    : MeasurementSelector({{GeometryIdentifier(), cuts}}) {}
+
+MeasurementSelector::MeasurementSelector(const Config& config) {
+  std::vector<InternalConfig::InputElement> tmp;
+  tmp.reserve(config.size());
+  for (std::size_t i = 0; i < config.size(); ++i) {
+    GeometryIdentifier geoID = config.idAt(i);
+    MeasurementSelectorCuts cuts = config.valueAt(i);
+
+    InternalCutBins internalCuts = convertCutBins(cuts);
+    tmp.emplace_back(geoID, std::move(internalCuts));
+  }
+  m_config = InternalConfig(std::move(tmp));
+}
+
+MeasurementSelector::InternalCutBins MeasurementSelector::convertCutBins(
+    const MeasurementSelectorCuts& config) {
+  InternalCutBins cutBins;
+
+  auto getEtaOrInf = [](const auto& vec, std::size_t bin) {
+    if (bin >= vec.size()) {
+      return std::numeric_limits<double>::infinity();
+    }
+    assert(vec[bin] >= 0 && "Eta bins must be positive");
+    return vec[bin];
+  };
+
+  auto getBinOrBackOrMax = [](const auto& vec, std::size_t bin) {
+    using Value = std::remove_reference_t<decltype(vec[0])>;
+    static constexpr Value max = std::numeric_limits<Value>::max();
+    return vec.empty() ? max : (bin < vec.size() ? vec[bin] : vec.back());
+  };
 
-MeasurementSelector::MeasurementSelector(Config config)
-    : m_config(std::move(config)) {}
+  for (std::size_t bin = 0; bin < config.etaBins.size() + 1; ++bin) {
+    InternalCutBin cuts;
+    cuts.maxTheta = getEtaOrInf(config.etaBins, bin);
+    cuts.maxNumMeasurements =
+        getBinOrBackOrMax(config.numMeasurementsCutOff, bin);
+    cuts.maxChi2Measurement = getBinOrBackOrMax(config.chi2CutOff, bin);
+    cuts.maxChi2Outlier = getBinOrBackOrMax(config.chi2CutOffOutlier, bin);
+    cutBins.push_back(cuts);
+  }
+
+  return cutBins;
+}
 
 double MeasurementSelector::calculateChi2(
     const double* fullCalibrated, const double* fullCalibratedCovariance,
@@ -64,4 +109,31 @@ double MeasurementSelector::calculateChi2(
       });
 }
 
+MeasurementSelector::Cuts MeasurementSelector::getCutsByTheta(
+    const InternalCutBins& config, double theta) {
+  // since theta is in [0, pi] and we have a symmetric cut in eta, we can just
+  // look at the positive half of the Z axis
+  const double constrainedTheta = std::min(theta, M_PI - theta);
+
+  auto it = std::find_if(config.begin(), config.end(),
+                         [constrainedTheta](const InternalCutBin& cuts) {
+                           return constrainedTheta < cuts.maxTheta;
+                         });
+  assert(it != config.end());
+  return {it->maxNumMeasurements, it->maxChi2Measurement, it->maxChi2Outlier};
+}
+
+Result<MeasurementSelector::Cuts> MeasurementSelector::getCuts(
+    const GeometryIdentifier& geoID, double theta) const {
+  // Find the appropriate cuts
+  auto cuts = m_config.find(geoID);
+  if (cuts == m_config.end()) {
+    // for now we consider missing cuts an unrecoverable error
+    // TODO consider other options e.g. do not add measurements at all (not
+    // even as outliers)
+    return CombinatorialKalmanFilterError::MeasurementSelectionFailed;
+  }
+  return getCutsByTheta(*cuts, theta);
+}
+
 }  // namespace Acts

Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp

@@ -241,8 +241,10 @@ class BranchStopper {
     bool tooManyHoles =
         track.nHoles() > singleConfig->maxHoles || tooManyHolesPS;
     bool tooManyOutliers = track.nOutliers() > singleConfig->maxOutliers;
+    bool tooManyHolesAndOutliers = (track.nHoles() + track.nOutliers()) >
+                                   singleConfig->maxHolesAndOutliers;
 
-    if (tooManyHoles || tooManyOutliers) {
+    if (tooManyHoles || tooManyOutliers || tooManyHolesAndOutliers) {
       ++m_nStoppedBranches;
       return enoughMeasurements ? BranchStopperResult::StopAndKeep
                                 : BranchStopperResult::StopAndDrop;