Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-12-08 09:44:50

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