Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-04-11 07:46:41

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008 
0009 #pragma once
0010 
0011 #include "Acts/EventData/TrackProxyConcept.hpp"
0012 #include "Acts/Geometry/GeometryHierarchyMap.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Utilities/AngleHelpers.hpp"
0015 
0016 #include <algorithm>
0017 #include <cmath>
0018 #include <functional>
0019 #include <limits>
0020 #include <ostream>
0021 #include <vector>
0022 
0023 #include <boost/container/small_vector.hpp>
0024 
0025 namespace Acts {
0026 
0027 /// Class which performs filtering of tracks. It accepts an input and an output
0028 /// track container and uses the built-in copy facility to copy tracks into the
0029 /// output container.
0030 class TrackSelector {
0031   static constexpr double inf = std::numeric_limits<double>::infinity();
0032 
0033  public:
0034   /// Measurement counting rules for geometry-specific requirements.
0035   struct MeasurementCounter {
0036     /// Type combining a geometry hierarchy map and a minimum hit count
0037     using CounterElement =
0038         std::pair<GeometryHierarchyMap<unsigned int>, unsigned int>;
0039 
0040     /// Collection of counter elements
0041     boost::container::small_vector<CounterElement, 4> counters;
0042 
0043     /// Check if track satisfies all measurement requirements
0044     /// @param track The track to check
0045     /// @return True if track satisfies all counter thresholds
0046     template <TrackProxyConcept track_proxy_t>
0047     bool isValidTrack(const track_proxy_t& track) const;
0048 
0049     /// Add a new counter with threshold for specified geometry
0050     /// @param identifiers Geometry identifiers to count measurements in
0051     /// @param threshold Minimum number of required measurements
0052     void addCounter(const std::vector<GeometryIdentifier>& identifiers,
0053                     unsigned int threshold) {
0054       std::vector<GeometryHierarchyMap<unsigned int>::InputElement> elements;
0055       for (const auto& id : identifiers) {
0056         elements.emplace_back(id, 0);
0057       }
0058       counters.emplace_back(std::move(elements), threshold);
0059     }
0060   };
0061 
0062   /// Configuration of a set of cuts for a single eta bin
0063   /// Default construction yields a set of cuts that accepts everything.
0064   struct Config {
0065     // Minimum/maximum local positions.
0066     /// Minimum local position in first coordinate
0067     double loc0Min = -inf;
0068     /// Maximum local position in first coordinate
0069     double loc0Max = inf;
0070     /// Minimum local position in second coordinate
0071     double loc1Min = -inf;
0072     /// Maximum local position in second coordinate
0073     double loc1Max = inf;
0074     // Minimum/maximum track time.
0075     /// Minimum track time cut
0076     double timeMin = -inf;
0077     /// Maximum track time cut
0078     double timeMax = inf;
0079     // Direction cuts.
0080     /// Minimum phi cut for track selection
0081     double phiMin = -inf;
0082     /// Maximum phi cut for track selection
0083     double phiMax = inf;
0084     /// Minimum eta cut for track selection
0085     double etaMin = -inf;
0086     /// Maximum eta cut for track selection
0087     double etaMax = inf;
0088     /// Minimum absolute eta cut for track selection
0089     double absEtaMin = 0.0;
0090     /// Maximum absolute eta cut for track selection
0091     double absEtaMax = inf;
0092     // Momentum cuts.
0093     /// Minimum transverse momentum cut
0094     double ptMin = 0.0;
0095     /// Maximum transverse momentum cut
0096     double ptMax = inf;
0097 
0098     /// Minimum number of measurements required
0099     std::size_t minMeasurements = 0;
0100     /// Maximum number of holes allowed
0101     std::size_t maxHoles = std::numeric_limits<std::size_t>::max();
0102     /// Maximum number of outliers allowed
0103     std::size_t maxOutliers = std::numeric_limits<std::size_t>::max();
0104     /// Maximum number of holes and outliers combined
0105     std::size_t maxHolesAndOutliers = std::numeric_limits<std::size_t>::max();
0106     /// Maximum number of shared hits allowed
0107     std::size_t maxSharedHits = std::numeric_limits<std::size_t>::max();
0108     /// Maximum chi-squared cut for track selection
0109     double maxChi2 = inf;
0110 
0111     /// Whether a reference surface is required for the track
0112     /// If false, the parameter cuts are not evaluated
0113     bool requireReferenceSurface = true;
0114 
0115     // Defaults to: no cut
0116     /// Counter for geometry-specific measurement requirements
0117     MeasurementCounter measurementCounter;
0118 
0119     // Helper factory functions to produce a populated config object more
0120     // conveniently
0121 
0122     /// Set loc0 acceptance range
0123     /// @param min Minimum value
0124     /// @param max Maximum value
0125     /// @return Reference to this object
0126     Config& loc0(double min, double max);
0127 
0128     /// Set loc1 acceptance range
0129     /// @param min Minimum value
0130     /// @param max Maximum value
0131     /// @return Reference to this object
0132     Config& loc1(double min, double max);
0133 
0134     /// Set time acceptance range
0135     /// @param min Minimum value
0136     /// @param max Maximum value
0137     /// @return Reference to this object
0138     Config& time(double min, double max);
0139 
0140     /// Set phi acceptance range
0141     /// @param min Minimum value
0142     /// @param max Maximum value
0143     /// @return Reference to this object
0144     Config& phi(double min, double max);
0145 
0146     /// Set the eta acceptance range
0147     /// @param min Minimum value
0148     /// @param max Maximum value
0149     /// @return Reference to this object
0150     Config& eta(double min, double max);
0151 
0152     /// Set the absolute eta acceptance range
0153     /// @param min Minimum value
0154     /// @param max Maximum value
0155     /// @return Reference to this object
0156     Config& absEta(double min, double max);
0157 
0158     /// Set the pt acceptance range
0159     /// @param min Minimum value
0160     /// @param max Maximum value
0161     /// @return Reference to this object
0162     Config& pt(double min, double max);
0163 
0164     /// Print this set of cuts to an output stream
0165     /// @param os Output stream
0166     /// @param cuts Cuts to print
0167     /// @return Reference to the output stream
0168     friend std::ostream& operator<<(std::ostream& os, const Config& cuts);
0169   };
0170 
0171   /// Main config object for the track selector. Combines a set of cut
0172   /// configurations and corresponding eta bins
0173   struct EtaBinnedConfig {
0174     /// Cut sets for each eta bin
0175     std::vector<Config> cutSets = {};
0176 
0177     /// Eta bin edges for varying cuts by eta
0178     std::vector<double> absEtaEdges = {0, inf};
0179 
0180     /// Get the number of eta bins
0181     /// @return Number of eta bins
0182     std::size_t nEtaBins() const { return absEtaEdges.size() - 1; }
0183 
0184     /// Construct an empty (accepts everything) configuration.
0185     /// Results in a single cut set and one abs eta bin from 0 to infinity.
0186     EtaBinnedConfig() : cutSets{{}} {};
0187 
0188     /// Constructor to create a config object that is not upper-bounded.
0189     /// This is useful to use the "fluent" API to populate the configuration.
0190     /// @param etaMin Minimum eta bin edge
0191     explicit EtaBinnedConfig(double etaMin) : cutSets{}, absEtaEdges{etaMin} {}
0192 
0193     /// Constructor from a vector of eta bin edges. This automatically
0194     /// initializes all the cuts to be the same for all eta and be essentially
0195     /// no-op.
0196     /// @param absEtaEdgesIn is the vector of eta bin edges
0197     explicit EtaBinnedConfig(std::vector<double> absEtaEdgesIn)
0198         : absEtaEdges{std::move(absEtaEdgesIn)} {
0199       cutSets.resize(nEtaBins());
0200     }
0201 
0202     /// Auto-converting constructor from a single cut configuration.
0203     /// Results in a single absolute eta bin from 0 to infinity.
0204     /// @param cutSet Single cut configuration to use for all eta values
0205     explicit EtaBinnedConfig(Config cutSet) : cutSets{std::move(cutSet)} {}
0206 
0207     /// Add a new eta bin with the given upper bound.
0208     /// @param etaMax Upper bound of the new eta bin
0209     /// @param callback Callback to configure the cuts for this eta bin
0210     /// @return Reference to this object
0211     EtaBinnedConfig& addCuts(double etaMax,
0212                              const std::function<void(Config&)>& callback = {});
0213 
0214     /// Add a new eta bin with an upper bound of +infinity.
0215     /// @param callback Callback to configure the cuts for this eta bin
0216     /// @return Reference to this object
0217     EtaBinnedConfig& addCuts(const std::function<void(Config&)>& callback = {});
0218 
0219     /// Print this configuration to an output stream
0220     /// @param os Output stream
0221     /// @param cfg Configuration to print
0222     /// @return Reference to the output stream
0223     friend std::ostream& operator<<(std::ostream& os,
0224                                     const EtaBinnedConfig& cfg);
0225 
0226     /// Check if the configuration has a bin for a given eta
0227     /// @param eta Eta value
0228     /// @return True if the configuration has a bin for the given eta
0229     bool hasCuts(double eta) const;
0230 
0231     /// Get the index of the eta bin for a given eta.
0232     /// throws an exception if Eta is outside the abs eta bin edges.
0233     /// @param eta Eta value
0234     /// @return Index of the eta bin
0235     std::size_t binIndex(double eta) const;
0236 
0237     /// Get the index of the eta bin for a given eta
0238     /// @param eta Eta value
0239     /// @return Index of the eta bin, or >= nEtaBins() if Eta is outside the abs eta bin edges.
0240     std::size_t binIndexNoCheck(double eta) const;
0241 
0242     /// Get the cuts for a given eta
0243     /// @param eta Eta value
0244     /// @return Cuts for the given eta
0245     const Config& getCuts(double eta) const;
0246   };
0247 
0248   /// Constructor from a single cut config object
0249   /// @param config is the configuration object
0250   explicit TrackSelector(const Config& config);
0251 
0252   /// Constructor from a multi-eta
0253   /// @param config is the configuration object
0254   explicit TrackSelector(const EtaBinnedConfig& config);
0255 
0256   /// Select tracks from an input container and copy them into an output
0257   /// container
0258   /// @tparam input_tracks_t is the type of the input track container
0259   /// @tparam output_tracks_t is the type of the output track container
0260   /// @param inputTracks is the input track container
0261   /// @param outputTracks is the output track container
0262   template <typename input_tracks_t, typename output_tracks_t>
0263   void selectTracks(const input_tracks_t& inputTracks,
0264                     output_tracks_t& outputTracks) const;
0265 
0266   /// Helper function to check if a track is valid
0267   /// @tparam track_proxy_t is the type of the track proxy
0268   /// @param track is the track proxy
0269   /// @return true if the track is valid
0270   template <TrackProxyConcept track_proxy_t>
0271   bool isValidTrack(const track_proxy_t& track) const;
0272 
0273   /// Get readonly access to the config parameters
0274   /// @return the config object
0275   const EtaBinnedConfig& config() const { return m_cfg; }
0276 
0277  private:
0278   EtaBinnedConfig m_cfg;
0279   bool m_isUnbinned = false;
0280 };
0281 
0282 inline TrackSelector::Config& TrackSelector::Config::loc0(double min,
0283                                                           double max) {
0284   loc0Min = min;
0285   loc0Max = max;
0286   return *this;
0287 }
0288 
0289 inline TrackSelector::Config& TrackSelector::Config::loc1(double min,
0290                                                           double max) {
0291   loc1Min = min;
0292   loc1Max = max;
0293   return *this;
0294 }
0295 
0296 inline TrackSelector::Config& TrackSelector::Config::time(double min,
0297                                                           double max) {
0298   timeMin = min;
0299   timeMax = max;
0300   return *this;
0301 }
0302 
0303 inline TrackSelector::Config& TrackSelector::Config::phi(double min,
0304                                                          double max) {
0305   phiMin = min;
0306   phiMax = max;
0307   return *this;
0308 }
0309 
0310 inline TrackSelector::Config& TrackSelector::Config::eta(double min,
0311                                                          double max) {
0312   if (absEtaMin != 0.0 || absEtaMax != inf) {
0313     throw std::invalid_argument(
0314         "Cannot set both eta and absEta cuts in the same cut set");
0315   }
0316   etaMin = min;
0317   etaMax = max;
0318   return *this;
0319 }
0320 
0321 inline TrackSelector::Config& TrackSelector::Config::absEta(double min,
0322                                                             double max) {
0323   if (etaMin != -inf || etaMax != inf) {
0324     throw std::invalid_argument(
0325         "Cannot set both eta and absEta cuts in the same cut set");
0326   }
0327   absEtaMin = min;
0328   absEtaMax = max;
0329   return *this;
0330 }
0331 
0332 inline TrackSelector::Config& TrackSelector::Config::pt(double min,
0333                                                         double max) {
0334   ptMin = min;
0335   ptMax = max;
0336   return *this;
0337 }
0338 
0339 inline std::ostream& operator<<(std::ostream& os,
0340                                 const TrackSelector::Config& cuts) {
0341   // for printing cuts set up with `within`
0342   auto printMinMax = [&](const char* name, const auto& min, const auto& max) {
0343     os << " - " << min << " <= " << name << " < " << max << "\n";
0344   };
0345   // for printing cuts set up with `checkMin`
0346   auto printMin = [&](const char* name, const auto& min) {
0347     os << " - " << min << " <= " << name << "\n";
0348   };
0349   // for printing cuts set up with `checkMax`
0350   auto printMax = [&](const char* name, const auto& max) {
0351     os << " - " << name << " <= " << max << "\n";
0352   };
0353 
0354   printMinMax("loc0", cuts.loc0Min, cuts.loc0Max);
0355   printMinMax("loc1", cuts.loc1Min, cuts.loc1Max);
0356   printMinMax("time", cuts.timeMin, cuts.timeMax);
0357   printMinMax("phi", cuts.phiMin, cuts.phiMax);
0358   printMinMax("eta", cuts.etaMin, cuts.etaMax);
0359   printMinMax("absEta", cuts.absEtaMin, cuts.absEtaMax);
0360   printMinMax("pt", cuts.ptMin, cuts.ptMax);
0361   printMax("nHoles", cuts.maxHoles);
0362   printMax("nOutliers", cuts.maxOutliers);
0363   printMax("nHoles + nOutliers", cuts.maxHolesAndOutliers);
0364   printMax("nSharedHits", cuts.maxSharedHits);
0365   printMax("chi2", cuts.maxChi2);
0366   printMin("nMeasurements", cuts.minMeasurements);
0367   return os;
0368 }
0369 
0370 inline TrackSelector::EtaBinnedConfig& TrackSelector::EtaBinnedConfig::addCuts(
0371     double etaMax, const std::function<void(Config&)>& callback) {
0372   if (etaMax <= absEtaEdges.back()) {
0373     throw std::invalid_argument{
0374         "Abs Eta bin edges must be in increasing order"};
0375   }
0376 
0377   if (etaMax < 0.0) {
0378     throw std::invalid_argument{"Abs Eta bin edges must be positive"};
0379   }
0380 
0381   absEtaEdges.push_back(etaMax);
0382   cutSets.emplace_back();
0383   if (callback) {
0384     callback(cutSets.back());
0385   }
0386   return *this;
0387 }
0388 
0389 inline TrackSelector::EtaBinnedConfig& TrackSelector::EtaBinnedConfig::addCuts(
0390     const std::function<void(Config&)>& callback) {
0391   return addCuts(inf, callback);
0392 }
0393 
0394 inline bool TrackSelector::EtaBinnedConfig::hasCuts(double eta) const {
0395   return std::abs(eta) < absEtaEdges.back();
0396 }
0397 
0398 inline std::size_t TrackSelector::EtaBinnedConfig::binIndex(double eta) const {
0399   std::size_t index = binIndexNoCheck(eta);
0400   if (!(index < nEtaBins())) {
0401     throw std::invalid_argument{"Eta is outside the abs eta bin edges"};
0402   }
0403   return index;
0404 }
0405 
0406 inline std::size_t TrackSelector::EtaBinnedConfig::binIndexNoCheck(
0407     double eta) const {
0408   auto binIt = std::ranges::upper_bound(absEtaEdges, std::abs(eta));
0409   std::size_t index = std::ranges::distance(absEtaEdges.begin(), binIt);
0410   if (index == 0) {
0411     index = absEtaEdges.size() + 1;  // positive value to check for underflow
0412   }
0413   return index - 1;
0414 }
0415 
0416 inline const TrackSelector::Config& TrackSelector::EtaBinnedConfig::getCuts(
0417     double eta) const {
0418   return nEtaBins() == 1 ? cutSets[0] : cutSets[binIndex(eta)];
0419 }
0420 
0421 inline std::ostream& operator<<(std::ostream& os,
0422                                 const TrackSelector::EtaBinnedConfig& cfg) {
0423   os << "TrackSelector::EtaBinnedConfig:\n";
0424 
0425   for (std::size_t i = 1; i < cfg.absEtaEdges.size(); i++) {
0426     os << cfg.absEtaEdges[i - 1] << " <= eta < " << cfg.absEtaEdges[i] << "\n";
0427     os << cfg.cutSets[i - 1];
0428   }
0429 
0430   return os;
0431 }
0432 
0433 template <typename input_tracks_t, typename output_tracks_t>
0434 void TrackSelector::selectTracks(const input_tracks_t& inputTracks,
0435                                  output_tracks_t& outputTracks) const {
0436   for (auto track : inputTracks) {
0437     if (!isValidTrack(track)) {
0438       continue;
0439     }
0440     auto destProxy = outputTracks.makeTrack();
0441     destProxy.copyFromShallow(track);
0442   }
0443 }
0444 
0445 template <TrackProxyConcept track_proxy_t>
0446 bool TrackSelector::isValidTrack(const track_proxy_t& track) const {
0447   auto checkMin = [](auto x, auto min) { return min <= x; };
0448   auto checkMax = [](auto x, auto max) { return x <= max; };
0449   auto within = [](double x, double min, double max) {
0450     return (min <= x) && (x < max);
0451   };
0452 
0453   const auto theta = track.theta();
0454 
0455   constexpr double kUnset = -std::numeric_limits<double>::infinity();
0456 
0457   double _eta = kUnset;
0458   double _absEta = kUnset;
0459 
0460   auto absEta = [&]() {
0461     if (_absEta == kUnset) {
0462       _eta = AngleHelpers::etaFromTheta(theta);
0463       _absEta = std::abs(_eta);
0464     }
0465     return _absEta;
0466   };
0467 
0468   const Config* cutsPtr{nullptr};
0469   if (!m_isUnbinned) {
0470     // return false if |eta| is outside its range, or nan.
0471     if (!(absEta() >= m_cfg.absEtaEdges.front() &&
0472           _absEta < m_cfg.absEtaEdges.back())) {
0473       return false;
0474     }
0475     cutsPtr = &m_cfg.getCuts(_eta);
0476   } else {
0477     cutsPtr = &m_cfg.cutSets.front();
0478   }
0479 
0480   const Config& cuts = *cutsPtr;
0481 
0482   auto parameterCuts = [&]() {
0483     return within(track.transverseMomentum(), cuts.ptMin, cuts.ptMax) &&
0484            (!m_isUnbinned ||
0485             (within(absEta(), cuts.absEtaMin, cuts.absEtaMax) &&
0486              within(_eta, cuts.etaMin, cuts.etaMax))) &&
0487            within(track.phi(), cuts.phiMin, cuts.phiMax) &&
0488            within(track.loc0(), cuts.loc0Min, cuts.loc0Max) &&
0489            within(track.loc1(), cuts.loc1Min, cuts.loc1Max) &&
0490            within(track.time(), cuts.timeMin, cuts.timeMax);
0491   };
0492 
0493   auto trackCuts = [&]() {
0494     return checkMin(track.nMeasurements(), cuts.minMeasurements) &&
0495            checkMax(track.nHoles(), cuts.maxHoles) &&
0496            checkMax(track.nOutliers(), cuts.maxOutliers) &&
0497            checkMax(track.nHoles() + track.nOutliers(),
0498                     cuts.maxHolesAndOutliers) &&
0499            checkMax(track.nSharedHits(), cuts.maxSharedHits) &&
0500            checkMax(track.chi2(), cuts.maxChi2) &&
0501            cuts.measurementCounter.isValidTrack(track);
0502   };
0503 
0504   if (cuts.requireReferenceSurface) {
0505     return track.hasReferenceSurface() && parameterCuts() && trackCuts();
0506   } else {
0507     return trackCuts();
0508   }
0509 }
0510 
0511 inline TrackSelector::TrackSelector(
0512     const TrackSelector::EtaBinnedConfig& config)
0513     : m_cfg(config) {
0514   if (m_cfg.cutSets.size() != m_cfg.nEtaBins()) {
0515     throw std::invalid_argument{
0516         "TrackSelector cut / eta bin configuration is inconsistent"};
0517   }
0518 
0519   if (m_cfg.nEtaBins() == 1) {
0520     static const std::vector<double> infVec = {0, inf};
0521     m_isUnbinned =
0522         m_cfg.absEtaEdges == infVec;  // single bin, no eta edges given
0523   }
0524 
0525   if (!m_isUnbinned) {
0526     for (const auto& cuts : m_cfg.cutSets) {
0527       if (cuts.etaMin != -inf || cuts.etaMax != inf || cuts.absEtaMin != 0.0 ||
0528           cuts.absEtaMax != inf) {
0529         throw std::invalid_argument{
0530             "Explicit eta cuts are only valid for single eta bin"};
0531       }
0532     }
0533   }
0534 }
0535 
0536 inline TrackSelector::TrackSelector(const Config& config)
0537     : TrackSelector{EtaBinnedConfig{config}} {}
0538 
0539 template <TrackProxyConcept track_proxy_t>
0540 bool TrackSelector::MeasurementCounter::isValidTrack(
0541     const track_proxy_t& track) const {
0542   // No hit cuts, accept everything
0543   if (counters.empty()) {
0544     return true;
0545   }
0546 
0547   boost::container::small_vector<unsigned int, 4> counterValues;
0548   counterValues.resize(counters.size(), 0);
0549 
0550   for (const auto& ts : track.trackStatesReversed()) {
0551     if (!ts.typeFlags().isMeasurement()) {
0552       continue;
0553     }
0554 
0555     const auto geoId = ts.referenceSurface().geometryId();
0556 
0557     for (std::size_t i = 0; i < counters.size(); i++) {
0558       const auto& [counterMap, threshold] = counters[i];
0559       if (const auto it = counterMap.find(geoId); it != counterMap.end()) {
0560         counterValues[i]++;
0561       }
0562     }
0563   }
0564 
0565   for (std::size_t i = 0; i < counters.size(); i++) {
0566     const auto& [counterMap, threshold] = counters[i];
0567     const unsigned int value = counterValues[i];
0568     if (value < threshold) {
0569       return false;
0570     }
0571   }
0572 
0573   return true;
0574 }
0575 }  // namespace Acts