From 01d32b0e3919d6eb2a118877a5973bcf271e0bba Mon Sep 17 00:00:00 2001 From: Joana Niermann Date: Fri, 21 Nov 2025 16:58:16 +0100 Subject: [PATCH 1/2] Skip the Kalman fitter on the first CKF step, since it juyst retrieves the first measurement of the seed --- .../details/combinatorial_kalman_filter.hpp | 49 ++++++++++++-- .../finding/device/impl/find_tracks.ipp | 64 +++++++++++++++---- 2 files changed, 97 insertions(+), 16 deletions(-) diff --git a/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp b/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp index 0c27892bb9..fa3dd9af5d 100644 --- a/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp +++ b/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp @@ -260,12 +260,51 @@ combinatorial_kalman_filter( auto trk_state = edm::make_track_state(measurements, meas_id); - const bool is_line = detail::is_line(sf); + kalman_fitter_status res{kalman_fitter_status::ERROR_OTHER}; + // Check if the measurement is even remotely close to the track + if (step == 0 && !sf.has_material()) { + detray::dmatrix meas_local; + edm::get_measurement_local(meas, meas_local); + + const subspace subs( + meas.subspace()); + detray::dmatrix H = + subs.template projector<2>(); + const auto residual = meas_local - H * in_param.vector(); + + const scalar_type res_0{getter::element(residual, 0, 0)}; + const scalar_type res_1{ + meas.dimensions() == 1 + ? 0.f + : getter::element(residual, 1, 0)}; + const scalar_type search_rad = + math::sqrt(res_0 * res_0 + res_1 * res_1); + if (search_rad > 0.f) { + TRACCC_VERBOSE_HOST( + "Measurement outside search radius: " + << search_rad << " mm for meas. " << meas_id); + continue; + } + res = kalman_fitter_status::SUCCESS; + trk_state.filtered_params() = in_param; + trk_state.filtered_chi2() = 0.f; + + // Set the covariance to the measurement variance + detray::dmatrix V; + edm::get_measurement_covariance(meas, V); + auto& filtered_cov = + trk_state.filtered_params().covariance(); + getter::element(filtered_cov, e_bound_loc0, e_bound_loc0) = + getter::element(V, 0, 0); + getter::element(filtered_cov, e_bound_loc1, e_bound_loc1) = + getter::element(V, 1, 1); + } else { + const bool is_line = detail::is_line(sf); - // Run the Kalman update on a copy of the track parameters - const kalman_fitter_status res = - gain_matrix_updater{}(trk_state, measurements, - in_param, is_line); + // Run the Kalman update on a copy of the track parameters + res = gain_matrix_updater{}( + trk_state, measurements, in_param, is_line); + } const traccc::scalar chi2 = trk_state.filtered_chi2(); diff --git a/device/common/include/traccc/finding/device/impl/find_tracks.ipp b/device/common/include/traccc/finding/device/impl/find_tracks.ipp index 09ec69239f..e8f8397ebe 100644 --- a/device/common/include/traccc/finding/device/impl/find_tracks.ipp +++ b/device/common/include/traccc/finding/device/impl/find_tracks.ipp @@ -44,6 +44,8 @@ TRACCC_HOST_DEVICE inline void find_tracks( const finding_config& cfg, const find_tracks_payload& payload, const find_tracks_shared_payload& shared_payload) { + using algebra_t = typename detector_t::algebra_type; + const unsigned int in_param_id = thread_id.getGlobalThreadIdX(); const bool last_step = payload.step == cfg.max_track_candidates_per_track - 1; @@ -165,10 +167,9 @@ TRACCC_HOST_DEVICE inline void find_tracks( barrier.blockBarrier(); - std::optional::device::object_type, - unsigned int, unsigned int>> + std::optional::device::object_type, + unsigned int, unsigned int>> result = std::nullopt; /* @@ -207,17 +208,58 @@ TRACCC_HOST_DEVICE inline void find_tracks( if (use_measurement) { auto trk_state = - edm::make_track_state( - measurements, meas_idx); + edm::make_track_state(measurements, meas_idx); const detray::tracking_surface sf{det, in_par.surface_link()}; - const bool is_line = detail::is_line(sf); - - // Run the Kalman update - const kalman_fitter_status res = - gain_matrix_updater{}( + kalman_fitter_status res{kalman_fitter_status::ERROR_OTHER}; + // Check if the measurement is even remotely close to the track + if (payload.step == 0 && !sf.has_material()) { + const auto& meas = measurements.at(meas_idx); + + detray::dmatrix meas_local; + edm::get_measurement_local(meas, meas_local); + + const subspace subs( + meas.subspace()); + detray::dmatrix H = + subs.template projector<2>(); + + const auto residual = + meas_local - + H * in_params.at(owner_global_thread_id).vector(); + + const traccc::scalar res_0{getter::element(residual, 0, 0)}; + const traccc::scalar res_1{ + meas.dimensions() == 1 + ? 0.f + : getter::element(residual, 1, 0)}; + const traccc::scalar search_rad = + math::sqrt(res_0 * res_0 + res_1 * res_1); + if (search_rad > 0.f) { + use_measurement = false; + } + res = kalman_fitter_status::SUCCESS; + trk_state.filtered_params() = + in_params.at(owner_global_thread_id); + trk_state.filtered_chi2() = 0.f; + + // Set the covariance to the measurement variance + detray::dmatrix V; + edm::get_measurement_covariance(meas, V); + auto& filtered_cov = + trk_state.filtered_params().covariance(); + getter::element(filtered_cov, e_bound_loc0, e_bound_loc0) = + getter::element(V, 0, 0); + getter::element(filtered_cov, e_bound_loc1, e_bound_loc1) = + getter::element(V, 1, 1); + } else { + const bool is_line = detail::is_line(sf); + + // Run the Kalman update on a copy of the track parameters + res = gain_matrix_updater{}( trk_state, measurements, in_par, is_line); + } TRACCC_DEBUG_DEVICE("KF status: %d", res); From 2d959e08a7d137db5bce7caa0fb432a112012d26 Mon Sep 17 00:00:00 2001 From: Joana Niermann Date: Thu, 4 Dec 2025 17:13:35 +0100 Subject: [PATCH 2/2] Try general radius criterion --- core/CMakeLists.txt | 1 + .../details/combinatorial_kalman_filter.hpp | 64 +++--- .../include/traccc/fitting/fitting_config.hpp | 4 + .../kalman_filter/gain_matrix_updater.hpp | 21 +- .../fitting/kalman_filter/kalman_actor.hpp | 8 +- .../kalman_filter/measurement_selector.hpp | 206 ++++++++++++++++++ .../kalman_filter/two_filters_smoother.hpp | 19 +- .../finding/device/impl/find_tracks.ipp | 99 ++++----- 8 files changed, 304 insertions(+), 118 deletions(-) create mode 100644 core/include/traccc/fitting/kalman_filter/measurement_selector.hpp diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index a970dfb48c..642d3e6ca8 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -92,6 +92,7 @@ traccc_add_library( traccc_core core TYPE SHARED "include/traccc/fitting/kalman_filter/kalman_actor.hpp" "include/traccc/fitting/kalman_filter/kalman_fitter.hpp" "include/traccc/fitting/kalman_filter/kalman_step_aborter.hpp" + "include/traccc/fitting/kalman_filter/measurement_selector.hpp" "include/traccc/fitting/kalman_filter/statistics_updater.hpp" "include/traccc/fitting/kalman_filter/two_filters_smoother.hpp" "include/traccc/fitting/details/kalman_fitting_types.hpp" diff --git a/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp b/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp index fa3dd9af5d..029ba7af9c 100644 --- a/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp +++ b/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp @@ -17,6 +17,7 @@ #include "traccc/finding/finding_config.hpp" #include "traccc/fitting/kalman_filter/gain_matrix_updater.hpp" #include "traccc/fitting/kalman_filter/is_line_visitor.hpp" +#include "traccc/fitting/kalman_filter/measurement_selector.hpp" #include "traccc/fitting/status_codes.hpp" #include "traccc/sanity/contiguous_on.hpp" #include "traccc/utils/logging.hpp" @@ -248,6 +249,9 @@ combinatorial_kalman_filter( bound_track_parameters>> best_links; + const bool is_line = detail::is_line(sf); + measurement_selector::config calib_cfg{}; + // Iterate over the measurements TRACCC_VERBOSE_HOST("No. measurements: " << (up - lo)); for (unsigned int meas_id = lo; meas_id < up; meas_id++) { @@ -256,42 +260,34 @@ combinatorial_kalman_filter( // The measurement on surface to handle. const auto meas = measurements.at(meas_id); + const scalar_type chi2 = measurement_selector::predicted_chi2( + meas, in_param, calib_cfg, is_line); + + // If the measurement is outside the chi2 cut, skip it + if (chi2 > config.chi2_max) { + continue; + } + // Create a standalone track state object. auto trk_state = edm::make_track_state(measurements, meas_id); + // Kalman filter status code kalman_fitter_status res{kalman_fitter_status::ERROR_OTHER}; - // Check if the measurement is even remotely close to the track + + // Don't run the filter on the first measurement if (step == 0 && !sf.has_material()) { - detray::dmatrix meas_local; - edm::get_measurement_local(meas, meas_local); - - const subspace subs( - meas.subspace()); - detray::dmatrix H = - subs.template projector<2>(); - const auto residual = meas_local - H * in_param.vector(); - - const scalar_type res_0{getter::element(residual, 0, 0)}; - const scalar_type res_1{ - meas.dimensions() == 1 - ? 0.f - : getter::element(residual, 1, 0)}; - const scalar_type search_rad = - math::sqrt(res_0 * res_0 + res_1 * res_1); - if (search_rad > 0.f) { - TRACCC_VERBOSE_HOST( - "Measurement outside search radius: " - << search_rad << " mm for meas. " << meas_id); - continue; - } res = kalman_fitter_status::SUCCESS; + + // Copy the full track parameters + // TODO: Apply calibration ? trk_state.filtered_params() = in_param; - trk_state.filtered_chi2() = 0.f; - // Set the covariance to the measurement variance - detray::dmatrix V; - edm::get_measurement_covariance(meas, V); + // Update measurement covariance + const auto V = + measurement_selector::calibrated_measurement_covariance< + algebra_type, 2>(meas, calib_cfg); + auto& filtered_cov = trk_state.filtered_params().covariance(); getter::element(filtered_cov, e_bound_loc0, e_bound_loc0) = @@ -299,20 +295,18 @@ combinatorial_kalman_filter( getter::element(filtered_cov, e_bound_loc1, e_bound_loc1) = getter::element(V, 1, 1); } else { - const bool is_line = detail::is_line(sf); - - // Run the Kalman update on a copy of the track parameters - res = gain_matrix_updater{}( - trk_state, measurements, in_param, is_line); + // Run the Kalman update on the track state + constexpr gain_matrix_updater + kalman_updater{}; + res = kalman_updater(trk_state, meas, in_param, is_line); } - const traccc::scalar chi2 = trk_state.filtered_chi2(); + trk_state.filtered_chi2() = chi2; TRACCC_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)) { + if (res == kalman_fitter_status::SUCCESS) { TRACCC_VERBOSE_HOST("Found measurement: " << meas_id); diff --git a/core/include/traccc/fitting/fitting_config.hpp b/core/include/traccc/fitting/fitting_config.hpp index e00e983af1..8d6f720008 100644 --- a/core/include/traccc/fitting/fitting_config.hpp +++ b/core/include/traccc/fitting/fitting_config.hpp @@ -29,6 +29,10 @@ struct fitting_config { float min_p = 100.f * traccc::unit::MeV; float min_pT = 600.f * traccc::unit::MeV; + // Maximal local radius a measurement is allowed to be away from the + // predicted track position in order to be considered as a condidate + float max_measurement_radius{5.f * traccc::unit::mm}; + /// Particle hypothesis traccc::pdg_particle ptc_hypothesis = traccc::muon(); diff --git a/core/include/traccc/fitting/kalman_filter/gain_matrix_updater.hpp b/core/include/traccc/fitting/kalman_filter/gain_matrix_updater.hpp index 413573d576..68e0581222 100644 --- a/core/include/traccc/fitting/kalman_filter/gain_matrix_updater.hpp +++ b/core/include/traccc/fitting/kalman_filter/gain_matrix_updater.hpp @@ -11,12 +11,10 @@ #include "traccc/definitions/qualifiers.hpp" #include "traccc/definitions/track_parametrization.hpp" #include "traccc/edm/measurement_collection.hpp" -#include "traccc/edm/measurement_helpers.hpp" -#include "traccc/edm/track_state_collection.hpp" #include "traccc/fitting/details/regularize_covariance.hpp" +#include "traccc/fitting/kalman_filter/measurement_selector.hpp" #include "traccc/fitting/status_codes.hpp" #include "traccc/utils/logging.hpp" -#include "traccc/utils/subspace.hpp" namespace traccc { @@ -42,18 +40,16 @@ struct gain_matrix_updater { /// @param bound_params bound parameter /// /// @return true if the update succeeds - template + template [[nodiscard]] TRACCC_HOST_DEVICE inline kalman_fitter_status operator()( typename edm::track_state& trk_state, - const edm::measurement_collection::const_device& - measurements, + const edm::measurement& measurement, const bound_track_parameters& bound_params, const bool is_line) const { static constexpr unsigned int D = 2; - [[maybe_unused]] const unsigned int dim{ - measurements.at(trk_state.measurement_index()).dimensions()}; + [[maybe_unused]] const unsigned int dim{measurement.dimensions()}; TRACCC_VERBOSE_HOST_DEVICE("In gain-matrix-updater..."); TRACCC_VERBOSE_HOST_DEVICE("Measurement dim: %d", dim); @@ -70,8 +66,7 @@ struct gain_matrix_updater { // Measurement data on surface matrix_type meas_local; - edm::get_measurement_local( - measurements.at(trk_state.measurement_index()), meas_local); + edm::get_measurement_local(measurement, meas_local); assert((dim > 1) || (getter::element(meas_local, 1u, 0u) == 0.f)); @@ -83,8 +78,7 @@ struct gain_matrix_updater { // Predicted covaraince of bound track parameters const bound_matrix_type& predicted_cov = bound_params.covariance(); - const subspace subs( - measurements.at(trk_state.measurement_index()).subspace()); + const subspace subs(measurement.subspace()); matrix_type H = subs.template projector(); // Flip the sign of projector matrix element in case the first element @@ -101,8 +95,7 @@ struct gain_matrix_updater { // Spatial resolution (Measurement covariance) matrix_type V; - edm::get_measurement_covariance( - measurements.at(trk_state.measurement_index()), V); + edm::get_measurement_covariance(measurement, V); // @TODO: Fix properly if (/*dim == 1*/ getter::element(meas_local, 1u, 0u) == 0.f) { getter::element(V, 1u, 1u) = 1000.f; diff --git a/core/include/traccc/fitting/kalman_filter/kalman_actor.hpp b/core/include/traccc/fitting/kalman_filter/kalman_actor.hpp index affbb4f9ed..92dc96350b 100644 --- a/core/include/traccc/fitting/kalman_filter/kalman_actor.hpp +++ b/core/include/traccc/fitting/kalman_filter/kalman_actor.hpp @@ -406,6 +406,8 @@ struct kalman_actor : detray::actor { const auto sf = navigation.current_surface(); const bool is_line = detail::is_line(sf); + const auto measurement = + actor_state.m_measurements.at(trk_state.measurement_index()); if (!actor_state.backward_mode) { if constexpr (direction_e == kalman_actor_direction::FORWARD_ONLY || @@ -417,8 +419,7 @@ struct kalman_actor : detray::actor { // Forward filter TRACCC_DEBUG_HOST_DEVICE("Run filtering..."); actor_state.fit_result = gain_matrix_updater{}( - trk_state, actor_state.m_measurements, bound_param, - is_line); + trk_state, measurement, bound_param, is_line); // Update the propagation flow bound_param = trk_state.filtered_params(); @@ -447,8 +448,7 @@ struct kalman_actor : detray::actor { } else { actor_state.fit_result = two_filters_smoother{}( - trk_state, actor_state.m_measurements, - bound_param, is_line); + trk_state, measurement, bound_param, is_line); } } else { assert(false); diff --git a/core/include/traccc/fitting/kalman_filter/measurement_selector.hpp b/core/include/traccc/fitting/kalman_filter/measurement_selector.hpp new file mode 100644 index 0000000000..bfb5fcdeee --- /dev/null +++ b/core/include/traccc/fitting/kalman_filter/measurement_selector.hpp @@ -0,0 +1,206 @@ +/** TRACCC library, part of the ACTS project (R&D line) + * + * (c) 2025 CERN for the benefit of the ACTS project + * + * Mozilla Public License Version 2.0 + */ + +#pragma once + +// Project include(s). +#include "traccc/definitions/qualifiers.hpp" +#include "traccc/definitions/track_parametrization.hpp" +#include "traccc/edm/measurement_collection.hpp" +#include "traccc/edm/measurement_helpers.hpp" +#include "traccc/utils/logging.hpp" +#include "traccc/utils/subspace.hpp" + +namespace traccc { + +/// Associate a measurement to a candidate track +struct measurement_selector { + + template + using matrix_t = detray::dmatrix; + + // Where to get the calibration from + struct config { /*TODO: implement calibration handling*/ + }; + + /// Get the obersavtion model for a given measurement + /// + /// @brief Based on "Application of Kalman filtering to track and vertex + /// fitting", R.Frühwirth, NIM A + /// + /// @param measurement the measurement + /// @param is_line whether the measurement belong to a line surface + /// + /// @returns the projection matrix H + template + TRACCC_HOST_DEVICE static constexpr detray::dmatrix + observation_model( + const edm::measurement& measurement, + const bound_track_parameters& bound_params, + const bool is_line) { + + // Oservation model: Subspace of measurment space for this measurement + const subspace subs(measurement.subspace()); + detray::dmatrix H = + subs.template projector(); + + // Flip the sign of projector matrix element in case the first element + // of line measurement is negative + if (is_line && bound_params.bound_local()[e_bound_loc0] < 0) { + getter::element(H, 0u, e_bound_loc0) = -1; + } + + if (measurement.dimensions() == 1) { + getter::element(H, 1u, 0u) = 0.f; + getter::element(H, 1u, 1u) = 0.f; + } + + TRACCC_DEBUG_HOST("Observation model (H):\n" << H); + + return H; + } + + /// Get the calibrated measurement position + /// + /// @brief Based on "Application of Kalman filtering to track and vertex + /// fitting", R.Frühwirth, NIM A + /// + /// @param measurement the measurement + /// @param is_line whether the measurement belong to a line surface + /// @param cfg how to apply calibrations + /// + /// @returns the projection matrix H + template + TRACCC_HOST_DEVICE static constexpr detray::dmatrix + calibrated_measurement_position( + const edm::measurement& measurement, + const config& /*cfg*/) { + + // Measurement data on surface + detray::dmatrix meas_local; + edm::get_measurement_local(measurement, meas_local); + + TRACCC_DEBUG_HOST( + "Measurement position (uncalibrated): " << meas_local); + + assert((measurement.dimensions() > 1) || + (getter::element(meas_local, 1u, 0u) == 0.f)); + + return meas_local; + } + + /// Get the calibrated measurement covariance + /// + /// @brief Based on "Application of Kalman filtering to track and vertex + /// fitting", R.Frühwirth, NIM A + /// + /// @param measurement the measurement + /// @param is_line whether the measurement belong to a line surface + /// @param cfg how to apply calibrations + /// + /// @returns the projection matrix H + template + TRACCC_HOST_DEVICE static constexpr detray::dmatrix + calibrated_measurement_covariance( + const edm::measurement& measurement, + const config& /*cfg*/) { + + // Measurement covariance + detray::dmatrix V; + edm::get_measurement_covariance(measurement, V); + + if (measurement.dimensions() == 1) { + getter::element(V, 1u, 1u) = + std::numeric_limits>::max(); + } + + TRACCC_DEBUG_HOST("Measurement covariance (uncalibrated):\n" << V); + + return V; + } + + /// Caluculate the predicted chi2 + /// + /// @brief Based on "Application of Kalman filtering to track and vertex + /// fitting", R.Frühwirth, NIM A + /// + /// @param measurement the measurement + /// @param bound_params bound track parameters (state vector, covariance) + /// @param is_line whether the measurement belong to a line surface + /// + /// @returns the predicted chi2 + template + TRACCC_HOST_DEVICE static constexpr detray::dscalar + predicted_chi2(const edm::measurement& measurement, + const bound_track_parameters& bound_params, + const config& cfg, const bool is_line) { + + // Measurement maximal dimension + constexpr unsigned int D = 2; + + TRACCC_VERBOSE_HOST_DEVICE("In measurement selector..."); + TRACCC_VERBOSE_HOST_DEVICE("Measurement dim: %d", + measurement.dimensions()); + + assert(measurement.dimensions() == 1u || + measurement.dimensions() == 2u); + + TRACCC_DEBUG_HOST("Predicted param.: " << bound_params); + + assert(!bound_params.is_invalid()); + assert(!bound_params.surface_link().is_invalid()); + + // Predicted vector and covariance of bound track parameters + const traccc::bound_vector& predicted_vec = + bound_params.vector(); + const traccc::bound_matrix& predicted_cov = + bound_params.covariance(); + + const matrix_t meas_local = + calibrated_measurement_position(measurement, cfg); + + // Measurement covariance + const matrix_t V = + calibrated_measurement_covariance(measurement, cfg); + + // Get the projection matrix for this measurement (obeservation model) + const matrix_t H = + observation_model(measurement, bound_params, is_line); + + // Project the predicted covariance to the observation + const matrix_t R = + H * algebra::matrix::transposed_product(predicted_cov, + H) + + V; + + TRACCC_DEBUG_HOST("R:\n" << R); + TRACCC_DEBUG_HOST_DEVICE("det(R): %f", matrix::determinant(R)); + TRACCC_DEBUG_HOST("R_inv:\n" << matrix::inverse(R)); + + // Residual between measurement and (projected) vector (innovation) + const matrix_t residual = + meas_local - H * predicted_vec; + + TRACCC_VERBOSE_HOST("Predicted residual: " << residual); + + const matrix_t chi2_mat = + algebra::matrix::transposed_product( + residual, matrix::inverse(R)) * + residual; + + TRACCC_VERBOSE_HOST_DEVICE("Chi2: %f", getter::element(chi2_mat, 0, 0)); + + return getter::element(chi2_mat, 0, 0); + } +}; + +} // namespace traccc diff --git a/core/include/traccc/fitting/kalman_filter/two_filters_smoother.hpp b/core/include/traccc/fitting/kalman_filter/two_filters_smoother.hpp index 634a812ae7..eb45104e73 100644 --- a/core/include/traccc/fitting/kalman_filter/two_filters_smoother.hpp +++ b/core/include/traccc/fitting/kalman_filter/two_filters_smoother.hpp @@ -36,18 +36,16 @@ struct two_filters_smoother { /// @param bound_params bound parameter /// /// @return true if the update succeeds + template [[nodiscard]] TRACCC_HOST_DEVICE inline kalman_fitter_status operator()( - typename edm::track_state_collection::device::proxy_type& - trk_state, - const typename edm::measurement_collection::const_device& - measurements, + typename edm::track_state& trk_state, + const edm::measurement& measurement, bound_track_parameters& bound_params, const bool is_line) const { static constexpr unsigned int D = 2; - [[maybe_unused]] const unsigned int dim{ - measurements.at(trk_state.measurement_index()).dimensions()}; + [[maybe_unused]] const unsigned int dim{measurement.dimensions()}; assert(dim == 1u || dim == 2u); @@ -65,8 +63,7 @@ struct two_filters_smoother { // Measurement data on surface matrix_type meas_local; - edm::get_measurement_local( - measurements.at(trk_state.measurement_index()), meas_local); + edm::get_measurement_local(measurement, meas_local); assert((dim > 1) || (getter::element(meas_local, 1u, 0u) == 0.f)); @@ -132,8 +129,7 @@ struct two_filters_smoother { // Wrap the phi and theta angles in their valid ranges normalize_angles(trk_state.smoothed_params()); - const subspace subs( - measurements.at(trk_state.measurement_index()).subspace()); + const subspace subs(measurement.subspace()); matrix_type H = subs.template projector(); // @TODO: Fix properly if (getter::element(meas_local, 1u, 0u) == 0.f /*dim == 1*/) { @@ -145,8 +141,7 @@ struct two_filters_smoother { // Spatial resolution (Measurement covariance) matrix_type V; - edm::get_measurement_covariance( - measurements.at(trk_state.measurement_index()), V); + edm::get_measurement_covariance(measurement, V); // @TODO: Fix properly if (getter::element(meas_local, 1u, 0u) == 0.f /*dim == 1*/) { getter::element(V, 1u, 1u) = 1000.f; diff --git a/device/common/include/traccc/finding/device/impl/find_tracks.ipp b/device/common/include/traccc/finding/device/impl/find_tracks.ipp index e8f8397ebe..594117e166 100644 --- a/device/common/include/traccc/finding/device/impl/find_tracks.ipp +++ b/device/common/include/traccc/finding/device/impl/find_tracks.ipp @@ -207,58 +207,57 @@ TRACCC_HOST_DEVICE inline void find_tracks( if (use_measurement) { + const detray::tracking_surface sf{det, in_par.surface_link()}; + const bool is_line = detail::is_line(sf); + + measurement_selector::config calib_cfg{}; + + const edm::measurement meas = measurements.at(meas_idx); + + // The filtered track state auto trk_state = edm::make_track_state(measurements, meas_idx); - const detray::tracking_surface sf{det, in_par.surface_link()}; + const traccc::scalar chi2 = + measurement_selector::predicted_chi2(meas, in_par, + calib_cfg, is_line); + + trk_state.filtered_chi2() = chi2; + + // If the measurement is outside the chi2 cut, skip it + if (chi2 > cfg.chi2_max) { + use_measurement = false; + } + // Kalman filter status code kalman_fitter_status res{kalman_fitter_status::ERROR_OTHER}; - // Check if the measurement is even remotely close to the track - if (payload.step == 0 && !sf.has_material()) { - const auto& meas = measurements.at(meas_idx); - - detray::dmatrix meas_local; - edm::get_measurement_local(meas, meas_local); - - const subspace subs( - meas.subspace()); - detray::dmatrix H = - subs.template projector<2>(); - - const auto residual = - meas_local - - H * in_params.at(owner_global_thread_id).vector(); - - const traccc::scalar res_0{getter::element(residual, 0, 0)}; - const traccc::scalar res_1{ - meas.dimensions() == 1 - ? 0.f - : getter::element(residual, 1, 0)}; - const traccc::scalar search_rad = - math::sqrt(res_0 * res_0 + res_1 * res_1); - if (search_rad > 0.f) { - use_measurement = false; + + if (use_measurement) { + if (payload.step == 0 && chi2 == 0.f && + !sf.has_material()) { + res = kalman_fitter_status::SUCCESS; + + trk_state.filtered_params() = in_par; + + // Update measurement covariance + const auto V = measurement_selector:: + calibrated_measurement_covariance( + meas, calib_cfg); + + auto& filtered_cov = + trk_state.filtered_params().covariance(); + getter::element(filtered_cov, e_bound_loc0, + e_bound_loc0) = + getter::element(V, 0, 0); + getter::element(filtered_cov, e_bound_loc1, + e_bound_loc1) = + getter::element(V, 1, 1); + } else { + // Run the Kalman update on a copy of the track + // parameters + res = gain_matrix_updater{}(trk_state, meas, + in_par, is_line); } - res = kalman_fitter_status::SUCCESS; - trk_state.filtered_params() = - in_params.at(owner_global_thread_id); - trk_state.filtered_chi2() = 0.f; - - // Set the covariance to the measurement variance - detray::dmatrix V; - edm::get_measurement_covariance(meas, V); - auto& filtered_cov = - trk_state.filtered_params().covariance(); - getter::element(filtered_cov, e_bound_loc0, e_bound_loc0) = - getter::element(V, 0, 0); - getter::element(filtered_cov, e_bound_loc1, e_bound_loc1) = - getter::element(V, 1, 1); - } else { - const bool is_line = detail::is_line(sf); - - // Run the Kalman update on a copy of the track parameters - res = gain_matrix_updater{}( - trk_state, measurements, in_par, is_line); } TRACCC_DEBUG_DEVICE("KF status: %d", res); @@ -274,13 +273,7 @@ TRACCC_HOST_DEVICE inline void find_tracks( * but, more importantly, allows us to more easily use * block-wide synchronization primitives. */ - if (const traccc::scalar chi2 = trk_state.filtered_chi2(); - res != kalman_fitter_status::SUCCESS || - chi2 >= cfg.chi2_max) { - use_measurement = false; - } - - if (use_measurement) { + if (use_measurement && res == kalman_fitter_status::SUCCESS) { TRACCC_VERBOSE_DEVICE("Found measurement: %d", meas_idx); result.emplace(trk_state, meas_idx, owner_local_thread_id); }