traccc/finding/measurement_selector.hpp

0001 /** TRACCC library, part of the ACTS project (R&D line)
0002  *
0003  * (c) 2026 CERN for the benefit of the ACTS project
0004  *
0005  * Mozilla Public License Version 2.0
0006  */
0007
0008 #pragma once
0009
0010 // Project include(s).
0011 #include "traccc/definitions/primitives.hpp"
0012 #include "traccc/definitions/qualifiers.hpp"
0013 #include "traccc/definitions/track_parametrization.hpp"
0014 #include "traccc/edm/measurement_collection.hpp"
0015 #include "traccc/edm/measurement_helpers.hpp"
0016 #include "traccc/edm/track_parameters.hpp"
0017 #include "traccc/utils/logging.hpp"
0018 #include "traccc/utils/matrix_helpers.hpp"
0019 #include "traccc/utils/subspace.hpp"
0020
0021 // System include(s)
0022 #include <limits>
0023
0024 namespace traccc {
0025
0026 /// Potential next measurement to be added to track
0027 struct candidate_measurement {
0028   unsigned int meas_idx{std::numeric_limits<unsigned int>::max()};
0029   float chi2{std::numeric_limits<float>::max()};
0030
0031   /// Define comparisons
0032   constexpr bool operator<=>(const candidate_measurement& other) const =
0033       default;
0034 };
0035
0036 /// Associate a measurement to a candidate track
0037 struct measurement_selector {
0038   template <detray::concepts::algebra A, unsigned int ROWS, unsigned int COLS>
0039   using matrix_t = detray::dmatrix<A, ROWS, COLS>;
0040
0041   // Where to get the calibration from
0042   struct config { /*TODO: implement calibration handling*/
0043   };
0044
0045   /// Get the observation model for a given measurement
0046   ///
0047   /// @param measurement the measurement
0048   /// @param bound_params predicted bound track parameters
0049   /// @param is_line whether the measurement belong to a line surface
0050   ///
0051   /// @returns the projection matrix H
0052   template <detray::concepts::algebra algebra_t, unsigned int D,
0053             typename measurement_backend_t>
0054   TRACCC_HOST_DEVICE static detray::dmatrix<algebra_t, D, e_bound_size>
0055   observation_model(const edm::measurement<measurement_backend_t>& measurement,
0056                     const bound_track_parameters<algebra_t>& bound_params,
0057                     const bool is_line) {
0058     // Oservation model: Subspace of measurement space for this measurement
0059     subspace<algebra_t, e_bound_size> subs(measurement.subspace());
0060
0061     // Flip the sign of projector matrix element in case the first element
0062     // of a line measurement is negative
0063     if (is_line && bound_params.bound_local()[e_bound_loc0] < 0) {
0064       subs.set_sign(0, true);
0065     }
0066
0067     detray::dmatrix<algebra_t, D, 1> meas_local;
0068     edm::get_measurement_local<algebra_t>(measurement, meas_local);
0069     if (measurement.dimensions() == 1) {
0070       subs.set_invalid(1);
0071     }
0072
0073     const detray::dmatrix<algebra_t, D, e_bound_size> H =
0074         subs.template projector<D>();
0075
0076     TRACCC_DEBUG_HOST("--> Observation model (H):\n" << H);
0077
0078     return H;
0079   }
0080
0081   /// Get the calibrated measurement position
0082   ///
0083   /// @param measurement the measurement
0084   /// @param cfg how to apply calibrations
0085   ///
0086   /// @returns the projection matrix H
0087   template <detray::concepts::algebra algebra_t, unsigned int D,
0088             typename measurement_backend_t>
0089   TRACCC_HOST_DEVICE static detray::dmatrix<algebra_t, D, 1>
0090   calibrated_measurement_position(
0091       const edm::measurement<measurement_backend_t>& measurement,
0092       const config& /*cfg*/) {
0093     // Measurement local position on surface
0094     detray::dmatrix<algebra_t, D, 1> meas_local;
0095     edm::get_measurement_local<algebra_t>(measurement, meas_local);
0096
0097     TRACCC_DEBUG_HOST("--> Measurement position (uncalibrated):\n"
0098                       << meas_local);
0099
0100     assert((measurement.dimensions() > 1) ||
0101            (getter::element(meas_local, 1u, 0u) == 0.f));
0102
0103     return meas_local;
0104   }
0105
0106   /// Get the calibrated measurement covariance
0107   ///
0108   /// @param measurement the measurement
0109   /// @param cfg how to apply calibrations
0110   ///
0111   /// @returns the projection matrix H
0112   template <detray::concepts::algebra algebra_t, unsigned int D,
0113             typename measurement_backend_t>
0114   TRACCC_HOST_DEVICE static detray::dmatrix<algebra_t, D, D>
0115   calibrated_measurement_covariance(
0116       const edm::measurement<measurement_backend_t>& measurement,
0117       const config& /*cfg*/) {
0118     // Measurement covariance
0119     detray::dmatrix<algebra_t, D, D> V;
0120     edm::get_measurement_covariance<algebra_t>(measurement, V);
0121
0122     detray::dmatrix<algebra_t, D, 1> meas_local;
0123     edm::get_measurement_local<algebra_t>(measurement, meas_local);
0124     if (measurement.dimensions() == 1) {
0125       getter::element(V, 1u, 1u) =
0126           std::numeric_limits<detray::dscalar<algebra_t>>::max();
0127     }
0128
0129     TRACCC_DEBUG_HOST("--> Measurement covariance (uncalibrated):\n" << V);
0130
0131     return V;
0132   }
0133
0134   /// Calculate the predicted chi2
0135   ///
0136   /// @brief Based on "Application of Kalman filtering to track and vertex
0137   /// fitting", R.Frühwirth, NIM A
0138   ///
0139   /// @param measurement the measurement
0140   /// @param bound_params predicted bound track parameters
0141   /// @param cfg the calibration configuration
0142   /// @param is_line whether the measurement belong to a line surface
0143   ///
0144   /// @returns the predicted chi2 of the calibrated measurement
0145   template <typename measurement_backend_t, detray::concepts::algebra algebra_t>
0146   TRACCC_HOST_DEVICE static detray::dscalar<algebra_t> predicted_chi2(
0147       const edm::measurement<measurement_backend_t>& measurement,
0148       const bound_track_parameters<algebra_t>& bound_params, const config& cfg,
0149       const bool is_line) {
0150     using scalar_t = detray::dscalar<algebra_t>;
0151
0152     // Measurement maximal dimension
0153     constexpr unsigned int D = 2;
0154
0155     TRACCC_VERBOSE_HOST_DEVICE("--> dim: %d", measurement.dimensions());
0156
0157     assert(measurement.dimensions() == 1u || measurement.dimensions() == 2u);
0158
0159     assert(!bound_params.is_invalid());
0160     assert(!bound_params.surface_link().is_invalid());
0161
0162     // Get calibrated measurement and covariance
0163     const matrix_t<algebra_t, D, 1> meas_local =
0164         calibrated_measurement_position<algebra_t, D>(measurement, cfg);
0165
0166     const matrix_t<algebra_t, D, D> V =
0167         calibrated_measurement_covariance<algebra_t, D>(measurement, cfg);
0168
0169     // Project the predicted covariance to the observation
0170     const matrix_t<algebra_t, D, e_bound_size> H =
0171         observation_model<algebra_t, D>(measurement, bound_params, is_line);
0172
0173     const matrix_t<algebra_t, D, D> R =
0174         H * matrix::transposed_product<false, true>(bound_params.covariance(),
0175                                                     H) +
0176         V;
0177
0178     const matrix_t<algebra_t, D, D> R_inv =
0179         masked_inverse<algebra_t>(R, measurement.dimensions());
0180
0181     TRACCC_DEBUG_HOST("--> R:\n" << R);
0182     TRACCC_DEBUG_HOST("--> R_inv:\n" << R_inv);
0183
0184     // Residual between measurement and (projected) vector (innovation)
0185     const matrix_t<algebra_t, D, 1> residual =
0186         meas_local - H * bound_params.vector();
0187
0188     TRACCC_DEBUG_HOST("--> Predicted residual:\n" << residual);
0189
0190     const matrix_t<algebra_t, 1, 1> pred_chi2 =
0191         matrix::transposed_product<true, false>(residual, R_inv) * residual;
0192
0193     const scalar_t pred_chi2_val{getter::element(pred_chi2, 0, 0)};
0194
0195     if (!std::isfinite(pred_chi2_val)) {
0196       TRACCC_WARNING_HOST_DEVICE("Infinite predicted chi2 value!");
0197     } else if (pred_chi2_val < 0.f) {
0198       TRACCC_WARNING_HOST_DEVICE("Negative predicted chi2 value!");
0199     } else {
0200       TRACCC_VERBOSE_HOST_DEVICE("--> chi2: %.10e", pred_chi2_val);
0201     }
0202
0203     return pred_chi2_val;
0204   }
0205
0206   /// Measurement selection (optimal)
0207   ///
0208   /// @param bound_params predicted bound track parameters
0209   /// @param measurements the measurement container
0210   /// @param meas_range contains the index ranges into the measurements
0211   /// @param cfg the calibration configuration
0212   /// @param is_line whether the measurement belong to a line surface
0213   ///
0214   /// @returns the optimal candidate measurement for the input params
0215   template <detray::concepts::algebra algebra_t>
0216   TRACCC_HOST_DEVICE static candidate_measurement find_optimal_measurement(
0217       const bound_track_parameters<algebra_t>& bound_params,
0218       const typename edm::measurement_collection::const_device& measurements,
0219       vecmem::device_vector<unsigned int> meas_ranges, const config& cfg,
0220       const bool is_line) {
0221     using scalar_t = detray::dscalar<algebra_t>;
0222
0223     // The optimal candidate
0224     candidate_measurement cand{};
0225
0226     // Iterate over the measurements for this surface
0227     const unsigned int sf_idx{bound_params.surface_link().index()};
0228     const unsigned int lo{sf_idx == 0u ? 0u : meas_ranges[sf_idx - 1]};
0229     const unsigned int up{meas_ranges[sf_idx]};
0230
0231     TRACCC_VERBOSE_HOST_DEVICE("Have %d measurement(s) on surface %d...",
0232                                up - lo, sf_idx);
0233
0234     // Find the best fitting measurement by prediced chi2
0235     // TODO: Load balancing
0236     for (unsigned int meas_idx = lo; meas_idx < up; meas_idx++) {
0237       TRACCC_VERBOSE_HOST_DEVICE("-> measurement %d:", meas_idx);
0238
0239       // Predicted chi2
0240       const scalar_t chi2 = measurement_selector::predicted_chi2(
0241           measurements.at(meas_idx), bound_params, cfg, is_line);
0242
0243       // Check predicted chi2 cut
0244       if (chi2 < cand.chi2 && chi2 >= 0.f) {
0245         cand = {meas_idx, static_cast<float>(chi2)};
0246         // Found optimal
0247         if (cand.chi2 <= std::numeric_limits<scalar_t>::epsilon()) {
0248           return cand;
0249         }
0250       }
0251     }
0252
0253     return cand;
0254   }
0255
0256   /// Measurement selection (collection of compatible measurements)
0257   ///
0258   /// @param bound_params predicted bound track parameters
0259   /// @param measurements the measurement container
0260   /// @param meas_range contains the index ranges into the measurements
0261   /// @param cfg the calibration configuration
0262   /// @param is_line whether the measurement belong to a line surface
0263   ///
0264   /// @returns a collection of compatible measurements, sorted by pred.
0265   /// chi2
0266   template <detray::concepts::algebra algebra_t>
0267   TRACCC_HOST_DEVICE static vecmem::vector<candidate_measurement>
0268   find_compatible_measurements(
0269       const bound_track_parameters<algebra_t>& /*bound_params*/,
0270       const typename edm::measurement_collection::const_device&
0271       /*measurements*/,
0272       vecmem::device_vector<unsigned int> /*meas_ranges*/,
0273       const config& /*cfg*/, const bool /*is_line*/) {
0274     /* TODO: Implement*/
0275     assert(false);
0276     return {};
0277   }
0278 };
0279
0280 }  // namespace traccc