cpu/simulation/track_generators.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 // Project include(s)
0010 #include "detray/test/common/track_generators.hpp"
0011
0012 #include "detray/definitions/units.hpp"
0013 #include "detray/tracks/tracks.hpp"
0014 #include "detray/tracks/trajectories.hpp"
0015
0016 // Detray test include(s)
0017 #include "detray/test/framework/types.hpp"
0018 #include "detray/test/utils/statistics.hpp"
0019
0020 // GTest include(s)
0021 #include <gtest/gtest.h>
0022
0023 using namespace detray;
0024
0025 using test_algebra = test::algebra;
0026 using scalar = test::scalar;
0027 using vector3 = test::vector3;
0028 using point3 = test::point3;
0029
0030 GTEST_TEST(detray_simulation, uniform_track_generator) {
0031   using generator_t =
0032       uniform_track_generator<free_track_parameters<test_algebra>>;
0033
0034   constexpr const scalar tol{1e-5f};
0035   constexpr const scalar max_pi{generator_t::configuration::k_max_pi};
0036
0037   constexpr std::size_t phi_steps{50u};
0038   constexpr std::size_t theta_steps{50u};
0039   constexpr const scalar charge{-1.f};
0040
0041   std::array<vector3, phi_steps * theta_steps> momenta{};
0042
0043   // Loop over theta values [0,pi)
0044   for (std::size_t itheta{0u}; itheta < theta_steps; ++itheta) {
0045     const scalar theta{static_cast<scalar>(itheta) * max_pi /
0046                        static_cast<scalar>(theta_steps - 1u)};
0047
0048     // Loop over phi values [-pi, pi)
0049     for (std::size_t iphi{0u}; iphi < phi_steps; ++iphi) {
0050       // The direction
0051       const scalar phi{-constant<scalar>::pi +
0052                        static_cast<scalar>(iphi) *
0053                            (constant<scalar>::pi + max_pi) /
0054                            static_cast<scalar>(phi_steps)};
0055
0056       // initialize a track
0057       vector3 mom{std::cos(phi) * std::sin(theta),
0058                   std::sin(phi) * std::sin(theta), std::cos(theta)};
0059       mom = vector::normalize(mom);
0060       free_track_parameters<test_algebra> traj({0.f, 0.f, 0.f}, 0.f, mom, -1.f);
0061
0062       momenta[itheta * phi_steps + iphi] = traj.mom(-1.f);
0063     }
0064   }
0065
0066   // Now run the track generator and compare
0067   std::size_t n_tracks{0u};
0068   generator_t::configuration trk_gen_cfg{};
0069   trk_gen_cfg.charge(charge);
0070   trk_gen_cfg.phi_steps(phi_steps).theta_steps(theta_steps);
0071
0072   for (const auto track : generator_t(trk_gen_cfg)) {
0073     vector3& expected = momenta[n_tracks];
0074     vector3 result = track.mom(trk_gen_cfg.charge());
0075
0076     // Compare with for loop
0077     EXPECT_NEAR(vector::norm(expected - result), 0.f, tol)
0078         << "Track: \n"
0079         << expected[0] << "\t" << result[0] << "\n"
0080         << expected[1] << "\t" << result[1] << "\n"
0081         << expected[2] << "\t" << result[2] << std::endl;
0082
0083     ++n_tracks;
0084   }
0085   ASSERT_EQ(momenta.size(), n_tracks);
0086
0087   // Generate rays
0088   n_tracks = 0u;
0089   for (const auto r : generator_t(phi_steps, theta_steps)) {
0090     vector3& expected = momenta[n_tracks];
0091     vector3 result = r.dir();
0092
0093     // Compare with for loop
0094     EXPECT_NEAR(vector::norm(expected - result), 0.f, tol)
0095         << "Ray: \n"
0096         << expected[0] << "\t" << result[0] << "\n"
0097         << expected[1] << "\t" << result[1] << "\n"
0098         << expected[2] << "\t" << result[2] << std::endl;
0099
0100     ++n_tracks;
0101   }
0102   ASSERT_EQ(momenta.size(), n_tracks);
0103
0104   // Generate helical trajectories
0105   const vector3 B{0.f * unit<scalar>::T, 0.f * unit<scalar>::T,
0106                   2.f * unit<scalar>::T};
0107   n_tracks = 0u;
0108   for (const auto track : generator_t(phi_steps, theta_steps)) {
0109     detail::helix<test_algebra> helix_traj(track, B);
0110     vector3& expected = momenta[n_tracks];
0111     vector3 result = helix_traj.dir(0.f);
0112
0113     // Compare with for loop
0114     EXPECT_NEAR(vector::norm(expected - result), 0.f, tol)
0115         << "Helix: \n"
0116         << expected[0] << "\t" << result[0] << "\n"
0117         << expected[1] << "\t" << result[1] << "\n"
0118         << expected[2] << "\t" << result[2] << std::endl;
0119
0120     ++n_tracks;
0121   }
0122   ASSERT_EQ(momenta.size(), n_tracks);
0123 }
0124
0125 GTEST_TEST(detray_simulation, uniform_track_generator_eta) {
0126   using generator_t =
0127       uniform_track_generator<free_track_parameters<test_algebra>>;
0128
0129   constexpr const scalar tol{1e-5f};
0130   constexpr const scalar max_pi{generator_t::configuration::k_max_pi};
0131   constexpr const scalar charge{-1.f};
0132
0133   constexpr std::size_t phi_steps{50u};
0134   constexpr std::size_t eta_steps{50u};
0135
0136   std::array<vector3, phi_steps * eta_steps> momenta{};
0137
0138   // Loop over eta values [-5, 5]
0139   for (std::size_t ieta{0u}; ieta < eta_steps; ++ieta) {
0140     const scalar eta{-5.f + static_cast<scalar>(ieta) * 10.f /
0141                                 static_cast<scalar>(eta_steps - 1u)};
0142     const scalar theta{2.f * std::atan(std::exp(-eta))};
0143
0144     // Loop over phi values [-pi, pi)
0145     for (std::size_t iphi{0u}; iphi < phi_steps; ++iphi) {
0146       // The direction
0147       const scalar phi{-constant<scalar>::pi +
0148                        static_cast<scalar>(iphi) *
0149                            (constant<scalar>::pi + max_pi) /
0150                            static_cast<scalar>(phi_steps)};
0151
0152       // initialize a track
0153       vector3 mom{std::cos(phi) * std::sin(theta),
0154                   std::sin(phi) * std::sin(theta), std::cos(theta)};
0155       mom = vector::normalize(mom);
0156       free_track_parameters<test_algebra> traj({0.f, 0.f, 0.f}, 0.f, mom, -1.f);
0157
0158       momenta[ieta * phi_steps + iphi] = traj.mom(charge);
0159     }
0160   }
0161
0162   // Now run the track generator and compare
0163   std::size_t n_tracks{0u};
0164
0165   auto trk_generator = generator_t{};
0166   trk_generator.config().phi_steps(phi_steps).eta_steps(eta_steps);
0167   trk_generator.config().charge(charge);
0168
0169   for (const auto track : trk_generator) {
0170     vector3& expected = momenta[n_tracks];
0171     vector3 result = track.mom(trk_generator.config().charge());
0172
0173     // Compare with for loop
0174     EXPECT_NEAR(vector::norm(expected - result), 0.f, tol)
0175         << "Expected\tResult: \n"
0176         << expected[0] << "\t" << result[0] << "\n"
0177         << expected[1] << "\t" << result[1] << "\n"
0178         << expected[2] << "\t" << result[2] << std::endl;
0179
0180     ++n_tracks;
0181   }
0182   ASSERT_EQ(momenta.size(), n_tracks);
0183 }
0184
0185 GTEST_TEST(detray_simulation, uniform_track_generator_with_range) {
0186   using generator_t =
0187       uniform_track_generator<free_track_parameters<test_algebra>>;
0188
0189   constexpr const scalar tol{1e-5f};
0190
0191   std::vector<std::array<scalar, 2>> theta_phi;
0192
0193   auto trk_gen_cfg = generator_t::configuration{};
0194   trk_gen_cfg.phi_range(-2.f, 2.f).phi_steps(4u);
0195   trk_gen_cfg.theta_range(1.f, 2.f).theta_steps(2u);
0196
0197   for (const auto track : generator_t{trk_gen_cfg}) {
0198     const auto dir = track.dir();
0199     theta_phi.push_back({vector::theta(dir), vector::phi(dir)});
0200   }
0201
0202   EXPECT_EQ(theta_phi.size(), 8u);
0203   EXPECT_NEAR(theta_phi[0][0], 1.f, tol);
0204   EXPECT_NEAR(theta_phi[0][1], -2.f, tol);
0205   EXPECT_NEAR(theta_phi[1][0], 1.f, tol);
0206   EXPECT_NEAR(theta_phi[1][1], -1.f, tol);
0207   EXPECT_NEAR(theta_phi[2][0], 1.f, tol);
0208   EXPECT_NEAR(theta_phi[2][1], 0.f, tol);
0209   EXPECT_NEAR(theta_phi[3][0], 1.f, tol);
0210   EXPECT_NEAR(theta_phi[3][1], 1.f, tol);
0211   EXPECT_NEAR(theta_phi[4][0], 2.f, tol);
0212   EXPECT_NEAR(theta_phi[4][1], -2.f, tol);
0213   EXPECT_NEAR(theta_phi[5][0], 2.f, tol);
0214   EXPECT_NEAR(theta_phi[5][1], -1.f, tol);
0215   EXPECT_NEAR(theta_phi[6][0], 2.f, tol);
0216   EXPECT_NEAR(theta_phi[6][1], 0.f, tol);
0217   EXPECT_NEAR(theta_phi[7][0], 2.f, tol);
0218   EXPECT_NEAR(theta_phi[7][1], 1.f, tol);
0219 }
0220
0221 /// Tests a random number based track state generator - uniform distribution
0222 GTEST_TEST(detray_simulation, random_track_generator_uniform) {
0223   // Use deterministic random number generator for testing
0224   using uniform_gen_t =
0225       detail::random_numbers<scalar, std::uniform_real_distribution<scalar>>;
0226   using trk_generator_t =
0227       random_track_generator<free_track_parameters<test_algebra>,
0228                              uniform_gen_t>;
0229
0230   // Tolerance depends on sample size
0231   constexpr scalar tol{0.02f};
0232
0233   // Track counter
0234   std::size_t n_tracks{0u};
0235   constexpr std::size_t n_gen_tracks{10000u};
0236
0237   // Other params
0238   trk_generator_t::configuration trk_gen_cfg{};
0239   trk_gen_cfg.n_tracks(n_gen_tracks);
0240   trk_gen_cfg.seed(42u);
0241   trk_gen_cfg.phi_range(-0.9f * constant<scalar>::pi,
0242                         0.8f * constant<scalar>::pi);
0243   trk_gen_cfg.eta_range(-4.f, 4.f);
0244   trk_gen_cfg.mom_range(1.f * unit<scalar>::GeV, 2.f * unit<scalar>::GeV);
0245   trk_gen_cfg.origin_stddev(0.1f * unit<scalar>::mm, 0.f * unit<scalar>::mm,
0246                             0.2f * unit<scalar>::mm);
0247
0248   // Catch the results
0249   std::array<scalar, n_gen_tracks> x{};
0250   std::array<scalar, n_gen_tracks> y{};
0251   std::array<scalar, n_gen_tracks> z{};
0252   std::array<scalar, n_gen_tracks> mom{};
0253   std::array<scalar, n_gen_tracks> phi{};
0254   std::array<scalar, n_gen_tracks> theta{};
0255
0256   for (const auto track : trk_generator_t{trk_gen_cfg}) {
0257     const auto pos = track.pos();
0258     x[n_tracks] = pos[0];
0259     y[n_tracks] = pos[1];
0260     z[n_tracks] = pos[2];
0261     mom[n_tracks] = track.p(trk_gen_cfg.charge());
0262     phi[n_tracks] = vector::phi(track.dir());
0263     theta[n_tracks] = vector::theta(track.dir());
0264
0265     ++n_tracks;
0266   }
0267
0268   ASSERT_EQ(n_gen_tracks, n_tracks);
0269
0270   // Check uniform distribution
0271   const auto& ori = trk_gen_cfg.origin();
0272   const auto& ori_stddev = trk_gen_cfg.origin_stddev();
0273   const auto& phi_range = trk_gen_cfg.phi_range();
0274   const auto& theta_range = trk_gen_cfg.theta_range();
0275   ASSERT_NEAR(theta_range[0], 0.0366f, tol);
0276   ASSERT_NEAR(theta_range[1], 3.105f, tol);
0277   const auto& mom_range = trk_gen_cfg.mom_range();
0278
0279   // Mean
0280   EXPECT_NEAR(statistics::mean(x), ori[0], tol);
0281   EXPECT_NEAR(statistics::mean(y), ori[1], tol);
0282   EXPECT_NEAR(statistics::mean(z), ori[2], tol);
0283   EXPECT_NEAR(statistics::mean(mom), 0.5f * (mom_range[0] + mom_range[1]), tol);
0284   EXPECT_NEAR(statistics::mean(phi), 0.5f * (phi_range[0] + phi_range[1]), tol);
0285   EXPECT_NEAR(statistics::mean(theta), 0.5f * (theta_range[0] + theta_range[1]),
0286               tol);
0287
0288   // variance
0289   EXPECT_NEAR(statistics::variance(x), ori_stddev[0] * ori_stddev[0], tol);
0290   EXPECT_NEAR(statistics::variance(y), ori_stddev[1] * ori_stddev[1], tol);
0291   EXPECT_NEAR(statistics::variance(z), ori_stddev[2] * ori_stddev[2], tol);
0292   EXPECT_NEAR(statistics::variance(mom),
0293               1.0f / 12.0f * std::pow(mom_range[1] - mom_range[0], 2.f), tol);
0294   EXPECT_NEAR(statistics::variance(phi),
0295               1.0f / 12.0f * std::pow(phi_range[1] - phi_range[0], 2.f), tol);
0296   EXPECT_NEAR(statistics::variance(theta),
0297               1.0f / 12.0f * std::pow(theta_range[1] - theta_range[0], 2.f),
0298               tol);
0299 }
0300
0301 /// Tests a random number based track state generator - normal distribution
0302 GTEST_TEST(detray_simulation, random_track_generator_normal) {
0303   // Use deterministic random number generator for testing
0304   using normal_gen_t =
0305       detail::random_numbers<scalar, std::normal_distribution<scalar>>;
0306   using trk_generator_t =
0307       random_track_generator<free_track_parameters<test_algebra>, normal_gen_t>;
0308
0309   // Tolerance depends on sample size
0310   constexpr scalar tol{0.02f};
0311
0312   // Track counter
0313   std::size_t n_tracks{0u};
0314   constexpr std::size_t n_gen_tracks{10000u};
0315
0316   // Other params
0317   trk_generator_t::configuration trk_gen_cfg{};
0318   trk_gen_cfg.n_tracks(n_gen_tracks);
0319   trk_gen_cfg.seed(42u);
0320   trk_gen_cfg.phi_range(-0.9f * constant<scalar>::pi,
0321                         0.8f * constant<scalar>::pi);
0322   trk_gen_cfg.eta_range(-4.f, 4.f);
0323   trk_gen_cfg.mom_range(1.f * unit<scalar>::GeV, 2.f * unit<scalar>::GeV);
0324   trk_gen_cfg.origin(0.f, 0.f, 0.f);
0325   trk_gen_cfg.origin_stddev(0.1f * unit<scalar>::mm, 0.5f * unit<scalar>::mm,
0326                             0.3f * unit<scalar>::mm);
0327
0328   // Catch the results
0329   std::array<scalar, n_gen_tracks> x{};
0330   std::array<scalar, n_gen_tracks> y{};
0331   std::array<scalar, n_gen_tracks> z{};
0332   std::array<scalar, n_gen_tracks> mom{};
0333   std::array<scalar, n_gen_tracks> phi{};
0334   std::array<scalar, n_gen_tracks> theta{};
0335
0336   for (const auto track : trk_generator_t{trk_gen_cfg}) {
0337     const auto pos = track.pos();
0338     x[n_tracks] = pos[0];
0339     y[n_tracks] = pos[1];
0340     z[n_tracks] = pos[2];
0341     mom[n_tracks] = track.p(trk_gen_cfg.charge());
0342     phi[n_tracks] = vector::phi(track.dir());
0343     theta[n_tracks] = vector::theta(track.dir());
0344
0345     ++n_tracks;
0346   }
0347
0348   ASSERT_EQ(n_gen_tracks, n_tracks);
0349
0350   // check gaussian distribution - values are clamped to phi/theta range
0351   const auto& ori = trk_gen_cfg.origin();
0352   const auto& ori_stddev = trk_gen_cfg.origin_stddev();
0353   const auto& phi_range = trk_gen_cfg.phi_range();
0354   const auto& theta_range = trk_gen_cfg.theta_range();
0355   ASSERT_NEAR(theta_range[0], 0.0366f, tol);
0356   ASSERT_NEAR(theta_range[1], 3.105f, tol);
0357   const auto& mom_range = trk_gen_cfg.mom_range();
0358
0359   // Mean
0360   EXPECT_NEAR(statistics::mean(x), ori[0], tol);
0361   EXPECT_NEAR(statistics::mean(y), ori[1], tol);
0362   EXPECT_NEAR(statistics::mean(z), ori[2], tol);
0363   EXPECT_NEAR(statistics::mean(mom),
0364               mom_range[0] + 0.5f * (mom_range[1] - mom_range[0]), tol);
0365   EXPECT_NEAR(statistics::mean(phi),
0366               phi_range[0] + 0.5f * (phi_range[1] - phi_range[0]), tol);
0367   EXPECT_NEAR(statistics::mean(theta),
0368               theta_range[0] + 0.5f * (theta_range[1] - theta_range[0]), tol);
0369
0370   // variance
0371   EXPECT_NEAR(statistics::variance(x), ori_stddev[0] * ori_stddev[0], tol);
0372   EXPECT_NEAR(statistics::variance(y), ori_stddev[1] * ori_stddev[1], tol);
0373   EXPECT_NEAR(statistics::variance(z), ori_stddev[2] * ori_stddev[2], tol);
0374   EXPECT_NEAR(statistics::variance(mom),
0375               std::pow(0.5f / 3.0f * (mom_range[1] - mom_range[0]), 2.f), tol);
0376   EXPECT_NEAR(statistics::variance(phi),
0377               std::pow(0.5f / 3.0f * (phi_range[1] - phi_range[0]), 2.f), tol);
0378   EXPECT_NEAR(statistics::variance(theta),
0379               std::pow(0.5f / 3.0f * (theta_range[1] - theta_range[0]), 2.f),
0380               tol);
0381 }