Add tests for momentum resolution validation

performance/include/traccc/utils/ranges.hpp

@@ -16,21 +16,29 @@
 
 namespace traccc {
 
+template <typename scalar_t>
+TRACCC_HOST_DEVICE inline scalar_t eta_to_theta(const scalar_t eta) {
+    return 2.f * math::atan(std::exp(-eta));
+}
+
+template <typename scalar_t>
+TRACCC_HOST_DEVICE inline scalar_t theta_to_eta(const scalar_t theta) {
+    return -math::log(math::tan(theta * 0.5f));
+}
+
 template <typename range_t>
 TRACCC_HOST_DEVICE inline range_t eta_to_theta_range(const range_t& eta_range) {
     // @NOTE: eta_range[0] is converted to theta_range[1] and eta_range[1]
     // to theta_range[0] because theta(minEta) > theta(maxEta)
-    return {2.f * math::atan(std::exp(-eta_range[1])),
-            2.f * math::atan(std::exp(-eta_range[0]))};
+    return {eta_to_theta(eta_range[1]), eta_to_theta(eta_range[0])};
 }
 
 template <typename range_t>
 TRACCC_HOST_DEVICE inline range_t theta_to_eta_range(
     const range_t& theta_range) {
     // @NOTE: theta_range[0] is converted to eta_range[1] and theta_range[1]
     // to eta_range[0]
-    return {-math::log(math::tan(theta_range[1] * 0.5f)),
-            -math::log(math::tan(theta_range[0] * 0.5f))};
+    return {theta_to_eta(theta_range[1]), theta_to_eta(theta_range[0])};
 }
 
 }  // namespace traccc

tests/CMakeLists.txt

@@ -13,10 +13,12 @@ add_library( traccc_tests_common STATIC
     "common/tests/cca_test.hpp"
     "common/tests/ckf_telescope_test.hpp"
     "common/tests/data_test.hpp"
+    "common/tests/kalman_fitting_momentum_resolution_test.hpp"
     "common/tests/kalman_fitting_test.hpp"
     "common/tests/kalman_fitting_telescope_test.hpp"
     "common/tests/kalman_fitting_toy_detector_test.hpp"
     "common/tests/kalman_fitting_wire_chamber_test.hpp"
+    "common/tests/kalman_fitting_momentum_resolution_test.cpp"
     "common/tests/kalman_fitting_test.cpp" )
 target_include_directories( traccc_tests_common
     PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/common )

tests/common/tests/kalman_fitting_momentum_resolution_test.cpp

@@ -0,0 +1,160 @@
+/** 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
+ */
+
+#include <cmath>
+
+// Local include(s).
+#include "kalman_fitting_momentum_resolution_test.hpp"
+
+// ROOT include(s).
+#ifdef TRACCC_HAVE_ROOT
+#include <TF1.h>
+#include <TFile.h>
+#include <TFitResult.h>
+#include <TFitResultPtr.h>
+#include <TH1.h>
+#endif  // TRACCC_HAVE_ROOT
+
+// System include(s).
+#include <iostream>
+#include <memory>
+#include <stdexcept>
+
+namespace traccc {
+
+void KalmanFittingMomentumResolutionTests::consistency_tests(
+    const track_state_collection_types::host& track_states_per_track) const {
+
+    // The nubmer of track states is supposed be equal to the number
+    // of planes
+    ASSERT_EQ(track_states_per_track.size(), std::get<11>(GetParam()));
+}
+
+void KalmanFittingMomentumResolutionTests::momentum_resolution_tests([
+    [maybe_unused]] std::string_view file_name) const {
+
+#ifdef TRACCC_HAVE_ROOT
+
+    // Open the file with the histograms.
+    std::unique_ptr<TFile> ifile(TFile::Open(file_name.data(), "READ"));
+    if ((!ifile) || ifile->IsZombie()) {
+        throw std::runtime_error(std::string("Could not open file \"") +
+                                 file_name.data() + "\"");
+    }
+
+    // Access the histogram.
+    TH1* residual_qopT_hist = dynamic_cast<TH1*>(ifile->Get("res_qopT"));
+    if (!residual_qopT_hist) {
+        throw std::runtime_error(
+            std::string("Could not access histogram residual_qopT in file \"") +
+            file_name.data() + "\"");
+    }
+
+    // Function used for the fit.
+    TF1 gaus{"gaus", "gaus", -5, 5};
+    double fit_par[3];
+
+    // Set the mean seed to 0
+    gaus.SetParameters(1, 0.);
+    gaus.SetParLimits(1, -1., 1.);
+
+    // Set the standard deviation seed to 0.01
+    gaus.SetParameters(2, 0.01f);
+    gaus.SetParLimits(2, 0.f, 0.1f);
+
+    auto res = residual_qopT_hist->Fit("gaus", "Q0S");
+
+    gaus.GetParameters(&fit_par[0]);
+
+    // Mean check
+    EXPECT_NEAR(fit_par[1], 0, 0.05) << " Residual qopT mean value error";
+
+    // Expected momentum resolution
+    const std::size_t N = std::get<11>(GetParam());
+    const auto smearing = std::get<15>(GetParam());
+    const scalar epsilon = smearing[0u];
+    const scalar Bz = std::get<13>(GetParam())[2u];
+    const scalar spacing = std::get<12>(GetParam());
+    const scalar L = static_cast<scalar>(N - 1u) * spacing;
+
+    // Curvature (1/helix_radius) resolution from detector noise Eq. (35.61) of
+    // PDG 2024
+    const scalar dkres =
+        epsilon / (L * L) * math::sqrt(720.f / static_cast<scalar>(N + 4));
+
+    // Curvature (1/helix_radius) resolution from multiple scattering Eq.
+    // (35.63) of PDG 2024
+    // @NOTE: The calculation of the multiple scattering term is in work in
+    // progress. Currently we are validating the momentum resolution without any
+    // material on the modules, therefore, the multiple scattering contribution
+    // to the momentum resolution is set to zero.
+    const scalar dkms = 0.f;
+
+    // Eq. (35.60) of PDG 2024
+    const scalar dk = math::sqrt(dkres * dkres + dkms * dkms);
+
+    // σ(q/pT) = σ(k/B) = σ(k)/B
+    const scalar expected_sigma_qopT = dk / Bz;
+
+    // Sigma check
+    EXPECT_NEAR(fit_par[2], expected_sigma_qopT, expected_sigma_qopT * 0.05f)
+        << " Residual qopT sigma value error";
+
+#else
+    std::cout << "Pull value tests not performed without ROOT" << std::endl;
+#endif  // TRACCC_HAVE_ROOT
+
+    return;
+}
+
+void KalmanFittingMomentumResolutionTests::SetUp() {
+
+    vecmem::host_memory_resource host_mr;
+
+    const scalar offset = std::get<10>(GetParam());
+    const unsigned int n_planes = std::get<11>(GetParam());
+    const scalar spacing = std::get<12>(GetParam());
+
+    std::vector<scalar> plane_positions;
+    for (unsigned int i = 0; i < n_planes; i++) {
+        plane_positions.push_back(offset * unit<scalar>::mm +
+                                  static_cast<scalar>(i) * spacing *
+                                      unit<scalar>::mm);
+    }
+
+    /// Plane alignment direction (aligned to x-axis)
+    const vector3 bfield = std::get<13>(GetParam());
+
+    const auto p = std::get<3>(GetParam());
+    const detray::detail::helix<traccc::default_algebra> hlx(
+        {0, 0, 0}, 0, {1, 0, 0}, -1.f / p, &bfield);
+
+    constexpr scalar rect_half_length = 500.f;
+    constexpr detray::mask<detray::rectangle2D, traccc::default_algebra>
+        rectangle{0u, rect_half_length, rect_half_length};
+
+    detray::tel_det_config<default_algebra, detray::rectangle2D,
+                           detray::detail::helix>
+        tel_cfg(rectangle, hlx);
+    tel_cfg.positions(plane_positions);
+    tel_cfg.module_material(std::get<14>(GetParam()));
+    tel_cfg.mat_thickness(thickness);
+    tel_cfg.envelope(rect_half_length);
+
+    // Create telescope detector
+    auto [det, name_map] = build_telescope_detector(host_mr, tel_cfg);
+
+    // Write detector file
+    auto writer_cfg = detray::io::detector_writer_config{}
+                          .format(detray::io::format::json)
+                          .replace_files(true)
+                          .write_material(true)
+                          .path(std::get<0>(GetParam()));
+    detray::io::write_detector(det, name_map, writer_cfg);
+}
+
+}  // namespace traccc

