Skip to content

Commit d8261ad

Browse files
niermann999niermann999
authored
Regularize the track covariances in the Kalman filter and patch strip measurement handling (#1184)
Add the check for negative diagonal elements to the KF. If the variance is just slightly negative, take the absolute value, otherwise throw an error. The error is a new Kalman fitter status "ERROR_UPDATER_INVALID_COVARIANCE" or "ERROR_SMOOTHER_INVALID_COVARIANCE". The requirements for the smoothed chi2 have been relaxed, so that the fit will only be aborted if it is below -10. Edit: Also contains a hotfix for the strip measurement projector matrix, which will get properly fixed in a later PR
1 parent 1cf94b9 commit d8261ad

File tree

4 files changed

+117
-50
lines changed

4 files changed

+117
-50
lines changed

‎core/include/traccc/fitting/details/regularize_covariance.hpp‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
/** TRACCC library, part of the ACTS project (R&D line)
2+
*
3+
* (c) 2025 CERN for the benefit of the ACTS project
4+
*
5+
* Mozilla Public License Version 2.0
6+
*/
7+
8+
#pragma once
9+
10+
// Project include(s).
11+
#include "traccc/edm/track_parameters.hpp"
12+
#include "traccc/utils/logging.hpp"
13+
14+
namespace traccc::details {
15+
16+
/// Check the covariance matirx and try to make it positive semi-definite
17+
///
18+
/// @param[out] cov covariance matrix
19+
/// @param[in] min_var variance threshold below which to flag an error
20+
template <detray::concepts::algebra algebra_t>
21+
TRACCC_HOST_DEVICE constexpr bool regularize_covariance(
22+
traccc::bound_matrix<algebra_t>& cov,
23+
const detray::dscalar<algebra_t> min_var) {
24+
25+
if (getter::element(cov, 0, 0) < min_var ||
26+
getter::element(cov, 1, 1) < min_var ||
27+
getter::element(cov, 2, 2) < min_var ||
28+
getter::element(cov, 3, 3) < min_var ||
29+
getter::element(cov, 4, 4) < min_var ||
30+
getter::element(cov, 5, 5) < min_var) {
31+
TRACCC_ERROR_HOST_DEVICE("Negative variance");
32+
return false;
33+
} else if (getter::element(cov, 0, 0) < 0.f ||
34+
getter::element(cov, 1, 1) < 0.f ||
35+
getter::element(cov, 2, 2) < 0.f ||
36+
getter::element(cov, 3, 3) < 0.f ||
37+
getter::element(cov, 4, 4) < 0.f ||
38+
getter::element(cov, 5, 5) < 0.f) {
39+
TRACCC_WARNING_HOST_DEVICE("Negative variance: Regularize...");
40+
}
41+
42+
getter::element(cov, 0, 0) = math::fabs(getter::element(cov, 0, 0));
43+
getter::element(cov, 1, 1) = math::fabs(getter::element(cov, 1, 1));
44+
getter::element(cov, 2, 2) = math::fabs(getter::element(cov, 2, 2));
45+
getter::element(cov, 3, 3) = math::fabs(getter::element(cov, 3, 3));
46+
getter::element(cov, 4, 4) = math::fabs(getter::element(cov, 4, 4));
47+
getter::element(cov, 5, 5) = math::fabs(getter::element(cov, 5, 5));
48+
49+
return true;
50+
}
51+
52+
} // namespace traccc::details

‎core/include/traccc/fitting/kalman_filter/gain_matrix_updater.hpp‎

Lines changed: 19 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313
#include "traccc/edm/measurement_collection.hpp"
1414
#include "traccc/edm/measurement_helpers.hpp"
1515
#include "traccc/edm/track_state_collection.hpp"
16+
#include "traccc/fitting/details/regularize_covariance.hpp"
1617
#include "traccc/fitting/status_codes.hpp"
1718
#include "traccc/utils/logging.hpp"
1819
#include "traccc/utils/subspace.hpp"
@@ -49,24 +50,11 @@ struct gain_matrix_updater {
4950
const bound_track_parameters<algebra_t>& bound_params,
5051
const bool is_line) const {
5152

52-
const auto D =
53-
measurements.at(trk_state.measurement_index()).dimensions();
54-
55-
assert(D == 1u || D == 2u);
56-
57-
return update(trk_state, measurements, bound_params, D, is_line);
58-
}
59-
60-
template <typename track_state_backend_t>
61-
[[nodiscard]] TRACCC_HOST_DEVICE inline kalman_fitter_status update(
62-
typename edm::track_state<track_state_backend_t>& trk_state,
63-
const edm::measurement_collection<default_algebra>::const_device&
64-
measurements,
65-
const bound_track_parameters<algebra_t>& bound_params,
66-
const unsigned int dim, const bool is_line) const {
67-
6853
static constexpr unsigned int D = 2;
6954

55+
[[maybe_unused]] const unsigned int dim{
56+
measurements.at(trk_state.measurement_index()).dimensions()};
57+
7058
TRACCC_VERBOSE_HOST_DEVICE("In gain-matrix-updater...");
7159
TRACCC_VERBOSE_HOST_DEVICE("Measurement dim: %d", dim);
7260

@@ -105,7 +93,8 @@ struct gain_matrix_updater {
10593
getter::element(H, 0u, e_bound_loc0) = -1;
10694
}
10795

108-
if (dim == 1) {
96+
// @TODO: Fix properly
97+
if (/*dim == 1*/ getter::element(meas_local, 1u, 0u) == 0.f) {
10998
getter::element(H, 1u, 0u) = 0.f;
11099
getter::element(H, 1u, 1u) = 0.f;
111100
}
@@ -114,9 +103,9 @@ struct gain_matrix_updater {
114103
matrix_type<D, D> V;
115104
edm::get_measurement_covariance<algebra_t>(
116105
measurements.at(trk_state.measurement_index()), V);
117-
118-
if (dim == 1) {
119-
getter::element(V, 1u, 1u) = 1.f;
106+
// @TODO: Fix properly
107+
if (/*dim == 1*/ getter::element(meas_local, 1u, 0u) == 0.f) {
108+
getter::element(V, 1u, 1u) = 1000.f;
120109
}
121110

122111
TRACCC_DEBUG_HOST("Measurement position: " << meas_local);
@@ -140,13 +129,21 @@ struct gain_matrix_updater {
140129
const matrix_type<6, 1> filtered_vec =
141130
predicted_vec + K * (meas_local - H * predicted_vec);
142131
const matrix_type<6, 6> i_minus_kh = I66 - K * H;
143-
const matrix_type<6, 6> filtered_cov =
132+
matrix_type<6, 6> filtered_cov =
144133
i_minus_kh * predicted_cov * matrix::transpose(i_minus_kh) +
145134
K * V * matrix::transpose(K);
146135

147136
TRACCC_DEBUG_HOST("Filtered param:\n" << filtered_vec);
148137
TRACCC_DEBUG_HOST("Filtered cov:\n" << filtered_cov);
149138

139+
// Check the covariance for consistency
140+
// @TODO: Need to understand why negative variance happens
141+
if (constexpr traccc::scalar min_var{-0.01f};
142+
!details::regularize_covariance<algebra_t>(filtered_cov, min_var)) {
143+
TRACCC_ERROR_HOST_DEVICE("Negative variance after filtering");
144+
return kalman_fitter_status::ERROR_UPDATER_INVALID_COVARIANCE;
145+
}
146+
150147
// Return false if track is parallel to z-axis or phi is not finite
151148
if (!std::isfinite(getter::element(filtered_vec, e_bound_theta, 0))) {
152149
TRACCC_ERROR_HOST_DEVICE(
@@ -204,7 +201,7 @@ struct gain_matrix_updater {
204201
TRACCC_ERROR_HOST_DEVICE(
205202
"Hit theta pole after filtering : %f (unrecoverable error "
206203
"pre-normalization)",
207-
theta);
204+
trk_state.filtered_params().theta());
208205
return kalman_fitter_status::ERROR_THETA_POLE;
209206
}
210207

‎core/include/traccc/fitting/kalman_filter/two_filters_smoother.hpp‎

Lines changed: 36 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313
#include "traccc/edm/measurement_collection.hpp"
1414
#include "traccc/edm/measurement_helpers.hpp"
1515
#include "traccc/edm/track_state_collection.hpp"
16+
#include "traccc/fitting/details/regularize_covariance.hpp"
1617
#include "traccc/fitting/status_codes.hpp"
1718
#include "traccc/utils/logging.hpp"
1819

@@ -43,25 +44,11 @@ struct two_filters_smoother {
4344
bound_track_parameters<algebra_t>& bound_params,
4445
const bool is_line) const {
4546

46-
const auto D =
47-
measurements.at(trk_state.measurement_index()).dimensions();
48-
assert(D == 1u || D == 2u);
49-
50-
return smoothe(trk_state, measurements, bound_params, D, is_line);
51-
}
52-
53-
// Reference: The Optimun Linear Smoother as a Combination of Two Optimum
54-
// Linear Filters
55-
[[nodiscard]] TRACCC_HOST_DEVICE inline kalman_fitter_status smoothe(
56-
typename edm::track_state_collection<algebra_t>::device::proxy_type&
57-
trk_state,
58-
const typename edm::measurement_collection<algebra_t>::const_device&
59-
measurements,
60-
bound_track_parameters<algebra_t>& bound_params, const unsigned int dim,
61-
const bool is_line) const {
62-
6347
static constexpr unsigned int D = 2;
6448

49+
[[maybe_unused]] const unsigned int dim{
50+
measurements.at(trk_state.measurement_index()).dimensions()};
51+
6552
assert(dim == 1u || dim == 2u);
6653

6754
assert(!bound_params.is_invalid());
@@ -102,9 +89,17 @@ struct two_filters_smoother {
10289
predicted_cov_inv + filtered_cov_inv;
10390

10491
assert(matrix::determinant(smoothed_cov_inv) != 0.f);
105-
const matrix_type<e_bound_size, e_bound_size> smoothed_cov =
92+
matrix_type<e_bound_size, e_bound_size> smoothed_cov =
10693
matrix::inverse(smoothed_cov_inv);
10794

95+
// Check the covariance for consistency
96+
// @TODO: Need to understand why negative variance happens
97+
if (constexpr traccc::scalar min_var{-0.01f};
98+
!details::regularize_covariance<algebra_t>(smoothed_cov, min_var)) {
99+
TRACCC_ERROR_HOST_DEVICE("Negative variance after smoothing");
100+
return kalman_fitter_status::ERROR_SMOOTHER_INVALID_COVARIANCE;
101+
}
102+
108103
// Eq (3.38) of "Pattern Recognition, Tracking and Vertex
109104
// Reconstruction in Particle Detectors"
110105
const matrix_type<e_bound_size, 1u> smoothed_vec =
@@ -140,7 +135,8 @@ struct two_filters_smoother {
140135
const subspace<algebra_t, e_bound_size> subs(
141136
measurements.at(trk_state.measurement_index()).subspace());
142137
matrix_type<D, e_bound_size> H = subs.template projector<D>();
143-
if (dim == 1) {
138+
// @TODO: Fix properly
139+
if (getter::element(meas_local, 1u, 0u) == 0.f /*dim == 1*/) {
144140
getter::element(H, 1u, 0u) = 0.f;
145141
getter::element(H, 1u, 1u) = 0.f;
146142
}
@@ -151,8 +147,9 @@ struct two_filters_smoother {
151147
matrix_type<D, D> V;
152148
edm::get_measurement_covariance<algebra_t>(
153149
measurements.at(trk_state.measurement_index()), V);
154-
if (dim == 1) {
155-
getter::element(V, 1u, 1u) = 1.f;
150+
// @TODO: Fix properly
151+
if (getter::element(meas_local, 1u, 0u) == 0.f /*dim == 1*/) {
152+
getter::element(V, 1u, 1u) = 1000.f;
156153
}
157154

158155
TRACCC_DEBUG_HOST("Measurement position: " << meas_local);
@@ -180,12 +177,16 @@ struct two_filters_smoother {
180177
TRACCC_DEBUG_HOST("R_smt:\n" << R_smt);
181178
TRACCC_DEBUG_HOST_DEVICE("det(R_smt): %f", matrix::determinant(R_smt));
182179
TRACCC_DEBUG_HOST("R_smt_inv:\n" << matrix::inverse(R_smt));
183-
TRACCC_VERBOSE_HOST_DEVICE("Chi2: %f", chi2_smt_value);
180+
TRACCC_VERBOSE_HOST_DEVICE("Smoothed chi2: %f", chi2_smt_value);
184181

185182
if (chi2_smt_value < 0.f) {
186183
TRACCC_ERROR_HOST_DEVICE("Smoothed chi2 negative: %f",
187184
chi2_smt_value);
188-
return kalman_fitter_status::ERROR_SMOOTHER_CHI2_NEGATIVE;
185+
186+
// @TODO: Need to understand why negative chi2 happens
187+
if (chi2_smt_value < -10.f) {
188+
return kalman_fitter_status::ERROR_SMOOTHER_CHI2_NEGATIVE;
189+
}
189190
}
190191

191192
if (!std::isfinite(chi2_smt_value)) {
@@ -226,10 +227,18 @@ struct two_filters_smoother {
226227
const matrix_type<6, 1> filtered_vec =
227228
predicted_vec + K * (meas_local - H * predicted_vec);
228229
const matrix_type<6, 6> i_minus_kh = I66 - K * H;
229-
const matrix_type<6, 6> filtered_cov =
230+
matrix_type<6, 6> filtered_cov =
230231
i_minus_kh * predicted_cov * matrix::transpose(i_minus_kh) +
231232
K * V * matrix::transpose(K);
232233

234+
// Check the covariance for consistency
235+
// @TODO: Need to understand why negative variance happens
236+
if (constexpr traccc::scalar min_var{-0.01f};
237+
!details::regularize_covariance<algebra_t>(filtered_cov, min_var)) {
238+
TRACCC_ERROR_HOST_DEVICE("Negative variance after filtering");
239+
return kalman_fitter_status::ERROR_SMOOTHER_INVALID_COVARIANCE;
240+
}
241+
233242
// Update the bound track parameters
234243
bound_params.set_vector(filtered_vec);
235244

@@ -273,15 +282,15 @@ struct two_filters_smoother {
273282
TRACCC_DEBUG_HOST("R:\n" << R);
274283
TRACCC_DEBUG_HOST_DEVICE("det(R): %f", matrix::determinant(R));
275284
TRACCC_DEBUG_HOST("R_inv:\n" << matrix::inverse(R));
276-
TRACCC_VERBOSE_HOST_DEVICE("Chi2: %f", chi2_val);
285+
TRACCC_VERBOSE_HOST_DEVICE("Filtered chi2: %f", chi2_val);
277286

278287
if (chi2_val < 0.f) {
279-
TRACCC_ERROR_HOST_DEVICE("Chi2 negative: %f", chi2_val);
288+
TRACCC_ERROR_HOST_DEVICE("Filtered chi2 negative: %f", chi2_val);
280289
return kalman_fitter_status::ERROR_SMOOTHER_CHI2_NEGATIVE;
281290
}
282291

283292
if (!std::isfinite(chi2_val)) {
284-
TRACCC_ERROR_HOST_DEVICE("Chi2 infinite");
293+
TRACCC_ERROR_HOST_DEVICE("Filtered chi2 infinite");
285294
return kalman_fitter_status::ERROR_SMOOTHER_CHI2_NOT_FINITE;
286295
}
287296

‎core/include/traccc/fitting/status_codes.hpp‎

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,8 +18,10 @@ enum class kalman_fitter_status : uint32_t {
1818
ERROR_INVERSION,
1919
ERROR_SMOOTHER_CHI2_NEGATIVE,
2020
ERROR_SMOOTHER_CHI2_NOT_FINITE,
21+
ERROR_SMOOTHER_INVALID_COVARIANCE,
2122
ERROR_UPDATER_CHI2_NEGATIVE,
2223
ERROR_UPDATER_CHI2_NOT_FINITE,
24+
ERROR_UPDATER_INVALID_COVARIANCE,
2325
ERROR_BARCODE_SEQUENCE_OVERFLOW,
2426
ERROR_INVALID_TRACK_STATE,
2527
ERROR_OTHER,
@@ -48,12 +50,18 @@ struct fitter_debug_msg {
4850
case ERROR_SMOOTHER_CHI2_NOT_FINITE: {
4951
return msg + "Invalid chi2 in smoother";
5052
}
53+
case ERROR_SMOOTHER_INVALID_COVARIANCE: {
54+
return msg + "Invalid track covariance during smoothing";
55+
}
5156
case ERROR_UPDATER_CHI2_NEGATIVE: {
5257
return msg + "Negative chi2 in gain matrix updater";
5358
}
5459
case ERROR_UPDATER_CHI2_NOT_FINITE: {
5560
return msg + "Invalid chi2 in gain matrix updater";
5661
}
62+
case ERROR_UPDATER_INVALID_COVARIANCE: {
63+
return msg + "Invalid track covariance in forward fit";
64+
}
5765
case ERROR_BARCODE_SEQUENCE_OVERFLOW: {
5866
return msg + "Barcode sequence overflow in direct navigator";
5967
}
@@ -65,7 +73,8 @@ struct fitter_debug_msg {
6573
return msg + "Unspecified error";
6674
}
6775
default: {
68-
return "";
76+
return "Not a defined error code: " +
77+
std::to_string(static_cast<int>(m_error_code));
6978
}
7079
}
7180
}

0 commit comments

Comments
 (0)