Improve kalman fitter error handling (#1191)

Expand on the number of Kalman fit status codes and track fit outcomes to make a more fine-grained distinction between different error types in the Kalman fitter possible. Takes up the improvements in #1177.

Edit: The wire chamber Kalman fitter test had to be adapted, because more error states are recognized

Edit2: Also fixed the condition of is_complete for the backward fit to point one position below the first state

Author: niermann999

core/include/traccc/edm/impl/track_state_helpers.ipp

@@ -19,8 +19,9 @@ TRACCC_HOST_DEVICE
     // Create the result object.
     typename track_state_collection<algebra_t>::device::object_type state;
 
-    // Set it to be a hole by default, with the appropriate (measurement) index.
-    state.set_hole();
+    // Set it not to be a hole by default, with the appropriate (measurement)
+    // index.
+    state.set_hole(false);
     state.measurement_index() = mindex;
 
     // Set the correct surface link for the track parameters.

core/include/traccc/edm/track_fit_outcome.hpp

@@ -17,7 +17,12 @@ enum class track_fit_outcome : std::uint16_t {
     UNKNOWN,
     SUCCESS,
     FAILURE_NON_POSITIVE_NDF,
+    FAILURE_NOT_ALL_FITTED,
     FAILURE_NOT_ALL_SMOOTHED,
+    FAILURE_FITTER,
+    FAILURE_SMOOTHER,
+    FAILURE_FORWARD_PROPAGATION,
+    FAILURE_BACKWARD_PROPAGATION,
     MAX_OUTCOME
 };
 

core/include/traccc/finding/actors/ckf_aborter.hpp

@@ -51,14 +51,10 @@ struct ckf_aborter : detray::actor {
             abrt_state.path_from_surface > abrt_state.min_step_length) {
             prop_state._heartbeat &= navigation.pause();
             abrt_state.success = true;
-        }
-
-        TRACCC_VERBOSE_HOST_DEVICE("-> Found sensitive surface: %d",
-                                   navigation.barcode().index());
-
-        // Reset path from surface
-        if (navigation.is_on_sensitive()) {
             abrt_state.path_from_surface = 0.f;
+
+            TRACCC_VERBOSE_HOST_DEVICE("-> Found sensitive surface: %d",
+                                       navigation.barcode().index());
         }
 
         if (abrt_state.count > abrt_state.max_count) {

core/include/traccc/finding/details/combinatorial_kalman_filter.hpp

@@ -522,15 +522,16 @@ combinatorial_kalman_filter(
             TRACCC_DEBUG_HOST("Propagating... ");
             propagator.propagate(propagation,
                                  detray::tie(s0, s1, s2, s3, s4, s5));
-            TRACCC_DEBUG_HOST("Finished propagation: On surface "
-                              << propagation._navigation.barcode());
+            TRACCC_DEBUG_HOST("Finished propagation");
 
             // If a surface found, add the parameter for the next
             // step
             bool valid_track{s5.success};
             if (valid_track) {
                 assert(propagation._navigation.is_on_sensitive());
                 assert(!propagation._stepping.bound_params().is_invalid());
+                TRACCC_DEBUG_HOST(
+                    "On surface: " << propagation._navigation.barcode());
 
                 const auto& out_param = propagation._stepping.bound_params();
 

core/include/traccc/finding/details/combinatorial_kalman_filter_types.hpp

@@ -14,15 +14,7 @@
 #include "traccc/finding/actors/interaction_register.hpp"
 
 // Detray include(s).
-#include <detray/navigation/navigator.hpp>
-#include <detray/propagator/actor_chain.hpp>
-#include <detray/propagator/actors/aborters.hpp>
-#include <detray/propagator/actors/parameter_resetter.hpp>
-#include <detray/propagator/actors/parameter_transporter.hpp>
-#include <detray/propagator/actors/pointwise_material_interactor.hpp>
-#include <detray/propagator/constrained_step.hpp>
-#include <detray/propagator/propagator.hpp>
-#include <detray/propagator/rk_stepper.hpp>
+#include "traccc/utils/propagation.hpp"
 
 // System include(s).
 #include <type_traits>

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

@@ -210,8 +210,7 @@ struct gain_matrix_updater {
 
         const scalar theta = trk_state.filtered_params().theta();
         if (theta <= 0.f || theta >= 2.f * constant<traccc::scalar>::pi) {
-            TRACCC_ERROR_HOST_DEVICE("Hit theta pole after filtering : %f",
-                                     theta);
+            TRACCC_ERROR_HOST_DEVICE("Hit theta pole in filtering : %f", theta);
             return kalman_fitter_status::ERROR_THETA_POLE;
         }
 

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

@@ -52,13 +52,29 @@ struct kalman_actor_state {
         reset();
     }
 
+    /// Get the track state at a given position along the track
+    TRACCC_HOST_DEVICE
+    typename edm::track_state_collection<algebra_t>::device::proxy_type at(
+        unsigned int i) {
+        assert(m_track.constituent_links().at(i).type ==
+               edm::track_constituent_link::track_state);
+        return m_track_states.at(m_track.constituent_links().at(i).index);
+    }
+
+    /// Get the track state at a given position along the track
+    TRACCC_HOST_DEVICE
+    auto at(unsigned int i) const {
+        assert(m_track.constituent_links().at(i).type ==
+               edm::track_constituent_link::track_state);
+        return m_track_states.at(m_track.constituent_links().at(i).index);
+    }
+
     /// @return the reference of track state pointed by the iterator
     TRACCC_HOST_DEVICE
     typename edm::track_state_collection<algebra_t>::device::proxy_type
     operator()() {
-        assert(m_track.constituent_links().at(m_idx).type ==
-               edm::track_constituent_link::track_state);
-        return m_track_states.at(m_track.constituent_links().at(m_idx).index);
+        assert(m_idx >= 0);
+        return at(static_cast<unsigned int>(m_idx));
     }
 
     /// Reset the iterator
@@ -67,8 +83,9 @@ struct kalman_actor_state {
         if (!backward_mode) {
             m_idx = 0;
         } else {
-            m_idx = m_track.constituent_links().size() - 1;
+            m_idx = static_cast<int>(size()) - 1;
         }
+        n_holes = 0u;
     }
 
     /// Advance the iterator
@@ -81,22 +98,65 @@ struct kalman_actor_state {
         }
     }
 
+    /// @TODO: Const-correctness broken due to a vecmem bug
+    /// @returns the number of track states
+    TRACCC_HOST_DEVICE
+    unsigned int size() /*const*/ { return m_track.constituent_links().size(); }
+
     /// @return true if the iterator reaches the end of vector
     TRACCC_HOST_DEVICE
-    bool is_complete() {
-        if (!backward_mode && m_idx == m_track.constituent_links().size()) {
-            return true;
-        } else if (backward_mode &&
-                   m_idx > m_track.constituent_links().size()) {
-            return true;
+    bool is_complete() /*const*/ {
+        return (!backward_mode && m_idx == static_cast<int>(size())) ||
+               (backward_mode && m_idx == -1);
+    }
+
+    /// @TODO: Const-correctness broken due to a vecmem bug
+    TRACCC_HOST_DEVICE
+    bool is_state() /* const*/ {
+        assert(m_idx >= 0);
+        return (m_track.constituent_links()
+                    .at(static_cast<unsigned int>(m_idx))
+                    .type == edm::track_constituent_link::track_state);
+    }
+
+    /// @TODO: Const-correctness broken due to a vecmem bug
+    /// @returns the current number of missed states during forward fit
+    TRACCC_HOST_DEVICE
+    unsigned int count_missed_fit() /*const*/ {
+        unsigned int n_missed{0u};
+
+        for (unsigned int i = 0u; i < size(); ++i) {
+            const auto trk_state = at(i);
+            if (!trk_state.is_hole() &&
+                trk_state.filtered_params().is_invalid()) {
+                TRACCC_DEBUG_HOST_DEVICE(
+                    "Missed track state %d/%d on surface %d during forward fit",
+                    i, size(), at(i).filtered_params().surface_link().index());
+                ++n_missed;
+            }
         }
-        return false;
+
+        return n_missed;
     }
 
+    /// @TODO: Const-correctness broken due to a vecmem bug
+    /// @returns the current number of missed states during smoothing
     TRACCC_HOST_DEVICE
-    bool is_state() {
-        return (m_track.constituent_links().at(m_idx).type ==
-                edm::track_constituent_link::track_state);
+    unsigned int count_missed_smoother() /*const*/ {
+        unsigned int n_missed{0u};
+
+        for (unsigned int i = 0u; i < size(); ++i) {
+            const auto trk_state = at(i);
+            if (!trk_state.is_hole() &&
+                trk_state.smoothed_params().is_invalid()) {
+                TRACCC_DEBUG_HOST_DEVICE(
+                    "Missed track state %d/%d on surface %d during smoothing",
+                    i, size(), at(i).smoothed_params().surface_link().index());
+                ++n_missed;
+            }
+        }
+
+        return n_missed;
     }
 
     /// Object describing the track fit
@@ -108,14 +168,17 @@ struct kalman_actor_state {
         m_measurements;
 
     /// Index of the current track state
-    unsigned int m_idx;
+    int m_idx;
 
-    // The number of holes (The number of sensitive surfaces which do not
-    // have a measurement for the track pattern)
+    /// The number of holes (The number of sensitive surfaces which do not
+    /// have a measurement for the track pattern)
     unsigned int n_holes{0u};
 
-    // Run back filtering for smoothing, if true
+    /// Run back filtering for smoothing, if true
     bool backward_mode = false;
+
+    /// Result of the fitter pass
+    kalman_fitter_status fit_result = kalman_fitter_status::SUCCESS;
 };
 
 /// Detray actor for Kalman filtering
@@ -142,20 +205,18 @@ struct kalman_actor : detray::actor {
             return;
         }
 
+        TRACCC_VERBOSE_HOST_DEVICE("In Kalman actor...");
+        TRACCC_VERBOSE_HOST(
+            "Expected: " << actor_state().filtered_params().surface_link());
+
         // triggered only for sensitive surfaces
         if (navigation.is_on_sensitive()) {
 
-            TRACCC_VERBOSE_HOST_DEVICE("In Kalman actor...");
             TRACCC_DEBUG_HOST("-> on surface: " << navigation.get_surface());
 
             typename edm::track_state_collection<algebra_t>::device::proxy_type
                 trk_state = actor_state();
 
-            TRACCC_DEBUG_HOST(
-                "-> expecting: " << actor_state.m_measurements
-                                        .at(trk_state.measurement_index())
-                                        .surface_link);
-
             // Increase the hole counts if the propagator fails to find the next
             // measurement
             if (navigation.barcode() !=
@@ -167,36 +228,28 @@ struct kalman_actor : detray::actor {
                 return;
             }
 
-            TRACCC_VERBOSE_HOST_DEVICE("Found next track state to fit");
-
-            // This track state is not a hole
-            if (!actor_state.backward_mode) {
-                trk_state.set_hole(false);
-            }
+            auto& bound_param = stepping.bound_params();
 
             // Run Kalman Gain Updater
             const auto sf = navigation.get_surface();
-
             const bool is_line = detail::is_line(sf);
 
-            kalman_fitter_status res = kalman_fitter_status::SUCCESS;
-
             if (!actor_state.backward_mode) {
                 if constexpr (direction_e ==
                                   kalman_actor_direction::FORWARD_ONLY ||
                               direction_e ==
                                   kalman_actor_direction::BIDIRECTIONAL) {
                     // Wrap the phi and theta angles in their valid ranges
-                    normalize_angles(propagation._stepping.bound_params());
+                    normalize_angles(bound_param);
 
                     // Forward filter
-                    TRACCC_DEBUG_HOST_DEVICE("Run filtering");
-                    res = gain_matrix_updater<algebra_t>{}(
-                        trk_state, actor_state.m_measurements,
-                        propagation._stepping.bound_params(), is_line);
+                    TRACCC_DEBUG_HOST_DEVICE("Run filtering...");
+                    actor_state.fit_result = gain_matrix_updater<algebra_t>{}(
+                        trk_state, actor_state.m_measurements, bound_param,
+                        is_line);
 
                     // Update the propagation flow
-                    stepping.bound_params() = trk_state.filtered_params();
+                    bound_param = trk_state.filtered_params();
                 } else {
                     assert(false);
                 }
@@ -205,39 +258,50 @@ struct kalman_actor : detray::actor {
                                   kalman_actor_direction::BACKWARD_ONLY ||
                               direction_e ==
                                   kalman_actor_direction::BIDIRECTIONAL) {
-                    // Backward filter for smoothing
-                    TRACCC_DEBUG_HOST_DEVICE("Run smoothing");
-                    res = two_filters_smoother<algebra_t>{}(
-                        trk_state, actor_state.m_measurements,
-                        propagation._stepping.bound_params(), is_line);
+                    TRACCC_DEBUG_HOST_DEVICE("Run smoothing...");
+
+                    // Forward filter did not find this state: cannot smoothe
+                    if (trk_state.filtered_params().is_invalid()) {
+                        TRACCC_ERROR_HOST_DEVICE(
+                            "Track state not filtered by forward fit. "
+                            "Skipping");
+                        actor_state.fit_result =
+                            kalman_fitter_status::ERROR_UPDATER_SKIPPED_STATE;
+                    } else {
+                        actor_state.fit_result =
+                            two_filters_smoother<algebra_t>{}(
+                                trk_state, actor_state.m_measurements,
+                                bound_param, is_line);
+                    }
                 } else {
                     assert(false);
                 }
             }
 
             // Abort if the Kalman update fails
-            if (res != kalman_fitter_status::SUCCESS) {
+            if (actor_state.fit_result != kalman_fitter_status::SUCCESS) {
                 if (actor_state.backward_mode) {
                     TRACCC_ERROR_DEVICE("Abort backward fit: KF status %d",
-                                        res);
+                                        actor_state.fit_result);
                     TRACCC_ERROR_HOST(
-                        "Abort backward fit: " << fitter_debug_msg{res}());
+                        "Abort backward fit: "
+                        << fitter_debug_msg{actor_state.fit_result}());
                 } else {
-                    TRACCC_ERROR_DEVICE("Abort forward fit: KF status %d", res);
-                    TRACCC_ERROR_HOST(
-                        "Abort forward fit: " << fitter_debug_msg{res}());
+                    TRACCC_ERROR_DEVICE("Abort forward fit: KF status %d",
+                                        actor_state.fit_result);
+                    TRACCC_ERROR_HOST("Abort forward fit: " << fitter_debug_msg{
+                                          actor_state.fit_result}());
                 }
                 propagation._heartbeat &=
-                    navigation.abort(fitter_debug_msg{res});
+                    navigation.abort(fitter_debug_msg{actor_state.fit_result});
                 return;
             }
 
             // Change the charge of hypothesized particles when the sign of qop
             // is changed (This rarely happens when qop is set with a poor seed
             // resolution)
             propagation.set_particle(detail::correct_particle_hypothesis(
-                stepping.particle_hypothesis(),
-                propagation._stepping.bound_params()));
+                stepping.particle_hypothesis(), bound_param));
 
             // Update iterator
             actor_state.next();