diff --git a/core/include/traccc/edm/track_parameters.hpp b/core/include/traccc/edm/track_parameters.hpp index 0a33dc3ffe..0f5ec78187 100644 --- a/core/include/traccc/edm/track_parameters.hpp +++ b/core/include/traccc/edm/track_parameters.hpp @@ -14,6 +14,7 @@ #include "traccc/definitions/qualifiers.hpp" #include "traccc/definitions/track_parametrization.hpp" #include "traccc/edm/container.hpp" +#include "traccc/utils/trigonometric_helpers.hpp" // detray include(s). #include @@ -48,19 +49,15 @@ using bound_track_parameters_collection_types = // Wrap the phi of track parameters to [-pi,pi] template -TRACCC_HOST_DEVICE inline void wrap_phi( +TRACCC_HOST_DEVICE constexpr void normalize_angles( bound_track_parameters& param) { + traccc::scalar phi; + traccc::scalar theta; + + std::tie(phi, theta) = detail::wrap_phi_theta(param.phi(), param.theta()); - traccc::scalar phi = param.phi(); - static constexpr traccc::scalar TWOPI = - 2.f * traccc::constant::pi; - phi = math::fmod(phi, TWOPI); - if (phi > traccc::constant::pi) { - phi -= TWOPI; - } else if (phi < -traccc::constant::pi) { - phi += TWOPI; - } param.set_phi(phi); + param.set_theta(theta); } /// Covariance inflation used for track fitting diff --git a/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp b/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp index 27586ba155..aaf0e2cd66 100644 --- a/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp +++ b/core/include/traccc/finding/details/combinatorial_kalman_filter.hpp @@ -126,8 +126,10 @@ combinatorial_kalman_filter( std::vector> tips; // Create propagator + auto prop_cfg{config.propagation}; + prop_cfg.navigation.estimate_scattering_noise = false; traccc::details::ckf_propagator_t propagator( - config.propagation); + prop_cfg); // Create the input seeds container. bound_track_parameters_collection_types::const_device seeds{seeds_view}; @@ -240,7 +242,6 @@ combinatorial_kalman_filter( // Iterate over the measurements for (unsigned int meas_id = lo; meas_id < up; meas_id++) { - // The measurement on surface to handle. const measurement& meas = measurements.at(meas_id); @@ -248,7 +249,7 @@ combinatorial_kalman_filter( auto trk_state = edm::make_track_state(measurements, meas_id); - const bool is_line = sf.template visit_mask(); + const bool is_line = detail::is_line(sf); // Run the Kalman update on a copy of the track parameters const kalman_fitter_status res = @@ -484,8 +485,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_cfg}; typename detray::momentum_aborter::state s4{}; typename ckf_aborter::state s5; // Update the actor config @@ -495,21 +496,41 @@ combinatorial_kalman_filter( s5.max_count = config.max_step_counts_for_next_surface; // 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)); // If a surface found, add the parameter for the next // step - if (s5.success) { + bool valid_track{s5.success}; + if (valid_track) { assert(propagation._navigation.is_on_sensitive()); assert(!propagation._stepping.bound_params().is_invalid()); - out_params.push_back(propagation._stepping.bound_params()); - param_to_link[step].push_back(link_id); + const auto& out_param = propagation._stepping.bound_params(); + + const scalar theta = out_param.theta(); + if (theta <= 0.f || + theta >= 2.f * constant::pi) { + valid_track = false; + } + + if (!std::isfinite(out_param.phi())) { + valid_track = false; + } + + if (math::fabs(out_param.qop()) == 0.f) { + valid_track = false; + } + + if (valid_track) { + out_params.push_back(out_param); + 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))) { + if (!valid_track && + (step >= (config.min_track_candidates_per_track - 1u))) { tips.push_back({step, link_id}); } diff --git a/core/include/traccc/finding/finding_config.hpp b/core/include/traccc/finding/finding_config.hpp index 97f7ef9639..4c1a7d8337 100644 --- a/core/include/traccc/finding/finding_config.hpp +++ b/core/include/traccc/finding/finding_config.hpp @@ -43,7 +43,7 @@ struct finding_config { unsigned int max_step_counts_for_next_surface = 100; /// Maximum Chi-square that is allowed for branching - float chi2_max = 100.f; + float chi2_max = 30.f; /// Propagation configuration detray::propagation::config propagation{}; diff --git a/core/include/traccc/fitting/details/kalman_fitting.hpp b/core/include/traccc/fitting/details/kalman_fitting.hpp index 2c0eab17bb..435f5b7ab1 100644 --- a/core/include/traccc/fitting/details/kalman_fitting.hpp +++ b/core/include/traccc/fitting/details/kalman_fitting.hpp @@ -88,7 +88,7 @@ typename edm::track_container::host kalman_fitting( result_tracks_device.at(result_tracks_device.size() - 1), typename edm::track_state_collection::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 diff --git a/core/include/traccc/fitting/fitting_config.hpp b/core/include/traccc/fitting/fitting_config.hpp index 9d956203a9..1e30ca0298 100644 --- a/core/include/traccc/fitting/fitting_config.hpp +++ b/core/include/traccc/fitting/fitting_config.hpp @@ -34,7 +34,7 @@ struct fitting_config { std::size_t barcode_sequence_size_factor = 5; std::size_t min_barcode_sequence_capacity = 100; traccc::scalar backward_filter_mask_tolerance = - 5.f * traccc::unit::mm; + 10.f * traccc::unit::m; }; } // namespace traccc 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 c96b2d4f76..16631f0b0a 100644 --- a/core/include/traccc/fitting/kalman_filter/gain_matrix_updater.hpp +++ b/core/include/traccc/fitting/kalman_filter/gain_matrix_updater.hpp @@ -123,6 +123,30 @@ struct gain_matrix_updater { predicted_vec + K * (meas_local - H * predicted_vec); const matrix_type<6, 6> filtered_cov = (I66 - K * H) * predicted_cov; + // Check the covariance for consistency + if (getter::element(filtered_cov, 0, 0) <= 0.f || + getter::element(filtered_cov, 1, 1) <= 0.f || + getter::element(filtered_cov, 2, 2) <= 0.f || + getter::element(filtered_cov, 3, 3) <= 0.f || + getter::element(filtered_cov, 4, 4) <= 0.f || + getter::element(filtered_cov, 5, 5) <= 0.f) { + return kalman_fitter_status::ERROR_UPDATER_INVALID_COVARIANCE; + } + + // Return false if track is parallel to z-axis or phi is not finite + if (!std::isfinite(getter::element(filtered_vec, e_bound_theta, 0))) { + return kalman_fitter_status::ERROR_INVERSION; + } + + if (!std::isfinite(getter::element(filtered_vec, e_bound_phi, 0))) { + return kalman_fitter_status::ERROR_INVERSION; + } + + if (math::fabs(getter::element(filtered_vec, e_bound_qoverp, 0)) == + 0.f) { + return kalman_fitter_status::ERROR_QOP_ZERO; + } + // Residual between measurement and (projected) filtered vector const matrix_type residual = meas_local - H * filtered_vec; @@ -133,36 +157,28 @@ struct gain_matrix_updater { residual, matrix::inverse(R)) * residual; - // 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::pi) { - return kalman_fitter_status::ERROR_THETA_ZERO; - } - - if (!std::isfinite(bound_params.phi())) { - return kalman_fitter_status::ERROR_INVERSION; - } + const scalar chi2_val{getter::element(chi2, 0, 0)}; - if (std::abs(bound_params.qop()) == 0.f) { - return kalman_fitter_status::ERROR_QOP_ZERO; - } - - if (getter::element(chi2, 0, 0) < 0.f) { + if (chi2_val < 0.f) { return kalman_fitter_status::ERROR_UPDATER_CHI2_NEGATIVE; } - if (!std::isfinite(getter::element(chi2, 0, 0))) { + if (!std::isfinite(chi2_val)) { return kalman_fitter_status::ERROR_UPDATER_CHI2_NOT_FINITE; } - // Set the track state parameters + // Set the chi2 for this track and measurement trk_state.filtered_params().set_vector(filtered_vec); trk_state.filtered_params().set_covariance(filtered_cov); - trk_state.filtered_chi2() = getter::element(chi2, 0, 0); + trk_state.filtered_chi2() = chi2_val; - // Wrap the phi in the range of [-pi, pi] - wrap_phi(trk_state.filtered_params()); + // Wrap the phi and theta angles in their valid ranges + normalize_angles(trk_state.filtered_params()); + + const scalar theta = trk_state.filtered_params().theta(); + if (theta <= 0.f || theta >= 2.f * constant::pi) { + return kalman_fitter_status::ERROR_THETA_POLE; + } assert(!trk_state.filtered_params().is_invalid()); diff --git a/core/include/traccc/fitting/kalman_filter/is_line_visitor.hpp b/core/include/traccc/fitting/kalman_filter/is_line_visitor.hpp index 34c3d5a14a..a4846c5774 100644 --- a/core/include/traccc/fitting/kalman_filter/is_line_visitor.hpp +++ b/core/include/traccc/fitting/kalman_filter/is_line_visitor.hpp @@ -7,20 +7,31 @@ #pragma once +#include #include +#include #include "traccc/definitions/qualifiers.hpp" -namespace traccc { +namespace traccc::detail { -struct is_line_visitor { - template - [[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> || - std::is_same_v>; +/// @returns true if the surface has "line" shape +template +[[nodiscard]] TRACCC_HOST_DEVICE bool constexpr is_line( + const detray::geometry::surface sf) { + using algebra_t = typename detector_t::algebra_type; + using straw_tube = detray::mask, algebra_t>; + using wire_cell = detray::mask, algebra_t>; + + if constexpr (detector_t::masks::template is_defined() || + detector_t::masks::template is_defined()) { + return (sf.shape_id() == + detector_t::masks::template get_id()) || + (sf.shape_id() == + detector_t::masks::template get_id()); + } else { + return false; } }; -} // namespace traccc +} // namespace traccc::detail diff --git a/core/include/traccc/fitting/kalman_filter/kalman_actor.hpp b/core/include/traccc/fitting/kalman_filter/kalman_actor.hpp index 94b70c1f91..9d5f995192 100644 --- a/core/include/traccc/fitting/kalman_filter/kalman_actor.hpp +++ b/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(); + const bool is_line = detail::is_line(sf); kalman_fitter_status res = kalman_fitter_status::SUCCESS; @@ -173,6 +173,9 @@ struct kalman_actor : detray::actor { 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()); + // Forward filter res = gain_matrix_updater{}( trk_state, actor_state.m_measurements, @@ -215,7 +218,9 @@ struct kalman_actor : detray::actor { actor_state.next(); // Flag renavigation of the current candidate - navigation.set_high_trust(); + if (math::fabs(navigation()) > 1.f * unit::um) { + navigation.set_high_trust(); + } } } }; diff --git a/core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp b/core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp index 06eddb4d85..066e8c6f19 100644 --- a/core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp +++ b/core/include/traccc/fitting/kalman_filter/kalman_fitter.hpp @@ -55,6 +55,7 @@ class kalman_fitter { // Actor types using aborter = detray::pathlimit_aborter; + using momentum_aborter = detray::momentum_aborter; using transporter = detray::parameter_transporter; using interactor = detray::pointwise_material_interactor; using forward_fit_actor = @@ -71,11 +72,13 @@ class kalman_fitter { using forward_actor_chain_type = detray::actor_chain; + resetter, momentum_aborter, barcode_sequencer, + kalman_step_aborter>; using backward_actor_chain_type = detray::actor_chain; + interactor, resetter, momentum_aborter, + kalman_step_aborter>; // Navigator type for backward propagator using direct_navigator_type = detray::direct_navigator; @@ -110,11 +113,13 @@ class kalman_fitter { track_states, const measurement_collection_types::const_device& measurements, vecmem::data::vector_view - sequence_buffer) + sequence_buffer, + const detray::propagation::config& prop_cfg) : m_fit_actor_state{track, track_states, measurements}, m_sequencer_state( vecmem::device_vector( sequence_buffer)), + m_parameter_resetter{prop_cfg}, m_fit_res{track}, m_sequence_buffer(sequence_buffer) {} @@ -122,7 +127,8 @@ class kalman_fitter { 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_fit_actor_state, m_parameter_resetter, + m_montum_aborter_state, m_sequencer_state, m_step_aborter_state); } @@ -131,7 +137,8 @@ class kalman_fitter { 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_montum_aborter_state, m_step_aborter_state); } /// Individual actor states @@ -140,6 +147,8 @@ 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{}; + typename momentum_aborter::state m_montum_aborter_state{}; /// Fitting result per track typename edm::track_collection::device::proxy_type @@ -277,16 +286,20 @@ class kalman_fitter { } auto last = fitter_state.m_fit_actor_state(); - const scalar theta = last.filtered_params().theta(); - if (theta <= 0.f || theta >= constant::pi) { - return kalman_fitter_status::ERROR_THETA_ZERO; + + if (!std::isfinite(last.filtered_params().theta())) { + return kalman_fitter_status::ERROR_INVERSION; } if (!std::isfinite(last.filtered_params().phi())) { return kalman_fitter_status::ERROR_INVERSION; } + if (theta <= 0.f || theta >= 2.f * constant::pi) { + return kalman_fitter_status::ERROR_THETA_POLE; + } + last.smoothed_params().set_parameter_vector(last.filtered_params()); last.smoothed_params().set_covariance( last.filtered_params().covariance()); @@ -324,6 +337,8 @@ class kalman_fitter { propagation._navigation.set_direction( detray::navigation::direction::e_backward); + propagation._navigation.safe_step_size = + 0.1f * traccc::unit::mm; // Synchronize the current barcode with the input track parameter while (propagation._navigation.get_target_barcode() != 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 81dcbe44d8..5c9c1941e9 100644 --- a/core/include/traccc/fitting/kalman_filter/two_filters_smoother.hpp +++ b/core/include/traccc/fitting/kalman_filter/two_filters_smoother.hpp @@ -99,6 +99,16 @@ struct two_filters_smoother { const matrix_type smoothed_cov = matrix::inverse(smoothed_cov_inv); + // Check the covariance for consistency + if (getter::element(smoothed_cov, 0, 0) <= 0.f || + getter::element(smoothed_cov, 1, 1) <= 0.f || + getter::element(smoothed_cov, 2, 2) <= 0.f || + getter::element(smoothed_cov, 3, 3) <= 0.f || + getter::element(smoothed_cov, 4, 4) <= 0.f || + getter::element(smoothed_cov, 5, 5) <= 0.f) { + return kalman_fitter_status::ERROR_SMOOTHER_INVALID_COVARIANCE; + } + // Eq (3.38) of "Pattern Recognition, Tracking and Vertex // Reconstruction in Particle Detectors" const matrix_type smoothed_vec = @@ -107,8 +117,25 @@ struct two_filters_smoother { predicted_cov_inv * predicted_vec); trk_state.smoothed_params().set_vector(smoothed_vec); + + // Return false if track is parallel to z-axis or phi is not finite + if (!std::isfinite(trk_state.smoothed_params().theta())) { + return kalman_fitter_status::ERROR_INVERSION; + } + + if (!std::isfinite(trk_state.smoothed_params().phi())) { + return kalman_fitter_status::ERROR_INVERSION; + } + + if (math::fabs(trk_state.smoothed_params().qop()) == 0.f) { + return kalman_fitter_status::ERROR_QOP_ZERO; + } + trk_state.smoothed_params().set_covariance(smoothed_cov); + // Wrap the phi and theta angles in their valid ranges + normalize_angles(trk_state.smoothed_params()); + matrix_type H = measurements.at(trk_state.measurement_index()) .subs.template projector(); @@ -141,11 +168,13 @@ struct two_filters_smoother { residual_smt, matrix::inverse(R_smt)) * residual_smt; - if (getter::element(chi2_smt, 0, 0) < 0.f) { + const scalar chi2_smt_value{getter::element(chi2_smt, 0, 0)}; + + if (chi2_smt_value < 0.f) { return kalman_fitter_status::ERROR_SMOOTHER_CHI2_NEGATIVE; } - if (!std::isfinite(getter::element(chi2_smt, 0, 0))) { + if (!std::isfinite(chi2_smt_value)) { return kalman_fitter_status::ERROR_SMOOTHER_CHI2_NOT_FINITE; } @@ -180,49 +209,66 @@ struct two_filters_smoother { predicted_vec + K * (meas_local - H * predicted_vec); const matrix_type<6, 6> filtered_cov = (I66 - K * H) * predicted_cov; - // Residual between measurement and (projected) filtered vector - const matrix_type residual = meas_local - H * filtered_vec; - - // Calculate backward chi2 - const matrix_type R = (I_m - H * K) * V; - // assert(matrix::determinant(R) != 0.f); - assert(std::isfinite(matrix::determinant(R))); - const matrix_type<1, 1> chi2 = - algebra::matrix::transposed_product( - residual, matrix::inverse(R)) * - residual; + // Check the covariance for consistency + if (getter::element(filtered_cov, 0, 0) <= 0.f || + getter::element(filtered_cov, 1, 1) <= 0.f || + getter::element(filtered_cov, 2, 2) <= 0.f || + getter::element(filtered_cov, 3, 3) <= 0.f || + getter::element(filtered_cov, 4, 4) <= 0.f || + getter::element(filtered_cov, 5, 5) <= 0.f) { + return kalman_fitter_status::ERROR_SMOOTHER_INVALID_COVARIANCE; + } // 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::pi) { - return kalman_fitter_status::ERROR_THETA_ZERO; + if (!std::isfinite(bound_params.theta())) { + return kalman_fitter_status::ERROR_INVERSION; } if (!std::isfinite(bound_params.phi())) { return kalman_fitter_status::ERROR_INVERSION; } - if (std::abs(bound_params.qop()) == 0.f) { + if (math::fabs(bound_params.qop()) == 0.f) { return kalman_fitter_status::ERROR_QOP_ZERO; } - if (getter::element(chi2, 0, 0) < 0.f) { - return kalman_fitter_status::ERROR_UPDATER_CHI2_NEGATIVE; + bound_params.set_covariance(filtered_cov); + + // Residual between measurement and (projected) filtered vector + const matrix_type residual = meas_local - H * filtered_vec; + + // Calculate backward chi2 + const matrix_type R = (I_m - H * K) * V; + // assert(matrix::determinant(R) != 0.f); // @TODO: This fails + assert(std::isfinite(matrix::determinant(R))); + const matrix_type<1, 1> chi2 = + algebra::matrix::transposed_product( + residual, matrix::inverse(R)) * + residual; + + const scalar chi2_val{getter::element(chi2, 0, 0)}; + + if (chi2_val < 0.f) { + return kalman_fitter_status::ERROR_SMOOTHER_CHI2_NEGATIVE; } - if (!std::isfinite(getter::element(chi2, 0, 0))) { - return kalman_fitter_status::ERROR_UPDATER_CHI2_NOT_FINITE; + if (!std::isfinite(chi2_val)) { + return kalman_fitter_status::ERROR_SMOOTHER_CHI2_NOT_FINITE; } // Set backward chi2 - trk_state.backward_chi2() = getter::element(chi2, 0, 0); + trk_state.backward_chi2() = chi2_val; - // Wrap the phi in the range of [-pi, pi] - wrap_phi(bound_params); + // Wrap the phi and theta angles in their valid ranges + normalize_angles(bound_params); + + const scalar theta = bound_params.theta(); + if (theta <= 0.f || theta >= 2.f * constant::pi) { + return kalman_fitter_status::ERROR_THETA_POLE; + } trk_state.set_smoothed(); diff --git a/core/include/traccc/fitting/status_codes.hpp b/core/include/traccc/fitting/status_codes.hpp index ae7ce138af..05962e5add 100644 --- a/core/include/traccc/fitting/status_codes.hpp +++ b/core/include/traccc/fitting/status_codes.hpp @@ -14,12 +14,14 @@ namespace traccc { enum class kalman_fitter_status : uint32_t { SUCCESS, ERROR_QOP_ZERO, - ERROR_THETA_ZERO, + ERROR_THETA_POLE, ERROR_INVERSION, ERROR_SMOOTHER_CHI2_NEGATIVE, ERROR_SMOOTHER_CHI2_NOT_FINITE, + ERROR_SMOOTHER_INVALID_COVARIANCE, ERROR_UPDATER_CHI2_NEGATIVE, ERROR_UPDATER_CHI2_NOT_FINITE, + ERROR_UPDATER_INVALID_COVARIANCE, ERROR_BARCODE_SEQUENCE_OVERFLOW, ERROR_INVALID_TRACK_STATE, ERROR_OTHER, @@ -36,8 +38,8 @@ struct fitter_debug_msg { case ERROR_QOP_ZERO: { return msg + "Track qop is zero"; } - case ERROR_THETA_ZERO: { - return msg + "Track theta is zero"; + case ERROR_THETA_POLE: { + return msg + "Track theta hit pole"; } case ERROR_INVERSION: { return msg + "Failed matrix inversion"; @@ -49,10 +51,10 @@ struct fitter_debug_msg { return msg + "Invalid chi2 in smoother"; } case ERROR_UPDATER_CHI2_NEGATIVE: { - return msg + "Negative chi2 in gain matrix update"; + return msg + "Negative chi2 in gain matrix updater"; } case ERROR_UPDATER_CHI2_NOT_FINITE: { - return msg + "Invalid chi2 in gain matrix update"; + return msg + "Invalid chi2 in gain matrix updater"; } case ERROR_BARCODE_SEQUENCE_OVERFLOW: { return msg + "Barcode sequence overflow in direct navigator"; @@ -61,6 +63,12 @@ struct fitter_debug_msg { return msg + "Invalid track state in forward pass (skipped or error)"; } + case ERROR_UPDATER_INVALID_COVARIANCE: { + return msg + "Invalid track covariance after filtering"; + } + case ERROR_SMOOTHER_INVALID_COVARIANCE: { + return msg + "Invalid track covariance after smoothing"; + } case ERROR_OTHER: { return msg + "Unspecified error"; } diff --git a/core/include/traccc/seeding/doublet_finding_helper.hpp b/core/include/traccc/seeding/doublet_finding_helper.hpp index 2135cb2323..089955ed92 100644 --- a/core/include/traccc/seeding/doublet_finding_helper.hpp +++ b/core/include/traccc/seeding/doublet_finding_helper.hpp @@ -77,7 +77,7 @@ bool TRACCC_HOST_DEVICE doublet_finding_helper::isCompatible( (math::fabs(cotTheta) >= config.cotThetaMax * deltaR) || (zOrigin <= config.collisionRegionMin * deltaR) || (zOrigin >= config.collisionRegionMax * deltaR) || - std::abs(cotTheta) >= config.deltaZMax) { + math::fabs(cotTheta) >= config.deltaZMax) { return false; } diff --git a/core/include/traccc/utils/trigonometric_helpers.hpp b/core/include/traccc/utils/trigonometric_helpers.hpp new file mode 100644 index 0000000000..e9a15b4a27 --- /dev/null +++ b/core/include/traccc/utils/trigonometric_helpers.hpp @@ -0,0 +1,83 @@ +/** 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/common.hpp" +#include "traccc/definitions/math.hpp" +#include "traccc/definitions/primitives.hpp" +#include "traccc/definitions/qualifiers.hpp" + +/// @see +/// https://github.com/acts-project/acts/blob/8098e6953ac35771c34a4e3b13dbfee50869a0c1/Core/include/Acts/Utilities/detail/periodic.hpp +namespace traccc::detail { + +/// Wrap a periodic value back into the nominal range. +TRACCC_HOST_DEVICE +constexpr traccc::scalar wrap_periodic(traccc::scalar value, + traccc::scalar start, + traccc::scalar range) { + // only wrap if really necessary + const traccc::scalar diff{value - start}; + return ((0 <= diff) && (diff < range)) + ? value + : (value - range * math::floor(diff / range)); +} + +/// Calculate the equivalent angle in the [-pi, pi) range. +TRACCC_HOST_DEVICE +constexpr traccc::scalar wrap_phi(traccc::scalar phi) { + constexpr traccc::scalar PI{traccc::constant::pi}; + constexpr traccc::scalar TWOPI{2.f * traccc::constant::pi}; + return wrap_periodic(phi, -PI, TWOPI); +} + +/// Calculate the equivalent angle in the [0, 2*pi) range. +TRACCC_HOST_DEVICE +constexpr traccc::scalar wrap_theta(traccc::scalar theta) { + constexpr traccc::scalar TWOPI{2.f * traccc::constant::pi}; + return wrap_periodic(theta, 0.f, TWOPI); +} + +/// Ensure both phi and theta direction angles are within the allowed range. +/// +/// @param[in] phi Transverse direction angle +/// @param[in] theta Longitudinal direction angle +/// @return pair containing the updated angles +/// +/// The phi angle is truly cyclic, i.e. all values outside the nominal range +/// [-pi,pi) have a corresponding value inside nominal range, independent from +/// the theta angle. The theta angle is more complicated. Imagine that the two +/// angles describe a position on the unit sphere. If theta moves outside its +/// nominal range [0,pi], we are moving over one of the two poles of the unit +/// sphere along the great circle defined by phi. The angles still describe a +/// valid position on the unit sphere, but to describe it with angles within +/// their nominal range, both phi and theta need to be updated; when moving over +/// the poles, phi needs to be flipped by 180degree to allow theta to remain +/// within its nominal range. +constexpr std::pair wrap_phi_theta( + traccc::scalar phi, traccc::scalar theta) { + constexpr traccc::scalar PI{traccc::constant::pi}; + constexpr traccc::scalar TWOPI{2.f * traccc::constant::pi}; + + // wrap to [0,2pi). while the nominal range of theta is [0,pi], it is + // periodic, i.e. describes identical positions, in the full [0,2pi) range. + // moving it first to the periodic range simplifies further steps as the + // possible range of theta becomes fixed. + theta = wrap_theta(theta); + if (PI < theta) { + // theta is in the second half of the great circle and outside its + // nominal range. need to change both phi and theta to be within range. + phi += PI; + theta = TWOPI - theta; + } + + return {wrap_phi(phi), theta}; +} + +} // namespace traccc::detail diff --git a/device/alpaka/src/finding/combinatorial_kalman_filter.hpp b/device/alpaka/src/finding/combinatorial_kalman_filter.hpp index b5176c0379..68bc6fe93d 100644 --- a/device/alpaka/src/finding/combinatorial_kalman_filter.hpp +++ b/device/alpaka/src/finding/combinatorial_kalman_filter.hpp @@ -28,7 +28,6 @@ #include "traccc/finding/device/propagate_to_next_surface.hpp" #include "traccc/finding/device/remove_duplicates.hpp" #include "traccc/finding/finding_config.hpp" -#include "traccc/utils/logging.hpp" #include "traccc/utils/memory_resource.hpp" #include "traccc/utils/projections.hpp" #include "traccc/utils/propagation.hpp" 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 52ca4aec37..c4ea2f0d10 100644 --- a/device/common/include/traccc/finding/device/impl/find_tracks.ipp +++ b/device/common/include/traccc/finding/device/impl/find_tracks.ipp @@ -209,7 +209,7 @@ TRACCC_HOST_DEVICE inline void find_tracks( const detray::tracking_surface sf{det, in_par.surface_link()}; - const bool is_line = sf.template visit_mask(); + const bool is_line = detail::is_line(sf); // Run the Kalman update const kalman_fitter_status res = @@ -527,8 +527,8 @@ TRACCC_HOST_DEVICE inline void find_tracks( } /* - * Compute the offset at which this block will write, as well as the index - * at which this block will write. + * Compute the offset at which this block will write, as well as the + * index at which this block will write. */ if (in_param_is_live) { local_num_params = std::get<1>(decode_insertion_mutex( diff --git a/device/common/include/traccc/finding/device/impl/propagate_to_next_surface.ipp b/device/common/include/traccc/finding/device/impl/propagate_to_next_surface.ipp index eb4ddc91db..e01fa37476 100644 --- a/device/common/include/traccc/finding/device/impl/propagate_to_next_surface.ipp +++ b/device/common/include/traccc/finding/device/impl/propagate_to_next_surface.ipp @@ -62,7 +62,9 @@ TRACCC_HOST_DEVICE inline void propagate_to_next_surface( const bound_track_parameters<> in_par = params.at(param_id); // Create propagator - propagator_t propagator(cfg.propagation); + auto prop_cfg{cfg.propagation}; + prop_cfg.navigation.estimate_scattering_noise = false; + propagator_t propagator(prop_cfg); // Create propagator state typename propagator_t::state propagation(in_par, payload.field_data, det); @@ -87,9 +89,8 @@ TRACCC_HOST_DEVICE inline void propagate_to_next_surface( typename detray::detail::tuple_element<2, actor_tuple_type>::type::state s2{ s3}; // Parameter resetter - // typename detray::detail::tuple_element<4, actor_tuple_type>::type::state - // s4{ - // cfg.propagation}; + typename detray::detail::tuple_element<4, actor_tuple_type>::type::state s4{ + prop_cfg}; // Momentum aborter typename detray::detail::tuple_element<5, actor_tuple_type>::type::state s5; // CKF aborter @@ -101,7 +102,7 @@ TRACCC_HOST_DEVICE inline void propagate_to_next_surface( s5.min_p(static_cast(cfg.min_p)); // Propagate to the next surface - propagator.propagate(propagation, detray::tie(s0, s2, s3, s5, s6)); + propagator.propagate(propagation, detray::tie(s0, s2, s3, s4, s5, s6)); // If a surface found, add the parameter for the next step if (s6.success) { @@ -110,13 +111,27 @@ TRACCC_HOST_DEVICE inline void propagate_to_next_surface( params[param_id] = propagation._stepping.bound_params(); params_liveness[param_id] = 1u; + + const scalar theta = params[param_id].theta(); + if (theta <= 0.f || theta >= 2.f * constant::pi) { + params_liveness[param_id] = 0u; + } + + if (!std::isfinite(params[param_id].phi())) { + params_liveness[param_id] = 0u; + } + + if (math::fabs(params[param_id].qop()) == 0.f) { + params_liveness[param_id] = 0u; + } } else { params_liveness[param_id] = 0u; + } - if (n_cands >= cfg.min_track_candidates_per_track) { - auto tip_pos = tips.push_back(link_idx); - tip_lengths.at(tip_pos) = n_cands; - } + if (params_liveness[param_id] == 0 && + n_cands >= cfg.min_track_candidates_per_track) { + auto tip_pos = tips.push_back(link_idx); + tip_lengths.at(tip_pos) = n_cands; } } diff --git a/device/common/include/traccc/fitting/device/fit_backward.hpp b/device/common/include/traccc/fitting/device/fit_backward.hpp index 70305decf9..83d2c90078 100644 --- a/device/common/include/traccc/fitting/device/fit_backward.hpp +++ b/device/common/include/traccc/fitting/device/fit_backward.hpp @@ -39,7 +39,8 @@ TRACCC_HOST_DEVICE inline void fit_backward( if (param_liveness.at(param_id) > 0u) { typename fitter_t::state fitter_state( track, tracks.states, tracks.measurements, - *(payload.barcodes_view.ptr() + param_id)); + *(payload.barcodes_view.ptr() + param_id), + fitter.config().propagation); kalman_fitter_status fit_status = fitter.smooth(fitter_state); diff --git a/device/common/include/traccc/fitting/device/fit_forward.hpp b/device/common/include/traccc/fitting/device/fit_forward.hpp index 5b3495e05c..c67c035fb8 100644 --- a/device/common/include/traccc/fitting/device/fit_forward.hpp +++ b/device/common/include/traccc/fitting/device/fit_forward.hpp @@ -42,7 +42,7 @@ TRACCC_HOST_DEVICE inline void fit_forward( typename fitter_t::state fitter_state( track, tracks.states, tracks.measurements, - *(payload.barcodes_view.ptr() + param_id)); + *(payload.barcodes_view.ptr() + param_id), fitter.config().propagation); kalman_fitter_status fit_status = fitter.filter(params, fitter_state); diff --git a/examples/options/src/track_propagation.cpp b/examples/options/src/track_propagation.cpp index c1221517f5..881d544754 100644 --- a/examples/options/src/track_propagation.cpp +++ b/examples/options/src/track_propagation.cpp @@ -52,6 +52,25 @@ track_propagation::track_propagation() "search-window", po::value(&m_search_window)->default_value(m_search_window), "Size of the grid surface search window"); + m_desc.add_options()( + "mask-tolerance-scaling", + po::value(&(m_config.navigation.mask_tolerance_scalor)) + ->default_value(m_config.navigation.mask_tolerance_scalor), + "Scale factor between min. and max. mask tolerance with surface " + "distance"); + + m_desc.add_options()( + "accumulated-noise-factor", + po::value(&(m_config.navigation.accumulated_error)) + ->default_value(m_config.navigation.accumulated_error), + "Scale factor on the total track path length to model accumualted " + "noise [%]"); + + m_desc.add_options()( + "scattering-stddevs", + po::value(&(m_config.navigation.n_scattering_stddev)) + ->default_value(m_config.navigation.n_scattering_stddev), + "Number of angle standard deviations to use for the noise modelling"); m_desc.add_options()("rk-tolerance-mm", po::value(&(m_config.stepping.rk_error_tol)) ->default_value(m_config.stepping.rk_error_tol / @@ -101,6 +120,7 @@ void track_propagation::read(const po::variables_map &) { m_config.navigation.min_mask_tolerance *= traccc::unit::mm; m_config.navigation.max_mask_tolerance *= traccc::unit::mm; m_config.navigation.search_window = m_search_window; + m_config.navigation.accumulated_error /= 100.f; m_config.stepping.min_stepsize *= traccc::unit::mm; m_config.stepping.path_limit *= traccc::unit::m; @@ -125,7 +145,7 @@ std::unique_ptr track_propagation::as_printable() traccc::unit::mm) + " mm")); cat_nav->add_child(std::make_unique( - "Mask tolerance scalar", + "Mask tolerance scaling", std::to_string(m_config.navigation.mask_tolerance_scalor))); cat_nav->add_child(std::make_unique( "Path tolerance", std::to_string(m_config.navigation.path_tolerance / @@ -140,6 +160,12 @@ std::unique_ptr track_propagation::as_printable() "Search window", std::to_string(m_config.navigation.search_window[0]) + " x " + std::to_string(m_config.navigation.search_window[1]))); + cat_nav->add_child(std::make_unique( + "Scale factor for accumulated noise", + std::to_string(m_config.navigation.accumulated_error * 100.f) + " %")); + cat_nav->add_child(std::make_unique( + "# scattering stddevs to assume", + std::to_string(m_config.navigation.n_scattering_stddev))); auto cat_tsp = std::make_unique("Transport"); diff --git a/extern/algebra-plugins/CMakeLists.txt b/extern/algebra-plugins/CMakeLists.txt index 15869a9522..7707510a75 100644 --- a/extern/algebra-plugins/CMakeLists.txt +++ b/extern/algebra-plugins/CMakeLists.txt @@ -13,7 +13,7 @@ message( STATUS "Building Algebra Plugins as part of the TRACCC project" ) # Declare where to get Algebra Plugins from. set( TRACCC_ALGEBRA_PLUGINS_SOURCE - "URL;https://github.com/acts-project/algebra-plugins/archive/refs/tags/v0.28.0.tar.gz;URL_MD5;24fa671f564a332858599df60fb2228f" + "URL;https://github.com/acts-project/algebra-plugins/archive/refs/tags/v0.30.0.tar.gz;URL_MD5;d1ec191c6fea5bd0e73f62d4b0319b8c" CACHE STRING "Source for Algebra Plugins, when built as part of this project" ) mark_as_advanced( TRACCC_ALGEBRA_PLUGINS_SOURCE ) FetchContent_Declare( AlgebraPlugins SYSTEM ${TRACCC_ALGEBRA_PLUGINS_SOURCE} ) diff --git a/extern/detray/CMakeLists.txt b/extern/detray/CMakeLists.txt index c3ce19a867..19c60b9218 100644 --- a/extern/detray/CMakeLists.txt +++ b/extern/detray/CMakeLists.txt @@ -13,7 +13,7 @@ message( STATUS "Building Detray as part of the TRACCC project" ) # Declare where to get Detray from. set( TRACCC_DETRAY_SOURCE - "URL;https://github.com/acts-project/detray/archive/refs/tags/v0.103.0.tar.gz;URL_MD5;7df38072d676b63ee3f0f88bd84c0106" + "URL;https://github.com/acts-project/detray/archive/refs/tags/v0.104.1.tar.gz;URL_MD5;2233a3b9945122e1b3f16ac973c7d1be" CACHE STRING "Source for Detray, when built as part of this project" ) mark_as_advanced( TRACCC_DETRAY_SOURCE ) FetchContent_Declare( Detray SYSTEM ${TRACCC_DETRAY_SOURCE} ) diff --git a/simulation/include/traccc/simulation/simulator.hpp b/simulation/include/traccc/simulation/simulator.hpp index d6ba2f8183..bbc9043fe0 100644 --- a/simulation/include/traccc/simulation/simulator.hpp +++ b/simulation/include/traccc/simulation/simulator.hpp @@ -82,6 +82,15 @@ struct simulator { m_cfg.ptc_type = ptc_type; m_track_generator->config().charge(ptc_type.charge()); + + // Turn off tracking features + m_cfg.propagation.stepping.do_covariance_transport = false; + m_cfg.propagation.stepping.use_eloss_gradient = false; + m_cfg.propagation.stepping.use_field_gradient = false; + m_cfg.propagation.navigation.estimate_scattering_noise = false; + + m_resetter = typename detray::parameter_resetter::state{ + m_cfg.propagation}; } config& get_config() { return m_cfg; } @@ -108,8 +117,8 @@ struct simulator { m_scatterer.set_seed(event_id); writer_state.set_seed(event_id); - auto actor_states = - detray::tie(m_aborter_state, m_scatterer, writer_state); + auto actor_states = detray::tie(m_aborter_state, m_scatterer, + m_resetter, writer_state); for (auto track : *m_track_generator.get()) { @@ -149,6 +158,7 @@ struct simulator { /// Actor states typename detray::momentum_aborter::state m_aborter_state{}; typename detray::random_scatterer::state m_scatterer{}; + typename detray::parameter_resetter::state m_resetter{}; }; } // namespace traccc diff --git a/tests/cpu/test_kalman_fitter_momentum_resolution.cpp b/tests/cpu/test_kalman_fitter_momentum_resolution.cpp index 314653cb78..24c609bb97 100644 --- a/tests/cpu/test_kalman_fitter_momentum_resolution.cpp +++ b/tests/cpu/test_kalman_fitter_momentum_resolution.cpp @@ -187,9 +187,9 @@ TEST_P(KalmanFittingMomentumResolutionTests, Run) { // n_trakcs = 100 ASSERT_GE(static_cast(n_tracks), - 0.95 * static_cast(n_truth_tracks)); + 0.94 * static_cast(n_truth_tracks)); ASSERT_GE(static_cast(n_fitted_tracks), - 0.95 * static_cast(n_truth_tracks)); + 0.94 * static_cast(n_truth_tracks)); for (std::size_t i_trk = 0; i_trk < n_tracks; i_trk++) { @@ -242,7 +242,7 @@ TEST_P(KalmanFittingMomentumResolutionTests, Run) { float success_rate = static_cast(n_success) / static_cast(n_truth_tracks * n_events); - ASSERT_GE(success_rate, 0.98f); + ASSERT_GE(success_rate, 0.97f); } // Muon with 1, 10, 100 GeV/c, no materials diff --git a/tests/cpu/test_kalman_fitter_wire_chamber.cpp b/tests/cpu/test_kalman_fitter_wire_chamber.cpp index 4fa6e077d4..38ac359b18 100644 --- a/tests/cpu/test_kalman_fitter_wire_chamber.cpp +++ b/tests/cpu/test_kalman_fitter_wire_chamber.cpp @@ -132,6 +132,8 @@ TEST_P(KalmanFittingWireChamberTests, Run) { fit_cfg.propagation.navigation.min_mask_tolerance = static_cast(mask_tolerance); fit_cfg.propagation.navigation.search_window = search_window; + // TODO: Disable until overlaps are handled correctly + fit_cfg.propagation.navigation.estimate_scattering_noise = false; fit_cfg.ptc_hypothesis = ptc; traccc::host::kalman_fitting_algorithm fitting(fit_cfg, host_mr, copy);