Protect against negative chi-square in updater

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

@@ -147,6 +147,10 @@ struct gain_matrix_smoother {
         const matrix_type<1, 1> chi2 =
             matrix::transpose(residual) * matrix::inverse(R) * residual;
 
+        if (getter::element(chi2, 0, 0) < 0.f) {
+            return kalman_fitter_status::ERROR_SMOOTHER_CHI2_NEGATIVE;
+        }
+
         cur_state.smoothed_chi2() = getter::element(chi2, 0, 0);
         cur_state.is_smoothed = true;
 

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

@@ -142,6 +142,10 @@ struct gain_matrix_updater {
             return kalman_fitter_status::ERROR_QOP_ZERO;
         }
 
+        if (getter::element(chi2, 0, 0) < 0.f) {
+            return kalman_fitter_status::ERROR_UPDATER_CHI2_NEGATIVE;
+        }
+
         // Set the track state parameters
         trk_state.filtered().set_vector(filtered_vec);
         trk_state.filtered().set_covariance(filtered_cov);

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

@@ -119,6 +119,10 @@ struct two_filters_smoother {
                                            matrix::inverse(R_smt) *
                                            residual_smt;
 
+        if (getter::element(chi2_smt, 0, 0) < 0.f) {
+            return kalman_fitter_status::ERROR_SMOOTHER_CHI2_NEGATIVE;
+        }
+
         trk_state.smoothed_chi2() = getter::element(chi2_smt, 0, 0);
 
         /*************************************
@@ -167,6 +171,10 @@ struct two_filters_smoother {
             return kalman_fitter_status::ERROR_QOP_ZERO;
         }
 
+        if (getter::element(chi2, 0, 0) < 0.f) {
+            return kalman_fitter_status::ERROR_UPDATER_CHI2_NEGATIVE;
+        }
+
         // Set backward chi2
         trk_state.backward_chi2() = getter::element(chi2, 0, 0);
 

core/include/traccc/fitting/status_codes.hpp

@@ -14,6 +14,8 @@ enum class kalman_fitter_status : uint32_t {
     ERROR_QOP_ZERO,
     ERROR_THETA_ZERO,
     ERROR_INVERSION,
+    ERROR_SMOOTHER_CHI2_NEGATIVE,
+    ERROR_UPDATER_CHI2_NEGATIVE,
     ERROR_OTHER,
     MAX_STATUS
 };

device/common/include/traccc/finding/device/impl/find_tracks.ipp

@@ -219,6 +219,8 @@ TRACCC_DEVICE inline void find_tracks(
                 // Add measurement candidates to link
                 const unsigned int l_pos = num_total_candidates.fetch_add(1);
 
+                assert(trk_state.filtered_chi2() >= 0.f);
+
                 if (l_pos >= payload.n_max_candidates) {
                     *payload.n_total_candidates = payload.n_max_candidates;
                 } else {

tests/cpu/test_kalman_fitter_momentum_resolution.cpp

@@ -171,8 +171,10 @@ TEST_P(KalmanFittingMomentumResolutionTests, Run) {
         const std::size_t n_fitted_tracks = count_fitted_tracks(track_states);
 
         // n_trakcs = 100
-        ASSERT_EQ(n_tracks, n_truth_tracks);
-        ASSERT_EQ(n_tracks, n_fitted_tracks);
+        ASSERT_GE(static_cast<float>(n_tracks),
+                  0.98 * static_cast<float>(n_truth_tracks));
+        ASSERT_GE(static_cast<float>(n_tracks),
+                  0.98 * static_cast<float>(n_fitted_tracks));
 
         for (std::size_t i_trk = 0; i_trk < n_tracks; i_trk++) {
 
@@ -219,7 +221,7 @@ TEST_P(KalmanFittingMomentumResolutionTests, Run) {
     float success_rate = static_cast<float>(n_success) /
                          static_cast<float>(n_truth_tracks * n_events);
 
-    ASSERT_FLOAT_EQ(success_rate, 1.00f);
+    ASSERT_GE(success_rate, 0.98f);
 }
 
 // Muon with 1, 10, 100 GeV/c, no materials