Abort backward propagation when path length gets too negative

core/include/traccc/edm/track_state.hpp

@@ -200,6 +200,14 @@ struct track_state {
         return m_smoothed;
     }
 
+    /// @return the non-const path length
+    TRACCC_HOST_DEVICE
+    inline scalar_type& path_length() { return m_path_length; }
+
+    /// @return the const path length
+    TRACCC_HOST_DEVICE
+    inline const scalar_type& path_length() const { return m_path_length; }
+
     public:
     bool is_hole{true};
     bool is_smoothed{false};
@@ -214,6 +222,7 @@ struct track_state {
     scalar_type m_smoothed_chi2 = 0.f;
     bound_track_parameters_type m_smoothed;
     scalar_type m_backward_chi2 = 0.f;
+    scalar_type m_path_length = detray::detail::invalid_value<scalar_type>();
 };
 
 /// Declare all track_state collection types

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

@@ -94,6 +94,8 @@ struct kalman_actor : detray::actor {
 
         // Run back filtering for smoothing, if true
         bool backward_mode = false;
+        traccc::scalar min_path_length =
+            -100.f * traccc::unit<traccc::scalar>::mm;
     };
 
     /// Actor operation to perform the Kalman filtering
@@ -145,6 +147,9 @@ struct kalman_actor : detray::actor {
                 // Update the propagation flow
                 stepping.bound_params() = trk_state.filtered();
 
+                // Set path length
+                trk_state.path_length() = propagation._stepping.path_length();
+
                 // Set full jacobian
                 trk_state.jacobian() = stepping.full_jacobian();
             } else {
@@ -172,6 +177,15 @@ struct kalman_actor : detray::actor {
             // Flag renavigation of the current candidate
             navigation.set_high_trust();
         }
+
+        // Abort if the path length becomes too negative during the backward
+        // propagation
+        // TODO: Use configuration instead of hardcoded value
+        if (actor_state.backward_mode &&
+            propagation._stepping.path_length() < actor_state.min_path_length) {
+            propagation._heartbeat &= navigation.abort();
+            return;
+        }
     }
 };
 

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

@@ -270,13 +270,19 @@ class kalman_fitter {
 
             typename backward_propagator_type::state propagation(
                 last.smoothed(), m_field, m_detector);
+
             propagation.set_particle(detail::correct_particle_hypothesis(
                 m_cfg.ptc_hypothesis, last.smoothed()));
 
             assert(std::signbit(
                        propagation._stepping.particle_hypothesis().charge()) ==
                    std::signbit(propagation._stepping.bound_params().qop()));
 
+            // Update path langth from the forward filter
+            assert(last.path_length() !=
+                   detray::detail::invalid_value<scalar_type>());
+            propagation._stepping.update_path_lengths(last.path_length());
+
             inflate_covariance(propagation._stepping.bound_params(),
                                m_cfg.covariance_inflation_factor);
 

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

@@ -149,6 +149,10 @@ struct two_filters_smoother {
         const matrix_type<1, 1> chi2 =
             matrix::transpose(residual) * matrix::inverse(R) * residual;
 
+        // Update the bound track parameters
+        bound_params.set_vector(filtered_vec);
+        bound_params.set_covariance(filtered_cov);
+
         // Return false if track is parallel to z-axis or phi is not finite
         const scalar theta = bound_params.theta();
         if (theta <= 0.f || theta >= constant<traccc::scalar>::pi) {
@@ -159,9 +163,9 @@ struct two_filters_smoother {
             return kalman_fitter_status::ERROR_INVERSION;
         }
 
-        // Update the bound track parameters
-        bound_params.set_vector(filtered_vec);
-        bound_params.set_covariance(filtered_cov);
+        if (std::abs(bound_params.qop()) == 0.f) {
+            return kalman_fitter_status::ERROR_QOP_ZERO;
+        }
 
         // Set backward chi2
         trk_state.backward_chi2() = getter::element(chi2, 0, 0);

tests/common/tests/kalman_fitting_test.cpp

@@ -158,7 +158,7 @@ void KalmanFittingTests::ndf_tests(
 
     for (const auto& state : track_states_per_track) {
 
-        if (!state.is_hole) {
+        if (!state.is_hole && state.is_smoothed) {
 
             dim_sum += static_cast<scalar>(state.get_measurement().meas_dim);
             n_effective_states++;

tests/cpu/test_kalman_fitter_momentum_resolution.cpp

@@ -127,6 +127,8 @@ TEST_P(KalmanFittingMomentumResolutionTests, Run) {
         std::move(smearer_writer_cfg), full_path);
     sim.get_config().propagation.stepping.rk_error_tol =
         1e-8f * unit<float>::mm;
+    sim.get_config().propagation.navigation.overstep_tolerance =
+        -1.f * unit<float>::mm;
     sim.run();
 
     /***************
@@ -140,6 +142,7 @@ TEST_P(KalmanFittingMomentumResolutionTests, Run) {
     traccc::fitting_config fit_cfg;
     fit_cfg.ptc_hypothesis = ptc;
     fit_cfg.propagation.stepping.rk_error_tol = 1e-8f * unit<float>::mm;
+    fit_cfg.propagation.navigation.overstep_tolerance = -1.f * unit<float>::mm;
     fit_cfg.use_backward_filter = true;
     traccc::host::kalman_fitting_algorithm fitting(fit_cfg, host_mr);