update to detray v0.104.1

Author: niermann999

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

@@ -47,7 +47,7 @@ struct ckf_aborter : detray::actor {
         if (navigation.is_on_sensitive() &&
             abrt_state.path_from_surface > abrt_state.min_step_length) {
             prop_state._heartbeat &= navigation.pause();
-            abrt_state.success = true;
+            abrt_state.success = navigation.is_alive();
         }
 
         // Reset path from surface
@@ -59,6 +59,7 @@ struct ckf_aborter : detray::actor {
             prop_state._heartbeat &= navigation.abort(
                 "CKF: Maximum number of steps to reach next sensitive surface "
                 "exceeded");
+            abrt_state.success = false;
         }
     }
 };

core/include/traccc/finding/candidate_link.hpp

@@ -28,6 +28,9 @@ struct candidate_link {
     // Index to the initial seed
     unsigned int seed_idx;
 
+    // How many candidates were added to the track
+    unsigned int n_cand;
+
     // How many times it skipped a surface
     unsigned int n_skipped;
 

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

@@ -7,6 +7,8 @@
 #pragma once
 
 // Project include(s).
+#include <detray/utils/log.hpp>
+
 #include "traccc/edm/measurement.hpp"
 #include "traccc/edm/track_container.hpp"
 #include "traccc/edm/track_state_helpers.hpp"
@@ -126,8 +128,11 @@ combinatorial_kalman_filter(
     std::vector<std::pair<unsigned int, unsigned int>> tips;
 
     // Create propagator
+    // @TODO: Turn off the mask tolerance scaling for now
+    auto prop_config{config.propagation};
+    prop_config.navigation.estimate_scattering_noise = false;
     traccc::details::ckf_propagator_t<detector_t, bfield_t> propagator(
-        config.propagation);
+        prop_config);
 
     // Create the input seeds container.
     bound_track_parameters_collection_types::const_device seeds{seeds_view};
@@ -139,15 +144,20 @@ combinatorial_kalman_filter(
 
     std::vector<bound_track_parameters<algebra_type>> out_params;
 
+    std::vector<std::uint8_t> in_is_edge(seeds.size(), false);
+    std::vector<std::uint8_t> out_is_edge;
+
     for (unsigned int step = 0u; step < config.max_track_candidates_per_track;
          step++) {
 
+        DETRAY_DEBUG_HOST("PRPGAGATION STEP: " << step);
         TRACCC_VERBOSE("Starting step "
                        << step + 1 << " / "
                        << config.max_track_candidates_per_track);
 
         // Iterate over input parameters
         const std::size_t n_in_params = in_params.size();
+        DETRAY_DEBUG_HOST("# PARAMS: " << n_in_params);
 
         // Terminate if there is no parameter to proceed
         if (n_in_params == 0) {
@@ -156,6 +166,7 @@ combinatorial_kalman_filter(
 
         // Rough estimation on out parameters size
         out_params.reserve(n_in_params);
+        out_is_edge.reserve(n_in_params);
 
         // Previous step ID
         std::fill(n_trks_per_seed.begin(), n_trks_per_seed.end(), 0u);
@@ -169,13 +180,22 @@ combinatorial_kalman_filter(
             bound_track_parameters<algebra_type>& in_param =
                 in_params[in_param_id];
 
+            DETRAY_DEBUG_HOST("PARAMS: " << in_param);
+
+            const bool is_edge{in_is_edge[in_param_id] > 0u};
+
             assert(!in_param.is_invalid());
 
             const unsigned int orig_param_id =
                 (step == 0
                      ? in_param_id
                      : links[step - 1][param_to_link[step - 1][in_param_id]]
                            .seed_idx);
+            const unsigned int n_cand =
+                (step == 0
+                     ? 0
+                     : links[step - 1][param_to_link[step - 1][in_param_id]]
+                           .n_cand);
             const unsigned int skip_counter =
                 (step == 0
                      ? 0
@@ -239,7 +259,9 @@ combinatorial_kalman_filter(
                 best_links;
 
             // Iterate over the measurements
+            DETRAY_DEBUG_HOST("# MEAS: " << (up - lo));
             for (unsigned int meas_id = lo; meas_id < up; meas_id++) {
+                DETRAY_DEBUG_HOST("TESTING MEAS: " << meas_id);
 
                 // The measurement on surface to handle.
                 const measurement& meas = measurements.at(meas_id);
@@ -248,7 +270,7 @@ combinatorial_kalman_filter(
                 auto trk_state =
                     edm::make_track_state<algebra_type>(measurements, meas_id);
 
-                const bool is_line = sf.template visit_mask<is_line_visitor>();
+                const bool is_line = detail::is_line(sf);
 
                 // Run the Kalman update on a copy of the track parameters
                 const kalman_fitter_status res =
@@ -257,15 +279,20 @@ combinatorial_kalman_filter(
 
                 const traccc::scalar chi2 = trk_state.filtered_chi2();
 
+                DETRAY_DEBUG_HOST("KF STATUS: " << fitter_debug_msg{res}());
+
                 // The chi2 from Kalman update should be less than chi2_max
                 if (res == kalman_fitter_status::SUCCESS &&
                     chi2 < config.chi2_max) {
 
+                    DETRAY_DEBUG_HOST("FOUND MEAS: " << meas_id);
+
                     best_links.push_back(
                         {{.step = step,
                           .previous_candidate_idx = in_param_id,
                           .meas_idx = meas_id,
                           .seed_idx = orig_param_id,
+                          .n_cand = n_cand + 1,
                           .n_skipped = skip_counter,
                           .n_consecutive_skipped = 0,
                           .chi2 = chi2,
@@ -312,12 +339,15 @@ combinatorial_kalman_filter(
                      .previous_candidate_idx = in_param_id,
                      .meas_idx = std::numeric_limits<unsigned int>::max(),
                      .seed_idx = orig_param_id,
-                     .n_skipped = skip_counter + 1,
+                     .n_cand = n_cand,
+                     .n_skipped = is_edge ? skip_counter : skip_counter + 1,
                      .n_consecutive_skipped = consecutive_skipped + 1,
                      .chi2 = std::numeric_limits<traccc::scalar>::max(),
                      .chi2_sum = prev_chi2_sum,
                      .ndf_sum = prev_ndf_sum});
 
+                DETRAY_DEBUG_HOST("HOLE: " << std::boolalpha << !is_edge);
+
                 updated_params.push_back(in_param);
                 TRACCC_VERBOSE("updated_params["
                                << updated_params.size() - 1
@@ -363,7 +393,7 @@ combinatorial_kalman_filter(
                         prob(Lthisbase.chi2_sum,
                              static_cast<scalar>(Lthisbase.ndf_sum - 5));
 
-                    if (step + 1 - Lthisbase.n_skipped <=
+                    if (Lthisbase.n_cand <=
                             config.duplicate_removal_minimum_length ||
                         Lthisbase.ndf_sum <= 5) {
                         continue;
@@ -377,7 +407,7 @@ combinatorial_kalman_filter(
                         auto Lthis = Lthisbase;
                         auto Lthat = links.at(step).at(tracks.at(j));
 
-                        if (step + 1 - Lthat.n_skipped <=
+                        if (Lthat.n_cand <=
                                 config.duplicate_removal_minimum_length ||
                             Lthisbase.ndf_sum <= 5) {
                             continue;
@@ -438,6 +468,7 @@ combinatorial_kalman_filter(
                         }
 
                         if (this_is_dominated) {
+                            DETRAY_DEBUG_HOST("Track is DEAD!");
                             param_liveness.at(tracks.at(i)) = 0u;
                             break;
                         }
@@ -465,6 +496,8 @@ combinatorial_kalman_filter(
             // link to be a tip
             if (links.at(step).at(link_id).n_skipped >
                 config.max_num_skipping_per_cand) {
+
+                DETRAY_DEBUG_HOST("MAKE TIP: MAX NO. HOLES");
                 tips.push_back({step, link_id});
                 continue;
             }
@@ -484,8 +517,8 @@ combinatorial_kalman_filter(
             traccc::details::ckf_interactor_t::state s2;
             typename interaction_register<
                 traccc::details::ckf_interactor_t>::state s1{s2};
-            // typename detray::parameter_resetter<
-            //     typename detector_t::algebra_type>::state s3{};
+            typename detray::parameter_resetter<
+                typename detector_t::algebra_type>::state s3{prop_config};
             typename detray::momentum_aborter<scalar_type>::state s4{};
             typename ckf_aborter::state s5;
             // Update the actor config
@@ -494,8 +527,13 @@ combinatorial_kalman_filter(
             s5.min_step_length = config.min_step_length_for_next_surface;
             s5.max_count = config.max_step_counts_for_next_surface;
 
+            DETRAY_DEBUG_HOST("PROPAGATING... ");
             // Propagate to the next surface
-            propagator.propagate(propagation, detray::tie(s0, s1, s2, s4, s5));
+            propagator.propagate(propagation,
+                                 detray::tie(s0, s1, s2, s3, s4, s5));
+
+            DETRAY_DEBUG_HOST(
+                "FINISHED PROPAGATING. PATH:  " << s5.path_from_surface);
 
             // If a surface found, add the parameter for the next
             // step
@@ -504,25 +542,32 @@ combinatorial_kalman_filter(
                 assert(!propagation._stepping.bound_params().is_invalid());
 
                 out_params.push_back(propagation._stepping.bound_params());
+                out_is_edge.push_back(
+                    propagation._navigation.is_edge_candidate());
                 param_to_link[step].push_back(link_id);
             }
             // Unless the track found a surface, it is considered a
             // tip
             else if (!s5.success &&
                      (step >= (config.min_track_candidates_per_track - 1u))) {
+
+                DETRAY_DEBUG_HOST("MAKE TIP: NO SENSITIVE");
                 tips.push_back({step, link_id});
             }
 
             // If no more CKF step is expected, current candidate is
             // kept as a tip
             if (s5.success &&
                 (step == (config.max_track_candidates_per_track - 1u))) {
+                DETRAY_DEBUG_HOST("MAKE TIP: MAX CANDIDATES");
                 tips.push_back({step, link_id});
             }
         }
 
         in_params = std::move(out_params);
+        in_is_edge = std::move(out_is_edge);
         out_params.clear();
+        out_is_edge.clear();
     }
 
     /**********************
@@ -538,7 +583,7 @@ combinatorial_kalman_filter(
         // Get the link corresponding to tip
         auto L = links.at(tip.first).at(tip.second);
 
-        const unsigned int n_cands = tip.first + 1 - L.n_skipped;
+        const unsigned int n_cands = L.n_cand;
 
         // Skip if the number of tracks candidates is too small
         if (n_cands < config.min_track_candidates_per_track ||

core/include/traccc/fitting/details/kalman_fitting.hpp

@@ -88,7 +88,7 @@ typename edm::track_container<algebra_t>::host kalman_fitting(
             result_tracks_device.at(result_tracks_device.size() - 1),
             typename edm::track_state_collection<algebra_t>::device{
                 vecmem::get_data(result.states)},
-            measurements, seqs_buffer);
+            measurements, seqs_buffer, fitter.config().propagation);
 
         // Run the fitter. The status that it returns is not used here. The main
         // failure modes are saved onto the fitted track itself. Not sure what

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

@@ -7,20 +7,31 @@
 
 #pragma once
 
+#include <detray/geometry/mask.hpp>
 #include <detray/geometry/shapes/line.hpp>
+#include <detray/geometry/surface.hpp>
 
 #include "traccc/definitions/qualifiers.hpp"
 
-namespace traccc {
+namespace traccc::detail {
 
-struct is_line_visitor {
-    template <typename mask_group_t, typename index_t>
-    [[nodiscard]] TRACCC_HOST_DEVICE inline bool operator()(
-        const mask_group_t& /*mask_group*/, const index_t& /*index*/) const {
-        using shape_type = typename mask_group_t::value_type::shape;
-        return std::is_same_v<shape_type, detray::line<true>> ||
-               std::is_same_v<shape_type, detray::line<false>>;
+/// @returns true if the surface has "line" shape
+template <typename detector_t>
+[[nodiscard]] TRACCC_HOST_DEVICE bool constexpr is_line(
+    const detray::geometry::surface<detector_t> sf) {
+    using algebra_t = typename detector_t::algebra_type;
+    using straw_tube = detray::mask<detray::line<false>, algebra_t>;
+    using wire_cell = detray::mask<detray::line<true>, algebra_t>;
+
+    if constexpr (detector_t::masks::template is_defined<straw_tube>() ||
+                  detector_t::masks::template is_defined<wire_cell>()) {
+        return (sf.shape_id() ==
+                detector_t::masks::template get_id<straw_tube>()) ||
+               (sf.shape_id() ==
+                detector_t::masks::template get_id<wire_cell>());
+    } else {
+        return false;
     }
 };
 
-}  // namespace traccc
+}  // namespace traccc::detail

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

@@ -164,7 +164,7 @@ struct kalman_actor : detray::actor {
             // Run Kalman Gain Updater
             const auto sf = navigation.get_surface();
 
-            const bool is_line = sf.template visit_mask<is_line_visitor>();
+            const bool is_line = detail::is_line(sf);
 
             kalman_fitter_status res = kalman_fitter_status::SUCCESS;
 
@@ -215,7 +215,9 @@ struct kalman_actor : detray::actor {
             actor_state.next();
 
             // Flag renavigation of the current candidate
+            // if (math::fabs(navigation()) > 1.f * unit<float>::um) {
             navigation.set_high_trust();
+            //}
         }
     }
 };

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

@@ -110,28 +110,31 @@ class kalman_fitter {
                 track_states,
             const measurement_collection_types::const_device& measurements,
             vecmem::data::vector_view<detray::geometry::barcode>
-                sequence_buffer)
+                sequence_buffer,
+            const detray::propagation::config& prop_cfg)
             : m_fit_actor_state{track, track_states, measurements},
               m_sequencer_state(
                   vecmem::device_vector<detray::geometry::barcode>(
                       sequence_buffer)),
+              m_parameter_resetter{prop_cfg},
               m_fit_res{track},
               m_sequence_buffer(sequence_buffer) {}
 
         /// @return the actor chain state
         TRACCC_HOST_DEVICE
         typename forward_actor_chain_type::state_ref_tuple operator()() {
             return detray::tie(m_aborter_state, m_interactor_state,
-                               m_fit_actor_state, m_sequencer_state,
-                               m_step_aborter_state);
+                               m_fit_actor_state, m_parameter_resetter,
+                               m_sequencer_state, m_step_aborter_state);
         }
 
         /// @return the actor chain state
         TRACCC_HOST_DEVICE
         typename backward_actor_chain_type::state_ref_tuple
         backward_actor_state() {
             return detray::tie(m_aborter_state, m_fit_actor_state,
-                               m_interactor_state, m_step_aborter_state);
+                               m_interactor_state, m_parameter_resetter,
+                               m_step_aborter_state);
         }
 
         /// Individual actor states
@@ -140,6 +143,7 @@ class kalman_fitter {
         typename forward_fit_actor::state m_fit_actor_state;
         typename barcode_sequencer::state m_sequencer_state;
         kalman_step_aborter::state m_step_aborter_state{};
+        typename resetter::state m_parameter_resetter{};
 
         /// Fitting result per track
         typename edm::track_collection<algebra_type>::device::proxy_type

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

@@ -185,7 +185,7 @@ struct two_filters_smoother {
 
         // Calculate backward chi2
         const matrix_type<D, D> R = (I_m - H * K) * V;
-        // assert(matrix::determinant(R) != 0.f);
+        assert(matrix::determinant(R) != 0.f);
         assert(std::isfinite(matrix::determinant(R)));
         const matrix_type<1, 1> chi2 =
             algebra::matrix::transposed_product<true, false>(