tests/cpu/test_kalman_fitter_momentum_resolution.cpp

0001 /** TRACCC library, part of the ACTS project (R&D line)
0002  *
0003  * (c) 2025-2026 CERN for the benefit of the ACTS project
0004  *
0005  * Mozilla Public License Version 2.0
0006  */
0007
0008 // Project include(s).
0009 #include "traccc/bfield/construct_const_bfield.hpp"
0010 #include "traccc/fitting/kalman_fitting_algorithm.hpp"
0011 #include "traccc/io/read_detector.hpp"
0012 #include "traccc/io/utils.hpp"
0013 #include "traccc/resolution/fitting_performance_writer.hpp"
0014 #include "traccc/simulation/event_generators.hpp"
0015 #include "traccc/simulation/simulator.hpp"
0016 #include "traccc/utils/ranges.hpp"
0017 #include "traccc/utils/seed_generator.hpp"
0018
0019 // Test include(s).
0020 #include "tests/kalman_fitting_momentum_resolution_test.hpp"
0021
0022 // VecMem include(s).
0023 #include <vecmem/memory/host_memory_resource.hpp>
0024 #include <vecmem/utils/copy.hpp>
0025
0026 // GTest include(s).
0027 #include <gtest/gtest.h>
0028
0029 // System include(s).
0030 #include <filesystem>
0031 #include <string>
0032
0033 using namespace traccc;
0034
0035 TEST_P(KalmanFittingMomentumResolutionTests, Run) {
0036   // Get the parameters
0037   const std::string name = std::get<0>(GetParam());
0038   const std::array<scalar, 3u> origin = std::get<1>(GetParam());
0039   const std::array<scalar, 3u> origin_stddev = std::get<2>(GetParam());
0040   const scalar p = std::get<3>(GetParam());
0041   const scalar eta = std::get<4>(GetParam());
0042   const scalar theta = eta_to_theta(eta);
0043   const scalar phi = std::get<5>(GetParam());
0044   const traccc::pdg_particle<scalar> ptc = std::get<6>(GetParam());
0045   const unsigned int n_truth_tracks = std::get<7>(GetParam());
0046   const unsigned int n_events = std::get<8>(GetParam());
0047   const bool random_charge = std::get<9>(GetParam());
0048
0049   // Performance writer
0050   traccc::fitting_performance_writer::config fit_writer_cfg;
0051   fit_writer_cfg.res_config.var_binning["residual_qopT"] =
0052       plot_helpers::binning("r_{q/p_{T}} [c/GeV]", 1000, -0.1f, 0.1f);
0053   fit_writer_cfg.file_path = "performance_track_fitting_" + name + ".root";
0054   traccc::fitting_performance_writer fit_performance_writer(
0055       fit_writer_cfg, traccc::getDefaultLogger("FittingPerformanceWriter",
0056                                                traccc::Logging::Level::INFO));
0057
0058   // Set qop stddev to 10% of truth qop
0059   const scalar qop_stddev = 0.1f / p;
0060
0061   stddevs[4] = qop_stddev;
0062
0063   /*****************************
0064    * Build a telescope geometry
0065    *****************************/
0066
0067   // Memory resources used by the application.
0068   vecmem::host_memory_resource host_mr;
0069   // Copy obejct
0070   vecmem::copy copy;
0071
0072   // Read back detector file
0073   const std::string path = name + "/";
0074   traccc::host_detector detector;
0075   traccc::io::read_detector(
0076       detector, host_mr,
0077       std::filesystem::absolute(
0078           std::filesystem::path(path + "telescope_detector_geometry.json"))
0079           .native(),
0080       (std::get<14>(GetParam()) != detray::vacuum<scalar>()
0081            ? std::filesystem::absolute(
0082                  std::filesystem::path(
0083                      path + "telescope_detector_homogeneous_material.json"))
0084                  .native()
0085            : ""));
0086
0087   auto field = traccc::construct_const_bfield(std::get<13>(GetParam()));
0088
0089   const auto vol0 = detray::tracking_volume{detector.as<detector_traits>(), 0u};
0090
0091   // The number of sensitive surfaces = # of total surfaces - # of portals
0092   // (=6)
0093   const std::size_t n_sensitive_surfaces = vol0.surfaces().size() - 6u;
0094
0095   ASSERT_EQ(n_sensitive_surfaces, std::get<11>(GetParam()));
0096
0097   /***************************
0098    * Generate simulation data
0099    ***************************/
0100
0101   // Track generator
0102   using generator_type =
0103       detray::random_track_generator<traccc::free_track_parameters<>,
0104                                      uniform_gen_t>;
0105   generator_type::configuration gen_cfg{};
0106   gen_cfg.n_tracks(n_truth_tracks);
0107   gen_cfg.origin(origin);
0108   gen_cfg.origin_stddev(origin_stddev);
0109   gen_cfg.phi_range(phi, phi);
0110   gen_cfg.theta_range(theta, theta);
0111   gen_cfg.mom_range(p, p);
0112   gen_cfg.randomize_charge(random_charge);
0113   generator_type generator(gen_cfg);
0114
0115   // Smearing value for measurements
0116   const auto smearing{std::get<15>(GetParam())};
0117   traccc::measurement_smearer<traccc::default_algebra> meas_smearer(
0118       smearing[0], smearing[1]);
0119
0120   using writer_type = traccc::smearing_writer<
0121       traccc::measurement_smearer<traccc::default_algebra>>;
0122
0123   typename writer_type::config smearer_writer_cfg{meas_smearer};
0124   traccc::seed_generator<host_detector_type>::config seed_cfg{};
0125   seed_cfg.initial_sigmas = stddevs;
0126
0127   // Run simulator
0128   const std::string full_path = io::data_directory() + path;
0129   std::filesystem::create_directories(full_path);
0130   auto sim = traccc::simulator<host_detector_type, b_field_t, generator_type,
0131                                writer_type>(
0132       ptc, n_events, detector.as<detector_traits>(),
0133       field.as_field<traccc::const_bfield_backend_t<traccc::scalar>>(),
0134       std::move(generator), std::move(smearer_writer_cfg), full_path);
0135   sim.run();
0136
0137   /***************
0138    * Run fitting
0139    ***************/
0140
0141   // Seed generator
0142   seed_generator<host_detector_type> sg(detector.as<detector_traits>(),
0143                                         seed_cfg);
0144
0145   // Fitting algorithm object
0146   traccc::fitting_config fit_cfg;
0147   fit_cfg.ptc_hypothesis = ptc;
0148   traccc::host::kalman_fitting_algorithm fitting(fit_cfg, host_mr, copy);
0149
0150   // Iterate over events
0151   for (std::size_t i_evt = 0; i_evt < n_events; i_evt++) {
0152     // Event map
0153     traccc::event_data evt_data(path, i_evt, host_mr);
0154     // Truth Track Candidates
0155     traccc::edm::measurement_collection::host measurements(host_mr);
0156     traccc::edm::track_container<traccc::default_algebra>::host
0157         track_candidates{host_mr};
0158     evt_data.generate_truth_candidates(track_candidates, measurements, sg,
0159                                        host_mr);
0160     track_candidates.measurements = vecmem::get_data(measurements);
0161
0162     // n_trakcs = 100
0163     ASSERT_EQ(track_candidates.tracks.size(), n_truth_tracks);
0164
0165     // The nubmer of track candidates per track should be equal to the
0166     // number of planes
0167     for (std::size_t i_trk = 0; i_trk < n_truth_tracks; i_trk++) {
0168       ASSERT_EQ(track_candidates.tracks.at(i_trk).constituent_links().size(),
0169                 std::get<11>(GetParam()));
0170     }
0171
0172     // Run fitting
0173     auto track_states = fitting(
0174         detector, field,
0175         traccc::edm::track_container<traccc::default_algebra>::const_data(
0176             track_candidates));
0177
0178     // Iterator over tracks
0179     const std::size_t n_tracks = track_states.tracks.size();
0180     const std::size_t n_fitted_tracks =
0181         count_successfully_fitted_tracks(track_states.tracks);
0182
0183     // n_trakcs = 100
0184     ASSERT_GE(static_cast<float>(n_tracks),
0185               0.95 * static_cast<float>(n_truth_tracks));
0186     ASSERT_GE(static_cast<float>(n_fitted_tracks),
0187               0.95 * static_cast<float>(n_truth_tracks));
0188
0189     for (std::size_t i_trk = 0; i_trk < n_tracks; i_trk++) {
0190       // Some fits fail. The results of those cannot be reasonably tested.
0191       if (track_states.tracks.at(i_trk).fit_outcome() !=
0192           traccc::track_fit_outcome::SUCCESS) {
0193         continue;
0194       }
0195
0196       consistency_tests(track_states.tracks.at(i_trk), track_states.states);
0197
0198       ndf_tests(track_states.tracks.at(i_trk), track_states.states,
0199                 measurements);
0200
0201       ASSERT_EQ(track_states.tracks.at(i_trk).nholes(), 0u);
0202
0203       fit_performance_writer.write(track_states.tracks.at(i_trk),
0204                                    track_states.states, measurements,
0205                                    detector.as<detector_traits>(), evt_data);
0206     }
0207   }
0208
0209   fit_performance_writer.finalize();
0210
0211   /********************
0212    * Pull value test
0213    ********************/
0214
0215   static const std::vector<std::string> pull_names{
0216       "pull_d0", "pull_z0", "pull_phi", "pull_theta", "pull_qop"};
0217   pull_value_tests(fit_writer_cfg.file_path, pull_names);
0218
0219   /********************
0220    * P-value test
0221    ********************/
0222
0223   p_value_tests(fit_writer_cfg.file_path);
0224
0225   //**************************/
0226   // Momentum resolution test
0227   //**************************/
0228
0229   momentum_resolution_tests(fit_writer_cfg.file_path);
0230
0231   /********************
0232    * Success rate test
0233    ********************/
0234
0235   float success_rate = static_cast<float>(n_success) /
0236                        static_cast<float>(n_truth_tracks * n_events);
0237
0238   ASSERT_GE(success_rate, 0.98f);
0239 }
0240
0241 // Muon with 1, 10, 100 GeV/c, no materials
0242 INSTANTIATE_TEST_SUITE_P(
0243     KalmanFitMomentumResolutionValidation0,
0244     KalmanFittingMomentumResolutionTests,
0245     ::testing::Values(
0246         std::make_tuple(
0247             "mom_resolution_1_GeV_muon", std::array<scalar, 3u>{0.f, 0.f, 0.f},
0248             std::array<scalar, 3u>{0.f, 0.f, 0.f}, 1.f, 0.f, 0.f,
0249             traccc::muon<scalar>(), 100, 100, false, 20.f, 20u, 50.f,
0250             vector3{0, 0, 2 * traccc::unit<scalar>::T},
0251             detray::vacuum<scalar>(),
0252             std::array<scalar, 2u>{50.f * traccc::unit<scalar>::um,
0253                                    50.f * traccc::unit<scalar>::um}),
0254         std::make_tuple(
0255             "mom_resolution_10_GeV_muon", std::array<scalar, 3u>{0.f, 0.f, 0.f},
0256             std::array<scalar, 3u>{0.f, 0.f, 0.f}, 10.f, 0.f, 0.f,
0257             traccc::muon<scalar>(), 100, 100, false, 20.f, 20u, 50.f,
0258             vector3{0, 0, 2 * traccc::unit<scalar>::T},
0259             detray::vacuum<scalar>(),
0260             std::array<scalar, 2u>{50.f * traccc::unit<scalar>::um,
0261                                    50.f * traccc::unit<scalar>::um}),
0262         std::make_tuple("mom_resolution_100_GeV_muon",
0263                         std::array<scalar, 3u>{0.f, 0.f, 0.f},
0264                         std::array<scalar, 3u>{0.f, 0.f, 0.f}, 100.f, 0.f, 0.f,
0265                         traccc::muon<scalar>(), 100, 100, false, 20.f, 20u,
0266                         50.f, vector3{0, 0, 2 * traccc::unit<scalar>::T},
0267                         detray::vacuum<scalar>(),
0268                         std::array<scalar, 2u>{
0269                             50.f * traccc::unit<scalar>::um,
0270                             50.f * traccc::unit<scalar>::um})));
0271
0272 // Muon with 1 GeV/c and different smearing parameters, no materials
0273 INSTANTIATE_TEST_SUITE_P(
0274     KalmanFitMomentumResolutionValidation1,
0275     KalmanFittingMomentumResolutionTests,
0276     ::testing::Values(
0277         std::make_tuple("mom_resolution_1_GeV_muon_50_100_smearing",
0278                         std::array<scalar, 3u>{0.f, 0.f, 0.f},
0279                         std::array<scalar, 3u>{0.f, 0.f, 0.f}, 1.f, 0.f, 0.f,
0280                         traccc::muon<scalar>(), 100, 100, false, 20.f, 20u,
0281                         50.f, vector3{0, 0, 2 * traccc::unit<scalar>::T},
0282                         detray::vacuum<scalar>(),
0283                         std::array<scalar, 2u>{
0284                             50.f * traccc::unit<scalar>::um,
0285                             100.f * traccc::unit<scalar>::um}),
0286         std::make_tuple("mom_resolution_1_GeV_muon_100_50_smearing",
0287                         std::array<scalar, 3u>{0.f, 0.f, 0.f},
0288                         std::array<scalar, 3u>{0.f, 0.f, 0.f}, 1.f, 0.f, 0.f,
0289                         traccc::muon<scalar>(), 100, 100, false, 20.f, 20u,
0290                         50.f, vector3{0, 0, 2 * traccc::unit<scalar>::T},
0291                         detray::vacuum<scalar>(),
0292                         std::array<scalar, 2u>{
0293                             100.f * traccc::unit<scalar>::um,
0294                             50.f * traccc::unit<scalar>::um}),
0295         std::make_tuple("mom_resolution_1_GeV_muon_100_100_smearing",
0296                         std::array<scalar, 3u>{0.f, 0.f, 0.f},
0297                         std::array<scalar, 3u>{0.f, 0.f, 0.f}, 1.f, 0.f, 0.f,
0298                         traccc::muon<scalar>(), 100, 100, false, 20.f, 20u,
0299                         50.f, vector3{0, 0, 2 * traccc::unit<scalar>::T},
0300                         detray::vacuum<scalar>(),
0301                         std::array<scalar, 2u>{
0302                             100.f * traccc::unit<scalar>::um,
0303                             100.f * traccc::unit<scalar>::um})));