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