From 223d8f96b922b8ded51d2f623ac3932c51bdb60c Mon Sep 17 00:00:00 2001 From: Joana Niermann Date: Tue, 21 Oct 2025 17:01:42 +0200 Subject: [PATCH] Harmonize acts and traccc parameter smearing --- core/include/traccc/utils/seed_generator.hpp | 142 ++++++++++++++++-- .../cpu/misaligned_truth_fitting_example.cpp | 13 +- examples/run/cpu/truth_finding_example.cpp | 11 +- examples/run/cpu/truth_fitting_example.cpp | 11 +- .../run/cuda/truth_finding_example_cuda.cpp | 11 +- .../run/cuda/truth_fitting_example_cuda.cpp | 11 +- .../cpu/test_ckf_combinatorics_telescope.cpp | 3 +- .../cpu/test_ckf_sparse_tracks_telescope.cpp | 3 +- tests/cpu/test_kalman_fitter_hole_count.cpp | 3 +- ...test_kalman_fitter_momentum_resolution.cpp | 3 +- tests/cpu/test_kalman_fitter_telescope.cpp | 3 +- tests/cpu/test_kalman_fitter_wire_chamber.cpp | 3 +- .../cuda/test_ckf_combinatorics_telescope.cpp | 3 +- tests/cuda/test_ckf_toy_detector.cpp | 3 +- tests/cuda/test_kalman_fitter_telescope.cpp | 2 +- .../sycl/test_ckf_combinatorics_telescope.cpp | 3 +- tests/sycl/test_ckf_toy_detector.cpp | 3 +- tests/sycl/test_kalman_fitter_telescope.cpp | 2 +- 18 files changed, 144 insertions(+), 89 deletions(-) diff --git a/core/include/traccc/utils/seed_generator.hpp b/core/include/traccc/utils/seed_generator.hpp index 5c0c03e671..d3bccc3790 100644 --- a/core/include/traccc/utils/seed_generator.hpp +++ b/core/include/traccc/utils/seed_generator.hpp @@ -10,6 +10,7 @@ // Library include(s). #include "traccc/edm/track_parameters.hpp" #include "traccc/utils/particle.hpp" +#include "traccc/utils/trigonometric_helpers.hpp" // detray include(s). #include @@ -28,14 +29,53 @@ struct seed_generator { using algebra_type = typename detector_t::algebra_type; using ctx_t = typename detector_t::geometry_context; + /// Configure the seed generator + struct config { + // Smearing parameters + /// Constant term of the loc0 resolution. + double sigma_loc0 = 20 * traccc::unit::um; + /// Pt-dependent loc0 resolution of the form sigma_loc0 = + /// A*exp(-1.*abs(B)*pt). + double sigma_loc0_pT_a = 30 * traccc::unit::um; + double sigma_loc0_pT_b = 0.3 / traccc::unit::GeV; + /// Constant term of the loc1 resolution. + double sigma_loc1 = 20 * traccc::unit::um; + /// Pt-dependent loc1 resolution of the form sigma_loc1 = + /// A*exp(-1.*abs(B)*pt). + double sigma_loc1_pT_a = 30 * traccc::unit::um; + double sigma_loc1_pT_b = 0.3 / traccc::unit::GeV; + /// Time resolution. + double sigma_time = 1 * traccc::unit::ns; + /// Phi angular resolution. + double sigma_phi = 1 * traccc::unit::degree; + /// Theta angular resolution. + double sigma_theta = 1 * traccc::unit::degree; + /// Relative transverse momentum resolution. + double sigma_pT_rel = 0.05; + + /// Optional. Initial sigmas for the track parameters which overwrites + /// the smearing params if set. + std::optional> initial_sigmas; + /// Initial sigma(q/pt) for the track parameters. + /// @note The resulting q/p sigma is added to the one in `initialSigmas` + double initialsigma_qopt = + 0.1 * traccc::unit::e / traccc::unit::GeV; + /// Initial sigma(pt)/pt for the track parameters. + /// @note The resulting q/p sigma is added to the one in `initialSigmas` + double initialsigma_pT_rel = 0.1; + + /// Inflation factors for the variances + std::array cov_inflation{1., 1., 1., 1., 1., 1.}; + }; + /// Constructor with detector /// /// @param det input detector /// @param stddevs standard deviations for parameter smearing - seed_generator(const detector_t& det, - const std::array& stddevs, - const std::size_t sd = 0, ctx_t ctx = {}) - : m_detector(det), m_stddevs(stddevs), m_ctx(ctx) { + explicit seed_generator(const detector_t& det, const config& cfg = {}, + const std::size_t sd = std::mt19937::default_seed, + ctx_t ctx = {}) + : m_detector(det), m_ctx(ctx), m_cfg(cfg) { m_generator.seed(static_cast(sd)); } @@ -73,32 +113,102 @@ struct seed_generator { sf); } - for (std::size_t i = 0; i < e_bound_size; i++) { + // Call the smearing operator + (*this)(bound_param, ptc_type); - if (m_stddevs[i] != scalar{0}) { - bound_param[i] = std::normal_distribution( - bound_param[i], m_stddevs[i])(m_generator); - } + return bound_param; + } + + /// Seed generator operation + /// + /// @param vertex vertex of particle + /// @param stddevs standard deviations for track parameter smearing + void operator()(bound_track_parameters& bound_param, + const traccc::pdg_particle& ptc_type) { + + const scalar q{ptc_type.charge()}; + const auto pos = bound_param.bound_local(); + const auto time = bound_param.time(); + const auto phi = bound_param.phi(); + const auto theta = bound_param.theta(); + const auto pt = bound_param.pT(q); + const auto qop = bound_param.qop(); + + // Compute momentum-dependent resolutions + const double sigma_loc_0 = + m_cfg.sigma_loc0 + + m_cfg.sigma_loc0_pT_a * + math::exp(-1.0 * math::fabs(m_cfg.sigma_loc0_pT_b) * pt); + const double sigma_loc_1 = + m_cfg.sigma_loc1 + + m_cfg.sigma_loc1_pT_a * + math::exp(-1.0 * math::fabs(m_cfg.sigma_loc1_pT_b) * pt); + // Shortcuts for other resolutions + const double sigma_qop = math::sqrt( + math::pow(m_cfg.sigma_pT_rel * qop, 2) + + math::pow(m_cfg.sigma_theta * (qop * math::tan(theta)), 2)); + + // Smear the position/time + // Note that we smear d0 and z0 in the perigee frame + const scalar smeared_loc0 = std::normal_distribution( + pos[0], static_cast(sigma_loc_0))(m_generator); + const scalar smeared_loc1 = std::normal_distribution( + pos[1], static_cast(sigma_loc_1))(m_generator); + bound_param.set_bound_local({smeared_loc0, smeared_loc1}); + + // Time + bound_param.set_time(std::normal_distribution( + time, static_cast(m_cfg.sigma_time))(m_generator)); + + // Smear direction angles phi,theta ensuring correct bounds + const scalar smeared_phi = std::normal_distribution( + phi, static_cast(m_cfg.sigma_phi))(m_generator); + const scalar smeared_theta = std::normal_distribution( + theta, static_cast(m_cfg.sigma_theta))(m_generator); + const auto [new_phi, new_theta] = + detail::wrap_phi_theta(smeared_phi, smeared_theta); + bound_param.set_phi(new_phi); + bound_param.set_theta(new_theta); + + // Compute smeared q/p + bound_param.set_qop(std::normal_distribution( + qop, static_cast(sigma_qop))(m_generator)); + + // Build the track covariance matrix using the smearing sigmas + std::array sigmas{ + static_cast(sigma_loc_0), + static_cast(sigma_loc_1), + static_cast(m_cfg.sigma_phi), + static_cast(m_cfg.sigma_theta), + static_cast(sigma_qop), + static_cast(m_cfg.sigma_time)}; + + for (std::size_t i = e_bound_loc0; i < e_bound_size; ++i) { + double sigma = sigmas[i]; + double variance = sigma * sigma; + + // Inflate the initial covariance + variance *= m_cfg.cov_inflation[i]; getter::element(bound_param.covariance(), i, i) = - m_stddevs[i] * m_stddevs[i]; + static_cast(variance); } assert(!bound_param.is_invalid()); - - return bound_param; } private: - // Random generator + /// Random generator std::random_device m_rd{}; std::mt19937 m_generator{m_rd()}; - // Detector object + /// Detector object const detector_t& m_detector; - /// Standard deviations for parameter smearing - std::array m_stddevs; + /// Geometry context ctx_t m_ctx; + + /// Seed generator configuration + config m_cfg{}; }; } // namespace traccc diff --git a/examples/run/cpu/misaligned_truth_fitting_example.cpp b/examples/run/cpu/misaligned_truth_fitting_example.cpp index fa99388861..f2dc5d2a05 100644 --- a/examples/run/cpu/misaligned_truth_fitting_example.cpp +++ b/examples/run/cpu/misaligned_truth_fitting_example.cpp @@ -110,15 +110,6 @@ int main(int argc, char* argv[]) { * Do the reconstruction *****************************/ - /// Standard deviations for seed track parameters - static constexpr std::array stddevs = { - 0.03f * traccc::unit::mm, - 0.03f * traccc::unit::mm, - 0.017f, - 0.017f, - 0.001f / traccc::unit::GeV, - 1.f * traccc::unit::ns}; - // Fitting algorithm objects // Alg0 traccc::fitting_config fit_cfg0(fitting_opts); @@ -135,9 +126,9 @@ int main(int argc, char* argv[]) { // Seed generators traccc::seed_generator sg0( - host_det, stddevs, 0, fit_cfg0.propagation.context); + host_det, {}, 0, fit_cfg0.propagation.context); traccc::seed_generator sg1( - host_det, stddevs, 0, fit_cfg1.propagation.context); + host_det, {}, 0, fit_cfg1.propagation.context); // Iterate over events for (auto event = input_opts.skip; diff --git a/examples/run/cpu/truth_finding_example.cpp b/examples/run/cpu/truth_finding_example.cpp index 1561850edb..c72bd5fed5 100644 --- a/examples/run/cpu/truth_finding_example.cpp +++ b/examples/run/cpu/truth_finding_example.cpp @@ -90,15 +90,6 @@ int seq_run(const traccc::opts::track_finding& finding_opts, * Do the reconstruction *****************************/ - // Standard deviations for seed track parameters - static constexpr std::array stddevs = - {1e-4f * traccc::unit::mm, - 1e-4f * traccc::unit::mm, - 1e-3f, - 1e-3f, - 1e-4f / traccc::unit::GeV, - 1e-4f * traccc::unit::ns}; - // Propagation configuration detray::propagation::config propagation_config(propagation_opts); @@ -137,7 +128,7 @@ int seq_run(const traccc::opts::track_finding& finding_opts, [&]( const typename detector_traits_t::host& det) { traccc::seed_generator sg( - det, stddevs); + det); evt_data.generate_truth_candidates( truth_track_candidates, truth_measurements, sg, host_mr, truth_finding_opts.m_pT_min); diff --git a/examples/run/cpu/truth_fitting_example.cpp b/examples/run/cpu/truth_fitting_example.cpp index 7b3bb535a0..83abad866f 100644 --- a/examples/run/cpu/truth_fitting_example.cpp +++ b/examples/run/cpu/truth_fitting_example.cpp @@ -91,15 +91,6 @@ int main(int argc, char* argv[]) { * Do the reconstruction *****************************/ - /// Standard deviations for seed track parameters - static constexpr std::array stddevs = { - 0.03f * traccc::unit::mm, - 0.03f * traccc::unit::mm, - 0.017f, - 0.017f, - 0.001f / traccc::unit::GeV, - 1.f * traccc::unit::ns}; - // Fitting algorithm object traccc::fitting_config fit_cfg(fitting_opts); fit_cfg.propagation = propagation_opts; @@ -127,7 +118,7 @@ int main(int argc, char* argv[]) { const typename detector_traits_t::host& det) { // Seed generator traccc::seed_generator sg( - det, stddevs); + det); evt_data.generate_truth_candidates( truth_track_candidates, truth_measurements, sg, host_mr); }); diff --git a/examples/run/cuda/truth_finding_example_cuda.cpp b/examples/run/cuda/truth_finding_example_cuda.cpp index 1e5248dfa7..789cde6ed6 100644 --- a/examples/run/cuda/truth_finding_example_cuda.cpp +++ b/examples/run/cuda/truth_finding_example_cuda.cpp @@ -122,15 +122,6 @@ int seq_run(const traccc::opts::track_finding& finding_opts, traccc::buffer_from_host_detector(polymorphic_detector, device_mr, async_copy); - // Standard deviations for seed track parameters - static constexpr std::array stddevs = - {1e-4f * traccc::unit::mm, - 1e-4f * traccc::unit::mm, - 1e-3f, - 1e-3f, - 1e-4f / traccc::unit::GeV, - 1e-4f * traccc::unit::ns}; - // Propagation configuration detray::propagation::config propagation_config(propagation_opts); @@ -175,7 +166,7 @@ int seq_run(const traccc::opts::track_finding& finding_opts, const typename detector_traits_t::host& det) { // Seed generator traccc::seed_generator sg( - det, stddevs); + det); evt_data.generate_truth_candidates( truth_track_candidates, truth_measurements, sg, host_mr, truth_finding_opts.m_pT_min); diff --git a/examples/run/cuda/truth_fitting_example_cuda.cpp b/examples/run/cuda/truth_fitting_example_cuda.cpp index e7a1ed733d..8f2d66d08d 100644 --- a/examples/run/cuda/truth_fitting_example_cuda.cpp +++ b/examples/run/cuda/truth_fitting_example_cuda.cpp @@ -120,15 +120,6 @@ int main(int argc, char* argv[]) { traccc::buffer_from_host_detector(polymorphic_detector, device_mr, async_copy); - /// Standard deviations for seed track parameters - static constexpr std::array stddevs = { - 0.03f * traccc::unit::mm, - 0.03f * traccc::unit::mm, - 0.017f, - 0.017f, - 0.001f / traccc::unit::GeV, - 1.f * traccc::unit::ns}; - // Fitting algorithm object traccc::fitting_config fit_cfg(fitting_opts); fit_cfg.propagation = propagation_opts; @@ -160,7 +151,7 @@ int main(int argc, char* argv[]) { const typename detector_traits_t::host& det) { // Seed generator traccc::seed_generator sg( - det, stddevs); + det); evt_data.generate_truth_candidates( truth_track_candidates, truth_measurements, sg, host_mr); }); diff --git a/tests/cpu/test_ckf_combinatorics_telescope.cpp b/tests/cpu/test_ckf_combinatorics_telescope.cpp index f02c0584e4..bfb8f6adb1 100644 --- a/tests/cpu/test_ckf_combinatorics_telescope.cpp +++ b/tests/cpu/test_ckf_combinatorics_telescope.cpp @@ -114,8 +114,7 @@ TEST_P(CpuCkfCombinatoricsTelescopeTests, Run) { *****************************/ // Seed generator - seed_generator sg(detector.as(), - stddevs); + seed_generator sg(detector.as()); // Finding algorithm configuration traccc::finding_config cfg_no_limit; diff --git a/tests/cpu/test_ckf_sparse_tracks_telescope.cpp b/tests/cpu/test_ckf_sparse_tracks_telescope.cpp index 8bd7e04e39..2821020094 100644 --- a/tests/cpu/test_ckf_sparse_tracks_telescope.cpp +++ b/tests/cpu/test_ckf_sparse_tracks_telescope.cpp @@ -126,8 +126,7 @@ TEST_P(CkfSparseTrackTelescopeTests, Run) { *****************************/ // Seed generator - seed_generator sg(detector.as(), - stddevs); + seed_generator sg(detector.as()); // Finding algorithm configuration typename traccc::finding_config cfg; diff --git a/tests/cpu/test_kalman_fitter_hole_count.cpp b/tests/cpu/test_kalman_fitter_hole_count.cpp index 4c673d1b49..57a7afdd74 100644 --- a/tests/cpu/test_kalman_fitter_hole_count.cpp +++ b/tests/cpu/test_kalman_fitter_hole_count.cpp @@ -121,8 +121,7 @@ TEST_P(KalmanFittingHoleCountTests, Run) { ***************/ // Seed generator - seed_generator sg(detector.as(), - stddevs); + seed_generator sg(detector.as()); // Fitting algorithm object traccc::fitting_config fit_cfg; diff --git a/tests/cpu/test_kalman_fitter_momentum_resolution.cpp b/tests/cpu/test_kalman_fitter_momentum_resolution.cpp index 314653cb78..467b9e6250 100644 --- a/tests/cpu/test_kalman_fitter_momentum_resolution.cpp +++ b/tests/cpu/test_kalman_fitter_momentum_resolution.cpp @@ -142,8 +142,7 @@ TEST_P(KalmanFittingMomentumResolutionTests, Run) { ***************/ // Seed generator - seed_generator sg(detector.as(), - stddevs); + seed_generator sg(detector.as()); // Fitting algorithm object traccc::fitting_config fit_cfg; diff --git a/tests/cpu/test_kalman_fitter_telescope.cpp b/tests/cpu/test_kalman_fitter_telescope.cpp index a30be320a7..2a65f587f0 100644 --- a/tests/cpu/test_kalman_fitter_telescope.cpp +++ b/tests/cpu/test_kalman_fitter_telescope.cpp @@ -122,8 +122,7 @@ TEST_P(KalmanFittingTelescopeTests, Run) { ***************/ // Seed generator - seed_generator sg(detector.as(), - stddevs); + seed_generator sg(detector.as()); // Fitting algorithm object traccc::fitting_config fit_cfg; diff --git a/tests/cpu/test_kalman_fitter_wire_chamber.cpp b/tests/cpu/test_kalman_fitter_wire_chamber.cpp index 38ac359b18..1f68a2c570 100644 --- a/tests/cpu/test_kalman_fitter_wire_chamber.cpp +++ b/tests/cpu/test_kalman_fitter_wire_chamber.cpp @@ -124,8 +124,7 @@ TEST_P(KalmanFittingWireChamberTests, Run) { ***************/ // Seed generator - seed_generator sg(detector.as(), - stddevs); + seed_generator sg(detector.as()); // Fitting algorithm object traccc::fitting_config fit_cfg; diff --git a/tests/cuda/test_ckf_combinatorics_telescope.cpp b/tests/cuda/test_ckf_combinatorics_telescope.cpp index 4e3a85a2fd..5228ea2c7a 100644 --- a/tests/cuda/test_ckf_combinatorics_telescope.cpp +++ b/tests/cuda/test_ckf_combinatorics_telescope.cpp @@ -131,8 +131,7 @@ TEST_P(CudaCkfCombinatoricsTelescopeTests, Run) { vecmem::cuda::async_copy copy{stream.cudaStream()}; // Seed generator - seed_generator sg(detector.as(), - stddevs); + seed_generator sg(detector.as()); // Finding algorithm configuration typename traccc::cuda::combinatorial_kalman_filter_algorithm::config_type diff --git a/tests/cuda/test_ckf_toy_detector.cpp b/tests/cuda/test_ckf_toy_detector.cpp index fddcbbd837..9d222f9add 100644 --- a/tests/cuda/test_ckf_toy_detector.cpp +++ b/tests/cuda/test_ckf_toy_detector.cpp @@ -133,8 +133,7 @@ TEST_P(CkfToyDetectorTests, Run) { vecmem::cuda::async_copy copy{stream.cudaStream()}; // Seed generator - seed_generator sg(detector.as(), - stddevs); + seed_generator sg(detector.as()); // Finding algorithm configuration typename traccc::cuda::combinatorial_kalman_filter_algorithm::config_type diff --git a/tests/cuda/test_kalman_fitter_telescope.cpp b/tests/cuda/test_kalman_fitter_telescope.cpp index 51dc0966ff..a1d6538a8a 100644 --- a/tests/cuda/test_kalman_fitter_telescope.cpp +++ b/tests/cuda/test_kalman_fitter_telescope.cpp @@ -141,7 +141,7 @@ TEST_P(KalmanFittingTelescopeTests, Run) { // Seed generator seed_generator sg( - polymorphic_detector.as(), stddevs); + polymorphic_detector.as()); // Fitting algorithm object traccc::cuda::kalman_fitting_algorithm::config_type fit_cfg; diff --git a/tests/sycl/test_ckf_combinatorics_telescope.cpp b/tests/sycl/test_ckf_combinatorics_telescope.cpp index bc1113a984..0f40b8fe86 100644 --- a/tests/sycl/test_ckf_combinatorics_telescope.cpp +++ b/tests/sycl/test_ckf_combinatorics_telescope.cpp @@ -132,8 +132,7 @@ TEST_P(CkfCombinatoricsTelescopeTests, Run) { *****************************/ // Seed generator - seed_generator sg(host_detector.as(), - stddevs); + seed_generator sg(host_detector.as()); // Finding algorithm configuration traccc::sycl::combinatorial_kalman_filter_algorithm::config_type diff --git a/tests/sycl/test_ckf_toy_detector.cpp b/tests/sycl/test_ckf_toy_detector.cpp index 4164e8ba49..a453e4019f 100644 --- a/tests/sycl/test_ckf_toy_detector.cpp +++ b/tests/sycl/test_ckf_toy_detector.cpp @@ -131,8 +131,7 @@ TEST_P(CkfToyDetectorTests, Run) { *****************************/ // Seed generator - seed_generator sg(host_detector.as(), - stddevs); + seed_generator sg(host_detector.as()); // Finding algorithm configuration traccc::sycl::combinatorial_kalman_filter_algorithm::config_type cfg; diff --git a/tests/sycl/test_kalman_fitter_telescope.cpp b/tests/sycl/test_kalman_fitter_telescope.cpp index 1bc4fbf49b..96d76a4cab 100644 --- a/tests/sycl/test_kalman_fitter_telescope.cpp +++ b/tests/sycl/test_kalman_fitter_telescope.cpp @@ -146,7 +146,7 @@ TEST_P(KalmanFittingTelescopeTests, Run) { // Seed generator seed_generator sg( - polymorphic_detector.as(), stddevs); + polymorphic_detector.as()); // Fitting algorithm object typename traccc::sycl::kalman_fitting_algorithm::config_type fit_cfg;