Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-01 08:08:02

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