Skip to content

Commit be88ea1

Browse files
committed
Compute Kálmán filter chi2 from predicted state
This commit migrates the computation of the chi2 value in our gain matrix updater to use the predicted state rather than the filtered state. This is mathematically equivalent but potentially more stable.
1 parent e49d353 commit be88ea1

File tree

1 file changed

+25
-23
lines changed

1 file changed

+25
-23
lines changed

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

Lines changed: 25 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -66,7 +66,6 @@ struct gain_matrix_updater {
6666
// Some identity matrices
6767
// @TODO: Make constexpr work
6868
const auto I66 = matrix::identity<bound_matrix_type>();
69-
const auto I_m = matrix::identity<matrix_type<D, D>>();
7069

7170
// Measurement data on surface
7271
matrix_type<D, 1> meas_local;
@@ -113,6 +112,31 @@ struct gain_matrix_updater {
113112
TRACCC_DEBUG_HOST("Predicted residual: " << meas_local -
114113
H * predicted_vec);
115114

115+
scalar chi2_val;
116+
117+
{
118+
// Calculate the chi square
119+
const matrix_type<D, D> R =
120+
V + (H * predicted_cov * algebra::matrix::transpose(H));
121+
// Residual between measurement and predicted vector
122+
const matrix_type<D, 1> residual = meas_local - H * predicted_vec;
123+
const matrix_type<1, 1> chi2_mat =
124+
algebra::matrix::transposed_product<true, false>(
125+
residual, matrix::inverse(R)) *
126+
residual;
127+
chi2_val = getter::element(chi2_mat, 0, 0);
128+
129+
if (chi2_val < 0.f) {
130+
TRACCC_ERROR_HOST_DEVICE("Chi2 negative");
131+
return kalman_fitter_status::ERROR_UPDATER_CHI2_NEGATIVE;
132+
}
133+
134+
if (!std::isfinite(chi2_val)) {
135+
TRACCC_ERROR_HOST_DEVICE("Chi2 infinite");
136+
return kalman_fitter_status::ERROR_UPDATER_CHI2_NOT_FINITE;
137+
}
138+
}
139+
116140
const matrix_type<e_bound_size, D> projected_cov =
117141
algebra::matrix::transposed_product<false, true>(predicted_cov, H);
118142

@@ -163,34 +187,12 @@ struct gain_matrix_updater {
163187
return kalman_fitter_status::ERROR_QOP_ZERO;
164188
}
165189

166-
// Residual between measurement and (projected) filtered vector
167-
const matrix_type<D, 1> residual = meas_local - H * filtered_vec;
168-
169-
// Calculate the chi square
170-
const matrix_type<D, D> R = (I_m - H * K) * V;
171-
const matrix_type<1, 1> chi2 =
172-
algebra::matrix::transposed_product<true, false>(
173-
residual, matrix::inverse(R)) *
174-
residual;
175-
176-
const scalar chi2_val{getter::element(chi2, 0, 0)};
177-
178190
TRACCC_VERBOSE_HOST("Filtered residual: " << residual);
179191
TRACCC_DEBUG_HOST("R:\n" << R);
180192
TRACCC_DEBUG_HOST_DEVICE("det(R): %f", matrix::determinant(R));
181193
TRACCC_DEBUG_HOST("R_inv:\n" << matrix::inverse(R));
182194
TRACCC_VERBOSE_HOST_DEVICE("Chi2: %f", chi2_val);
183195

184-
if (chi2_val < 0.f) {
185-
TRACCC_ERROR_HOST_DEVICE("Chi2 negative");
186-
return kalman_fitter_status::ERROR_UPDATER_CHI2_NEGATIVE;
187-
}
188-
189-
if (!std::isfinite(chi2_val)) {
190-
TRACCC_ERROR_HOST_DEVICE("Chi2 infinite");
191-
return kalman_fitter_status::ERROR_UPDATER_CHI2_NOT_FINITE;
192-
}
193-
194196
// Set the chi2 for this track and measurement
195197
trk_state.filtered_params().set_vector(filtered_vec);
196198
trk_state.filtered_params().set_covariance(filtered_cov);

0 commit comments

Comments
 (0)