File indexing completed on 2026-04-11 07:46:41
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 <algorithm>
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
0028
0029
0030 class TrackSelector {
0031 static constexpr double inf = std::numeric_limits<double>::infinity();
0032
0033 public:
0034
0035 struct MeasurementCounter {
0036
0037 using CounterElement =
0038 std::pair<GeometryHierarchyMap<unsigned int>, unsigned int>;
0039
0040
0041 boost::container::small_vector<CounterElement, 4> counters;
0042
0043
0044
0045
0046 template <TrackProxyConcept track_proxy_t>
0047 bool isValidTrack(const track_proxy_t& track) const;
0048
0049
0050
0051
0052 void addCounter(const std::vector<GeometryIdentifier>& identifiers,
0053 unsigned int threshold) {
0054 std::vector<GeometryHierarchyMap<unsigned int>::InputElement> elements;
0055 for (const auto& id : identifiers) {
0056 elements.emplace_back(id, 0);
0057 }
0058 counters.emplace_back(std::move(elements), threshold);
0059 }
0060 };
0061
0062
0063
0064 struct Config {
0065
0066
0067 double loc0Min = -inf;
0068
0069 double loc0Max = inf;
0070
0071 double loc1Min = -inf;
0072
0073 double loc1Max = inf;
0074
0075
0076 double timeMin = -inf;
0077
0078 double timeMax = inf;
0079
0080
0081 double phiMin = -inf;
0082
0083 double phiMax = inf;
0084
0085 double etaMin = -inf;
0086
0087 double etaMax = inf;
0088
0089 double absEtaMin = 0.0;
0090
0091 double absEtaMax = inf;
0092
0093
0094 double ptMin = 0.0;
0095
0096 double ptMax = inf;
0097
0098
0099 std::size_t minMeasurements = 0;
0100
0101 std::size_t maxHoles = std::numeric_limits<std::size_t>::max();
0102
0103 std::size_t maxOutliers = std::numeric_limits<std::size_t>::max();
0104
0105 std::size_t maxHolesAndOutliers = std::numeric_limits<std::size_t>::max();
0106
0107 std::size_t maxSharedHits = std::numeric_limits<std::size_t>::max();
0108
0109 double maxChi2 = inf;
0110
0111
0112
0113 bool requireReferenceSurface = true;
0114
0115
0116
0117 MeasurementCounter measurementCounter;
0118
0119
0120
0121
0122
0123
0124
0125
0126 Config& loc0(double min, double max);
0127
0128
0129
0130
0131
0132 Config& loc1(double min, double max);
0133
0134
0135
0136
0137
0138 Config& time(double min, double max);
0139
0140
0141
0142
0143
0144 Config& phi(double min, double max);
0145
0146
0147
0148
0149
0150 Config& eta(double min, double max);
0151
0152
0153
0154
0155
0156 Config& absEta(double min, double max);
0157
0158
0159
0160
0161
0162 Config& pt(double min, double max);
0163
0164
0165
0166
0167
0168 friend std::ostream& operator<<(std::ostream& os, const Config& cuts);
0169 };
0170
0171
0172
0173 struct EtaBinnedConfig {
0174
0175 std::vector<Config> cutSets = {};
0176
0177
0178 std::vector<double> absEtaEdges = {0, inf};
0179
0180
0181
0182 std::size_t nEtaBins() const { return absEtaEdges.size() - 1; }
0183
0184
0185
0186 EtaBinnedConfig() : cutSets{{}} {};
0187
0188
0189
0190
0191 explicit EtaBinnedConfig(double etaMin) : cutSets{}, absEtaEdges{etaMin} {}
0192
0193
0194
0195
0196
0197 explicit EtaBinnedConfig(std::vector<double> absEtaEdgesIn)
0198 : absEtaEdges{std::move(absEtaEdgesIn)} {
0199 cutSets.resize(nEtaBins());
0200 }
0201
0202
0203
0204
0205 explicit EtaBinnedConfig(Config cutSet) : cutSets{std::move(cutSet)} {}
0206
0207
0208
0209
0210
0211 EtaBinnedConfig& addCuts(double etaMax,
0212 const std::function<void(Config&)>& callback = {});
0213
0214
0215
0216
0217 EtaBinnedConfig& addCuts(const std::function<void(Config&)>& callback = {});
0218
0219
0220
0221
0222
0223 friend std::ostream& operator<<(std::ostream& os,
0224 const EtaBinnedConfig& cfg);
0225
0226
0227
0228
0229 bool hasCuts(double eta) const;
0230
0231
0232
0233
0234
0235 std::size_t binIndex(double eta) const;
0236
0237
0238
0239
0240 std::size_t binIndexNoCheck(double eta) const;
0241
0242
0243
0244
0245 const Config& getCuts(double eta) const;
0246 };
0247
0248
0249
0250 explicit TrackSelector(const Config& config);
0251
0252
0253
0254 explicit TrackSelector(const EtaBinnedConfig& config);
0255
0256
0257
0258
0259
0260
0261
0262 template <typename input_tracks_t, typename output_tracks_t>
0263 void selectTracks(const input_tracks_t& inputTracks,
0264 output_tracks_t& outputTracks) const;
0265
0266
0267
0268
0269
0270 template <TrackProxyConcept track_proxy_t>
0271 bool isValidTrack(const track_proxy_t& track) const;
0272
0273
0274
0275 const EtaBinnedConfig& config() const { return m_cfg; }
0276
0277 private:
0278 EtaBinnedConfig m_cfg;
0279 bool m_isUnbinned = false;
0280 };
0281
0282 inline TrackSelector::Config& TrackSelector::Config::loc0(double min,
0283 double max) {
0284 loc0Min = min;
0285 loc0Max = max;
0286 return *this;
0287 }
0288
0289 inline TrackSelector::Config& TrackSelector::Config::loc1(double min,
0290 double max) {
0291 loc1Min = min;
0292 loc1Max = max;
0293 return *this;
0294 }
0295
0296 inline TrackSelector::Config& TrackSelector::Config::time(double min,
0297 double max) {
0298 timeMin = min;
0299 timeMax = max;
0300 return *this;
0301 }
0302
0303 inline TrackSelector::Config& TrackSelector::Config::phi(double min,
0304 double max) {
0305 phiMin = min;
0306 phiMax = max;
0307 return *this;
0308 }
0309
0310 inline TrackSelector::Config& TrackSelector::Config::eta(double min,
0311 double max) {
0312 if (absEtaMin != 0.0 || absEtaMax != inf) {
0313 throw std::invalid_argument(
0314 "Cannot set both eta and absEta cuts in the same cut set");
0315 }
0316 etaMin = min;
0317 etaMax = max;
0318 return *this;
0319 }
0320
0321 inline TrackSelector::Config& TrackSelector::Config::absEta(double min,
0322 double max) {
0323 if (etaMin != -inf || etaMax != inf) {
0324 throw std::invalid_argument(
0325 "Cannot set both eta and absEta cuts in the same cut set");
0326 }
0327 absEtaMin = min;
0328 absEtaMax = max;
0329 return *this;
0330 }
0331
0332 inline TrackSelector::Config& TrackSelector::Config::pt(double min,
0333 double max) {
0334 ptMin = min;
0335 ptMax = max;
0336 return *this;
0337 }
0338
0339 inline std::ostream& operator<<(std::ostream& os,
0340 const TrackSelector::Config& cuts) {
0341
0342 auto printMinMax = [&](const char* name, const auto& min, const auto& max) {
0343 os << " - " << min << " <= " << name << " < " << max << "\n";
0344 };
0345
0346 auto printMin = [&](const char* name, const auto& min) {
0347 os << " - " << min << " <= " << name << "\n";
0348 };
0349
0350 auto printMax = [&](const char* name, const auto& max) {
0351 os << " - " << name << " <= " << max << "\n";
0352 };
0353
0354 printMinMax("loc0", cuts.loc0Min, cuts.loc0Max);
0355 printMinMax("loc1", cuts.loc1Min, cuts.loc1Max);
0356 printMinMax("time", cuts.timeMin, cuts.timeMax);
0357 printMinMax("phi", cuts.phiMin, cuts.phiMax);
0358 printMinMax("eta", cuts.etaMin, cuts.etaMax);
0359 printMinMax("absEta", cuts.absEtaMin, cuts.absEtaMax);
0360 printMinMax("pt", cuts.ptMin, cuts.ptMax);
0361 printMax("nHoles", cuts.maxHoles);
0362 printMax("nOutliers", cuts.maxOutliers);
0363 printMax("nHoles + nOutliers", cuts.maxHolesAndOutliers);
0364 printMax("nSharedHits", cuts.maxSharedHits);
0365 printMax("chi2", cuts.maxChi2);
0366 printMin("nMeasurements", cuts.minMeasurements);
0367 return os;
0368 }
0369
0370 inline TrackSelector::EtaBinnedConfig& TrackSelector::EtaBinnedConfig::addCuts(
0371 double etaMax, const std::function<void(Config&)>& callback) {
0372 if (etaMax <= absEtaEdges.back()) {
0373 throw std::invalid_argument{
0374 "Abs Eta bin edges must be in increasing order"};
0375 }
0376
0377 if (etaMax < 0.0) {
0378 throw std::invalid_argument{"Abs Eta bin edges must be positive"};
0379 }
0380
0381 absEtaEdges.push_back(etaMax);
0382 cutSets.emplace_back();
0383 if (callback) {
0384 callback(cutSets.back());
0385 }
0386 return *this;
0387 }
0388
0389 inline TrackSelector::EtaBinnedConfig& TrackSelector::EtaBinnedConfig::addCuts(
0390 const std::function<void(Config&)>& callback) {
0391 return addCuts(inf, callback);
0392 }
0393
0394 inline bool TrackSelector::EtaBinnedConfig::hasCuts(double eta) const {
0395 return std::abs(eta) < absEtaEdges.back();
0396 }
0397
0398 inline std::size_t TrackSelector::EtaBinnedConfig::binIndex(double eta) const {
0399 std::size_t index = binIndexNoCheck(eta);
0400 if (!(index < nEtaBins())) {
0401 throw std::invalid_argument{"Eta is outside the abs eta bin edges"};
0402 }
0403 return index;
0404 }
0405
0406 inline std::size_t TrackSelector::EtaBinnedConfig::binIndexNoCheck(
0407 double eta) const {
0408 auto binIt = std::ranges::upper_bound(absEtaEdges, std::abs(eta));
0409 std::size_t index = std::ranges::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.copyFromShallow(track);
0442 }
0443 }
0444
0445 template <TrackProxyConcept track_proxy_t>
0446 bool TrackSelector::isValidTrack(const track_proxy_t& track) const {
0447 auto checkMin = [](auto x, auto min) { return min <= x; };
0448 auto checkMax = [](auto x, auto max) { return x <= max; };
0449 auto within = [](double x, double min, double max) {
0450 return (min <= x) && (x < max);
0451 };
0452
0453 const auto theta = track.theta();
0454
0455 constexpr double kUnset = -std::numeric_limits<double>::infinity();
0456
0457 double _eta = kUnset;
0458 double _absEta = kUnset;
0459
0460 auto absEta = [&]() {
0461 if (_absEta == kUnset) {
0462 _eta = AngleHelpers::etaFromTheta(theta);
0463 _absEta = std::abs(_eta);
0464 }
0465 return _absEta;
0466 };
0467
0468 const Config* cutsPtr{nullptr};
0469 if (!m_isUnbinned) {
0470
0471 if (!(absEta() >= m_cfg.absEtaEdges.front() &&
0472 _absEta < m_cfg.absEtaEdges.back())) {
0473 return false;
0474 }
0475 cutsPtr = &m_cfg.getCuts(_eta);
0476 } else {
0477 cutsPtr = &m_cfg.cutSets.front();
0478 }
0479
0480 const Config& cuts = *cutsPtr;
0481
0482 auto parameterCuts = [&]() {
0483 return within(track.transverseMomentum(), cuts.ptMin, cuts.ptMax) &&
0484 (!m_isUnbinned ||
0485 (within(absEta(), cuts.absEtaMin, cuts.absEtaMax) &&
0486 within(_eta, cuts.etaMin, cuts.etaMax))) &&
0487 within(track.phi(), cuts.phiMin, cuts.phiMax) &&
0488 within(track.loc0(), cuts.loc0Min, cuts.loc0Max) &&
0489 within(track.loc1(), cuts.loc1Min, cuts.loc1Max) &&
0490 within(track.time(), cuts.timeMin, cuts.timeMax);
0491 };
0492
0493 auto trackCuts = [&]() {
0494 return checkMin(track.nMeasurements(), cuts.minMeasurements) &&
0495 checkMax(track.nHoles(), cuts.maxHoles) &&
0496 checkMax(track.nOutliers(), cuts.maxOutliers) &&
0497 checkMax(track.nHoles() + track.nOutliers(),
0498 cuts.maxHolesAndOutliers) &&
0499 checkMax(track.nSharedHits(), cuts.maxSharedHits) &&
0500 checkMax(track.chi2(), cuts.maxChi2) &&
0501 cuts.measurementCounter.isValidTrack(track);
0502 };
0503
0504 if (cuts.requireReferenceSurface) {
0505 return track.hasReferenceSurface() && parameterCuts() && trackCuts();
0506 } else {
0507 return trackCuts();
0508 }
0509 }
0510
0511 inline TrackSelector::TrackSelector(
0512 const TrackSelector::EtaBinnedConfig& config)
0513 : m_cfg(config) {
0514 if (m_cfg.cutSets.size() != m_cfg.nEtaBins()) {
0515 throw std::invalid_argument{
0516 "TrackSelector cut / eta bin configuration is inconsistent"};
0517 }
0518
0519 if (m_cfg.nEtaBins() == 1) {
0520 static const std::vector<double> infVec = {0, inf};
0521 m_isUnbinned =
0522 m_cfg.absEtaEdges == infVec;
0523 }
0524
0525 if (!m_isUnbinned) {
0526 for (const auto& cuts : m_cfg.cutSets) {
0527 if (cuts.etaMin != -inf || cuts.etaMax != inf || cuts.absEtaMin != 0.0 ||
0528 cuts.absEtaMax != inf) {
0529 throw std::invalid_argument{
0530 "Explicit eta cuts are only valid for single eta bin"};
0531 }
0532 }
0533 }
0534 }
0535
0536 inline TrackSelector::TrackSelector(const Config& config)
0537 : TrackSelector{EtaBinnedConfig{config}} {}
0538
0539 template <TrackProxyConcept track_proxy_t>
0540 bool TrackSelector::MeasurementCounter::isValidTrack(
0541 const track_proxy_t& track) const {
0542
0543 if (counters.empty()) {
0544 return true;
0545 }
0546
0547 boost::container::small_vector<unsigned int, 4> counterValues;
0548 counterValues.resize(counters.size(), 0);
0549
0550 for (const auto& ts : track.trackStatesReversed()) {
0551 if (!ts.typeFlags().isMeasurement()) {
0552 continue;
0553 }
0554
0555 const auto geoId = ts.referenceSurface().geometryId();
0556
0557 for (std::size_t i = 0; i < counters.size(); i++) {
0558 const auto& [counterMap, threshold] = counters[i];
0559 if (const auto it = counterMap.find(geoId); it != counterMap.end()) {
0560 counterValues[i]++;
0561 }
0562 }
0563 }
0564
0565 for (std::size_t i = 0; i < counters.size(); i++) {
0566 const auto& [counterMap, threshold] = counters[i];
0567 const unsigned int value = counterValues[i];
0568 if (value < threshold) {
0569 return false;
0570 }
0571 }
0572
0573 return true;
0574 }
0575 }