File indexing completed on 2025-07-01 08:08:02
0001
0002
0003
0004
0005
0006
0007
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
0027
0028
0029 class TrackSelector {
0030 static constexpr double inf = std::numeric_limits<double>::infinity();
0031
0032 public:
0033 struct MeasurementCounter {
0034
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
0054
0055 struct Config {
0056
0057 double loc0Min = -inf;
0058 double loc0Max = inf;
0059 double loc1Min = -inf;
0060 double loc1Max = inf;
0061
0062 double timeMin = -inf;
0063 double timeMax = inf;
0064
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
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
0083 MeasurementCounter measurementCounter;
0084
0085
0086
0087
0088
0089
0090
0091
0092 Config& loc0(double min, double max);
0093
0094
0095
0096
0097
0098 Config& loc1(double min, double max);
0099
0100
0101
0102
0103
0104 Config& time(double min, double max);
0105
0106
0107
0108
0109
0110 Config& phi(double min, double max);
0111
0112
0113
0114
0115
0116 Config& eta(double min, double max);
0117
0118
0119
0120
0121
0122 Config& absEta(double min, double max);
0123
0124
0125
0126
0127
0128 Config& pt(double min, double max);
0129
0130
0131
0132
0133
0134 friend std::ostream& operator<<(std::ostream& os, const Config& cuts);
0135 };
0136
0137
0138
0139 struct EtaBinnedConfig {
0140
0141 std::vector<Config> cutSets = {};
0142
0143
0144 std::vector<double> absEtaEdges = {};
0145
0146
0147
0148 std::size_t nEtaBins() const { return absEtaEdges.size() - 1; }
0149
0150
0151
0152 EtaBinnedConfig() : cutSets{{}}, absEtaEdges{{0, inf}} {};
0153
0154
0155
0156
0157 EtaBinnedConfig(double etaMin) : cutSets{}, absEtaEdges{etaMin} {}
0158
0159
0160
0161
0162
0163 EtaBinnedConfig(std::vector<double> absEtaEdgesIn)
0164 : absEtaEdges{std::move(absEtaEdgesIn)} {
0165 cutSets.resize(absEtaEdges.size() - 1);
0166 }
0167
0168
0169
0170 EtaBinnedConfig(Config cutSet)
0171 : cutSets{std::move(cutSet)}, absEtaEdges{{0, inf}} {}
0172
0173
0174
0175
0176
0177 EtaBinnedConfig& addCuts(double etaMax,
0178 const std::function<void(Config&)>& callback = {});
0179
0180
0181
0182
0183 EtaBinnedConfig& addCuts(const std::function<void(Config&)>& callback = {});
0184
0185
0186
0187
0188
0189 friend std::ostream& operator<<(std::ostream& os,
0190 const EtaBinnedConfig& cfg);
0191
0192
0193
0194
0195 bool hasCuts(double eta) const;
0196
0197
0198
0199
0200 std::size_t binIndex(double eta) const;
0201
0202
0203
0204
0205 const Config& getCuts(double eta) const;
0206 };
0207
0208
0209
0210 TrackSelector(const Config& config);
0211
0212
0213
0214 TrackSelector(const EtaBinnedConfig& config);
0215
0216
0217
0218
0219
0220
0221
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
0227
0228
0229
0230 template <typename track_proxy_t>
0231 bool isValidTrack(const track_proxy_t& track) const;
0232
0233
0234
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;
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
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 }