1+ /* * TRACCC library, part of the ACTS project (R&D line)
2+ *
3+ * (c) 2025 CERN for the benefit of the ACTS project
4+ *
5+ * Mozilla Public License Version 2.0
6+ */
7+
8+ #include < cmath>
9+
10+ // Local include(s).
11+ #include " kalman_fitting_momentum_resolution_test.hpp"
12+
13+ // ROOT include(s).
14+ #ifdef TRACCC_HAVE_ROOT
15+ #include < TF1.h>
16+ #include < TFile.h>
17+ #include < TFitResult.h>
18+ #include < TFitResultPtr.h>
19+ #include < TH1.h>
20+ #endif // TRACCC_HAVE_ROOT
21+
22+ // System include(s).
23+ #include < iostream>
24+ #include < memory>
25+ #include < stdexcept>
26+
27+ namespace traccc {
28+
29+ void KalmanFittingMomentumResolutionTests::consistency_tests (
30+ const track_state_collection_types::host& track_states_per_track) const {
31+
32+ // The nubmer of track states is supposed be equal to the number
33+ // of planes
34+ ASSERT_EQ (track_states_per_track.size (), std::get<11 >(GetParam ()));
35+ }
36+
37+ void KalmanFittingMomentumResolutionTests::momentum_resolution_tests ([
38+ [maybe_unused]] std::string_view file_name) const {
39+
40+ #ifdef TRACCC_HAVE_ROOT
41+
42+ // Open the file with the histograms.
43+ std::unique_ptr<TFile> ifile (TFile::Open (file_name.data (), " READ"
44+ if ((!ifile) || ifile->IsZombie ()) {
45+ throw std::runtime_error (std::string (" Could not open file \" "
46+ file_name.data () + " \" "
47+ }
48+
49+ // Access the histogram.
50+ TH1* residual_qopT_hist = dynamic_cast <TH1*>(ifile->Get (" res_qopT"
51+ if (!residual_qopT_hist) {
52+ throw std::runtime_error (
53+ std::string (" Could not access histogram residual_qopT in file \" "
54+ file_name.data () + " \" "
55+ }
56+
57+ // Function used for the fit.
58+ TF1 gaus{" gaus" " gaus" 5 , 5 };
59+ double fit_par[3 ];
60+
61+ // Set the mean seed to 0
62+ gaus.SetParameters (1 , 0 .);
63+ gaus.SetParLimits (1 , -1 ., 1 .);
64+
65+ // Set the standard deviation seed to 0.01
66+ gaus.SetParameters (2 , 0 .01f );
67+ gaus.SetParLimits (2 , 0 .f , 0 .1f );
68+
69+ auto res = residual_qopT_hist->Fit (" gaus" " Q0S"
70+
71+ gaus.GetParameters (&fit_par[0 ]);
72+
73+ // Mean check
74+ EXPECT_NEAR (fit_par[1 ], 0 , 0.05 ) << " Residual qopT mean value error"
75+
76+ // Expected momentum resolution
77+ const std::size_t N = std::get<11 >(GetParam ());
78+ const auto smearing = std::get<15 >(GetParam ());
79+ const scalar epsilon = smearing[0u ];
80+ const scalar Bz = std::get<13 >(GetParam ())[2u ];
81+ const scalar spacing = std::get<12 >(GetParam ());
82+ const scalar L = static_cast <scalar>(N - 1u ) * spacing;
83+
84+ // Curvature (1/helix_radius) resolution from detector noise Eq. (35.61) of
85+ // PDG 2024
86+ const scalar dkres =
87+ epsilon / (L * L) * math::sqrt (720 .f / static_cast <scalar>(N + 4 ));
88+
89+ // Curvature (1/helix_radius) resolution from multiple scattering Eq.
90+ // (35.63) of PDG 2024
91+ // @NOTE: The calculation of the multiple scattering term is in work in
92+ // progress. Currently we are validating the momentum resolution without any
93+ // material on the modules, therefore, the multiple scattering contribution
94+ // to the momentum resolution is set to zero.
95+ const scalar dkms = 0 .f ;
96+
97+ // Eq. (35.60) of PDG 2024
98+ const scalar dk = math::sqrt (dkres * dkres + dkms * dkms);
99+
100+ // σ(q/pT) = σ(k/B) = σ(k)/B
101+ const scalar expected_sigma_qopT = dk / Bz;
102+
103+ // Check if the standard deviation of the qopT residuals is within the
104+ // theoretically expected range.
105+ EXPECT_NEAR (fit_par[2 ], expected_sigma_qopT, expected_sigma_qopT * 0 .05f )
106+ << " Residual qopT sigma value error"
107+
108+ #else
109+ std::cout << " Pull value tests not performed without ROOT"
110+ #endif // TRACCC_HAVE_ROOT
111+
112+ return ;
113+ }
114+
115+ void KalmanFittingMomentumResolutionTests::SetUp () {
116+
117+ vecmem::host_memory_resource host_mr;
118+
119+ const scalar offset = std::get<10 >(GetParam ());
120+ const unsigned int n_planes = std::get<11 >(GetParam ());
121+ const scalar spacing = std::get<12 >(GetParam ());
122+
123+ std::vector<scalar> plane_positions;
124+ for (unsigned int i = 0 ; i < n_planes; i++) {
125+ plane_positions.push_back (offset * unit<scalar>::mm +
126+ static_cast <scalar>(i) * spacing *
127+ unit<scalar>::mm);
128+ }
129+
130+ // / Plane alignment direction (aligned to x-axis)
131+ const vector3 bfield = std::get<13 >(GetParam ());
132+
133+ const auto p = std::get<3 >(GetParam ());
134+ const detray::detail::helix<traccc::default_algebra> hlx (
135+ {0 , 0 , 0 }, 0 , {1 , 0 , 0 }, -1 .f / p, &bfield);
136+
137+ constexpr scalar rect_half_length = 500 .f ;
138+ constexpr detray::mask<detray::rectangle2D, traccc::default_algebra>
139+ rectangle{0u , rect_half_length, rect_half_length};
140+
141+ detray::tel_det_config<default_algebra, detray::rectangle2D,
142+ detray::detail::helix>
143+ tel_cfg (rectangle, hlx);
144+ tel_cfg.positions (plane_positions);
145+ tel_cfg.module_material (std::get<14 >(GetParam ()));
146+ tel_cfg.mat_thickness (thickness);
147+ tel_cfg.envelope (rect_half_length);
148+
149+ // Create telescope detector
150+ auto [det, name_map] = build_telescope_detector (host_mr, tel_cfg);
151+
152+ // Write detector file
153+ auto writer_cfg = detray::io::detector_writer_config{}
154+ .format (detray::io::format::json)
155+ .replace_files (true )
156+ .write_material (true )
157+ .path (std::get<0 >(GetParam ()));
158+ detray::io::write_detector (det, name_map, writer_cfg);
159+ }
160+
161+ } // namespace traccc
0 commit comments