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