Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.ipp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2024 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include "Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Utilities/VectorHelpers.hpp"
0014
0015 #include <unordered_map>
0016
0017 namespace Acts {
0018
0019 inline const Logger& ScoreBasedAmbiguityResolution::logger() const {
0020   return *m_logger;
0021 }
0022
0023 template <typename track_container_t, typename traj_t,
0024           template <typename> class holder_t, typename source_link_hash_t,
0025           typename source_link_equality_t>
0026 std::vector<std::vector<ScoreBasedAmbiguityResolution::MeasurementInfo>>
0027 ScoreBasedAmbiguityResolution::computeInitialState(
0028     const TrackContainer<track_container_t, traj_t, holder_t>& tracks,
0029     source_link_hash_t sourceLinkHash,
0030     source_link_equality_t sourceLinkEquality,
0031     std::vector<std::vector<TrackFeatures>>& trackFeaturesVectors) const {
0032   auto MeasurementIndexMap =
0033       std::unordered_map<SourceLink, std::size_t, source_link_hash_t,
0034                          source_link_equality_t>(0, sourceLinkHash,
0035                                                  sourceLinkEquality);
0036
0037   std::vector<std::vector<MeasurementInfo>> measurementsPerTrack;
0038   measurementsPerTrack.reserve(tracks.size());
0039   ACTS_VERBOSE("Starting to compute initial state");
0040
0041   for (const auto& track : tracks) {
0042     int numberOfDetectors = m_cfg.detectorConfigs.size();
0043     int numberOfTrackStates = track.nTrackStates();
0044     std::vector<MeasurementInfo> measurements;
0045     measurements.reserve(numberOfTrackStates);
0046     std::vector<TrackFeatures> trackFeaturesVector(numberOfDetectors);
0047
0048     for (const auto& ts : track.trackStatesReversed()) {
0049       if (!ts.hasReferenceSurface()) {
0050         ACTS_ERROR("Track state has no reference surface");
0051         continue;
0052       }
0053       auto iVolume = ts.referenceSurface().geometryId().volume();
0054       auto volume_it = m_cfg.volumeMap.find(iVolume);
0055       if (volume_it == m_cfg.volumeMap.end()) {
0056         ACTS_ERROR("Volume " << iVolume << "not found in the volume map");
0057         continue;
0058       }
0059       auto detectorId = volume_it->second;
0060
0061       if (ts.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) {
0062         Acts::SourceLink sourceLink = ts.getUncalibratedSourceLink();
0063         ACTS_DEBUG("Track state type is MeasurementFlag");
0064
0065         if (ts.typeFlags().test(Acts::TrackStateFlag::SharedHitFlag)) {
0066           trackFeaturesVector[detectorId].nSharedHits++;
0067         }
0068         trackFeaturesVector[detectorId].nHits++;
0069
0070         // assign a new measurement index if the source link was not seen yet
0071         auto emplace = MeasurementIndexMap.try_emplace(
0072             sourceLink, MeasurementIndexMap.size());
0073
0074         bool isoutliner = false;
0075
0076         measurements.push_back({emplace.first->second, detectorId, isoutliner});
0077       } else if (ts.typeFlags().test(Acts::TrackStateFlag::OutlierFlag)) {
0078         Acts::SourceLink sourceLink = ts.getUncalibratedSourceLink();
0079         ACTS_DEBUG("Track state type is OutlierFlag");
0080         trackFeaturesVector[detectorId].nOutliers++;
0081
0082         // assign a new measurement index if the source link was not seen yet
0083         auto emplace = MeasurementIndexMap.try_emplace(
0084             sourceLink, MeasurementIndexMap.size());
0085
0086         bool isOutliner = true;
0087
0088         measurements.push_back({emplace.first->second, detectorId, isOutliner});
0089       } else if (ts.typeFlags().test(Acts::TrackStateFlag::HoleFlag)) {
0090         ACTS_DEBUG("Track state type is HoleFlag");
0091         trackFeaturesVector[detectorId].nHoles++;
0092       }
0093     }
0094     measurementsPerTrack.push_back(std::move(measurements));
0095     trackFeaturesVectors.push_back(std::move(trackFeaturesVector));
0096   }
0097
0098   return measurementsPerTrack;
0099 }
0100
0101 template <typename track_container_t, typename traj_t,
0102           template <typename> class holder_t, bool ReadOnly>
0103 std::vector<double> Acts::ScoreBasedAmbiguityResolution::simpleScore(
0104     const TrackContainer<track_container_t, traj_t, holder_t>& tracks,
0105     const std::vector<std::vector<TrackFeatures>>& trackFeaturesVectors,
0106     const OptionalCuts<track_container_t, traj_t, holder_t, ReadOnly>&
0107         optionalCuts) const {
0108   std::vector<double> trackScore;
0109   trackScore.reserve(tracks.size());
0110
0111   int iTrack = 0;
0112
0113   ACTS_VERBOSE("Number of detectors: " << m_cfg.detectorConfigs.size());
0114
0115   ACTS_INFO("Starting to score tracks");
0116
0117   // Loop over all the tracks in the container
0118   for (const auto& track : tracks) {
0119     // get the trackFeatures map for the track
0120     const auto& trackFeaturesVector = trackFeaturesVectors[iTrack];
0121     double score = 1;
0122     auto pT = Acts::VectorHelpers::perp(track.momentum());
0123     auto eta = Acts::VectorHelpers::eta(track.momentum());
0124     auto phi = Acts::VectorHelpers::phi(track.momentum());
0125     // cuts on pT
0126     if (pT < m_cfg.pTMin || pT > m_cfg.pTMax) {
0127       score = 0;
0128       iTrack++;
0129       trackScore.push_back(score);
0130       ACTS_DEBUG("Track: " << iTrack
0131                            << " has score = 0, due to pT cuts --- pT = " << pT);
0132       continue;
0133     }
0134
0135     // cuts on phi
0136     if (phi > m_cfg.phiMax || phi < m_cfg.phiMin) {
0137       score = 0;
0138       iTrack++;
0139       trackScore.push_back(score);
0140       ACTS_DEBUG("Track: " << iTrack
0141                            << " has score = 0, due to phi cuts --- phi =  "
0142                            << phi);
0143       continue;
0144     }
0145
0146     // cuts on eta
0147     if (eta > m_cfg.etaMax || eta < m_cfg.etaMin) {
0148       score = 0;
0149       iTrack++;
0150       trackScore.push_back(score);
0151       ACTS_DEBUG("Track: " << iTrack
0152                            << " has score = 0, due to eta cuts --- eta =  "
0153                            << eta);
0154       continue;
0155     }
0156
0157     // cuts on optional cuts
0158     for (const auto& cutFunction : optionalCuts.cuts) {
0159       if (cutFunction(track)) {
0160         score = 0;
0161         ACTS_DEBUG("Track: " << iTrack
0162                              << " has score = 0, due to optional cuts.");
0163         break;
0164       }
0165     }
0166
0167     if (score == 0) {
0168       iTrack++;
0169       trackScore.push_back(score);
0170       ACTS_DEBUG("Track: " << iTrack << " score : " << score);
0171       continue;
0172     }
0173     // Reject tracks which didn't pass the detector cuts.
0174     for (std::size_t detectorId = 0; detectorId < m_cfg.detectorConfigs.size();
0175          detectorId++) {
0176       const auto& detector = m_cfg.detectorConfigs.at(detectorId);
0177
0178       const auto& trackFeatures = trackFeaturesVector[detectorId];
0179
0180       ACTS_DEBUG("---> Found summary information");
0181       ACTS_DEBUG("---> Detector ID: " << detectorId);
0182       ACTS_DEBUG("---> Number of hits: " << trackFeatures.nHits);
0183       ACTS_DEBUG("---> Number of holes: " << trackFeatures.nHoles);
0184       ACTS_DEBUG("---> Number of outliers: " << trackFeatures.nOutliers);
0185
0186       if ((trackFeatures.nHits < detector.minHits) ||
0187           (trackFeatures.nHits > detector.maxHits) ||
0188           (trackFeatures.nHoles > detector.maxHoles) ||
0189           (trackFeatures.nOutliers > detector.maxOutliers)) {
0190         score = 0;
0191         ACTS_DEBUG("Track: " << iTrack
0192                              << " has score = 0, due to detector cuts");
0193         break;
0194       }
0195     }
0196
0197     if (score == 0) {
0198       iTrack++;
0199       trackScore.push_back(score);
0200       ACTS_DEBUG("Track: " << iTrack << " score : " << score);
0201       continue;
0202     }
0203
0204     // real scoring starts here
0205     // if the ambiguity scoring function is used, the score is processed with a
0206     // different algorithm than the simple score.
0207
0208     ACTS_VERBOSE("Using Simple Scoring function");
0209
0210     score = 100;
0211     // Adding the score for each detector.
0212     // detector score is determined by the number of hits/hole/outliers *
0213     // hit/hole/outlier scoreWeights in a detector.
0214     for (std::size_t detectorId = 0; detectorId < m_cfg.detectorConfigs.size();
0215          detectorId++) {
0216       const auto& detector = m_cfg.detectorConfigs.at(detectorId);
0217       const auto& trackFeatures = trackFeaturesVector[detectorId];
0218
0219       score += trackFeatures.nHits * detector.hitsScoreWeight;
0220       score += trackFeatures.nHoles * detector.holesScoreWeight;
0221       score += trackFeatures.nOutliers * detector.outliersScoreWeight;
0222       score += trackFeatures.nSharedHits * detector.otherScoreWeight;
0223     }
0224
0225     // Adding scores based on optional weights
0226     for (const auto& weightFunction : optionalCuts.weights) {
0227       weightFunction(track, score);
0228     }
0229
0230     // Adding the score based on the chi2/ndf
0231     if (track.chi2() > 0 && track.nDoF() > 0) {
0232       double p = 1. / std::log10(10. + 10. * track.chi2() / track.nDoF());
0233       if (p > 0) {
0234         score += p;
0235       } else {
0236         score -= 50;
0237       }
0238     }
0239
0240     iTrack++;
0241
0242     // Add the score to the vector
0243     trackScore.push_back(score);
0244     ACTS_VERBOSE("Track: " << iTrack << " score: " << score);
0245
0246   }  // end of loop over tracks
0247
0248   return trackScore;
0249 }
0250
0251 template <typename track_container_t, typename traj_t,
0252           template <typename> class holder_t, bool ReadOnly>
0253 std::vector<double> Acts::ScoreBasedAmbiguityResolution::ambiguityScore(
0254     const TrackContainer<track_container_t, traj_t, holder_t>& tracks,
0255     const std::vector<std::vector<TrackFeatures>>& trackFeaturesVectors,
0256     const OptionalCuts<track_container_t, traj_t, holder_t, ReadOnly>&
0257         optionalCuts) const {
0258   std::vector<double> trackScore;
0259   trackScore.reserve(tracks.size());
0260
0261   ACTS_VERBOSE("Using Ambiguity Scoring function");
0262
0263   int iTrack = 0;
0264
0265   ACTS_VERBOSE("Number of detectors: " << m_cfg.detectorConfigs.size());
0266
0267   ACTS_INFO("Starting to score tracks");
0268
0269   // Loop over all the tracks in the container
0270   for (const auto& track : tracks) {
0271     // get the trackFeatures map for the track
0272     const auto& trackFeaturesVector = trackFeaturesVectors[iTrack];
0273     double score = 1;
0274     auto pT = Acts::VectorHelpers::perp(track.momentum());
0275     auto eta = Acts::VectorHelpers::eta(track.momentum());
0276     auto phi = Acts::VectorHelpers::phi(track.momentum());
0277     // cuts on pT
0278     if (pT < m_cfg.pTMin || pT > m_cfg.pTMax) {
0279       score = 0;
0280       iTrack++;
0281       trackScore.push_back(score);
0282       ACTS_DEBUG("Track: " << iTrack
0283                            << " has score = 0, due to pT cuts --- pT = " << pT);
0284       continue;
0285     }
0286
0287     // cuts on phi
0288     if (phi > m_cfg.phiMax || phi < m_cfg.phiMin) {
0289       score = 0;
0290       iTrack++;
0291       trackScore.push_back(score);
0292       ACTS_DEBUG("Track: " << iTrack
0293                            << " has score = 0, due to phi cuts --- phi =  "
0294                            << phi);
0295       continue;
0296     }
0297
0298     // cuts on eta
0299     if (eta > m_cfg.etaMax || eta < m_cfg.etaMin) {
0300       score = 0;
0301       iTrack++;
0302       trackScore.push_back(score);
0303       ACTS_DEBUG("Track: " << iTrack
0304                            << " has score = 0, due to eta cuts --- eta =  "
0305                            << eta);
0306       continue;
0307     }
0308
0309     // cuts on optional cuts
0310     for (const auto& cutFunction : optionalCuts.cuts) {
0311       if (cutFunction(track)) {
0312         score = 0;
0313         ACTS_DEBUG("Track: " << iTrack
0314                              << " has score = 0, due to optional cuts.");
0315         break;
0316       }
0317     }
0318
0319     if (score == 0) {
0320       iTrack++;
0321       trackScore.push_back(score);
0322       ACTS_DEBUG("Track: " << iTrack << " score : " << score);
0323       continue;
0324     }
0325     // Reject tracks which didn't pass the detector cuts.
0326     for (std::size_t detectorId = 0; detectorId < m_cfg.detectorConfigs.size();
0327          detectorId++) {
0328       const auto& detector = m_cfg.detectorConfigs.at(detectorId);
0329
0330       const auto& trackFeatures = trackFeaturesVector[detectorId];
0331
0332       ACTS_DEBUG("---> Found summary information");
0333       ACTS_DEBUG("---> Detector ID: " << detectorId);
0334       ACTS_DEBUG("---> Number of hits: " << trackFeatures.nHits);
0335       ACTS_DEBUG("---> Number of holes: " << trackFeatures.nHoles);
0336       ACTS_DEBUG("---> Number of outliers: " << trackFeatures.nOutliers);
0337
0338       if ((trackFeatures.nHits < detector.minHits) ||
0339           (trackFeatures.nHits > detector.maxHits) ||
0340           (trackFeatures.nHoles > detector.maxHoles) ||
0341           (trackFeatures.nOutliers > detector.maxOutliers)) {
0342         score = 0;
0343         ACTS_DEBUG("Track: " << iTrack
0344                              << " has score = 0, due to detector cuts");
0345         break;
0346       }
0347     }
0348
0349     if (score == 0) {
0350       iTrack++;
0351       trackScore.push_back(score);
0352       ACTS_DEBUG("Track: " << iTrack << " score : " << score);
0353       continue;
0354     }
0355
0356     // start with larger score for tracks with higher pT.
0357     score = std::log10(pT / UnitConstants::MeV) - 1.;
0358     // pT in GeV, hence 100 MeV is minimum and gets score = 1
0359     ACTS_DEBUG("Modifier for pT = " << pT << " GeV is : " << score
0360                                     << "  New score now: " << score);
0361
0362     for (std::size_t detectorId = 0; detectorId < m_cfg.detectorConfigs.size();
0363          detectorId++) {
0364       const auto& detector = m_cfg.detectorConfigs.at(detectorId);
0365
0366       const auto& trackFeatures = trackFeaturesVector[detectorId];
0367
0368       // choosing a scaling factor based on the number of hits in a track per
0369       // detector.
0370       std::size_t nHits = trackFeatures.nHits;
0371       if (detector.factorHits.size() < nHits) {
0372         ACTS_WARNING("Detector " << detectorId
0373                                  << " has not enough factorhits in the "
0374                                     "detector.factorHits vector");
0375         continue;
0376       }
0377       if (nHits > detector.maxHits) {
0378         score = score * (detector.maxHits - nHits + 1);  // hits are good !
0379         nHits = detector.maxHits;
0380       }
0381       score = score * detector.factorHits[nHits];
0382       ACTS_DEBUG("Modifier for " << nHits
0383                                  << " hits: " << detector.factorHits[nHits]
0384                                  << "  New score now: " << score);
0385
0386       // choosing a scaling factor based on the number of holes in a track per
0387       // detector.
0388       std::size_t iHoles = trackFeatures.nHoles;
0389       if (detector.factorHoles.size() < iHoles) {
0390         ACTS_WARNING("Detector " << detectorId
0391                                  << " has not enough factorholes in the "
0392                                     "detector.factorHoles vector");
0393         continue;
0394       }
0395       if (iHoles > detector.maxHoles) {
0396         score /= (iHoles - detector.maxHoles + 1);  // holes are bad !
0397         iHoles = detector.maxHoles;
0398       }
0399       score = score * detector.factorHoles[iHoles];
0400       ACTS_DEBUG("Modifier for " << iHoles
0401                                  << " holes: " << detector.factorHoles[iHoles]
0402                                  << "  New score now: " << score);
0403     }
0404
0405     for (const auto& scoreFunction : optionalCuts.scores) {
0406       scoreFunction(track, score);
0407     }
0408
0409     if (track.chi2() > 0 && track.nDoF() > 0) {
0410       double chi2 = track.chi2();
0411       int indf = track.nDoF();
0412       double fac = 1. / std::log10(10. + 10. * chi2 / indf);
0413       score = score * fac;
0414       ACTS_DEBUG("Modifier for chi2 = " << chi2 << " and NDF = " << indf
0415                                         << " is : " << fac
0416                                         << "  New score now: " << score)
0417     }
0418     iTrack++;
0419
0420     // Add the score to the vector
0421     trackScore.push_back(score);
0422     ACTS_VERBOSE("Track: " << iTrack << " score: " << score);
0423
0424   }  // end of loop over tracks
0425
0426   return trackScore;
0427 }
0428 template <typename track_container_t, typename traj_t,
0429           template <typename> class holder_t, bool ReadOnly>
0430 std::vector<int> Acts::ScoreBasedAmbiguityResolution::solveAmbiguity(
0431     const TrackContainer<track_container_t, traj_t, holder_t>& tracks,
0432     const std::vector<std::vector<MeasurementInfo>>& measurementsPerTrack,
0433     const std::vector<std::vector<TrackFeatures>>& trackFeaturesVectors,
0434     const OptionalCuts<track_container_t, traj_t, holder_t, ReadOnly>&
0435         optionalCuts) const {
0436   ACTS_INFO("Number of tracks before Ambiguty Resolution: " << tracks.size());
0437   // vector of trackFeaturesVectors. where each trackFeaturesVector contains the
0438   // number of hits/hole/outliers for each detector in a track.
0439
0440   std::vector<double> trackScore;
0441   trackScore.reserve(tracks.size());
0442   if (m_cfg.useAmbiguityFunction) {
0443     trackScore = ambiguityScore(tracks, trackFeaturesVectors, optionalCuts);
0444   } else {
0445     trackScore = simpleScore(tracks, trackFeaturesVectors, optionalCuts);
0446   }
0447
0448   std::vector<bool> cleanTracks = getCleanedOutTracks(
0449       trackScore, trackFeaturesVectors, measurementsPerTrack);
0450
0451   std::vector<int> goodTracks;
0452   int cleanTrackIndex = 0;
0453   std::size_t iTrack = 0;
0454   for (const auto& track : tracks) {
0455     if (cleanTracks[iTrack]) {
0456       cleanTrackIndex++;
0457       if (trackScore[iTrack] >= m_cfg.minScore) {
0458         goodTracks.push_back(track.index());
0459       }
0460     }
0461     iTrack++;
0462   }
0463   ACTS_VERBOSE("Number of clean tracks: " << cleanTrackIndex);
0464   ACTS_VERBOSE("Min score: " << m_cfg.minScore);
0465   ACTS_INFO("Number of Good tracks: " << goodTracks.size());
0466   return goodTracks;
0467 }
0468
0469 }  // namespace Acts