1010// Library include(s).
1111#include " traccc/edm/track_parameters.hpp"
1212#include " traccc/utils/particle.hpp"
13+ #include " traccc/utils/trigonometric_helpers.hpp"
1314
1415// detray include(s).
1516#include < detray/geometry/barcode.hpp>
@@ -28,14 +29,53 @@ struct seed_generator {
2829 using algebra_type = typename detector_t ::algebra_type;
2930 using ctx_t = typename detector_t ::geometry_context;
3031
32+ // / Configure the seed generator
33+ struct config {
34+ // Smearing parameters
35+ // / Constant term of the loc0 resolution.
36+ double sigma_loc0 = 20 * traccc::unit<scalar>::um;
37+ // / Pt-dependent loc0 resolution of the form sigma_loc0 =
38+ // / A*exp(-1.*abs(B)*pt).
39+ double sigma_loc0_pT_a = 30 * traccc::unit<scalar>::um;
40+ double sigma_loc0_pT_b = 0.3 / traccc::unit<scalar>::GeV;
41+ // / Constant term of the loc1 resolution.
42+ double sigma_loc1 = 20 * traccc::unit<scalar>::um;
43+ // / Pt-dependent loc1 resolution of the form sigma_loc1 =
44+ // / A*exp(-1.*abs(B)*pt).
45+ double sigma_loc1_pT_a = 30 * traccc::unit<scalar>::um;
46+ double sigma_loc1_pT_b = 0.3 / traccc::unit<scalar>::GeV;
47+ // / Time resolution.
48+ double sigma_time = 1 * traccc::unit<scalar>::ns;
49+ // / Phi angular resolution.
50+ double sigma_phi = 1 * traccc::unit<scalar>::degree;
51+ // / Theta angular resolution.
52+ double sigma_theta = 1 * traccc::unit<scalar>::degree;
53+ // / Relative transverse momentum resolution.
54+ double sigma_pT_rel = 0.05 ;
55+
56+ // / Optional. Initial sigmas for the track parameters which overwrites
57+ // / the smearing params if set.
58+ std::optional<std::array<double , e_bound_size>> initial_sigmas;
59+ // / Initial sigma(q/pt) for the track parameters.
60+ // / @note The resulting q/p sigma is added to the one in `initialSigmas`
61+ double initialsigma_qopt =
62+ 0.1 * traccc::unit<scalar>::e / traccc::unit<scalar>::GeV;
63+ // / Initial sigma(pt)/pt for the track parameters.
64+ // / @note The resulting q/p sigma is added to the one in `initialSigmas`
65+ double initialsigma_pT_rel = 0.1 ;
66+
67+ // / Inflation factors for the variances
68+ std::array<scalar, e_bound_size> cov_inflation{1 ., 1 ., 1 ., 1 ., 1 ., 1 .};
69+ };
70+
3171 // / Constructor with detector
3272 // /
3373 // / @param det input detector
3474 // / @param stddevs standard deviations for parameter smearing
35- seed_generator (const detector_t & det,
36- const std::array<scalar, e_bound_size>& stddevs ,
37- const std:: size_t sd = 0 , ctx_t ctx = {})
38- : m_detector(det), m_stddevs(stddevs ), m_ctx(ctx ) {
75+ seed_generator (const detector_t & det, const config& cfg = {},
76+ const std::size_t sd = std::mt19937::default_seed ,
77+ ctx_t ctx = {})
78+ : m_detector(det), m_ctx(ctx ), m_cfg(cfg ) {
3979 m_generator.seed (static_cast <std::mt19937::result_type>(sd));
4080 }
4181
@@ -73,32 +113,97 @@ struct seed_generator {
73113 sf);
74114 }
75115
76- for (std::size_t i = 0 ; i < e_bound_size; i++) {
116+ // Call the smearing operator
117+ (*this )(bound_param, ptc_type);
77118
78- if (m_stddevs[i] != scalar{0 }) {
79- bound_param[i] = std::normal_distribution<scalar>(
80- bound_param[i], m_stddevs[i])(m_generator);
81- }
119+ return bound_param;
120+ }
82121
83- getter::element (bound_param.covariance (), i, i) =
84- m_stddevs[i] * m_stddevs[i];
122+ // / Seed generator operation
123+ // /
124+ // / @param vertex vertex of particle
125+ // / @param stddevs standard deviations for track parameter smearing
126+ void operator ()(bound_track_parameters<algebra_type>& bound_param,
127+ const traccc::pdg_particle<scalar>& ptc_type) {
128+
129+ const scalar q{ptc_type.charge ()};
130+ const auto pos = bound_param.bound_local ();
131+ const auto time = bound_param.time ();
132+ const auto phi = bound_param.phi ();
133+ const auto theta = bound_param.theta ();
134+ const auto pt = bound_param.pT (q);
135+ const auto qop = bound_param.qop ();
136+
137+ // Compute momentum-dependent resolutions
138+ const double sigma_loc_0 =
139+ m_cfg.sigma_loc0 +
140+ m_cfg.sigma_loc0_pT_a *
141+ math::exp (-1.0 * math::fabs (m_cfg.sigma_loc0_pT_b ) * pt);
142+ const double sigma_loc_1 =
143+ m_cfg.sigma_loc1 +
144+ m_cfg.sigma_loc1_pT_a *
145+ math::exp (-1.0 * math::fabs (m_cfg.sigma_loc1_pT_b ) * pt);
146+ // Shortcuts for other resolutions
147+ const double sigma_qop = math::sqrt (
148+ math::pow (m_cfg.sigma_pT_rel * qop, 2 ) +
149+ math::pow (m_cfg.sigma_theta * (qop * math::tan (theta)), 2 ));
150+
151+ // Smear the position/time
152+ // Note that we smear d0 and z0 in the perigee frame
153+ const scalar smeared_loc0 =
154+ std::normal_distribution<scalar>(pos[0 ], sigma_loc_0)(m_generator);
155+ const scalar smeared_loc1 =
156+ std::normal_distribution<scalar>(pos[1 ], sigma_loc_1)(m_generator);
157+ bound_param.set_bound_local ({smeared_loc0, smeared_loc1});
158+
159+ // Time
160+ bound_param.set_time (std::normal_distribution<scalar>(
161+ time, m_cfg.sigma_time )(m_generator));
162+
163+ // Smear direction angles phi,theta ensuring correct bounds
164+ const scalar smeared_phi =
165+ std::normal_distribution<scalar>(phi, m_cfg.sigma_phi )(m_generator);
166+ const scalar smeared_theta = std::normal_distribution<scalar>(
167+ theta, m_cfg.sigma_theta )(m_generator);
168+ const auto [new_phi, new_theta] =
169+ detail::wrap_phi_theta (smeared_phi, smeared_theta);
170+ bound_param.set_phi (new_phi);
171+ bound_param.set_theta (new_theta);
172+
173+ // Compute smeared q/p
174+ bound_param.set_qop (
175+ std::normal_distribution<scalar>(qop, sigma_qop)(m_generator));
176+
177+ // Build the track covariance matrix using the smearing sigmas
178+ std::array<scalar, e_bound_size> sigmas{
179+ sigma_loc_0, sigma_loc_1, m_cfg.sigma_phi ,
180+ m_cfg.sigma_theta , sigma_qop, m_cfg.sigma_time };
181+
182+ for (std::size_t i = e_bound_loc0; i < e_bound_size; ++i) {
183+ double sigma = sigmas[i];
184+ double variance = sigma * sigma;
185+
186+ // Inflate the initial covariance
187+ variance *= m_cfg.cov_inflation [i];
188+
189+ getter::element (bound_param.covariance (), i, i) = variance;
85190 }
86191
87192 assert (!bound_param.is_invalid ());
88-
89- return bound_param;
90193 }
91194
92195 private:
93- // Random generator
196+ // / Random generator
94197 std::random_device m_rd{};
95198 std::mt19937 m_generator{m_rd ()};
96199
97- // Detector object
200+ // / Detector object
98201 const detector_t & m_detector;
99- // / Standard deviations for parameter smearing
100- std::array<scalar, e_bound_size> m_stddevs;
202+ // / Geometry context
101203 ctx_t m_ctx;
204+
205+ // / Seed generator configuration
206+ config m_cfg{};
102207};
103208
104209} // namespace traccc
0 commit comments