Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /include/Acts/TrackFinding/TrackSelector.hpp was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

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