Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-07-06 08:03:15

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