diff --git a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp index b97cbb1a93b..40555986ce4 100644 --- a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp +++ b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp @@ -863,79 +863,57 @@ class CombinatorialKalmanFilter { auto currentBranch = result.activeBranches.back(); TrackIndexType prevTip = currentBranch.tipIndex(); - // Add a state only if there already is measurement on this branch - if (currentBranch.nMeasurements() > 0) { - ACTS_VERBOSE( - "Record hole or passive material state on surface since we " - "already have measurements on the branch"); - - // No source links on surface, add either hole or passive material - // TrackState. No storage allocation for uncalibrated/calibrated - // measurement and filtered parameter - auto stateMask = PM::Predicted | PM::Jacobian; - - // Transport the covariance to a curvilinear surface - stepper.transportCovarianceToCurvilinear(state.stepping); + // No source links on surface, add either hole or passive material + // TrackState. No storage allocation for uncalibrated/calibrated + // measurement and filtered parameter + auto stateMask = PM::Predicted | PM::Jacobian; - // Update state and stepper with pre material effects - materialInteractor(surface, state, stepper, navigator, - MaterialUpdateStage::PreUpdate); + // Transport the covariance to a curvilinear surface + stepper.transportCovarianceToCurvilinear(state.stepping); - // Transport & bind the state to the current surface - auto boundStateRes = - stepper.boundState(state.stepping, *surface, false); - if (!boundStateRes.ok()) { - return boundStateRes.error(); - } - auto& boundState = *boundStateRes; - auto& [boundParams, jacobian, pathLength] = boundState; - boundParams.covariance() = state.stepping.cov; - - // Add a hole or material track state - TrackIndexType currentTip = addNonSourcelinkState( - stateMask, boundState, result, isSensitive, prevTip); - auto nonSourcelinkState = - result.trackStates->getTrackState(currentTip); - currentBranch.tipIndex() = currentTip; - if (isSensitive) { - currentBranch.nHoles()++; - } + // Update state and stepper with pre material effects + materialInteractor(surface, state, stepper, navigator, + MaterialUpdateStage::PreUpdate); - BranchStopperResult branchStopperResult = - m_extensions.branchStopper(currentBranch, nonSourcelinkState); + // Transport & bind the state to the current surface + auto boundStateRes = + stepper.boundState(state.stepping, *surface, false); + if (!boundStateRes.ok()) { + return boundStateRes.error(); + } + auto& boundState = *boundStateRes; + auto& [boundParams, jacobian, pathLength] = boundState; + boundParams.covariance() = state.stepping.cov; - // Check the branch - if (branchStopperResult == BranchStopperResult::Continue) { - // Remembered the active branch and its state - } else { - // No branch on this surface - nBranchesOnSurface = 0; - if (branchStopperResult == BranchStopperResult::StopAndKeep) { - storeLastActiveBranch(result); - } - // Remove the branch from list - result.activeBranches.pop_back(); - } + // Add a hole or material track state to the multitrajectory + TrackIndexType currentTip = addNonSourcelinkState( + stateMask, boundState, result, isSensitive, prevTip); + auto nonSourcelinkState = result.trackStates->getTrackState(currentTip); + currentBranch.tipIndex() = currentTip; - // Update state and stepper with post material effects - materialInteractor(surface, state, stepper, navigator, - MaterialUpdateStage::PostUpdate); - } else if (isMaterial) { - // Since we did not find a measurement yet we don't create a track - // state but still need to do the material update - ACTS_VERBOSE("Handle material without creating a track state"); + if (isSensitive) { + currentBranch.nHoles()++; + } - // Transport the covariance to a curvilinear surface - stepper.transportCovarianceToCurvilinear(state.stepping); + BranchStopperResult branchStopperResult = + m_extensions.branchStopper(currentBranch, nonSourcelinkState); - // Update state and stepper with full material effects - materialInteractor(surface, state, stepper, navigator, - MaterialUpdateStage::FullUpdate); + // Check the branch + if (branchStopperResult == BranchStopperResult::Continue) { + // Remembered the active branch and its state } else { - ACTS_VERBOSE( - "Skipping surface as there are no measurements on the branch yet " - "and it does not have material"); + // No branch on this surface + nBranchesOnSurface = 0; + if (branchStopperResult == BranchStopperResult::StopAndKeep) { + storeLastActiveBranch(result); + } + // Remove the branch from list + result.activeBranches.pop_back(); } + + // Update state and stepper with post material effects + materialInteractor(surface, state, stepper, navigator, + MaterialUpdateStage::PostUpdate); } else { // Neither measurement nor material on surface, this branch is still // valid. Count the branch on current surface @@ -1166,31 +1144,13 @@ class CombinatorialKalmanFilter { auto currentBranch = result.activeBranches.back(); TrackIndexType currentTip = currentBranch.tipIndex(); - ACTS_VERBOSE("Find track with entry index = " - << currentTip << " and there are nMeasurements = " - << currentBranch.nMeasurements() + ACTS_VERBOSE("Storing track " + << currentBranch.index() << " with tip index " << currentTip + << ". nMeasurements = " << currentBranch.nMeasurements() << ", nOutliers = " << currentBranch.nOutliers() - << ", nHoles = " << currentBranch.nHoles() << " on track"); - - std::optional lastMeasurement; - for (const auto& trackState : currentBranch.trackStatesReversed()) { - if (trackState.typeFlags().test(TrackStateFlag::MeasurementFlag) && - !trackState.typeFlags().test(TrackStateFlag::OutlierFlag)) { - lastMeasurement = trackState; - break; - } - } + << ", nHoles = " << currentBranch.nHoles()); - if (lastMeasurement.has_value()) { - currentBranch.tipIndex() = lastMeasurement->index(); - result.collectedTracks.push_back(currentBranch); - ACTS_VERBOSE("Last measurement found on track with entry index = " - << currentTip << " and measurement index = " - << lastMeasurement->index()); - } else { - ACTS_VERBOSE( - "No measurement found on track with entry index = " << currentTip); - } + result.collectedTracks.push_back(currentBranch); } CombinatorialKalmanFilterExtensions m_extensions; diff --git a/Core/include/Acts/Utilities/TrackHelpers.hpp b/Core/include/Acts/Utilities/TrackHelpers.hpp index 11589981b9b..69e91b5a738 100644 --- a/Core/include/Acts/Utilities/TrackHelpers.hpp +++ b/Core/include/Acts/Utilities/TrackHelpers.hpp @@ -368,13 +368,14 @@ Result extrapolateTracksToReferenceSurface( /// Helper function to calculate a number of track level quantities and store /// them on the track itself -/// @note The input track needs to be mutable, so @c ReadOnly=false -/// @tparam track_container_t the track container backend -/// @tparam track_state_container_t the track state container backend -/// @tparam holder_t the holder type for the track container backends +/// @tparam track_proxy_t The track proxy type /// @param track A mutable track proxy to operate on template -void calculateTrackQuantities(track_proxy_t track) { +void calculateTrackQuantities(track_proxy_t track) + requires(!track_proxy_t::ReadOnly) +{ + using ConstTrackStateProxy = typename track_proxy_t::ConstTrackStateProxy; + track.chi2() = 0; track.nDoF() = 0; @@ -383,8 +384,8 @@ void calculateTrackQuantities(track_proxy_t track) { track.nSharedHits() = 0; track.nOutliers() = 0; - for (const auto &trackState : track.trackStatesReversed()) { - auto typeFlags = trackState.typeFlags(); + for (ConstTrackStateProxy trackState : track.trackStatesReversed()) { + ConstTrackStateType typeFlags = trackState.typeFlags(); if (typeFlags.test(Acts::TrackStateFlag::HoleFlag)) { track.nHoles()++; @@ -401,6 +402,97 @@ void calculateTrackQuantities(track_proxy_t track) { } } +/// Helper function to trim track states from the front of a track +/// @tparam track_proxy_t the track proxy type +/// @param track the track to trim +/// @param trimHoles whether to trim holes +/// @param trimOutliers whether to trim outliers +/// @param trimMaterial whether to trim pure material states +template +void trimTrackFront(track_proxy_t track, bool trimHoles, bool trimOutliers, + bool trimMaterial) + requires(!track_proxy_t::ReadOnly) +{ + using TrackStateProxy = typename track_proxy_t::TrackStateProxy; + + // TODO specialize if track is forward linked + + std::optional front; + + for (TrackStateProxy trackState : track.trackStatesReversed()) { + TrackStateType typeFlags = trackState.typeFlags(); + if (trimHoles && typeFlags.test(TrackStateFlag::HoleFlag)) { + continue; + } + if (trimOutliers && typeFlags.test(TrackStateFlag::OutlierFlag)) { + continue; + } + if (trimMaterial && typeFlags.test(TrackStateFlag::MaterialFlag) && + !typeFlags.test(TrackStateFlag::MeasurementFlag)) { + continue; + } + + front = trackState; + } + + if (front.has_value()) { + front.value().previous() = TrackStateProxy::kInvalid; + } +} + +/// Helper function to trim track states from the back of a track +/// @tparam track_proxy_t the track proxy type +/// @param track the track to trim +/// @param trimHoles whether to trim holes +/// @param trimOutliers whether to trim outliers +/// @param trimMaterial whether to trim pure material states +template +void trimTrackBack(track_proxy_t track, bool trimHoles, bool trimOutliers, + bool trimMaterial) + requires(!track_proxy_t::ReadOnly) +{ + using TrackStateProxy = typename track_proxy_t::TrackStateProxy; + + std::optional back; + + for (TrackStateProxy trackState : track.trackStatesReversed()) { + back = trackState; + + TrackStateType typeFlags = trackState.typeFlags(); + if (trimHoles && typeFlags.test(TrackStateFlag::HoleFlag)) { + continue; + } + if (trimOutliers && typeFlags.test(TrackStateFlag::OutlierFlag)) { + continue; + } + if (trimMaterial && typeFlags.test(TrackStateFlag::MaterialFlag) && + !typeFlags.test(TrackStateFlag::MeasurementFlag)) { + continue; + } + + break; + } + + if (back.has_value()) { + track.tipIndex() = back.value().index(); + } +} + +/// Helper function to trim track states from the front and back of a track +/// @tparam track_proxy_t the track proxy type +/// @param track the track to trim +/// @param trimHoles whether to trim holes +/// @param trimOutliers whether to trim outliers +/// @param trimMaterial whether to trim pure material states +template +void trimTrack(track_proxy_t track, bool trimHoles, bool trimOutliers, + bool trimMaterial) + requires(!track_proxy_t::ReadOnly) +{ + trimTrackFront(track, trimHoles, trimOutliers, trimMaterial); + trimTrackBack(track, trimHoles, trimOutliers, trimMaterial); +} + } // namespace Acts namespace std { diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp index 056b2145270..74c13f9a5dd 100644 --- a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp @@ -126,6 +126,8 @@ class TrackFindingAlgorithm final : public IAlgorithm { bool stayOnSeed = false; /// Compute shared hit information bool computeSharedHits = false; + /// Whether to trim the tracks + bool trimTracks = true; // Pixel and strip volume ids to be used for maxPixel/StripHoles cuts std::set pixelVolumes; diff --git a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp index 7d0fbf5753a..81e850d297e 100644 --- a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp @@ -419,6 +419,10 @@ ProcessCode TrackFindingAlgorithm::execute(const AlgorithmContext& ctx) const { } }); + // trim the track if requested + if (m_cfg.trimTracks) { + Acts::trimTrack(track, true, true, true); + } Acts::calculateTrackQuantities(track); if (m_trackSelector.has_value() && !m_trackSelector->isValidTrack(track)) { diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index 11300d14d5d..e47b435e72b 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -146,8 +146,9 @@ "stripVolumes", "maxPixelHoles", "maxStripHoles", + "trimTracks", ], - defaults=[15.0, 25.0, 10, None, None, None, None, None, None, None], + defaults=[15.0, 25.0, 10, None, None, None, None, None, None, None, None], ) AmbiguityResolutionConfig = namedtuple( @@ -1534,6 +1535,7 @@ def addCKFTracks( stripVolumes=ckfConfig.stripVolumes, maxPixelHoles=ckfConfig.maxPixelHoles, maxStripHoles=ckfConfig.maxStripHoles, + trimTracks=ckfConfig.trimTracks, ), ) s.addAlgorithm(trackFinder) diff --git a/Examples/Python/src/TrackFinding.cpp b/Examples/Python/src/TrackFinding.cpp index bae0eb7f80a..1473ded385b 100644 --- a/Examples/Python/src/TrackFinding.cpp +++ b/Examples/Python/src/TrackFinding.cpp @@ -341,6 +341,7 @@ void addTrackFinding(Context& ctx) { ACTS_PYTHON_MEMBER(stripVolumes); ACTS_PYTHON_MEMBER(maxPixelHoles); ACTS_PYTHON_MEMBER(maxStripHoles); + ACTS_PYTHON_MEMBER(trimTracks); ACTS_PYTHON_STRUCT_END(); } diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt index 54ffad36b9d..7d86645f683 100644 --- a/Examples/Python/tests/root_file_hashes.txt +++ b/Examples/Python/tests/root_file_hashes.txt @@ -41,11 +41,11 @@ test_ckf_tracks_example[generic-truth_estimated]__tracksummary_ckf.root: 5a973ec test_ckf_tracks_example[generic-truth_estimated]__performance_seeding.root: 1facb05c066221f6361b61f015cdf0918e94d9f3fce2269ec7b6a4dffeb2bc7e test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: b10f61d3b68ecc3d1910a17aeadd01a1e23b31b6418935809d1f12e95eac607a test_ckf_tracks_example[generic-truth_smeared]__tracksummary_ckf.root: 69733ede1fc08370b5c0d0535f274b59bb51217239bd6645ff62ec9dadaa1f41 -test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: 6b9b97545e9c4c833547adcd54c78235165e7e184d52548ecb5f16c421f41b37 -test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: 5144c31f9f5281f1c287de978035c2584e9c7b1d46aece6afaecd5c1c86b0020 +test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: 33a68b20005a92b88bf7d01179a4b4fe3c1476a2019176904356924c0ad68f4b +test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: 749a53533fed3d550c2e37f504bac5a8ab5c5cdf60d5055dcd545aa44613601a test_ckf_tracks_example[odd-full_seeding]__performance_seeding_trees.root: 43c58577aafe07645e5660c4f43904efadf91d8cda45c5c04c248bbe0f59814f test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: 7750f58b970c79dc3c937482e2a6e4576bc5d495fd6e4576d63ac2882d8283d4 -test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: 3754817f9cf0972b92147e80bf36b2dc79a9c96e161b236a483a2104d2e33cd6 +test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: f5d5d5521e367dd2e26365b896b079279a5a8b97b7e11d38eb0eb317740ad4dd test_ckf_tracks_example[odd-truth_estimated]__performance_seeding.root: 1a36b7017e59f1c08602ef3c2cb0483c51df248f112e3780c66594110719c575 test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: 33398059bf968f7279d4cc706f3b914fcf6f010ae82c43df8875785bda6f7cbe test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: 3d01335a51fb03c78a174c8b11d473b38ba9c4ed3d6cba201b3925f463411708 diff --git a/Tests/UnitTests/Core/Utilities/TrackHelpersTests.cpp b/Tests/UnitTests/Core/Utilities/TrackHelpersTests.cpp index fcba710499c..ba6580b05f8 100644 --- a/Tests/UnitTests/Core/Utilities/TrackHelpersTests.cpp +++ b/Tests/UnitTests/Core/Utilities/TrackHelpersTests.cpp @@ -9,49 +9,101 @@ #include #include "Acts/EventData/TrackContainer.hpp" +#include "Acts/EventData/TrackStateType.hpp" #include "Acts/EventData/VectorMultiTrajectory.hpp" #include "Acts/EventData/VectorTrackContainer.hpp" #include "Acts/Utilities/TrackHelpers.hpp" +#include +#include + namespace Acts::Test { +namespace { + +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); + } + } + + return t; +} + +} // namespace + BOOST_AUTO_TEST_SUITE(Utilities) BOOST_AUTO_TEST_CASE(CalculateQuantities) { TrackContainer tc{VectorTrackContainer{}, VectorMultiTrajectory{}}; - auto t = tc.makeTrack(); + auto t = createTestTrack(tc, std::vector>{ + {MeasurementFlag}, + {OutlierFlag}, + {MeasurementFlag, SharedHitFlag}, + {HoleFlag}, + {OutlierFlag}, + {HoleFlag}, + {MeasurementFlag, SharedHitFlag}, + {OutlierFlag}, + }); + + calculateTrackQuantities(t); - auto ts = t.appendTrackState(); - ts.typeFlags().set(MeasurementFlag); + BOOST_CHECK_EQUAL(t.nHoles(), 2); + BOOST_CHECK_EQUAL(t.nMeasurements(), 3); + BOOST_CHECK_EQUAL(t.nOutliers(), 3); + BOOST_CHECK_EQUAL(t.nSharedHits(), 2); +} - ts = t.appendTrackState(); - ts.typeFlags().set(OutlierFlag); +BOOST_AUTO_TEST_CASE(TrimTrack) { + TrackContainer tc{VectorTrackContainer{}, VectorMultiTrajectory{}}; + auto t = createTestTrack(tc, std::vector>{ + {HoleFlag}, + {MeasurementFlag}, + {OutlierFlag}, + {MeasurementFlag, SharedHitFlag}, + {HoleFlag}, + {OutlierFlag}, + {HoleFlag}, + {MeasurementFlag}, + {OutlierFlag}, + }); - ts = t.appendTrackState(); - ts.typeFlags().set(MeasurementFlag); - ts.typeFlags().set(SharedHitFlag); + calculateTrackQuantities(t); - ts = t.appendTrackState(); - ts.typeFlags().set(HoleFlag); + BOOST_CHECK_EQUAL(t.nHoles(), 3); + BOOST_CHECK_EQUAL(t.nMeasurements(), 3); + BOOST_CHECK_EQUAL(t.nOutliers(), 3); + BOOST_CHECK_EQUAL(t.nSharedHits(), 1); - ts = t.appendTrackState(); - ts.typeFlags().set(OutlierFlag); + trimTrackFront(t, true, true, true); + calculateTrackQuantities(t); - ts = t.appendTrackState(); - ts.typeFlags().set(HoleFlag); + BOOST_CHECK_EQUAL(t.nHoles(), 2); + BOOST_CHECK_EQUAL(t.nMeasurements(), 3); + BOOST_CHECK_EQUAL(t.nOutliers(), 3); + BOOST_CHECK_EQUAL(t.nSharedHits(), 1); - ts = t.appendTrackState(); - ts.typeFlags().set(MeasurementFlag); - ts.typeFlags().set(SharedHitFlag); + trimTrackBack(t, true, true, true); + calculateTrackQuantities(t); - ts = t.appendTrackState(); - ts.typeFlags().set(OutlierFlag); + BOOST_CHECK_EQUAL(t.nHoles(), 2); + BOOST_CHECK_EQUAL(t.nMeasurements(), 3); + BOOST_CHECK_EQUAL(t.nOutliers(), 2); + BOOST_CHECK_EQUAL(t.nSharedHits(), 1); + trimTrack(t, true, true, true); calculateTrackQuantities(t); + BOOST_CHECK_EQUAL(t.nHoles(), 2); BOOST_CHECK_EQUAL(t.nMeasurements(), 3); - BOOST_CHECK_EQUAL(t.nOutliers(), 3); - BOOST_CHECK_EQUAL(t.nSharedHits(), 2); + BOOST_CHECK_EQUAL(t.nOutliers(), 2); + BOOST_CHECK_EQUAL(t.nSharedHits(), 1); } BOOST_AUTO_TEST_SUITE_END()