File indexing completed on 2026-07-06 08:03:15
0001
0002
0003
0004
0005
0006
0007
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
0027
0028
0029 class TrackSelector {
0030 static constexpr double inf = std::numeric_limits<double>::infinity();
0031
0032 public:
0033
0034 struct MeasurementCounter {
0035
0036 using CounterElement =
0037 std::pair<GeometryHierarchyMap<unsigned int>, unsigned int>;
0038
0039
0040 boost::container::small_vector<CounterElement, 4> counters;
0041
0042
0043
0044
0045 template <TrackProxyConcept track_proxy_t>
0046 bool isValidTrack(const track_proxy_t& track) const;
0047
0048
0049
0050
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
0062
0063 struct Config {
0064
0065
0066 double loc0Min = -inf;
0067
0068 double loc0Max = inf;
0069
0070 double loc1Min = -inf;
0071
0072 double loc1Max = inf;
0073
0074
0075 double timeMin = -inf;
0076
0077 double timeMax = inf;
0078
0079
0080 double phiMin = -inf;
0081
0082 double phiMax = inf;
0083
0084 double etaMin = -inf;
0085
0086 double etaMax = inf;
0087
0088 double absEtaMin = 0.0;
0089
0090 double absEtaMax = inf;
0091
0092
0093 double ptMin = 0.0;
0094
0095 double ptMax = inf;
0096
0097
0098 std::size_t minMeasurements = 0;
0099
0100 std::size_t maxHoles = std::numeric_limits<std::size_t>::max();
0101
0102 std::size_t maxOutliers = std::numeric_limits<std::size_t>::max();
0103
0104 std::size_t maxHolesAndOutliers = std::numeric_limits<std::size_t>::max();
0105
0106 std::size_t maxSharedHits = std::numeric_limits<std::size_t>::max();
0107
0108 double maxChi2 = inf;
0109
0110
0111
0112 bool requireReferenceSurface = true;
0113
0114
0115
0116 MeasurementCounter measurementCounter;
0117
0118
0119
0120
0121
0122
0123
0124
0125 Config& loc0(double min, double max);
0126
0127
0128
0129
0130
0131 Config& loc1(double min, double max);
0132
0133
0134
0135
0136
0137 Config& time(double min, double max);
0138
0139
0140
0141
0142
0143 Config& phi(double min, double max);
0144
0145
0146
0147
0148
0149 Config& eta(double min, double max);
0150
0151
0152
0153
0154
0155 Config& absEta(double min, double max);
0156
0157
0158
0159
0160
0161 Config& pt(double min, double max);
0162
0163
0164
0165
0166
0167 friend std::ostream& operator<<(std::ostream& os, const Config& cuts);
0168 };
0169
0170
0171
0172 struct EtaBinnedConfig {
0173
0174 std::vector<Config> cutSets = {};
0175
0176
0177 std::vector<double> absEtaEdges = {0, inf};
0178
0179
0180
0181 std::size_t nEtaBins() const { return absEtaEdges.size() - 1; }
0182
0183
0184
0185 EtaBinnedConfig() : cutSets{{}} {};
0186
0187
0188
0189
0190 explicit EtaBinnedConfig(double etaMin) : cutSets{}, absEtaEdges{etaMin} {}
0191
0192
0193
0194
0195
0196 explicit EtaBinnedConfig(std::vector<double> absEtaEdgesIn)
0197 : absEtaEdges{std::move(absEtaEdgesIn)} {
0198 cutSets.resize(nEtaBins());
0199 }
0200
0201
0202
0203
0204 explicit EtaBinnedConfig(Config cutSet) : cutSets{std::move(cutSet)} {}
0205
0206
0207
0208
0209
0210 EtaBinnedConfig& addCuts(double etaMax,
0211 const std::function<void(Config&)>& callback = {});
0212
0213
0214
0215
0216 EtaBinnedConfig& addCuts(const std::function<void(Config&)>& callback = {});
0217
0218
0219
0220
0221
0222 friend std::ostream& operator<<(std::ostream& os,
0223 const EtaBinnedConfig& cfg);
0224
0225
0226
0227
0228 bool hasCuts(double eta) const;
0229
0230
0231
0232
0233
0234 std::size_t binIndex(double eta) const;
0235
0236
0237
0238
0239 std::size_t binIndexNoCheck(double eta) const;
0240
0241
0242
0243
0244 const Config& getCuts(double eta) const;
0245 };
0246
0247
0248
0249 explicit TrackSelector(const Config& config);
0250
0251
0252
0253 explicit TrackSelector(const EtaBinnedConfig& config);
0254
0255
0256
0257
0258
0259
0260
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
0266
0267
0268
0269 template <TrackProxyConcept track_proxy_t>
0270 bool isValidTrack(const track_proxy_t& track) const;
0271
0272
0273
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
0341 auto printMinMax = [&](const char* name, const auto& min, const auto& max) {
0342 os << " - " << min << " <= " << name << " < " << max << "\n";
0343 };
0344
0345 auto printMin = [&](const char* name, const auto& min) {
0346 os << " - " << min << " <= " << name << "\n";
0347 };
0348
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;
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
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;
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
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 }