tests/common/tests/kalman_fitting_momentum_resolution_test.hpp

@@ -0,0 +1,71 @@
+/** 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 "kalman_fitting_test.hpp"
+
+// Detray include(s).
+#include "detray/geometry/mask.hpp"
+#include "detray/geometry/shapes/rectangle2D.hpp"
+#include "detray/io/frontend/detector_writer.hpp"
+#include "detray/navigation/detail/ray.hpp"
+#include "detray/test/utils/detectors/build_telescope_detector.hpp"
+
+namespace traccc {
+
+/// Momentum Resolution Test with Telescope Detector
+///
+/// Test parameters:
+/// (1) name
+/// (2) origin
+/// (3) origin stddev
+/// (4) momentum range
+/// (5) eta range
+/// (6) phi range
+/// (7) particle type
+/// (8) number of tracks per event
+/// (9) number of events
+/// (10) random charge
+/// (11) offset from origin of the first plane in mm
+/// (12) Number of planes
+/// (13) Spacing between planes in mm
+/// (14) Magnetic field
+/// (15) Module material
+/// (16) Measurement smearing
+class KalmanFittingMomentumResolutionTests
+    : public KalmanFittingTests,
+      public testing::WithParamInterface<std::tuple<
+          std::string, std::array<scalar, 3u>, std::array<scalar, 3u>, scalar,
+          scalar, scalar, detray::pdg_particle<scalar>, unsigned int,
+          unsigned int, bool, scalar, unsigned int, scalar, vector3,
+          detray::material<scalar>, std::array<scalar, 2u>>> {
+
+    public:
+    /// Plane thickness
+    static constexpr scalar thickness = 0.5f * detray::unit<scalar>::mm;
+
+    /// Standard deviations for seed track parameters
+    std::array<scalar, e_bound_size> stddevs = {
+        5.f * detray::unit<scalar>::mm,
+        5.f * detray::unit<scalar>::mm,
+        0.05f,
+        0.05f,
+        0.1f / detray::unit<scalar>::GeV,
+        1.f * detray::unit<scalar>::ns};
+
+    void consistency_tests(
+        const track_state_collection_types::host& track_states_per_track) const;
+
+    void momentum_resolution_tests(std::string_view file_name) const;
+
+    protected:
+    virtual void SetUp() override;
+};
+
+}  // namespace traccc

tests/cpu/CMakeLists.txt

@@ -14,6 +14,7 @@ traccc_add_test(cpu
     "test_clusterization_resolution.cpp"
     "test_copy.cpp"
     "test_kalman_fitter_hole_count.cpp"
+    "test_kalman_fitter_momentum_resolution.cpp"
     "test_kalman_fitter_telescope.cpp"
     "test_kalman_fitter_wire_chamber.cpp"
     "test_ranges.cpp"