From be67ec3dbc62d20ecd2e9fbc836299af60b5dc3c Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Tue, 19 Nov 2024 10:47:03 +0100 Subject: [PATCH 1/3] fix: Estimate entry size before sorting in ROOT readers (#3871) Extending the fix from here https://github.com/acts-project/acts/pull/3870 to other places. See also discussion here https://mattermost.web.cern.ch/acts/pl/gp9mw1y1njr4bjxmbhurk1b4me blocked by - https://github.com/acts-project/acts/pull/3870 --- Examples/Io/Root/src/RootMaterialTrackReader.cpp | 2 +- Examples/Io/Root/src/RootParticleReader.cpp | 4 ++++ Examples/Io/Root/src/RootTrackSummaryReader.cpp | 4 ++++ Examples/Io/Root/src/RootVertexReader.cpp | 12 +++++++----- Examples/Scripts/TrackingPerformance/TreeReader.h | 9 +++++++-- 5 files changed, 23 insertions(+), 8 deletions(-) diff --git a/Examples/Io/Root/src/RootMaterialTrackReader.cpp b/Examples/Io/Root/src/RootMaterialTrackReader.cpp index d1b1c443c92..70b2ab8d997 100644 --- a/Examples/Io/Root/src/RootMaterialTrackReader.cpp +++ b/Examples/Io/Root/src/RootMaterialTrackReader.cpp @@ -91,7 +91,7 @@ RootMaterialTrackReader::RootMaterialTrackReader(const Config& config, { // necessary to guarantee that m_inputChain->GetV1() is valid for the // entire range - m_inputChain->SetEstimate(nentries); + m_inputChain->SetEstimate(nentries + 1); m_entryNumbers.resize(nentries); m_inputChain->Draw("event_id", "", "goff"); diff --git a/Examples/Io/Root/src/RootParticleReader.cpp b/Examples/Io/Root/src/RootParticleReader.cpp index a76255dbb4a..631fc25ef3e 100644 --- a/Examples/Io/Root/src/RootParticleReader.cpp +++ b/Examples/Io/Root/src/RootParticleReader.cpp @@ -80,6 +80,10 @@ RootParticleReader::RootParticleReader(const RootParticleReader::Config& config, // Sort the entry numbers of the events { + // necessary to guarantee that m_inputChain->GetV1() is valid for the + // entire range + m_inputChain->SetEstimate(m_events + 1); + m_entryNumbers.resize(m_events); m_inputChain->Draw("event_id", "", "goff"); RootUtility::stableSort(m_inputChain->GetEntries(), m_inputChain->GetV1(), diff --git a/Examples/Io/Root/src/RootTrackSummaryReader.cpp b/Examples/Io/Root/src/RootTrackSummaryReader.cpp index 1827c26ad30..9ff3590271a 100644 --- a/Examples/Io/Root/src/RootTrackSummaryReader.cpp +++ b/Examples/Io/Root/src/RootTrackSummaryReader.cpp @@ -99,6 +99,10 @@ RootTrackSummaryReader::RootTrackSummaryReader( // Sort the entry numbers of the events { + // necessary to guarantee that m_inputChain->GetV1() is valid for the + // entire range + m_inputChain->SetEstimate(m_events + 1); + m_entryNumbers.resize(m_events); m_inputChain->Draw("event_nr", "", "goff"); RootUtility::stableSort(m_inputChain->GetEntries(), m_inputChain->GetV1(), diff --git a/Examples/Io/Root/src/RootVertexReader.cpp b/Examples/Io/Root/src/RootVertexReader.cpp index 336e0364f9b..31a8e6d5bec 100644 --- a/Examples/Io/Root/src/RootVertexReader.cpp +++ b/Examples/Io/Root/src/RootVertexReader.cpp @@ -8,13 +8,12 @@ #include "ActsExamples/Io/Root/RootVertexReader.hpp" -#include "Acts/Definitions/PdgParticle.hpp" #include "Acts/Utilities/Logger.hpp" #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/Framework/AlgorithmContext.hpp" +#include "ActsExamples/Io/Root/RootUtility.hpp" #include "ActsFatras/EventData/ProcessType.hpp" -#include #include #include #include @@ -64,11 +63,14 @@ RootVertexReader::RootVertexReader(const RootVertexReader::Config& config, // Sort the entry numbers of the events { + // necessary to guarantee that m_inputChain->GetV1() is valid for the + // entire range + m_inputChain->SetEstimate(m_events + 1); + m_entryNumbers.resize(m_events); m_inputChain->Draw("event_id", "", "goff"); - // Sort to get the entry numbers of the ordered events - TMath::Sort(m_inputChain->GetEntries(), m_inputChain->GetV1(), - m_entryNumbers.data(), false); + RootUtility::stableSort(m_inputChain->GetEntries(), m_inputChain->GetV1(), + m_entryNumbers.data(), false); } } diff --git a/Examples/Scripts/TrackingPerformance/TreeReader.h b/Examples/Scripts/TrackingPerformance/TreeReader.h index 4ab47746dc0..b12ce19d373 100644 --- a/Examples/Scripts/TrackingPerformance/TreeReader.h +++ b/Examples/Scripts/TrackingPerformance/TreeReader.h @@ -145,6 +145,7 @@ struct TrackStatesReader : public TreeReader { // It's not necessary if you just need to read one file, but please do it to // synchronize events if multiple root files are read if (sortEvents) { + tree->SetEstimate(tree->GetEntries() + 1); entryNumbers.resize(tree->GetEntries()); tree->Draw("event_nr", "", "goff"); // Sort to get the entry numbers of the ordered events @@ -338,6 +339,7 @@ struct TrackSummaryReader : public TreeReader { // It's not necessary if you just need to read one file, but please do it to // synchronize events if multiple root files are read if (sortEvents) { + tree->SetEstimate(tree->GetEntries() + 1); entryNumbers.resize(tree->GetEntries()); tree->Draw("event_nr", "", "goff"); // Sort to get the entry numbers of the ordered events @@ -368,7 +370,8 @@ struct TrackSummaryReader : public TreeReader { std::vector>* outlierLayer = new std::vector>; std::vector* nMajorityHits = new std::vector; - std::vector* majorityParticleId = new std::vector; + std::vector* majorityParticleId = + new std::vector; std::vector* hasFittedParams = new std::vector; @@ -427,6 +430,7 @@ struct ParticleReader : public TreeReader { // It's not necessary if you just need to read one file, but please do it to // synchronize events if multiple root files are read if (sortEvents) { + tree->SetEstimate(tree->GetEntries() + 1); entryNumbers.resize(tree->GetEntries()); tree->Draw("event_id", "", "goff"); // Sort to get the entry numbers of the ordered events @@ -436,7 +440,8 @@ struct ParticleReader : public TreeReader { } // Get all the particles with requested event id - std::vector getParticles(const std::uint32_t& eventNumber) const { + std::vector getParticles( + const std::uint32_t& eventNumber) const { // Find the start entry and the batch size for this event std::string eventNumberStr = std::to_string(eventNumber); std::string findStartEntry = "event_id<" + eventNumberStr; From f2fbbdc737f846a00002a5ebf4c7375681e9c72f Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Tue, 19 Nov 2024 13:38:46 +0100 Subject: [PATCH 2/3] fix: User `localtime_r` instead of `localtime` (#3872) This popped up as a warning in the ATLAS build: ``` /builds/acts/acts-athena-ci/acts-install/include/Acts/Utilities/Logger.hpp:485:61: warning: Non reentrant function 'localtime' called. For threadsafe applications it is recommended to use the reentrant replacement function 'localtime_r'. [localtimeCalled] std::strftime(buffer, sizeof(buffer), m_format.c_str(), localtime(&t)); ``` --- Core/include/Acts/Utilities/Logger.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Core/include/Acts/Utilities/Logger.hpp b/Core/include/Acts/Utilities/Logger.hpp index 48cadb49ec4..83dc989ca85 100644 --- a/Core/include/Acts/Utilities/Logger.hpp +++ b/Core/include/Acts/Utilities/Logger.hpp @@ -482,7 +482,9 @@ class TimedOutputDecorator final : public OutputDecorator { char buffer[20]; time_t t{}; std::time(&t); - std::strftime(buffer, sizeof(buffer), m_format.c_str(), localtime(&t)); + struct tm tbuf {}; + std::strftime(buffer, sizeof(buffer), m_format.c_str(), + localtime_r(&t, &tbuf)); return buffer; } From 2cf38f7dda859e624a7f889c0bbcf7e227c49115 Mon Sep 17 00:00:00 2001 From: Ragansu Chakkappai <66349236+Ragansu@users.noreply.github.com> Date: Tue, 19 Nov 2024 15:27:41 +0100 Subject: [PATCH 3/3] refactor!: Rewrote the ambiguity solver for clarity and added Optional Hits Selector (#3805) The Score based solver needed one more feature to be more compatible with Athena, that is optional hit selector. The optional selector enables users to add optional function to the hit selection stage which was previously NOT possible. In order to do this I had to template the existing getCleanedOutTracks. This gave me an opportunity to rewrite the whole solver again for efficiency and clarity. Ragansu Chakkappai. --- .../ScoreBasedAmbiguityResolution.hpp | 67 +++-- .../ScoreBasedAmbiguityResolution.ipp | 251 +++++++++++++----- .../ScoreBasedAmbiguityResolution.cpp | 183 ------------- ...ScoreBasedAmbiguityResolutionAlgorithm.cpp | 10 +- .../ScoreBasedAmbiguityResolutionTest.cpp | 163 ++++++++---- 5 files changed, 348 insertions(+), 326 deletions(-) diff --git a/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp b/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp index fbccb611d1a..6ed24a71be8 100644 --- a/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp +++ b/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp @@ -12,6 +12,7 @@ #include "Acts/EventData/TrackContainer.hpp" #include "Acts/EventData/TrackContainerFrontendConcept.hpp" #include "Acts/EventData/TrackProxyConcept.hpp" +#include "Acts/EventData/TrackStateProxy.hpp" #include "Acts/Utilities/Delegate.hpp" #include "Acts/Utilities/Logger.hpp" @@ -81,11 +82,17 @@ class ScoreBasedAmbiguityResolution { std::size_t nSharedHits = 0; }; - /// @brief MeasurementInfo : contains the measurement ID and the detector ID - struct MeasurementInfo { - std::size_t iMeasurement = 0; - std::size_t detectorId = 0; - bool isOutlier = false; + enum class TrackStateTypes : std::uint8_t { + // A measurement not yet used in any other track + UnsharedHit, + // A measurement shared with another track + SharedHit, + // A hit that needs to be removed from the track + RejectedHit, + // An outlier, to be copied in case + Outlier, + // Other trackstate types to be copied in case + OtherTrackStateType }; /// @brief Configuration struct : contains the configuration for the ambiguity resolution. @@ -125,11 +132,17 @@ class ScoreBasedAmbiguityResolution { using OptionalScoreModifier = std::function; + + using OptionalHitSelection = std::function; + std::vector cuts = {}; std::vector weights = {}; /// applied only if useAmbiguityFunction is true std::vector scores = {}; + std::vector hitSelections = {}; }; ScoreBasedAmbiguityResolution( @@ -141,16 +154,10 @@ class ScoreBasedAmbiguityResolution { /// Compute the initial state of the tracks. /// /// @param tracks is the input track container - /// @param sourceLinkHash is the source links - /// @param sourceLinkEquality is the equality function for the source links - /// @param trackFeaturesVectors is the trackFeatures map from detector ID to trackFeatures - /// @return a vector of the initial state of the tracks - template - std::vector> computeInitialState( - const track_container_t& tracks, source_link_hash_t sourceLinkHash, - source_link_equality_t sourceLinkEquality, - std::vector>& trackFeaturesVectors) const; + /// @return trackFeaturesVectors is the trackFeatures map from detector ID to trackFeatures + template + std::vector> computeInitialState( + const track_container_t& tracks) const; /// Compute the score of each track. /// @@ -182,29 +189,35 @@ class ScoreBasedAmbiguityResolution { /// that have a score below a certain threshold or not enough hits. /// /// @brief Remove tracks that are not good enough based on cuts + /// @param track is the input track /// @param trackScore is the score of each track - /// @param trackFeaturesVectors is the trackFeatures map for each track /// @param measurementsPerTrack is the list of measurements for each track + /// @param nTracksPerMeasurement is the number of tracks per measurement + /// @param optionalHitSelections is the optional hit selections to be applied /// @return a vector of IDs of the tracks we want to keep - std::vector getCleanedOutTracks( - const std::vector& trackScore, - const std::vector>& trackFeaturesVectors, - const std::vector>& measurementsPerTrack) - const; + template + bool getCleanedOutTracks( + const track_proxy_t& track, const double& trackScore, + const std::vector& measurementsPerTrack, + const std::map& nTracksPerMeasurement, + const std::vector>& optionalHitSelections = {}) const; /// Remove tracks that are bad based on cuts and weighted scores. /// /// @brief Remove tracks that are not good enough /// @param tracks is the input track container - /// @param measurementsPerTrack is the list of measurements for each track - /// @param trackFeaturesVectors is the map of detector id to trackFeatures for each track + /// @param sourceLinkHash is the source links + /// @param sourceLinkEquality is the equality function for the source links /// @param optionalCuts is the optional cuts to be applied /// @return a vector of IDs of the tracks we want to keep - template + template std::vector solveAmbiguity( - const track_container_t& tracks, - const std::vector>& measurementsPerTrack, - const std::vector>& trackFeaturesVectors, + const track_container_t& tracks, source_link_hash_t sourceLinkHash, + source_link_equality_t sourceLinkEquality, const OptionalCuts& optionalCuts = {}) const; diff --git a/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.ipp b/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.ipp index 58742d7d11a..247b767dd3d 100644 --- a/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.ipp +++ b/Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.ipp @@ -21,27 +21,17 @@ inline const Logger& ScoreBasedAmbiguityResolution::logger() const { return *m_logger; } -template -std::vector> +template +std::vector> ScoreBasedAmbiguityResolution::computeInitialState( - const track_container_t& tracks, source_link_hash_t sourceLinkHash, - source_link_equality_t sourceLinkEquality, - std::vector>& trackFeaturesVectors) const { - auto MeasurementIndexMap = - std::unordered_map(0, sourceLinkHash, - sourceLinkEquality); - - std::vector> measurementsPerTrack; - measurementsPerTrack.reserve(tracks.size()); + const track_container_t& tracks) const { ACTS_VERBOSE("Starting to compute initial state"); + std::vector> trackFeaturesVectors; + trackFeaturesVectors.reserve(tracks.size()); for (const auto& track : tracks) { int numberOfDetectors = m_cfg.detectorConfigs.size(); - int numberOfTrackStates = track.nTrackStates(); - std::vector measurements; - measurements.reserve(numberOfTrackStates); + std::vector trackFeaturesVector(numberOfDetectors); for (const auto& ts : track.trackStatesReversed()) { @@ -58,43 +48,25 @@ ScoreBasedAmbiguityResolution::computeInitialState( auto detectorId = volume_it->second; if (ts.typeFlags().test(Acts::TrackStateFlag::HoleFlag)) { - ACTS_DEBUG("Track state type is HoleFlag"); + ACTS_VERBOSE("Track state type is HoleFlag"); trackFeaturesVector[detectorId].nHoles++; } else if (ts.typeFlags().test(Acts::TrackStateFlag::OutlierFlag)) { - Acts::SourceLink sourceLink = ts.getUncalibratedSourceLink(); - ACTS_DEBUG("Track state type is OutlierFlag"); + ACTS_VERBOSE("Track state type is OutlierFlag"); trackFeaturesVector[detectorId].nOutliers++; - // assign a new measurement index if the source link was not seen yet - auto emplace = MeasurementIndexMap.try_emplace( - sourceLink, MeasurementIndexMap.size()); - - bool isOutliner = true; - - measurements.push_back({emplace.first->second, detectorId, isOutliner}); } else if (ts.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) { - Acts::SourceLink sourceLink = ts.getUncalibratedSourceLink(); - ACTS_DEBUG("Track state type is MeasurementFlag"); + ACTS_VERBOSE("Track state type is MeasurementFlag"); if (ts.typeFlags().test(Acts::TrackStateFlag::SharedHitFlag)) { trackFeaturesVector[detectorId].nSharedHits++; } trackFeaturesVector[detectorId].nHits++; - - // assign a new measurement index if the source link was not seen yet - auto emplace = MeasurementIndexMap.try_emplace( - sourceLink, MeasurementIndexMap.size()); - - bool isoutliner = false; - - measurements.push_back({emplace.first->second, detectorId, isoutliner}); } } - measurementsPerTrack.push_back(std::move(measurements)); trackFeaturesVectors.push_back(std::move(trackFeaturesVector)); } - return measurementsPerTrack; + return trackFeaturesVectors; } template @@ -363,12 +335,6 @@ std::vector Acts::ScoreBasedAmbiguityResolution::ambiguityScore( // choosing a scaling factor based on the number of hits in a track per // detector. std::size_t nHits = trackFeatures.nHits; - if (detector.factorHits.size() < nHits) { - ACTS_WARNING("Detector " << detectorId - << " has not enough factorhits in the " - "detector.factorHits vector"); - continue; - } if (nHits > detector.maxHits) { score = score * (detector.maxHits - nHits + 1); // hits are good ! nHits = detector.maxHits; @@ -381,12 +347,6 @@ std::vector Acts::ScoreBasedAmbiguityResolution::ambiguityScore( // choosing a scaling factor based on the number of holes in a track per // detector. std::size_t iHoles = trackFeatures.nHoles; - if (detector.factorHoles.size() < iHoles) { - ACTS_WARNING("Detector " << detectorId - << " has not enough factorholes in the " - "detector.factorHoles vector"); - continue; - } if (iHoles > detector.maxHoles) { score /= (iHoles - detector.maxHoles + 1); // holes are bad ! iHoles = detector.maxHoles; @@ -421,17 +381,20 @@ std::vector Acts::ScoreBasedAmbiguityResolution::ambiguityScore( return trackScore; } -template +template std::vector Acts::ScoreBasedAmbiguityResolution::solveAmbiguity( - const track_container_t& tracks, - const std::vector>& measurementsPerTrack, - const std::vector>& trackFeaturesVectors, + const track_container_t& tracks, source_link_hash_t sourceLinkHash, + source_link_equality_t sourceLinkEquality, const OptionalCuts& optionalCuts) const { ACTS_INFO("Number of tracks before Ambiguty Resolution: " << tracks.size()); // vector of trackFeaturesVectors. where each trackFeaturesVector contains the // number of hits/hole/outliers for each detector in a track. + const std::vector> trackFeaturesVectors = + computeInitialState(tracks); + std::vector trackScore; trackScore.reserve(tracks.size()); if (m_cfg.useAmbiguityFunction) { @@ -440,14 +403,72 @@ std::vector Acts::ScoreBasedAmbiguityResolution::solveAmbiguity( trackScore = simpleScore(tracks, trackFeaturesVectors, optionalCuts); } - std::vector cleanTracks = getCleanedOutTracks( - trackScore, trackFeaturesVectors, measurementsPerTrack); + auto MeasurementIndexMap = + std::unordered_map(0, sourceLinkHash, + sourceLinkEquality); + + std::vector> measurementsPerTrackVector; + std::map nTracksPerMeasurement; + + // Stores tracks measurement into a vector or vectors + // (measurementsPerTrackVector) and counts the number of tracks per + // measurement.(nTracksPerMeasurement) + + for (const auto& track : tracks) { + std::vector measurementsPerTrack; + for (const auto& ts : track.trackStatesReversed()) { + if (!ts.typeFlags().test(Acts::TrackStateFlag::OutlierFlag) && + !ts.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) { + continue; + } + Acts::SourceLink sourceLink = ts.getUncalibratedSourceLink(); + // assign a new measurement index if the source link was not seen yet + auto emplace = MeasurementIndexMap.try_emplace( + sourceLink, MeasurementIndexMap.size()); + std::size_t iMeasurement = emplace.first->second; + measurementsPerTrack.push_back(iMeasurement); + if (nTracksPerMeasurement.find(iMeasurement) == + nTracksPerMeasurement.end()) { + nTracksPerMeasurement[iMeasurement] = 0; + } + nTracksPerMeasurement[iMeasurement]++; + } + measurementsPerTrackVector.push_back(std::move(measurementsPerTrack)); + } std::vector goodTracks; int cleanTrackIndex = 0; - std::size_t iTrack = 0; - for (const auto& track : tracks) { - if (cleanTracks[iTrack]) { + + auto optionalHitSelections = optionalCuts.hitSelections; + + // Loop over all the tracks in the container + // For each track, check if the track has too many shared hits to be accepted. + // If the track is accepted, remove bad hits and check if the track has a + // score higher than the minimum score for shared tracks. + // If the track is good, add it to the goodTracks vector. + for (std::size_t iTrack = 0; const auto& track : tracks) { + // Check if the track has too many shared hits to be accepted. + auto trackFeaturesVector = trackFeaturesVectors.at(iTrack); + bool trkCouldBeAccepted = true; + for (std::size_t detectorId = 0; detectorId < m_cfg.detectorConfigs.size(); + detectorId++) { + auto detector = m_cfg.detectorConfigs.at(detectorId); + if (trackFeaturesVector[detectorId].nSharedHits > + detector.maxSharedHits) { + trkCouldBeAccepted = false; + break; + } + } + if (!trkCouldBeAccepted) { + iTrack++; + continue; + } + + trkCouldBeAccepted = getCleanedOutTracks( + track, trackScore[iTrack], measurementsPerTrackVector[iTrack], + nTracksPerMeasurement, optionalHitSelections); + if (trkCouldBeAccepted) { cleanTrackIndex++; if (trackScore[iTrack] >= m_cfg.minScore) { goodTracks.push_back(track.index()); @@ -455,10 +476,120 @@ std::vector Acts::ScoreBasedAmbiguityResolution::solveAmbiguity( } iTrack++; } - ACTS_VERBOSE("Number of clean tracks: " << cleanTrackIndex); + ACTS_INFO("Number of clean tracks: " << cleanTrackIndex); ACTS_VERBOSE("Min score: " << m_cfg.minScore); ACTS_INFO("Number of Good tracks: " << goodTracks.size()); return goodTracks; } +template +bool Acts::ScoreBasedAmbiguityResolution::getCleanedOutTracks( + const track_proxy_t& track, const double& trackScore, + const std::vector& measurementsPerTrack, + const std::map& nTracksPerMeasurement, + const std::vector< + std::function>& optionalHitSelections) const { + // For tracks with shared hits, we need to check and remove bad hits + + std::vector trackStateTypes; + // Loop over all measurements of the track and for each hit a + // trackStateTypes is assigned. + for (std::size_t index = 0; const auto& ts : track.trackStatesReversed()) { + if (ts.typeFlags().test(Acts::TrackStateFlag::OutlierFlag) || + ts.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) { + std::size_t iMeasurement = measurementsPerTrack[index]; + auto it = nTracksPerMeasurement.find(iMeasurement); + if (it == nTracksPerMeasurement.end()) { + trackStateTypes.push_back(TrackStateTypes::OtherTrackStateType); + index++; + continue; + } + + std::size_t nTracksShared = it->second; + auto isoutliner = ts.typeFlags().test(Acts::TrackStateFlag::OutlierFlag); + + if (isoutliner) { + ACTS_VERBOSE("Measurement is outlier on a fitter track, copy it over"); + trackStateTypes.push_back(TrackStateTypes::Outlier); + continue; + } + if (nTracksShared == 1) { + ACTS_VERBOSE("Measurement is not shared, copy it over"); + + trackStateTypes.push_back(TrackStateTypes::UnsharedHit); + continue; + } else if (nTracksShared > 1) { + ACTS_VERBOSE("Measurement is shared, copy it over"); + trackStateTypes.push_back(TrackStateTypes::SharedHit); + continue; + } + } + } + std::vector newMeasurementsPerTrack; + std::size_t measurement = 0; + std::size_t nshared = 0; + + // Loop over all measurements of the track and process them according to the + // trackStateTypes and other conditions. + // Good measurements are copied to the newMeasurementsPerTrack vector. + for (std::size_t index = 0; auto ts : track.trackStatesReversed()) { + if (ts.typeFlags().test(Acts::TrackStateFlag::OutlierFlag) || + ts.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) { + if (!ts.hasReferenceSurface()) { + ACTS_DEBUG("Track state has no reference surface"); + continue; + } + measurement = measurementsPerTrack[index]; + + auto it = nTracksPerMeasurement.find(measurement); + if (it == nTracksPerMeasurement.end()) { + index++; + continue; + } + auto nTracksShared = it->second; + + // Loop over all optionalHitSelections and apply them to trackStateType of + // the TrackState. + + for (const auto& hitSelection : optionalHitSelections) { + hitSelection(track, ts, trackStateTypes[index]); + } + + if (trackStateTypes[index] == TrackStateTypes::RejectedHit) { + ACTS_DEBUG("Dropping rejected hit"); + } else if (trackStateTypes[index] != TrackStateTypes::SharedHit) { + ACTS_DEBUG("Good TSOS, copy hit"); + newMeasurementsPerTrack.push_back(measurement); + + // a counter called nshared is used to keep track of the number of + // shared hits accepted. + } else if (nshared >= m_cfg.maxShared) { + ACTS_DEBUG("Too many shared hit, drop it"); + } + // If the track is shared, the hit is only accepted if the track has + // score higher than the minimum score for shared tracks. + else { + ACTS_DEBUG("Try to recover shared hit "); + if (nTracksShared <= m_cfg.maxSharedTracksPerMeasurement && + trackScore > m_cfg.minScoreSharedTracks) { + ACTS_DEBUG("Accepted hit shared with " << nTracksShared << " tracks"); + newMeasurementsPerTrack.push_back(measurement); + nshared++; + } else { + ACTS_DEBUG("Rejected hit shared with " << nTracksShared << " tracks"); + } + } + index++; + } + } + // Check if the track has enough hits to be accepted. + if (newMeasurementsPerTrack.size() < 3) { + return false; + } else { + return true; + } +} + } // namespace Acts diff --git a/Core/src/AmbiguityResolution/ScoreBasedAmbiguityResolution.cpp b/Core/src/AmbiguityResolution/ScoreBasedAmbiguityResolution.cpp index 4dcbb5071d4..c7e47f41046 100644 --- a/Core/src/AmbiguityResolution/ScoreBasedAmbiguityResolution.cpp +++ b/Core/src/AmbiguityResolution/ScoreBasedAmbiguityResolution.cpp @@ -7,186 +7,3 @@ // file, You can obtain one at https://mozilla.org/MPL/2.0/. #include "Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp" - -#include "Acts/EventData/SourceLink.hpp" - -#include - -std::vector Acts::ScoreBasedAmbiguityResolution::getCleanedOutTracks( - const std::vector& trackScore, - const std::vector>& trackFeaturesVectors, - const std::vector>& measurementsPerTrack) - const { - std::vector cleanTracks(measurementsPerTrack.size(), false); - - ACTS_VERBOSE("Cleaning tracks"); - - if (trackScore.size() != measurementsPerTrack.size()) { - throw std::invalid_argument( - "Track score and measurementsPerTrack size mismatch"); - } - - std::size_t numberOfTracks = measurementsPerTrack.size(); - ACTS_DEBUG("Number of tracks: " << numberOfTracks); - - boost::container::flat_map> - tracksPerMeasurement; - - // Removes bad tracks and counts computes the vector of tracks per - // measurement. - for (std::size_t iTrack = 0; iTrack < numberOfTracks; ++iTrack) { - if (trackScore[iTrack] <= 0) { - continue; - } - for (auto measurementObjects : measurementsPerTrack[iTrack]) { - auto iMeasurement = measurementObjects.iMeasurement; - tracksPerMeasurement[iMeasurement].insert(iTrack); - } - } - - enum TrackStateTypes { - // A measurement not yet used in any other track - UnsharedHit = 1, - // A measurement shared with another track - SharedHit = 2, - // A hit that needs to be removed from the track - RejectedHit = 3, - // an outlier, to be copied in case - Outlier = 4, - // other trackstate types to be copied in case - OtherTrackStateType = 5 - }; - - std::vector> newMeasurements; - // Loop over all tracks in the track container - for (std::size_t iTrack = 0; iTrack < numberOfTracks; ++iTrack) { - double track_score = trackScore[iTrack]; - ACTS_DEBUG("Track score: " << track_score); - - if (track_score <= 0) { - ACTS_DEBUG("Track " << iTrack << " could not be accepted - low score"); - continue; - } - - const auto& trackFeaturesVector = trackFeaturesVectors.at(iTrack); - - bool trkCouldBeAccepted = true; - - // For tracks with shared hits, we need to check and remove bad hits - - std::vector trackStateTypes(measurementsPerTrack[iTrack].size(), - OtherTrackStateType); - int index = 0; - - // Loop over all measurements of the track and for each hit a - // trackStateTypes is assigned. - for (const auto& measurementObjects : measurementsPerTrack[iTrack]) { - auto iMeasurement = measurementObjects.iMeasurement; - auto isoutliner = measurementObjects.isOutlier; - auto detectorId = measurementObjects.detectorId; - - auto detector = m_cfg.detectorConfigs.at(detectorId); - if (isoutliner) { - ACTS_VERBOSE("Measurement is outlier on a fitter track, copy it over"); - trackStateTypes[index] = Outlier; - index++; - continue; - } - if (tracksPerMeasurement[iMeasurement].size() == 1) { - ACTS_VERBOSE("Measurement is not shared, copy it over"); - - trackStateTypes[index] = UnsharedHit; - - index++; - continue; - } - if (tracksPerMeasurement[iMeasurement].size() > 1) { - ACTS_VERBOSE("Measurement is shared, copy it over"); - - if (detector.sharedHitsFlag == true) { - ACTS_VERBOSE("Measurement is shared, Reject it"); - trackStateTypes[index] = RejectedHit; - index++; - continue; - } - - trackStateTypes[index] = SharedHit; - - index++; - continue; - } - } - - std::vector newMeasurementsPerTrack; - std::size_t measurement = 0; - std::size_t nshared = 0; - - // Loop over all measurements of the track and process them according to the - // trackStateTypes and other conditions. - // Good measurements are copied to the newMeasurementsPerTrack vector. - for (std::size_t i = 0; i < trackStateTypes.size(); i++) { - auto& measurementObjects = measurementsPerTrack[iTrack][i]; - measurement = measurementObjects.iMeasurement; - - if (trackStateTypes[i] == RejectedHit) { - ACTS_DEBUG("Dropping rejected hit"); - } else if (trackStateTypes[i] != SharedHit) { - ACTS_DEBUG("Good TSOS, copy hit"); - newMeasurementsPerTrack.push_back(measurement); - - // a counter called nshared is used to keep track of the number of - // shared hits accepted. - } else if (nshared >= m_cfg.maxShared) { - ACTS_DEBUG("Too many shared hit, drop it"); - } - // If the track is shared, the hit is only accepted if the track has score - // higher than the minimum score for shared tracks. - else { - ACTS_DEBUG("Try to recover shared hit "); - if (tracksPerMeasurement[measurement].size() < - m_cfg.maxSharedTracksPerMeasurement && - track_score > m_cfg.minScoreSharedTracks) { - ACTS_DEBUG("Accepted hit shared with " - << tracksPerMeasurement[measurement].size() << " tracks"); - newMeasurementsPerTrack.push_back(measurement); - nshared++; - } else { - ACTS_DEBUG("Rejected hit shared with " - << tracksPerMeasurement[measurement].size() << " tracks"); - } - } - } - - // Check if the track has enough hits to be accepted. - if (newMeasurementsPerTrack.size() < 3) { - trkCouldBeAccepted = false; - ACTS_DEBUG(std::endl - << "Track " << iTrack - << " could not be accepted - not enough hits"); - ACTS_DEBUG("Number of hits: " << measurementsPerTrack[iTrack].size()); - ACTS_DEBUG("Number of good hits: " << newMeasurementsPerTrack.size()); - continue; - } - - // Check if the track has too many shared hits to be accepted. - for (std::size_t detectorId = 0; detectorId < m_cfg.detectorConfigs.size(); - detectorId++) { - auto detector = m_cfg.detectorConfigs.at(detectorId); - if (trackFeaturesVector[detectorId].nSharedHits > - detector.maxSharedHits) { - trkCouldBeAccepted = false; - break; - } - } - - if (trkCouldBeAccepted) { - cleanTracks[iTrack] = true; - newMeasurements.push_back(newMeasurementsPerTrack); - ACTS_VERBOSE("Track " << iTrack << " is accepted"); - continue; - } - } - - return cleanTracks; -} diff --git a/Examples/Algorithms/AmbiguityResolution/src/ScoreBasedAmbiguityResolutionAlgorithm.cpp b/Examples/Algorithms/AmbiguityResolution/src/ScoreBasedAmbiguityResolutionAlgorithm.cpp index bd3584fa4b2..b1a6e251dba 100644 --- a/Examples/Algorithms/AmbiguityResolution/src/ScoreBasedAmbiguityResolutionAlgorithm.cpp +++ b/Examples/Algorithms/AmbiguityResolution/src/ScoreBasedAmbiguityResolutionAlgorithm.cpp @@ -116,19 +116,11 @@ ActsExamples::ScoreBasedAmbiguityResolutionAlgorithm::execute( const auto& tracks = m_inputTracks(ctx); // Read input data ACTS_VERBOSE("Number of input tracks: " << tracks.size()); - std::vector> - measurementsPerTracks; - - std::vector> - trackFeaturesVectors; - measurementsPerTracks = m_ambi.computeInitialState( - tracks, &sourceLinkHash, &sourceLinkEquality, trackFeaturesVectors); - Acts::ScoreBasedAmbiguityResolution::OptionalCuts optionalCuts; optionalCuts.cuts.push_back(doubleHolesFilter); std::vector goodTracks = m_ambi.solveAmbiguity( - tracks, measurementsPerTracks, trackFeaturesVectors, optionalCuts); + tracks, &sourceLinkHash, &sourceLinkEquality, optionalCuts); // Prepare the output track collection from the IDs TrackContainer solvedTracks{std::make_shared(), std::make_shared()}; diff --git a/Tests/UnitTests/Core/AmbiguityResolution/ScoreBasedAmbiguityResolutionTest.cpp b/Tests/UnitTests/Core/AmbiguityResolution/ScoreBasedAmbiguityResolutionTest.cpp index 568e817031a..55beab4ca2f 100644 --- a/Tests/UnitTests/Core/AmbiguityResolution/ScoreBasedAmbiguityResolutionTest.cpp +++ b/Tests/UnitTests/Core/AmbiguityResolution/ScoreBasedAmbiguityResolutionTest.cpp @@ -18,10 +18,7 @@ #include "Acts/EventData/VectorTrackContainer.hpp" #include "Acts/EventData/detail/TestSourceLink.hpp" #include "Acts/EventData/detail/TestTrackState.hpp" -#include "Acts/Geometry/GeometryContext.hpp" -#include "Acts/Surfaces/PerigeeSurface.hpp" -#include "Acts/TrackFinding/TrackSelector.hpp" -#include "Acts/Utilities/CalibrationContext.hpp" +#include "Acts/Utilities/TrackHelpers.hpp" #include @@ -30,7 +27,6 @@ using Acts::MultiTrajectoryTraits::IndexType; namespace Acts::Test { BOOST_AUTO_TEST_SUITE(ScoreBasedAmbiguityResolutionTest) -using MeasurementInfo = ScoreBasedAmbiguityResolution::MeasurementInfo; // Test fixture for ScoreBasedAmbiguityResolution struct Fixture { @@ -65,64 +61,137 @@ struct Fixture { ~Fixture() = default; }; -// Helper function to create a sample input for getCleanedOutTracks -std::vector> createSampleInput() { - Fixture fixture; - std::vector>> trackVolumes = { - {0, {19, 18, 18, 18, 10, 10, 10, 10, 10, 10, 10, 10, 10}}, - {1, {19, 18, 18, 18, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10}}, - {2, {13, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8}}, - {3, {13, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8}}, - {4, {19, 18, 18, 18, 10, 10, 10, 10, 10, 10, 10, 10, 10}}}; - - std::vector> measurementsPerTrack; - // Add sample measurements for each track - - for (const auto& trackVolume : trackVolumes) { - std::vector measurements; - for (std::size_t i = 0; i < trackVolume.second.size(); ++i) { - std::size_t detectorID = fixture.config.volumeMap[trackVolume.second[i]]; - measurements.push_back({i + 2, detectorID, false}); +template +auto createTestTrack(TrackContainer& tc, const FlagsPerState& flagsPerState) { + auto t = tc.makeTrack(); + for (const auto& flags : flagsPerState) { + auto ts = t.appendTrackState(); + for (auto f : flags) { + ts.typeFlags().set(f); } - measurementsPerTrack.push_back(measurements); } - return measurementsPerTrack; + calculateTrackQuantities(t); + + return t; } -BOOST_FIXTURE_TEST_CASE(GetCleanedOutTracksTest, Fixture) { +BOOST_FIXTURE_TEST_CASE(ComputeInitialStateTest, Fixture) { Fixture fixture; // Create an instance of ScoreBasedAmbiguityResolution ScoreBasedAmbiguityResolution tester(fixture.config); - // Create sample input - std::vector> measurementsPerTrack = - createSampleInput(); + VectorTrackContainer mutVtc; + VectorMultiTrajectory mutMtj; - std::vector TrackSore; - for (std::size_t i = 0; i < measurementsPerTrack.size(); i++) { - TrackSore.push_back(60 + 40 * i); - } - std::vector> trackFeaturesVectors = { - {{0, 14, 0, 0}, {0, 2, 0, 0}}, - {{0, 15, 0, 0}, {0, 2, 0, 0}}, - {{0, 17, 0, 0}, {0, 2, 0, 0}}, - {{0, 18, 0, 0}, {0, 2, 0, 0}}, - {{0, 14, 0, 0}, {0, 3, 0, 0}}}; + TrackContainer mutTc{mutVtc, mutMtj}; + static_assert(!mutTc.ReadOnly, "Unexpectedly read only"); + + auto t = createTestTrack(mutTc, std::vector>{ + {MeasurementFlag}, + {OutlierFlag}, + {MeasurementFlag, SharedHitFlag}, + {HoleFlag}, + {OutlierFlag}, + {HoleFlag}, + {MeasurementFlag, SharedHitFlag}, + {OutlierFlag}, + }); + + BOOST_CHECK_EQUAL(t.nHoles(), 2); + BOOST_CHECK_EQUAL(t.nMeasurements(), 3); + BOOST_CHECK_EQUAL(t.nOutliers(), 3); + BOOST_CHECK_EQUAL(t.nSharedHits(), 2); + + ConstVectorTrackContainer vtc{std::move(mutVtc)}; + ConstVectorMultiTrajectory mtj{std::move(mutMtj)}; + + TrackContainer ctc{vtc, mtj}; + + std::vector> trackFeaturesVectors; + trackFeaturesVectors = tester.computeInitialState(ctc); + + BOOST_CHECK_EQUAL(trackFeaturesVectors.size(), ctc.size()); + + std::vector trackScores; + trackScores = tester.simpleScore(ctc, trackFeaturesVectors); - // Call the function under testBOOST_FIXTURE_TEST_CASE - std::vector cleanTracks = tester.getCleanedOutTracks( - TrackSore, trackFeaturesVectors, measurementsPerTrack); + BOOST_CHECK_EQUAL(trackScores.size(), ctc.size()); + + trackScores = tester.ambiguityScore(ctc, trackFeaturesVectors); + + BOOST_CHECK_EQUAL(trackScores.size(), ctc.size()); // Assert the expected results - BOOST_CHECK_EQUAL(measurementsPerTrack.size(), 5); +} - for (std::size_t i = 0; i < cleanTracks.size(); i++) { - auto score = TrackSore[i]; - if (cleanTracks[i]) { - BOOST_CHECK_GT(score, fixture.config.minScoreSharedTracks); +BOOST_FIXTURE_TEST_CASE(GetCleanedOutTracksTest, Fixture) { + Fixture fixture; + // Create an instance of ScoreBasedAmbiguityResolution + ScoreBasedAmbiguityResolution tester(fixture.config); + + VectorTrackContainer mutVtc; + VectorMultiTrajectory mutMtj; + + TrackContainer mutTc{mutVtc, mutMtj}; + static_assert(!mutTc.ReadOnly, "Unexpectedly read only"); + + auto t = createTestTrack(mutTc, std::vector>{ + {MeasurementFlag}, + {OutlierFlag}, + {MeasurementFlag, SharedHitFlag}, + {HoleFlag}, + {OutlierFlag}, + {HoleFlag}, + {MeasurementFlag, SharedHitFlag}, + {OutlierFlag}, + }); + + BOOST_CHECK_EQUAL(t.nHoles(), 2); + BOOST_CHECK_EQUAL(t.nMeasurements(), 3); + BOOST_CHECK_EQUAL(t.nOutliers(), 3); + BOOST_CHECK_EQUAL(t.nSharedHits(), 2); + + ConstVectorTrackContainer vtc{std::move(mutVtc)}; + ConstVectorMultiTrajectory mtj{std::move(mutMtj)}; + + TrackContainer ctc{vtc, mtj}; + + std::vector> trackFeaturesVectors; + std::vector> measurementsPerTrackVector; + std::map nTracksPerMeasurement; + + trackFeaturesVectors = tester.computeInitialState(ctc); + + BOOST_CHECK_EQUAL(trackFeaturesVectors.size(), ctc.size()); + + for (std::size_t iMeasurement = 1; const auto& track : ctc) { + std::vector measurementsPerTrack; + for (auto ts : track.trackStatesReversed()) { + if (ts.typeFlags().test(Acts::TrackStateFlag::OutlierFlag) || + ts.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) { + measurementsPerTrack.push_back(iMeasurement); + if (nTracksPerMeasurement.find(iMeasurement) == + nTracksPerMeasurement.end()) { + nTracksPerMeasurement[iMeasurement] = 0; + } + nTracksPerMeasurement[iMeasurement]++; + } + + iMeasurement++; } + measurementsPerTrackVector.push_back(std::move(measurementsPerTrack)); } + + auto trackScore = tester.ambiguityScore(ctc, trackFeaturesVectors); + + const auto& track = ctc.getTrack(0); + std::size_t iTrack = track.index(); + bool accepted = tester.getCleanedOutTracks(track, trackScore[iTrack], + measurementsPerTrackVector[iTrack], + nTracksPerMeasurement); + + BOOST_CHECK_EQUAL(accepted, false); } BOOST_AUTO_TEST_SUITE_END()