From cc014c584e13648c751847af15bc899cee6d227f Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Mon, 14 Oct 2024 19:08:00 +0200 Subject: [PATCH] refactor: Combine material, measurement and hole handling in Core CKF `filter` (#3723) Combines the handling of material, measurements and holes in the Core CKF. --- This pull request includes several changes to the `CombinatorialKalmanFilter` class in the `Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp` file to improve error handling, state management, and logging. Additionally, there is a minor change in the `MeasurementSelector` class to handle empty measurement candidates more gracefully. ### Error Handling and State Management: * Added checks and logging for empty active branches to ensure the filter stops correctly when no branches are active. (`Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp`) [[1]](diffhunk://#diff-55dd901875f5a87b6a9fee495de245296efbfd45aa9f9483036bcd299a277e09R583-R590) [[2]](diffhunk://#diff-55dd901875f5a87b6a9fee495de245296efbfd45aa9f9483036bcd299a277e09R660-R673) * Modified the filter to reset and stop properly when all branches are stopped or when certain conditions are met, improving the robustness of the filtering process. (`Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp`) [[1]](diffhunk://#diff-55dd901875f5a87b6a9fee495de245296efbfd45aa9f9483036bcd299a277e09L592-R600) [[2]](diffhunk://#diff-55dd901875f5a87b6a9fee495de245296efbfd45aa9f9483036bcd299a277e09L623-L642) ### Logging Improvements: * Enhanced logging to provide more detailed information about the state transitions and decisions made during the filtering process. (`Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp`) [[1]](diffhunk://#diff-55dd901875f5a87b6a9fee495de245296efbfd45aa9f9483036bcd299a277e09L716-R761) [[2]](diffhunk://#diff-55dd901875f5a87b6a9fee495de245296efbfd45aa9f9483036bcd299a277e09L749-L925) [[3]](diffhunk://#diff-55dd901875f5a87b6a9fee495de245296efbfd45aa9f9483036bcd299a277e09L950-R907) ### Code Simplification and Cleanup: * Simplified the handling of material and measurement surfaces by refactoring conditional checks and removing redundant code. (`Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp`) [[1]](diffhunk://#diff-55dd901875f5a87b6a9fee495de245296efbfd45aa9f9483036bcd299a277e09L716-R761) [[2]](diffhunk://#diff-55dd901875f5a87b6a9fee495de245296efbfd45aa9f9483036bcd299a277e09L749-L925) ### Minor Change in MeasurementSelector: * Updated the `MeasurementSelector::select` method to return a success result when no measurement candidates are found, instead of returning an error. (`Core/include/Acts/TrackFinding/MeasurementSelector.ipp`) --- .../CombinatorialKalmanFilter.hpp | 328 ++++++++---------- .../Acts/TrackFinding/MeasurementSelector.ipp | 17 +- 2 files changed, 149 insertions(+), 196 deletions(-) diff --git a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp index afb74140b8d6..05e236e1e290 100644 --- a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp +++ b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp @@ -34,6 +34,7 @@ #include #include #include +#include #include namespace Acts { @@ -579,8 +580,14 @@ class CombinatorialKalmanFilter { ACTS_ERROR("Error in filter: " << res.error()); result.lastError = res.error(); } + + if (result.finished) { + return; + } } + assert(!result.activeBranches.empty() && "No active branches"); + const bool isEndOfWorldReached = endOfWorldReached.checkAbort(state, stepper, navigator, logger()); const bool isVolumeConstraintReached = volumeConstraintAborter.checkAbort( @@ -589,9 +596,8 @@ class CombinatorialKalmanFilter { state, stepper, navigator, logger()); const bool isTargetReached = targetReached.checkAbort(state, stepper, navigator, logger()); - const bool allBranchesStopped = result.activeBranches.empty(); - if (isEndOfWorldReached || isPathLimitReached || isTargetReached || - allBranchesStopped || isVolumeConstraintReached) { + if (isEndOfWorldReached || isVolumeConstraintReached || + isPathLimitReached || isTargetReached) { if (isEndOfWorldReached) { ACTS_VERBOSE("End of world reached"); } else if (isVolumeConstraintReached) { @@ -620,26 +626,12 @@ class CombinatorialKalmanFilter { stepper.releaseStepSize(state.stepping, ConstrainedStep::actor); } - if (!allBranchesStopped) { - // Record the active branch and remove it from the list - storeLastActiveBranch(result); - result.activeBranches.pop_back(); - } else { - // This can happen if we stopped all branches in the filter step - ACTS_VERBOSE("All branches stopped"); - } + // Record the active branch and remove it from the list + storeLastActiveBranch(result); + result.activeBranches.pop_back(); - // If no more active branches, done with filtering; Otherwise, reset - // propagation state to track state at next active branch - if (!result.activeBranches.empty()) { - ACTS_VERBOSE("Propagation jumps to branch with tip = " - << result.activeBranches.back().tipIndex()); - reset(state, stepper, navigator, result); - } else { - ACTS_VERBOSE("Stop Kalman filtering with " - << result.collectedTracks.size() << " found tracks"); - result.finished = true; - } + // Reset propagation state to track state at next active branch + reset(state, stepper, navigator, result); } } @@ -665,9 +657,20 @@ class CombinatorialKalmanFilter { typename navigator_t> void reset(propagator_state_t& state, const stepper_t& stepper, const navigator_t& navigator, result_type& result) const { + if (result.activeBranches.empty()) { + ACTS_VERBOSE("Stop CKF with " << result.collectedTracks.size() + << " found tracks"); + result.finished = true; + + return; + } + auto currentBranch = result.activeBranches.back(); auto currentState = currentBranch.outermostTrackState(); + ACTS_VERBOSE("Propagation jumps to branch with tip = " + << currentBranch.tipIndex()); + // Reset the stepping state stepper.resetState(state.stepping, currentState.filtered(), currentState.filteredCovariance(), @@ -713,40 +716,67 @@ class CombinatorialKalmanFilter { using PM = TrackStatePropMask; bool isSensitive = surface->associatedDetectorElement() != nullptr; - bool isMaterial = surface->surfaceMaterial() != nullptr; - - std::size_t nBranchesOnSurface = 0; + bool hasMaterial = surface->surfaceMaterial() != nullptr; + bool isMaterialOnly = hasMaterial && !isSensitive; + bool expectMeasurements = isSensitive; - if (auto [slBegin, slEnd] = m_sourceLinkAccessor(*surface); - slBegin != slEnd) { - // Screen output message + if (isSensitive) { ACTS_VERBOSE("Measurement surface " << surface->geometryId() << " detected."); + } else if (isMaterialOnly) { + ACTS_VERBOSE("Material surface " << surface->geometryId() + << " detected."); + } else { + ACTS_VERBOSE("Passive surface " << surface->geometryId() + << " detected."); + return Result::success(); + } - // Transport the covariance to the surface + using SourceLinkRange = decltype(m_sourceLinkAccessor(*surface)); + std::optional slRange = std::nullopt; + bool hasMeasurements = false; + if (isSensitive) { + slRange = m_sourceLinkAccessor(*surface); + hasMeasurements = slRange->first != slRange->second; + } + bool isHole = isSensitive && !hasMeasurements; + + if (isHole) { + ACTS_VERBOSE("Detected hole before measurement selection on surface " + << surface->geometryId()); + } + + // Transport the covariance to the surface + if (isHole || isMaterialOnly) { + stepper.transportCovarianceToCurvilinear(state.stepping); + } else { stepper.transportCovarianceToBound(state.stepping, *surface); + } - // Update state and stepper with pre material effects - materialInteractor(surface, state, stepper, navigator, - MaterialUpdateStage::PreUpdate); + // Update state and stepper with pre material effects + materialInteractor(surface, state, stepper, navigator, + MaterialUpdateStage::PreUpdate); - // Bind the transported 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; + // Bind the transported 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; + + auto currentBranch = result.activeBranches.back(); + TrackIndexType prevTip = currentBranch.tipIndex(); - auto currentBranch = result.activeBranches.back(); - TrackIndexType prevTip = currentBranch.tipIndex(); + // Create trackstates for all source links (will be filtered later) + using TrackStatesResult = + Acts::Result>; + TrackStatesResult tsRes = TrackStatesResult::success({}); + if (hasMeasurements) { + auto [slBegin, slEnd] = *slRange; - // Create trackstates for all source links (will be filtered later) - using TrackStatesResult = - Acts::Result>; - TrackStatesResult tsRes = trackStateCandidateCreator( + tsRes = trackStateCandidateCreator( state.geoContext, *calibrationContextPtr, *surface, boundState, slBegin, slEnd, prevTip, *result.trackStates, result.trackStateCandidates, *result.trackStates, logger()); @@ -755,175 +785,98 @@ class CombinatorialKalmanFilter { "Processing of selected track states failed: " << tsRes.error()); return tsRes.error(); } - const CkfTypes::BranchVector& newTrackStateList = - *tsRes; - - if (newTrackStateList.empty()) { - ACTS_VERBOSE("Detected hole after measurement selection on surface " - << surface->geometryId()); - - // Setting the number of branches on the surface to 1 as the hole - // still counts as a branch - nBranchesOnSurface = 1; - - auto stateMask = PM::Predicted | PM::Jacobian; - - // Add a hole track state to the multitrajectory - TrackIndexType currentTip = addNonSourcelinkState( - stateMask, boundState, result, true, prevTip); - auto nonSourcelinkState = - result.trackStates->getTrackState(currentTip); - currentBranch.tipIndex() = currentTip; - currentBranch.nHoles()++; + } + const CkfTypes::BranchVector& newTrackStateList = *tsRes; + + if (!newTrackStateList.empty()) { + Result procRes = + processNewTrackStates(state.geoContext, newTrackStateList, result); + if (!procRes.ok()) { + ACTS_ERROR("Processing of selected track states failed: " + << procRes.error()); + return procRes.error(); + } + unsigned int nBranchesOnSurface = *procRes; - BranchStopperResult branchStopperResult = - m_extensions.branchStopper(currentBranch, nonSourcelinkState); + if (nBranchesOnSurface == 0) { + ACTS_VERBOSE("All branches on surface " << surface->geometryId() + << " have been stopped"); - // 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(); - } - } else { - Result procRes = processNewTrackStates( - state.geoContext, newTrackStateList, result); - if (!procRes.ok()) { - ACTS_ERROR("Processing of selected track states failed: " - << procRes.error()); - return procRes.error(); - } - nBranchesOnSurface = *procRes; + reset(state, stepper, navigator, result); - if (nBranchesOnSurface == 0) { - // All branches on the surface have been stopped. Reset happens at - // the end of the function - ACTS_VERBOSE("All branches on surface " << surface->geometryId() - << " have been stopped"); - } else { - // `currentBranch` is invalidated after `processNewTrackStates` - currentBranch = result.activeBranches.back(); - prevTip = currentBranch.tipIndex(); - - auto currentState = currentBranch.outermostTrackState(); - - if (currentState.typeFlags().test(TrackStateFlag::OutlierFlag)) { - // We don't need to update the stepper given an outlier state - ACTS_VERBOSE("Outlier state detected on surface " - << surface->geometryId()); - } else { - // If there are measurement track states on this surface - ACTS_VERBOSE("Filtering step successful with " - << nBranchesOnSurface << " branches"); - // Update stepping state using filtered parameters of last track - // state on this surface - stepper.update(state.stepping, - MultiTrajectoryHelpers::freeFiltered( - state.options.geoContext, currentState), - currentState.filtered(), - currentState.filteredCovariance(), *surface); - ACTS_VERBOSE( - "Stepping state is updated with filtered parameter:"); - ACTS_VERBOSE("-> " << currentState.filtered().transpose() - << " of track state with tip = " - << currentState.index()); - } - } + return Result::success(); } - // Update state and stepper with post material effects - materialInteractor(surface, state, stepper, navigator, - MaterialUpdateStage::PostUpdate); - } else if (isSensitive || isMaterial) { - ACTS_VERBOSE("Handle " << (isSensitive ? "sensitive" : "passive") - << " surface: " << surface->geometryId() - << " without measurements"); - - // No splitting on the surface without source links. Set it to one - // first, but could be changed later - nBranchesOnSurface = 1; - - auto currentBranch = result.activeBranches.back(); - TrackIndexType prevTip = currentBranch.tipIndex(); + // `currentBranch` is invalidated after `processNewTrackStates` + currentBranch = result.activeBranches.back(); + prevTip = currentBranch.tipIndex(); + } else { + if (expectMeasurements) { + ACTS_VERBOSE("Detected hole after measurement selection on surface " + << surface->geometryId()); + } - // 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); - - // Update state and stepper with pre material effects - materialInteractor(surface, state, stepper, navigator, - MaterialUpdateStage::PreUpdate); - - // 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 to the multitrajectory TrackIndexType currentTip = addNonSourcelinkState( - stateMask, boundState, result, isSensitive, prevTip); - auto nonSourcelinkState = result.trackStates->getTrackState(currentTip); + stateMask, boundState, result, expectMeasurements, prevTip); currentBranch.tipIndex() = currentTip; - - if (isSensitive) { + auto currentState = currentBranch.outermostTrackState(); + if (expectMeasurements) { currentBranch.nHoles()++; } BranchStopperResult branchStopperResult = - m_extensions.branchStopper(currentBranch, nonSourcelinkState); + m_extensions.branchStopper(currentBranch, currentState); // 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(); - } - // 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 - nBranchesOnSurface = 1; - } + // Branch on the surface has been stopped - reset + ACTS_VERBOSE("Branch on surface " << surface->geometryId() + << " has been stopped"); - // Reset current branch if there is no branch on current surface - if (nBranchesOnSurface == 0) { - ACTS_DEBUG("Branch on surface " << surface->geometryId() - << " is stopped"); - if (!result.activeBranches.empty()) { - ACTS_VERBOSE("Propagation jumps to branch with tip = " - << result.activeBranches.back().tipIndex()); reset(state, stepper, navigator, result); - } else { - ACTS_VERBOSE("Stop Kalman filtering with " - << result.collectedTracks.size() << " found tracks"); - result.finished = true; + + return Result::success(); } } + auto currentState = currentBranch.outermostTrackState(); + + if (currentState.typeFlags().test(TrackStateFlag::OutlierFlag)) { + // We don't need to update the stepper given an outlier state + ACTS_VERBOSE("Outlier state detected on surface " + << surface->geometryId()); + } else if (currentState.typeFlags().test( + TrackStateFlag::MeasurementFlag)) { + // If there are measurement track states on this surface + // Update stepping state using filtered parameters of last track + // state on this surface + stepper.update(state.stepping, + MultiTrajectoryHelpers::freeFiltered( + state.options.geoContext, currentState), + currentState.filtered(), + currentState.filteredCovariance(), *surface); + ACTS_VERBOSE("Stepping state is updated with filtered parameter:"); + ACTS_VERBOSE("-> " << currentState.filtered().transpose() + << " of track state with tip = " + << currentState.index()); + } + + // Update state and stepper with post material effects + materialInteractor(surface, state, stepper, navigator, + MaterialUpdateStage::PostUpdate); + return Result::success(); } @@ -947,12 +900,13 @@ class CombinatorialKalmanFilter { auto rootBranch = result.activeBranches.back(); - // Build the new branches by forking the root branch. Reverse the order to - // process the best candidate first + // Build the new branches by forking the root branch. Reverse the order + // to process the best candidate first CkfTypes::BranchVector newBranches; for (auto it = newTrackStateList.rbegin(); it != newTrackStateList.rend(); ++it) { - // Keep the root branch as the first branch, make a copy for the others + // Keep the root branch as the first branch, make a copy for the + // others auto newBranch = (it == newTrackStateList.rbegin()) ? rootBranch : rootBranch.shallowCopy(); diff --git a/Core/include/Acts/TrackFinding/MeasurementSelector.ipp b/Core/include/Acts/TrackFinding/MeasurementSelector.ipp index 3e32c8845851..30d09c002be1 100644 --- a/Core/include/Acts/TrackFinding/MeasurementSelector.ipp +++ b/Core/include/Acts/TrackFinding/MeasurementSelector.ipp @@ -21,9 +21,9 @@ MeasurementSelector::select( ACTS_VERBOSE("Invoked MeasurementSelector"); - // Return error if no measurement + // Return if no measurement if (candidates.empty()) { - return CombinatorialKalmanFilterError::MeasurementSelectionFailed; + return Result::success(std::pair(candidates.begin(), candidates.end())); } // Get geoID of this surface @@ -92,20 +92,19 @@ MeasurementSelector::select( "No measurement candidate. Return an outlier measurement chi2=" << minChi2); isOutlier = true; - return Result::success(std::make_pair(candidates.begin() + minIndex, - candidates.begin() + minIndex + 1)); + return Result::success(std::pair(candidates.begin() + minIndex, + candidates.begin() + minIndex + 1)); } else { ACTS_VERBOSE("No measurement candidate. Return empty chi2=" << minChi2); - return Result::success( - std::make_pair(candidates.begin(), candidates.begin())); + return Result::success(std::pair(candidates.begin(), candidates.begin())); } } if (passedCandidates <= 1ul) { // return single item range, no sorting necessary ACTS_VERBOSE("Returning only 1 element chi2=" << minChi2); - return Result::success(std::make_pair(candidates.begin() + minIndex, - candidates.begin() + minIndex + 1)); + return Result::success(std::pair(candidates.begin() + minIndex, + candidates.begin() + minIndex + 1)); } std::sort( @@ -115,7 +114,7 @@ MeasurementSelector::select( ACTS_VERBOSE("Number of selected measurements: " << passedCandidates << ", max: " << cuts.numMeasurements); - return Result::success(std::make_pair( + return Result::success(std::pair( candidates.begin(), candidates.begin() + std::min(cuts.numMeasurements, passedCandidates))); }