tests/cpu/test_simulation.cpp

0001 /** TRACCC library, part of the ACTS project (R&D line)
0002  *
0003  * (c) 2023-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 "tests/test_detectors.hpp"
0010 #include "traccc/bfield/construct_const_bfield.hpp"
0011 #include "traccc/edm/track_parameters.hpp"
0012 #include "traccc/io/csv/make_hit_reader.hpp"
0013 #include "traccc/io/csv/make_measurement_hit_id_reader.hpp"
0014 #include "traccc/io/csv/make_measurement_reader.hpp"
0015 #include "traccc/io/csv/make_particle_reader.hpp"
0016 #include "traccc/simulation/event_generators.hpp"
0017 #include "traccc/simulation/simulator.hpp"
0018
0019 // Detray include(s).
0020 #include <detray/geometry/mask.hpp>
0021 #include <detray/geometry/shapes/line.hpp>
0022 #include <detray/geometry/shapes/rectangle2D.hpp>
0023 #include <detray/geometry/tracking_surface.hpp>
0024 #include <detray/test/utils/statistics.hpp>
0025
0026 // VecMem include(s).
0027 #include <vecmem/memory/host_memory_resource.hpp>
0028
0029 // GTest include(s).
0030 #include <gtest/gtest.h>
0031
0032 // System include(s).
0033 #include <filesystem>
0034
0035 using namespace traccc;
0036
0037 constexpr scalar tol{1e-7f};
0038
0039 TEST(traccc_simulation, simulation) {
0040   using line_t = detray::mask<detray::line<false>, traccc::default_algebra>;
0041   using rectangle_t =
0042       detray::mask<detray::rectangle2D, traccc::default_algebra>;
0043
0044   traccc::bound_track_parameters<traccc::default_algebra> bound_params{};
0045   bound_params.set_bound_local({1.f, 2.f});
0046
0047   measurement_smearer<traccc::default_algebra> smearer(0.f, 0.f);
0048
0049   traccc::io::csv::measurement iomeas1;
0050   smearer.template operator()<line_t>({-3.f, 2.f}, bound_params, iomeas1);
0051   ASSERT_NEAR(iomeas1.local0, 0.f, tol);
0052   ASSERT_NEAR(iomeas1.local1, 0.f, tol);
0053
0054   traccc::io::csv::measurement iomeas2;
0055   smearer.template operator()<line_t>({2.f, -5.f}, bound_params, iomeas2);
0056   ASSERT_NEAR(iomeas2.local0, 3.f, tol);
0057   ASSERT_NEAR(iomeas2.local1, 0.f, tol);
0058
0059   traccc::io::csv::measurement iomeas3;
0060   smearer.template operator()<rectangle_t>({2.f, -5.f}, bound_params, iomeas3);
0061   ASSERT_NEAR(iomeas3.local0, 3.f, tol);
0062   ASSERT_NEAR(iomeas3.local1, -3.f, tol);
0063 }
0064
0065 GTEST_TEST(traccc_simulation, toy_detector_simulation) {
0066   // Create geometry
0067   vecmem::host_memory_resource host_mr;
0068
0069   // Create B field
0070   using b_field_t = covfie::field<traccc::const_bfield_backend_t<scalar>>;
0071   const vector3 B{0.f, 0.f, 2.f * traccc::unit<scalar>::T};
0072   b_field_t field =
0073       traccc::construct_const_bfield(B)
0074           .as_field<traccc::const_bfield_backend_t<traccc::scalar>>();
0075
0076   // Create geometry
0077   detray::toy_det_config<scalar> toy_cfg{};
0078   const auto [detector, names] =
0079       detray::build_toy_detector<traccc::default_algebra>(host_mr, toy_cfg);
0080
0081   using geo_cxt_t = typename decltype(detector)::geometry_context;
0082   const geo_cxt_t ctx{};
0083
0084   // Create track generator
0085   using uniform_gen_t =
0086       detray::detail::random_numbers<scalar,
0087                                      std::uniform_real_distribution<scalar>>;
0088   using generator_type =
0089       detray::random_track_generator<traccc::free_track_parameters<>,
0090                                      uniform_gen_t>;
0091   generator_type::configuration gen_cfg{};
0092   constexpr unsigned int n_tracks{2500u};
0093   const vector3 ori{0.f, 0.f, 0.f};
0094   gen_cfg.n_tracks(n_tracks);
0095   gen_cfg.origin(ori);
0096   // @TODO The simulator sometimes gets stuck for lower momentum
0097   gen_cfg.p_tot(5.f * traccc::unit<scalar>::GeV);
0098   generator_type generator(gen_cfg);
0099
0100   // Create smearer
0101   measurement_smearer<traccc::default_algebra> smearer(
0102       67.f * traccc::unit<scalar>::um, 170.f * traccc::unit<scalar>::um);
0103
0104   std::size_t n_events{10u};
0105
0106   using detector_type = decltype(detector);
0107   using writer_type =
0108       smearing_writer<measurement_smearer<traccc::default_algebra>>;
0109
0110   typename writer_type::config writer_cfg{smearer};
0111
0112   auto sim = simulator<detector_type, b_field_t, generator_type, writer_type>(
0113       traccc::muon<scalar>(), n_events, detector, field, std::move(generator),
0114       std::move(writer_cfg));
0115
0116   // Lift step size constraints
0117   sim.get_config().propagation.stepping.step_constraint =
0118       std::numeric_limits<float>::max();
0119   sim.get_config().propagation.navigation.search_window = {3u, 3u};
0120
0121   // Do the simulation
0122   sim.run();
0123
0124   for (std::size_t i_event = 0u; i_event < n_events; i_event++) {
0125     std::vector<traccc::io::csv::particle> particles;
0126     auto particle_reader = traccc::io::csv::make_particle_reader(
0127         traccc::io::get_event_filename(i_event, "-particles_initial.csv"));
0128     traccc::io::csv::particle io_particle;
0129     while (particle_reader.read(io_particle)) {
0130       particles.push_back(io_particle);
0131     }
0132
0133     std::vector<traccc::io::csv::hit> hits;
0134     auto hit_reader = traccc::io::csv::make_hit_reader(
0135         traccc::io::get_event_filename(i_event, "-hits.csv"));
0136     traccc::io::csv::hit io_hit;
0137     while (hit_reader.read(io_hit)) {
0138       hits.push_back(io_hit);
0139     }
0140
0141     std::vector<traccc::io::csv::measurement> measurements;
0142     auto measurement_reader = traccc::io::csv::make_measurement_reader(
0143         traccc::io::get_event_filename(i_event, "-measurements.csv"));
0144     traccc::io::csv::measurement io_measurement;
0145     while (measurement_reader.read(io_measurement)) {
0146       measurements.push_back(io_measurement);
0147     }
0148
0149     std::vector<traccc::io::csv::measurement_hit_id> meas_hit_ids;
0150     auto measurement_hit_id_reader =
0151         traccc::io::csv::make_measurement_hit_id_reader(
0152             traccc::io::get_event_filename(i_event,
0153                                            "-measurement-simhit-map.csv"));
0154     traccc::io::csv::measurement_hit_id io_meas_hit_id;
0155     while (measurement_hit_id_reader.read(io_meas_hit_id)) {
0156       meas_hit_ids.push_back(io_meas_hit_id);
0157     }
0158
0159     ASSERT_EQ(particles.size(), n_tracks);
0160     ASSERT_TRUE(not measurements.empty());
0161     ASSERT_EQ(hits.size(), measurements.size());
0162     ASSERT_EQ(hits.size(), meas_hit_ids.size());
0163
0164     // Let's check if measurement smearing works correctly...
0165     std::vector<scalar> local0_diff;
0166     std::vector<scalar> local1_diff;
0167
0168     const std::size_t nhits = hits.size();
0169     for (std::size_t i = 0u; i < nhits; i++) {
0170       const point3 pos{hits[i].tx, hits[i].ty, hits[i].tz};
0171       const vector3 mom{hits[i].tpx, hits[i].tpy, hits[i].tpz};
0172       const auto truth_local =
0173           detray::tracking_surface{
0174               detector, detray::geometry::identifier(hits[i].geometry_id)}
0175               .global_to_local(ctx, pos, vector::normalize(mom));
0176
0177       local0_diff.push_back(truth_local[0] - measurements[i].local0);
0178       local1_diff.push_back(truth_local[1] - measurements[i].local1);
0179
0180       ASSERT_EQ(measurements[i].geometry_id, hits[i].geometry_id);
0181       ASSERT_EQ(meas_hit_ids[i].hit_id, i);
0182       ASSERT_EQ(meas_hit_ids[i].measurement_id, i);
0183     }
0184
0185     const auto var0 = detray::statistics::variance(local0_diff);
0186     const auto var1 = detray::statistics::variance(local1_diff);
0187
0188     EXPECT_NEAR((std::sqrt(var0) - smearer.stddev[0]) / smearer.stddev[0], 0.f,
0189                 0.1f);
0190     EXPECT_NEAR((std::sqrt(var1) - smearer.stddev[1]) / smearer.stddev[1], 0.f,
0191                 0.1f);
0192   }
0193 }
0194
0195 // Test parameters: <initial momentum, theta direction, charge>
0196 class TelescopeDetectorSimulation
0197     : public ::testing::TestWithParam<
0198           std::tuple<std::string, scalar, scalar, scalar>> {};
0199
0200 TEST_P(TelescopeDetectorSimulation, telescope_detector_simulation) {
0201   // Create geometry
0202   vecmem::host_memory_resource host_mr;
0203
0204   // Build from given module positions
0205   std::vector<scalar> positions = {0.f,   50.f,  100.f, 150.f, 200.f, 250.f,
0206                                    300.f, 350.f, 400.f, 450.f, 500.f};
0207
0208   // A thickness larger than 0.1 cm will flip the track direction of low
0209   // energy (or non-relativistic) particle due to the large scattering
0210   const scalar thickness = 0.005f * traccc::unit<scalar>::cm;
0211
0212   detray::tel_det_config<traccc::default_algebra, detray::rectangle2D> tel_cfg{
0213       1000.f * traccc::unit<scalar>::mm, 1000.f * traccc::unit<scalar>::mm};
0214   tel_cfg.positions(positions).mat_thickness(thickness);
0215
0216   const auto [detector, names] =
0217       detray::build_telescope_detector(host_mr, tel_cfg);
0218
0219   // Directory name
0220   const std::string directory = std::get<0>(GetParam()) + "/";
0221   std::filesystem::create_directory(directory);
0222
0223   // Field
0224   using b_field_t = covfie::field<traccc::const_bfield_backend_t<scalar>>;
0225   const vector3 B{0.f, 0.f, 2.f * traccc::unit<scalar>::T};
0226   b_field_t field =
0227       traccc::construct_const_bfield(B)
0228           .as_field<traccc::const_bfield_backend_t<traccc::scalar>>();
0229
0230   // Momentum
0231   const scalar mom = std::get<1>(GetParam());
0232
0233   // Create track generator
0234   constexpr unsigned int theta_steps{1u};
0235   constexpr unsigned int phi_steps{1u};
0236   const vector3 ori{0.f, 0.f, 0.f};
0237   const scalar theta = std::get<2>(GetParam());
0238
0239   const scalar charge = std::get<3>(GetParam());
0240
0241   // Track generator
0242   using generator_type =
0243       detray::uniform_track_generator<traccc::free_track_parameters<>>;
0244   generator_type::configuration gen_cfg{};
0245   gen_cfg.theta_steps(theta_steps);
0246   gen_cfg.phi_steps(phi_steps);
0247   gen_cfg.origin(ori);
0248   gen_cfg.theta_range(theta, theta);
0249   gen_cfg.p_tot(mom);
0250   gen_cfg.charge(charge);
0251   generator_type generator(gen_cfg);
0252
0253   // Create smearer
0254   measurement_smearer<traccc::default_algebra> smearer(
0255       50.f * traccc::unit<scalar>::um, 50.f * traccc::unit<scalar>::um);
0256
0257   std::size_t n_events{1000u};
0258
0259   using detector_type = decltype(detector);
0260   using generator_type = decltype(generator);
0261   using writer_type =
0262       smearing_writer<measurement_smearer<traccc::default_algebra>>;
0263
0264   typename writer_type::config writer_cfg{smearer};
0265
0266   auto sim = simulator<detector_type, b_field_t, generator_type, writer_type>(
0267       traccc::muon<scalar>(), n_events, detector, field, std::move(generator),
0268       std::move(writer_cfg), directory);
0269
0270   // Lift step size constraints
0271   sim.get_config().propagation.stepping.step_constraint =
0272       std::numeric_limits<float>::max();
0273
0274   // Run simulation
0275   sim.get_config().propagation.navigation.intersection.overstep_tolerance =
0276       -100.f * unit<float>::um;
0277   sim.get_config().propagation.navigation.intersection.max_mask_tolerance =
0278       1.f * unit<float>::mm;
0279   sim.run();
0280
0281   for (std::size_t i_event{0u}; i_event < n_events; i_event++) {
0282     std::vector<traccc::io::csv::measurement> measurements;
0283     auto measurement_reader = traccc::io::csv::make_measurement_reader(
0284         directory +
0285         traccc::io::get_event_filename(i_event, "-measurements.csv"));
0286     traccc::io::csv::measurement io_measurement;
0287     while (measurement_reader.read(io_measurement)) {
0288       measurements.push_back(io_measurement);
0289     }
0290
0291     // Make sure that number of measurements is equal to the number of
0292     // physical planes
0293     ASSERT_EQ(measurements.size(), positions.size());
0294   }
0295 }
0296
0297 INSTANTIATE_TEST_SUITE_P(
0298     Simulation, TelescopeDetectorSimulation,
0299     ::testing::Values(
0300         std::make_tuple("0", 0.1f * traccc::unit<scalar>::GeV, 0.01f, -1.f),
0301         std::make_tuple("1", 1.f * traccc::unit<scalar>::GeV, 0.01f, -1.f),
0302         std::make_tuple("2", 10.f * traccc::unit<scalar>::GeV, 0.01f, -1.f),
0303         std::make_tuple("3", 100.f * traccc::unit<scalar>::GeV, 0.01f, -1.f),
0304         std::make_tuple("4", 0.1f * traccc::unit<scalar>::GeV,
0305                         traccc::constant<scalar>::pi / 12.f, 1.f),
0306         std::make_tuple("5", 1.f * traccc::unit<scalar>::GeV,
0307                         traccc::constant<scalar>::pi / 12.f, 1.f),
0308         std::make_tuple("6", 10.f * traccc::unit<scalar>::GeV,
0309                         traccc::constant<scalar>::pi / 12.f, 1.f),
0310         std::make_tuple("7", 100.f * traccc::unit<scalar>::GeV,
0311                         traccc::constant<scalar>::pi / 12.f, 1.f),
0312         std::make_tuple("8", 0.1f * traccc::unit<scalar>::GeV,
0313                         traccc::constant<scalar>::pi / 8.f, -1.f),
0314         std::make_tuple("9", 1.f * traccc::unit<scalar>::GeV,
0315                         traccc::constant<scalar>::pi / 8.f, -1.f),
0316         std::make_tuple("10", 10.f * traccc::unit<scalar>::GeV,
0317                         traccc::constant<scalar>::pi / 8.f, -1.f),
0318         std::make_tuple("11", 100.f * traccc::unit<scalar>::GeV,
0319                         traccc::constant<scalar>::pi / 8.f, -1.f),
0320         std::make_tuple("12", 0.1f * traccc::unit<scalar>::GeV,
0321                         traccc::constant<scalar>::pi / 6.f, 1.f),
0322         std::make_tuple("13", 1.f * traccc::unit<scalar>::GeV,
0323                         traccc::constant<scalar>::pi / 6.f, 1.f),
0324         std::make_tuple("14", 10.f * traccc::unit<scalar>::GeV,
0325                         traccc::constant<scalar>::pi / 6.f, 1.f),
0326         std::make_tuple("15", 100.f * traccc::unit<scalar>::GeV,
0327                         traccc::constant<scalar>::pi / 6.f, 1.f)));