diff --git a/Core/include/Acts/Fitter/GainMatrixUpdater.hpp b/Core/include/Acts/Fitter/GainMatrixUpdater.hpp
index 5001fe3213980959694f6d4302ec76f3cdc565ee..e31adc5fc8a361e9bf5ad106a9e5d6d10a238aed 100644
--- a/Core/include/Acts/Fitter/GainMatrixUpdater.hpp
+++ b/Core/include/Acts/Fitter/GainMatrixUpdater.hpp
@@ -24,7 +24,7 @@ namespace Acts {
/// @brief Update step of Kalman Filter using gain matrix formalism
///
-/// @tparam parameters_t Type of the parameters to be updated
+/// @tparam parameters_t Type of the parameters to be filtered
/// @tparam jacobian_t Type of the Transport jacobian
///
template <typename parameters_t>
@@ -46,13 +46,15 @@ class GainMatrixUpdater {
///
/// @param gctx The current geometry context object, e.g. alignment
/// @param trackState the measured track state
+ /// @param direction the navigation direction
///
/// @return Bool indicating whether this update was 'successful'
/// @note Non-'successful' updates could be holes or outliers,
/// which need to be treated differently in calling code.
template <typename track_state_t>
- Result<void> operator()(const GeometryContext& /*gctx*/,
- track_state_t trackState) const {
+ Result<void> operator()(
+ const GeometryContext& /*gctx*/, track_state_t trackState,
+ const NavigationDirection& direction = forward) const {
ACTS_VERBOSE("Invoked GainMatrixUpdater");
// let's make sure the types are consistent
using SourceLink = typename track_state_t::SourceLink;
@@ -61,6 +63,9 @@ class GainMatrixUpdater {
static_assert(std::is_same_v<track_state_t, TrackStateProxy>,
"Given track state type is not a track state proxy");
+ using CovMatrix_t = typename parameters_t::CovMatrix_t;
+ using ParVector_t = typename parameters_t::ParVector_t;
+
// we should definitely have an uncalibrated measurement here
assert(trackState.hasUncalibrated());
// there should be a calibrated measurement
@@ -110,8 +115,11 @@ class GainMatrixUpdater {
ACTS_VERBOSE("Gain Matrix K:\n" << K);
if (K.hasNaN()) {
- error = KalmanFitterError::UpdateFailed; // set to error
- return false; // abort execution
+ error =
+ (direction == forward)
+ ? KalmanFitterError::ForwardUpdateFailed
+ : KalmanFitterError::BackwardUpdateFailed; // set to error
+ return false; // abort execution
}
filtered = predicted + K * (calibrated - H * predicted);
@@ -138,6 +146,11 @@ class GainMatrixUpdater {
return true; // continue execution
});
+ if (error) {
+ // error is set, return result
+ return *error;
+ }
+
// always succeed, no outlier logic yet
return Result<void>::success();
}
diff --git a/Core/include/Acts/Fitter/KalmanFitter.hpp b/Core/include/Acts/Fitter/KalmanFitter.hpp
index 832bbe63bc3b51e400e18e16bcaded12605f468b..fb89afbacdf608df4aab20a83c564a04a68c009a 100644
--- a/Core/include/Acts/Fitter/KalmanFitter.hpp
+++ b/Core/include/Acts/Fitter/KalmanFitter.hpp
@@ -23,6 +23,7 @@
#include "Acts/Propagator/ActionList.hpp"
#include "Acts/Propagator/ConstrainedStep.hpp"
#include "Acts/Propagator/DirectNavigator.hpp"
+#include "Acts/Propagator/Navigator.hpp"
#include "Acts/Propagator/Propagator.hpp"
#include "Acts/Propagator/detail/PointwiseMaterialInteraction.hpp"
#include "Acts/Propagator/detail/StandardAborters.hpp"
@@ -39,8 +40,9 @@ namespace Acts {
/// @brief Options struct how the Fitter is called
///
-/// It contains the context of the fitter call and the optional
-/// surface where to express the fit result
+/// It contains the context of the fitter call, the optional
+/// surface where to express the fit result and configurations for material
+/// effects and smoothing options
///
/// @note the context objects must be provided
struct KalmanFitterOptions {
@@ -53,17 +55,22 @@ struct KalmanFitterOptions {
/// @param mctx The magnetic context for this fit
/// @param cctx The calibration context for this fit
/// @param rSurface The reference surface for the fit to be expressed at
+ /// @param mScattering Whether to include multiple scattering
+ /// @param eLoss Whether to include energy loss
+ /// @param bwdFiltering Whether to run backward filtering as smoothing
KalmanFitterOptions(std::reference_wrapper<const GeometryContext> gctx,
std::reference_wrapper<const MagneticFieldContext> mctx,
std::reference_wrapper<const CalibrationContext> cctx,
const Surface* rSurface = nullptr,
- bool mScattering = true, bool eLoss = true)
+ bool mScattering = true, bool eLoss = true,
+ bool bwdFiltering = false)
: geoContext(gctx),
magFieldContext(mctx),
calibrationContext(cctx),
referenceSurface(rSurface),
multipleScattering(mScattering),
- energyLoss(eLoss) {}
+ energyLoss(eLoss),
+ backwardFiltering(bwdFiltering) {}
/// Context object for the geometry
std::reference_wrapper<const GeometryContext> geoContext;
@@ -75,11 +82,14 @@ struct KalmanFitterOptions {
/// The reference Surface
const Surface* referenceSurface = nullptr;
- /// Whether to consider multiple scattering.
+ /// Whether to consider multiple scattering
bool multipleScattering = true;
- /// Whether to consider energy loss.
+ /// Whether to consider energy loss
bool energyLoss = true;
+
+ /// Whether to run backward filtering.
+ bool backwardFiltering = false;
};
template <typename source_link_t>
@@ -108,7 +118,13 @@ struct KalmanFitterResult {
bool initialized = false;
// Measurement surfaces without hits
- std::vector<const Surface*> missedActiveSurfaces = {};
+ std::vector<const Surface*> missedActiveSurfaces;
+
+ // Indicator if forward filtering has been done
+ bool forwardFiltered = false;
+
+ // Measurement surfaces handled in both forward and backward filtering
+ std::vector<const Surface*> passedAgainSurfaces;
Result<void> result{Result<void>::success()};
};
@@ -230,6 +246,9 @@ class KalmanFitter {
/// Whether to consider energy loss.
bool energyLoss = true;
+ /// Whether run smoothing as backward filtering
+ bool backwardFiltering = false;
+
/// @brief Kalman actor operation
///
/// @tparam propagator_state_t is the type of Propagagor state
@@ -242,6 +261,7 @@ class KalmanFitter {
void operator()(propagator_state_t& state, const stepper_t& stepper,
result_type& result) const {
ACTS_VERBOSE("KalmanFitter step");
+
// Initialization:
// - Only when track states are not set
if (!result.initialized) {
@@ -256,41 +276,66 @@ class KalmanFitter {
// - Waiting for a current surface that has material
// -> a trackState will be created on surface with material
auto surface = state.navigation.currentSurface;
- if ((surface and surface->surfaceMaterial()) and not result.smoothed) {
+ if (surface and surface->surfaceMaterial()) {
// Check if the surface is in the measurement map
// -> Get the measurement / calibrate
// -> Create the predicted state
// -> Perform the kalman update
// -> Check outlier behavior (@todo)
// -> Fill strack state information & update stepper information
- ACTS_VERBOSE("Perform filter step");
- auto res = filter(surface, state, stepper, result);
- if (!res.ok()) {
- ACTS_ERROR("Error in filter: " << res.error());
- result.result = res.error();
+ if (state.stepping.navDir == forward and not result.smoothed and
+ not result.forwardFiltered) {
+ ACTS_VERBOSE("Perform forward filter step");
+ auto res = filter(surface, state, stepper, result);
+ if (!res.ok()) {
+ ACTS_ERROR("Error in forward filter: " << res.error());
+ result.result = res.error();
+ }
+ } else if (state.stepping.navDir == backward) {
+ ACTS_VERBOSE("Perform backward filter step");
+ auto res = backwardFilter(surface, state, stepper, result);
+ if (!res.ok()) {
+ ACTS_ERROR("Error in backward filter: " << res.error());
+ result.result = res.error();
+ }
}
}
// Finalization:
- // When all track states have been handled or the navigation is breaked
- if ((result.measurementStates == inputMeasurements.size() or
- (result.measurementStates > 0 and
- state.navigation.navigationBreak)) and
- not result.smoothed) {
- // -> Sort the track states (as now the path length is set)
- // -> Call the smoothing
- // -> Set a stop condition when all track states have been handled
- ACTS_VERBOSE("Finalize/run smoothing");
- auto res = finalize(state, stepper, result);
- if (!res.ok()) {
- ACTS_ERROR("Error in finalize: " << res.error());
- result.result = res.error();
+ // when all track states have been handled or the navigation is breaked,
+ // reset navigation&stepping before run backward filtering or
+ // proceed to run smoothing
+ if (state.stepping.navDir == forward) {
+ if (result.measurementStates == inputMeasurements.size() or
+ (result.measurementStates > 0 and
+ state.navigation.navigationBreak)) {
+ if (backwardFiltering and not result.forwardFiltered) {
+ ACTS_VERBOSE("Forward filtering done");
+ result.forwardFiltered = true;
+ // Run backward filtering
+ // Reverse navigation direction and reset navigation and stepping
+ // state to last measurement
+ ACTS_VERBOSE("Reverse navigation direction.");
+ reverse(state, stepper, result);
+ } else if (not result.smoothed) {
+ // -> Sort the track states (as now the path length is set)
+ // -> Call the smoothing
+ // -> Set a stop condition when all track states have been handled
+ ACTS_VERBOSE("Finalize/run smoothing");
+ auto res = finalize(state, stepper, result);
+ if (!res.ok()) {
+ ACTS_ERROR("Error in finalize: " << res.error());
+ result.result = res.error();
+ }
+ }
}
}
+
// Post-finalization:
// - Progress to target/reference surface and built the final track
// parameters
- if (result.smoothed and targetReached(state, stepper, *targetSurface)) {
+ if ((result.smoothed or state.stepping.navDir == backward) and
+ targetReached(state, stepper, *targetSurface)) {
ACTS_VERBOSE("Completing");
// Transport & bind the parameter to the final surface
auto fittedState =
@@ -299,6 +344,22 @@ class KalmanFitter {
result.fittedParameters = std::get<BoundParameters>(fittedState);
// Break the navigation for stopping the Propagation
state.navigation.navigationBreak = true;
+
+ // Reset smoothed status of states missed in backward filtering
+ if (backwardFiltering) {
+ result.fittedStates.applyBackwards(result.trackTip, [&](auto state) {
+ auto fSurface = &state.referenceSurface();
+ auto surface_it = std::find_if(
+ result.passedAgainSurfaces.begin(),
+ result.passedAgainSurfaces.end(),
+ [=](const Surface* surface) { return surface == fSurface; });
+ if (surface_it == result.passedAgainSurfaces.end()) {
+ // If backward filtering missed this surface, then there is no
+ // smoothed parameter
+ state.data().ismoothed = detail_lt::IndexData::kInvalid;
+ }
+ });
+ }
}
}
@@ -314,6 +375,69 @@ class KalmanFitter {
void initialize(propagator_state_t& /*state*/, const stepper_t& /*stepper*/,
result_type& /*result*/) const {}
+ /// @brief Kalman actor operation : reverse direction
+ ///
+ /// @tparam propagator_state_t is the type of Propagagor state
+ /// @tparam stepper_t Type of the stepper
+ ///
+ /// @param state is the mutable propagator state object
+ /// @param stepper The stepper in use
+ /// @param result is the mutable result state object
+ template <typename propagator_state_t, typename stepper_t>
+ void reverse(propagator_state_t& state, stepper_t& stepper,
+ result_type& result) const {
+ // Reset propagator options
+ state.options.direction = backward;
+ state.options.maxStepSize = -1.0 * state.options.maxStepSize;
+ // Not sure if reset of pathLimit during propagation makes any sense
+ state.options.pathLimit = -1.0 * state.options.pathLimit;
+
+ // Reset stepping&navigation state using last measurement track state on
+ // sensitive surface
+ state.navigation = typename propagator_t::NavigatorState();
+ result.fittedStates.applyBackwards(result.trackTip, [&](auto st) {
+ if (st.hasUncalibrated()) {
+ // Set the navigation state
+ state.navigation.startSurface = &st.referenceSurface();
+ state.navigation.startLayer =
+ state.navigation.startSurface->associatedLayer();
+ state.navigation.startVolume =
+ state.navigation.startLayer->trackingVolume();
+ state.navigation.targetSurface = targetSurface;
+ state.navigation.currentSurface = state.navigation.startSurface;
+ state.navigation.currentVolume = state.navigation.startVolume;
+
+ // Update the stepping state
+ stepper.update(state.stepping,
+ st.filteredParameters(state.options.geoContext));
+ // Reverse stepping direction
+ state.stepping.navDir = backward;
+ state.stepping.stepSize = ConstrainedStep(state.options.maxStepSize);
+ state.stepping.pathAccumulated = 0.;
+ // Reinitialize the stepping jacobian
+ st.referenceSurface().initJacobianToGlobal(
+ state.options.geoContext, state.stepping.jacToGlobal,
+ state.stepping.pos, state.stepping.dir,
+ st.filteredParameters(state.options.geoContext).parameters());
+ state.stepping.jacobian = BoundMatrix::Identity();
+ state.stepping.jacTransport = FreeMatrix::Identity();
+ state.stepping.derivative = FreeVector::Zero();
+
+ // For the last measurement state, smoothed is filtered
+ st.smoothed() = st.filtered();
+ st.smoothedCovariance() = st.filteredCovariance();
+ result.passedAgainSurfaces.push_back(&st.referenceSurface());
+
+ // Update material effects for last measurement state in backward
+ // direction
+ materialInteractor(state.navigation.currentSurface, state, stepper);
+
+ return false; // abort execution
+ }
+ return true; // continue execution
+ });
+ }
+
/// @brief Kalman actor operation : update
///
/// @tparam propagator_state_t is the type of Propagagor state
@@ -374,7 +498,7 @@ class KalmanFitter {
typeFlags.set(TrackStateFlag::ParameterFlag);
// If the update is successful, set covariance and
- auto updateRes = m_updater(state.geoContext, trackStateProxy);
+ auto updateRes = m_updater(state.geoContext, trackStateProxy, forward);
if (!updateRes.ok()) {
ACTS_ERROR("Update step failed: " << updateRes.error());
return updateRes.error();
@@ -396,71 +520,190 @@ class KalmanFitter {
// We count the processed state
++result.processedStates;
} else {
- // add a non-measurement TrackState entry multi trajectory
- // (this allocates storage for components except measurements, we will
- // set them later)
- result.trackTip = result.fittedStates.addTrackState(
- ~(TrackStatePropMask::Uncalibrated |
- TrackStatePropMask::Calibrated),
- result.trackTip);
+ if (result.measurementStates > 0) {
+ // No source links on surface, add either hole or passive material
+ // TrackState entry multi trajectory. No storage allocation for
+ // uncalibrated/calibrated measurement and filtered parameter
+ result.trackTip = result.fittedStates.addTrackState(
+ ~(TrackStatePropMask::Uncalibrated |
+ TrackStatePropMask::Calibrated | TrackStatePropMask::Filtered),
+ result.trackTip);
+
+ // now get track state proxy back
+ auto trackStateProxy =
+ result.fittedStates.getTrackState(result.trackTip);
+
+ // Set the surface
+ trackStateProxy.setReferenceSurface(surface->getSharedPtr());
+
+ // Set the track state flags
+ auto& typeFlags = trackStateProxy.typeFlags();
+ typeFlags.set(TrackStateFlag::MaterialFlag);
+ typeFlags.set(TrackStateFlag::ParameterFlag);
+
+ if (surface->associatedDetectorElement() != nullptr) {
+ ACTS_VERBOSE("Detected hole on " << surface->geoID());
+ // If the surface is sensitive, set the hole type flag
+ typeFlags.set(TrackStateFlag::HoleFlag);
+
+ // Count the missed surface
+ result.missedActiveSurfaces.push_back(surface);
+
+ // Transport & bind the state to the current surface
+ auto [boundParams, jacobian, pathLength] =
+ stepper.boundState(state.stepping, *surface, true);
+
+ // Fill the track state
+ trackStateProxy.predicted() = boundParams.parameters();
+ trackStateProxy.predictedCovariance() = *boundParams.covariance();
+ trackStateProxy.jacobian() = jacobian;
+ trackStateProxy.pathLength() = pathLength;
+ } else {
+ ACTS_VERBOSE("Detected in-sensitive surface " << surface->geoID());
+
+ // Transport & get curvilinear state instead of bound state
+ auto [curvilinearParams, jacobian, pathLength] =
+ stepper.curvilinearState(state.stepping, true);
+
+ // Fill the track state
+ trackStateProxy.predicted() = curvilinearParams.parameters();
+ trackStateProxy.predictedCovariance() =
+ *curvilinearParams.covariance();
+ trackStateProxy.jacobian() = jacobian;
+ trackStateProxy.pathLength() = pathLength;
+ }
+
+ // Set the filtered parameter index to be the same with predicted
+ // parameter
+ trackStateProxy.data().ifiltered = trackStateProxy.data().ipredicted;
+
+ // We count the processed state
+ ++result.processedStates;
+ }
- // now get track state proxy back
- auto trackStateProxy =
- result.fittedStates.getTrackState(result.trackTip);
+ // Update state and stepper with material effects
+ materialInteractor(surface, state, stepper, fullUpdate);
+ }
+ return Result<void>::success();
+ }
+
+ /// @brief Kalman actor operation : backward update
+ ///
+ /// @tparam propagator_state_t is the type of Propagagor state
+ /// @tparam stepper_t Type of the stepper
+ ///
+ /// @param surface The surface where the update happens
+ /// @param state The mutable propagator state object
+ /// @param stepper The stepper in use
+ /// @param result The mutable result state object
+ template <typename propagator_state_t, typename stepper_t>
+ Result<void> backwardFilter(const Surface* surface,
+ propagator_state_t& state,
+ const stepper_t& stepper,
+ result_type& result) const {
+ // Try to find the surface in the measurement surfaces
+ auto sourcelink_it = inputMeasurements.find(surface);
+ if (sourcelink_it != inputMeasurements.end()) {
+ // Screen output message
+ ACTS_VERBOSE("Measurement surface "
+ << surface->geoID()
+ << " detected in backward propagation.");
+
+ // No backward filtering for last measurement state, but still update
+ // with material effects
+ if (state.stepping.navDir == backward and
+ surface == state.navigation.startSurface) {
+ materialInteractor(surface, state, stepper);
+ return Result<void>::success();
+ }
- // Set the surface
- trackStateProxy.setReferenceSurface(surface->getSharedPtr());
+ // Update state and stepper with pre material effects
+ materialInteractor(surface, state, stepper, preUpdate);
- // Set the track state flags
- auto& typeFlags = trackStateProxy.typeFlags();
- typeFlags.set(TrackStateFlag::MaterialFlag);
- typeFlags.set(TrackStateFlag::ParameterFlag);
+ // Transport & bind the state to the current surface
+ auto [boundParams, jacobian, pathLength] =
+ stepper.boundState(state.stepping, *surface, true);
+
+ // Create a detached track state proxy
+ auto tempTrackTip =
+ result.fittedStates.addTrackState(TrackStatePropMask::All);
+
+ // Get the detached track state proxy back
+ auto trackStateProxy = result.fittedStates.getTrackState(tempTrackTip);
+
+ // Assign the source link to the detached track state
+ trackStateProxy.uncalibrated() = sourcelink_it->second;
+
+ // Fill the track state
+ trackStateProxy.predicted() = boundParams.parameters();
+ trackStateProxy.predictedCovariance() = *boundParams.covariance();
+ trackStateProxy.jacobian() = jacobian;
+ trackStateProxy.pathLength() = pathLength;
+
+ // We have predicted parameters, so calibrate the uncalibrated input
+ // measuerement
+ std::visit(
+ [&](const auto& calibrated) {
+ trackStateProxy.setCalibrated(calibrated);
+ },
+ m_calibrator(trackStateProxy.uncalibrated(),
+ trackStateProxy.predicted()));
+
+ // If the update is successful, set covariance and
+ auto updateRes = m_updater(state.geoContext, trackStateProxy, backward);
+ if (!updateRes.ok()) {
+ ACTS_ERROR("Backward update step failed: " << updateRes.error());
+ return updateRes.error();
+ } else {
+ // Update the stepping state with filtered parameters
+ ACTS_VERBOSE(
+ "Backward Filtering step successful, updated parameters are : \n"
+ << trackStateProxy.filtered().transpose());
+
+ // Fill the smoothed parameter for the existing track state
+ result.fittedStates.applyBackwards(result.trackTip, [&](auto state) {
+ auto fSurface = &state.referenceSurface();
+ if (fSurface == surface) {
+ result.passedAgainSurfaces.push_back(surface);
+ state.smoothed() = trackStateProxy.filtered();
+ state.smoothedCovariance() = trackStateProxy.filteredCovariance();
+ return false;
+ }
+ return true;
+ });
+ // update stepping state using filtered parameters after kalman update
+ // We need to (re-)construct a BoundParameters instance here, which is
+ // a bit awkward.
+ stepper.update(state.stepping, trackStateProxy.filteredParameters(
+ state.options.geoContext));
+
+ // Update state and stepper with post material effects
+ materialInteractor(surface, state, stepper, postUpdate);
+ }
+ } else {
+ // Transport covariance
if (surface->associatedDetectorElement() != nullptr) {
- ACTS_VERBOSE("Detected hole on " << surface->geoID());
- // If the surface is sensitive, set the hole type flag
- typeFlags.set(TrackStateFlag::HoleFlag);
-
- // Count the missed surface
- result.missedActiveSurfaces.push_back(surface);
-
- // Transport & bind the state to the current surface
- auto [boundParams, jacobian, pathLength] =
- stepper.boundState(state.stepping, *surface, true);
-
- // Fill the track state
- trackStateProxy.predicted() = boundParams.parameters();
- trackStateProxy.predictedCovariance() = *boundParams.covariance();
- trackStateProxy.jacobian() = jacobian;
- trackStateProxy.pathLength() = pathLength;
+ ACTS_VERBOSE("Detected hole on " << surface->geoID()
+ << " in backward filtering");
+ if (state.stepping.covTransport) {
+ stepper.covarianceTransport(state.stepping, *surface, true);
+ }
} else {
- ACTS_VERBOSE("Detected in-sensitive surface " << surface->geoID());
-
- // Transport & get curvilinear state instead of bound state
- auto [curvilinearParams, jacobian, pathLength] =
- stepper.curvilinearState(state.stepping, true);
-
- // Fill the track state
- trackStateProxy.predicted() = curvilinearParams.parameters();
- trackStateProxy.predictedCovariance() =
- *curvilinearParams.covariance();
- trackStateProxy.jacobian() = jacobian;
- trackStateProxy.pathLength() = pathLength;
+ ACTS_VERBOSE("Detected in-sensitive surface "
+ << surface->geoID() << " in backward filtering");
+ if (state.stepping.covTransport) {
+ stepper.covarianceTransport(state.stepping, true);
+ }
}
- // Update state and stepper with material effects
- materialInteractor(surface, state, stepper, fullUpdate);
-
- // Set the filtered parameter to be the same with predicted parameter
- // @Todo: shall we update the filterd parameter with material effects?
- // But it seems that the smoothing does not like this
- trackStateProxy.filtered() = trackStateProxy.predicted();
- trackStateProxy.filteredCovariance() =
- trackStateProxy.predictedCovariance();
+ // Not creating bound state here, so need manually reinitialize jacobian
+ state.stepping.jacobian = BoundMatrix::Identity();
- // We count the processed state
- ++result.processedStates;
+ // Update state and stepper with material effects
+ materialInteractor(surface, state, stepper);
}
+
return Result<void>::success();
}
@@ -475,9 +718,9 @@ class KalmanFitter {
/// @param updateStage The materal update stage
///
template <typename propagator_state_t, typename stepper_t>
- void materialInteractor(const Surface* surface, propagator_state_t& state,
- stepper_t& stepper,
- const MaterialUpdateStage& updateStage) const {
+ void materialInteractor(
+ const Surface* surface, propagator_state_t& state, stepper_t& stepper,
+ const MaterialUpdateStage& updateStage = fullUpdate) const {
// Prepare relevant input particle properties
detail::PointwiseMaterialInteraction interaction(surface, state, stepper);
@@ -519,23 +762,55 @@ class KalmanFitter {
// Remember you smoothed the track states
result.smoothed = true;
+ // Get the index of measurement states;
+ std::vector<size_t> measurementIndices;
+ auto lastState = result.fittedStates.getTrackState(result.trackTip);
+ if (lastState.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) {
+ measurementIndices.push_back(result.trackTip);
+ }
+ // Count track states to be smoothed
+ size_t nStates = 0;
+ result.fittedStates.applyBackwards(result.trackTip, [&](auto st) {
+ // Smoothing will start from the last measurement state
+ if (measurementIndices.empty()) {
+ // No smoothed parameter for the last few non-measurment states
+ st.data().ismoothed = detail_lt::IndexData::kInvalid;
+ } else {
+ nStates++;
+ }
+ size_t iprevious = st.previous();
+ if (iprevious != Acts::detail_lt::IndexData::kInvalid) {
+ auto previousState = result.fittedStates.getTrackState(iprevious);
+ if (previousState.typeFlags().test(
+ Acts::TrackStateFlag::MeasurementFlag)) {
+ measurementIndices.push_back(iprevious);
+ }
+ }
+ });
+ // Return error if the track has no measurement states (but this should
+ // not happen)
+ if (measurementIndices.empty()) {
+ return KalmanFitterError::SmoothFailed;
+ }
// Screen output for debugging
if (logger().doPrint(Logging::VERBOSE)) {
- // need to count track states
- size_t nStates = 0;
- result.fittedStates.visitBackwards(result.trackTip,
- [&](const auto) { nStates++; });
ACTS_VERBOSE("Apply smoothing on " << nStates
<< " filtered track states.");
}
- // Smooth the track states and obtain the last smoothed track parameters
- auto smoothRes =
- m_smoother(state.geoContext, result.fittedStates, result.trackTip);
+ // Smooth the track states
+ auto smoothRes = m_smoother(state.geoContext, result.fittedStates,
+ measurementIndices.front());
+
if (!smoothRes.ok()) {
ACTS_ERROR("Smoothing step failed: " << smoothRes.error());
return smoothRes.error();
}
- parameters_t smoothedPars = *smoothRes;
+ // Obtain the smoothed parameters at first measurement state
+ auto firstMeasurement =
+ result.fittedStates.getTrackState(measurementIndices.back());
+ parameters_t smoothedPars =
+ firstMeasurement.smoothedParameters(state.options.geoContext);
+
// Update the stepping parameters - in order to progress to destination
ACTS_VERBOSE(
"Smoothing successful, updating stepping state, "
@@ -544,6 +819,10 @@ class KalmanFitter {
// Reverse the propagation direction
state.stepping.stepSize =
ConstrainedStep(-1. * state.options.maxStepSize);
+ state.stepping.navDir = backward;
+ // Set accumulatd path to zero before targeting surface
+ state.stepping.pathAccumulated = 0.;
+ // Not sure if the following line helps anything
state.options.direction = backward;
return Result<void>::success();
@@ -639,6 +918,7 @@ class KalmanFitter {
kalmanActor.targetSurface = kfOptions.referenceSurface;
kalmanActor.multipleScattering = kfOptions.multipleScattering;
kalmanActor.energyLoss = kfOptions.energyLoss;
+ kalmanActor.backwardFiltering = kfOptions.backwardFiltering;
// also set logger on updater and smoother
kalmanActor.m_updater.m_logger = m_logger;
@@ -658,6 +938,7 @@ class KalmanFitter {
/// It could happen that the fit ends in zero processed states.
/// The result gets meaningless so such case is regarded as fit failure.
+ //@TODO: should we require the number of measurments >0 ?
if (kalmanResult.result.ok() and not kalmanResult.processedStates) {
kalmanResult.result = Result<void>(KalmanFitterError::PropagationInVain);
}
@@ -723,6 +1004,9 @@ class KalmanFitter {
kalmanActor.m_logger = m_logger.get();
kalmanActor.inputMeasurements = std::move(inputMeasurements);
kalmanActor.targetSurface = kfOptions.referenceSurface;
+ kalmanActor.multipleScattering = kfOptions.multipleScattering;
+ kalmanActor.energyLoss = kfOptions.energyLoss;
+ kalmanActor.backwardFiltering = kfOptions.backwardFiltering;
// also set logger on updater and smoother
kalmanActor.m_updater.m_logger = m_logger;
diff --git a/Core/include/Acts/Fitter/KalmanFitterError.hpp b/Core/include/Acts/Fitter/KalmanFitterError.hpp
index 56be3facb94a0e65c058f797182c7d0d7c30f1c4..78d599ff090eb8c4cfeb7d631c12e6a31cd6dddd 100644
--- a/Core/include/Acts/Fitter/KalmanFitterError.hpp
+++ b/Core/include/Acts/Fitter/KalmanFitterError.hpp
@@ -15,10 +15,11 @@
namespace Acts {
// This is the custom error code enum
enum class KalmanFitterError {
- UpdateFailed = 1,
- SmoothFailed = 2,
- OutputConversionFailed = 3,
- PropagationInVain = 4
+ ForwardUpdateFailed = 1,
+ BackwardUpdateFailed = 2,
+ SmoothFailed = 3,
+ OutputConversionFailed = 4,
+ PropagationInVain = 5
};
namespace detail {
@@ -30,8 +31,10 @@ class KalmanFitterErrorCategory : public std::error_category {
// Return what each enum means in text
std::string message(int c) const final {
switch (static_cast<KalmanFitterError>(c)) {
- case KalmanFitterError::UpdateFailed:
- return "Kalman update failed";
+ case KalmanFitterError::ForwardUpdateFailed:
+ return "Kalman forward update failed";
+ case KalmanFitterError::BackwardUpdateFailed:
+ return "Kalman backward update failed";
case KalmanFitterError::SmoothFailed:
return "Kalman smooth failed";
case KalmanFitterError::OutputConversionFailed:
diff --git a/Core/include/Acts/Propagator/detail/PointwiseMaterialInteraction.hpp b/Core/include/Acts/Propagator/detail/PointwiseMaterialInteraction.hpp
index 76990252591b314153655a9cf865feb3ef8ea0e6..766a6a166878818876bdcce86224e809616b97ea 100644
--- a/Core/include/Acts/Propagator/detail/PointwiseMaterialInteraction.hpp
+++ b/Core/include/Acts/Propagator/detail/PointwiseMaterialInteraction.hpp
@@ -144,9 +144,11 @@ struct PointwiseMaterialInteraction {
///
/// @param [in] variance A diagonal entry of the covariance matrix
/// @param [in] change The change that may be applied to it
+ /// @param [in] updateMode The noise update mode (in default: add noise)
///
/// @return The updated variance
- double updateVariance(double variance, double change) const;
+ double updateVariance(double variance, double change,
+ NoiseUpdateMode updateMode = addNoise) const;
};
} // namespace detail
} // end of namespace Acts
diff --git a/Core/include/Acts/Utilities/Definitions.hpp b/Core/include/Acts/Utilities/Definitions.hpp
index 562ea0321f606ef40474a97a719acc12e3a6d79c..7dbb61af0439c95dd26589225bbd78b2b40c3fb6 100644
--- a/Core/include/Acts/Utilities/Definitions.hpp
+++ b/Core/include/Acts/Utilities/Definitions.hpp
@@ -55,6 +55,12 @@ enum MaterialUpdateStage : int {
postUpdate = 1
};
+/// @enum NoiseUpdateMode to tell how to deal with noise term in covariance
+/// transport
+/// - removeNoise: subtract noise term
+/// - addNoise: add noise term
+enum NoiseUpdateMode : int { removeNoise = -1, addNoise = 1 };
+
// Eigen definitions
template <typename T, unsigned int rows, unsigned int cols>
using ActsMatrix = Eigen::Matrix<T, rows, cols>;
diff --git a/Core/src/Propagator/detail/PointwiseMaterialInteraction.cpp b/Core/src/Propagator/detail/PointwiseMaterialInteraction.cpp
index efb940d8d6627ed0213232786931c121d88e25ca..8669d5a197634cbd95ffd5a9674a480421eafc96 100644
--- a/Core/src/Propagator/detail/PointwiseMaterialInteraction.cpp
+++ b/Core/src/Propagator/detail/PointwiseMaterialInteraction.cpp
@@ -44,11 +44,11 @@ void PointwiseMaterialInteraction::covarianceContributions(
}
}
-double PointwiseMaterialInteraction::updateVariance(double variance,
- double change) const {
- // Add/Subtract the change (depending on the direction)
+double PointwiseMaterialInteraction::updateVariance(
+ double variance, double change, NoiseUpdateMode updateMode) const {
+ // Add/Subtract the change
// Protect the variance against becoming negative
- return std::max(0., variance + std::copysign(change, nav));
+ return std::max(0., variance + std::copysign(change, updateMode));
}
} // namespace detail
} // end of namespace Acts
diff --git a/Tests/UnitTests/Core/Fitter/KalmanFitterTests.cpp b/Tests/UnitTests/Core/Fitter/KalmanFitterTests.cpp
index 6a32c65f12311ebcdbb037fa8a2485f345e0b3d6..7023d3278d67f1e730c3424fff115b2ca30d0c6f 100644
--- a/Tests/UnitTests/Core/Fitter/KalmanFitterTests.cpp
+++ b/Tests/UnitTests/Core/Fitter/KalmanFitterTests.cpp
@@ -373,6 +373,25 @@ BOOST_AUTO_TEST_CASE(kalman_fitter_zero_field) {
BOOST_CHECK(!Acts::Test::checkCloseRel(fittedParameters.parameters(),
fittedWithHoleParameters.parameters(),
1e-6));
+ // Run KF fit in backward filtering mode
+ kfOptions.backwardFiltering = true;
+ // Fit the track
+ fitRes = kFitter.fit(sourcelinks, rStart, kfOptions);
+ BOOST_CHECK(fitRes.ok());
+ auto fittedWithBwdFiltering = *fitRes;
+ // Check the filtering and smoothing status flag
+ BOOST_CHECK(fittedWithBwdFiltering.forwardFiltered);
+ BOOST_CHECK(not fittedWithBwdFiltering.smoothed);
+
+ // Count the number of 'smoothed' states
+ auto& trackTip = fittedWithBwdFiltering.trackTip;
+ auto& mj = fittedWithBwdFiltering.fittedStates;
+ size_t nSmoothed = 0;
+ mj.visitBackwards(trackTip, [&](const auto& state) {
+ if (state.hasSmoothed())
+ nSmoothed++;
+ });
+ BOOST_CHECK_EQUAL(nSmoothed, 6u);
}
} // namespace Test