From ca84a1fe0400e49ab4087ab3924c64a43cb6e146 Mon Sep 17 00:00:00 2001 From: Joana Niermann Date: Mon, 2 Oct 2023 20:01:02 +0200 Subject: [PATCH] Enable regular steps in eta for uniform_track_generator (this is needed for the material scans). The stepping for theta/eta has been adjusted, so that it becomes symmetric for positive and negative eta. Theta and eta ranges are inclusive, but clamped to (0, pi) and [-0.001num_max, 0.001num_max], respectively. In order to facilitate the more complex configuration, the config structs of the track generators are used in more places and the possibility to pass complicated configurations to the parametrized constructors has been limited to the most common case. --- core/include/detray/propagator/navigator.hpp | 5 +- tests/benchmarks/cpu/intersect_all.cpp | 33 +-- tests/benchmarks/cpu/intersect_surfaces.cpp | 35 ++-- .../cuda/benchmark_propagator_cuda.cpp | 5 +- .../common/test_base/propagator_test.hpp | 7 +- .../helix_cylinder_intersector.hpp | 14 +- .../tests/common/tools/particle_gun.hpp | 3 +- tests/unit_tests/cpu/tools_particle_gun.cpp | 7 +- tests/unit_tests/cpu/tools_propagator.cpp | 7 +- tests/unit_tests/cpu/tools_stepper.cpp | 10 +- .../unit_tests/cpu/tools_track_generators.cpp | 173 +++++++++++----- .../unit_tests/device/cuda/navigator_cuda.cpp | 6 +- tests/unit_tests/svgtools/intersections.cpp | 11 +- .../detail/navigation_check_helper.hpp | 2 +- .../src/toy_detector_validation.cpp | 7 +- .../src/wire_chamber_validation.cpp | 6 +- .../src/cpu/detector/detector_ray_scan.cpp | 10 +- .../cpu/propagation/navigation_inspection.cpp | 5 +- tutorials/src/device/cuda/propagation.cpp | 5 +- .../random_track_generator.hpp | 102 ++++++---- .../uniform_track_generator.hpp | 189 +++++++++++++----- 21 files changed, 396 insertions(+), 246 deletions(-) diff --git a/core/include/detray/propagator/navigator.hpp b/core/include/detray/propagator/navigator.hpp index 244c3981e..f801d4702 100644 --- a/core/include/detray/propagator/navigator.hpp +++ b/core/include/detray/propagator/navigator.hpp @@ -126,7 +126,7 @@ class navigator { sf.template visit_mask( candidates, detail::ray(track), sf_descr, det.transform_store(), - tol); + sf.is_portal() ? 0.f : tol); } }; @@ -740,7 +740,8 @@ class navigator { // Check whether this candidate is reachable by the track return sf.template visit_mask( - detail::ray(track), candidate, det->transform_store(), tol); + detail::ray(track), candidate, det->transform_store(), + sf.is_portal() ? 0.f : tol); } /// Helper to evict all unreachable/invalid candidates from the cache: diff --git a/tests/benchmarks/cpu/intersect_all.cpp b/tests/benchmarks/cpu/intersect_all.cpp index 8e2fcc37a..702e12359 100644 --- a/tests/benchmarks/cpu/intersect_all.cpp +++ b/tests/benchmarks/cpu/intersect_all.cpp @@ -29,12 +29,15 @@ // Use the detray:: namespace implicitly. using namespace detray; +using trk_generator_t = + uniform_track_generator>; + +constexpr unsigned int theta_steps{100u}; +constexpr unsigned int phi_steps{100u}; + // This test runs intersection with all surfaces of the TrackML detector void BM_INTERSECT_ALL(benchmark::State &state) { - static const unsigned int theta_steps{100u}; - static const unsigned int phi_steps{100u}; - // Detector configuration vecmem::host_memory_resource host_mr; toy_det_config toy_cfg{}; @@ -48,18 +51,22 @@ void BM_INTERSECT_ALL(benchmark::State &state) { std::size_t hits{0u}; std::size_t missed{0u}; + std::size_t n_surfaces{0u}; + test::point3 origin{0.f, 0.f, 0.f}; + std::vector> + intersections{}; + + // Iterate through uniformly distributed momentum directions + auto trk_generator = trk_generator_t{}; + trk_generator.config() + .theta_steps(theta_steps) + .phi_steps(phi_steps) + .origin(origin); for (auto _ : state) { - test::point3 pos{0.f, 0.f, 0.f}; - std::vector> - intersections{}; - std::size_t n_surfaces{0u}; - - // Iterate through uniformly distributed momentum directions - for (const auto track : - uniform_track_generator>( - theta_steps, phi_steps, pos)) { + + for (const auto track : trk_generator) { // Loop over all surfaces in detector for (const auto &sf_desc : d.surface_lookup()) { diff --git a/tests/benchmarks/cpu/intersect_surfaces.cpp b/tests/benchmarks/cpu/intersect_surfaces.cpp index 0617b67fd..5637821ab 100644 --- a/tests/benchmarks/cpu/intersect_surfaces.cpp +++ b/tests/benchmarks/cpu/intersect_surfaces.cpp @@ -26,6 +26,8 @@ // Use the detray:: namespace implicitly. using namespace detray; +using ray_generator_t = uniform_track_generator>; + static const unsigned int theta_steps = 1000u; static const unsigned int phi_steps = 1000u; @@ -41,16 +43,17 @@ void BM_INTERSECT_PLANES(benchmark::State &state) { auto planes = test::planes_along_direction( dists, vector::normalize(test::vector3{1.f, 1.f, 1.f})); constexpr mask> rect{0u, 10.f, 20.f}; - test::point3 ori = {0.f, 0.f, 0.f}; + + // Iterate through uniformly distributed momentum directions + auto ray_generator = ray_generator_t{}; + ray_generator.config().theta_steps(theta_steps).phi_steps(phi_steps); for (auto _ : state) { benchmark::DoNotOptimize(sfhit); benchmark::DoNotOptimize(sfmiss); // Iterate through uniformly distributed momentum directions - for (const auto ray : - uniform_track_generator>( - theta_steps, phi_steps, ori, 1.f)) { + for (const auto ray : ray_generator) { for (const auto &plane : planes) { auto pi = rect.intersector>( - theta_steps, phi_steps, ori, 1.f)) { + for (const auto ray : ray_generator) { for (const auto &cylinder : cylinders) { auto ci = cylinder.intersector(); @@ -170,16 +173,16 @@ void BM_INTERSECT_PORTAL_CYLINDERS(benchmark::State &state) { plane_surface plane(test::transform3(), mask_link, material_link, 0u, false, surface_id::e_sensitive); - const test::point3 ori = {0.f, 0.f, 0.f}; + // Iterate through uniformly distributed momentum directions + auto ray_generator = ray_generator_t{}; + ray_generator.config().theta_steps(theta_steps).phi_steps(phi_steps); for (auto _ : state) { benchmark::DoNotOptimize(sfhit); benchmark::DoNotOptimize(sfmiss); // Iterate through uniformly distributed momentum directions - for (const auto ray : - uniform_track_generator>( - theta_steps, phi_steps, ori, 1.f)) { + for (const auto ray : ray_generator) { for (const auto &cylinder : cylinders) { auto cpi = cylinder.intersector(); @@ -222,14 +225,14 @@ void BM_INTERSECT_CONCETRIC_CYLINDERS(benchmark::State &state) { plane_surface plane(test::transform3(), mask_link, material_link, 0u, false, surface_id::e_sensitive); - const test::point3 ori = {0.f, 0.f, 0.f}; + // Iterate through uniformly distributed momentum directions + auto ray_generator = ray_generator_t{}; + ray_generator.config().theta_steps(theta_steps).phi_steps(phi_steps); for (auto _ : state) { // Iterate through uniformly distributed momentum directions - for (const auto ray : - uniform_track_generator>( - theta_steps, phi_steps, ori, 1.f)) { + for (const auto ray : ray_generator) { for (const auto &cylinder : cylinders) { auto cci = cylinder.intersector(); diff --git a/tests/benchmarks/cuda/benchmark_propagator_cuda.cpp b/tests/benchmarks/cuda/benchmark_propagator_cuda.cpp index 522a22119..d962310cd 100644 --- a/tests/benchmarks/cuda/benchmark_propagator_cuda.cpp +++ b/tests/benchmarks/cuda/benchmark_propagator_cuda.cpp @@ -31,13 +31,12 @@ toy_det_config toy_cfg{4u, 7u}; void fill_tracks(vecmem::vector> &tracks, const std::size_t theta_steps, const std::size_t phi_steps) { - // Set origin position of tracks - const point3 ori{0.f, 0.f, 0.f}; + // Set momentum of tracks const scalar mom_mag{10.f * unit::GeV}; // Iterate through uniformly distributed momentum directions for (auto traj : uniform_track_generator>( - theta_steps, phi_steps, ori, mom_mag)) { + phi_steps, theta_steps, mom_mag)) { tracks.push_back(traj); } } diff --git a/tests/common/include/tests/common/test_base/propagator_test.hpp b/tests/common/include/tests/common/test_base/propagator_test.hpp index 70be60af0..ec7096a53 100644 --- a/tests/common/include/tests/common/test_base/propagator_test.hpp +++ b/tests/common/include/tests/common/test_base/propagator_test.hpp @@ -139,14 +139,11 @@ inline vecmem::vector generate_tracks( // Track collection vecmem::vector tracks(mr); - // Set origin position of tracks - const point3 ori{0.f, 0.f, 0.f}; + // Set momentum of tracks const scalar p_mag{10.f * unit::GeV}; // Iterate through uniformly distributed momentum directions - for (auto track : uniform_track_generator( - ts, ps, ori, p_mag, {0.01f, constant::pi}, - {-constant::pi, constant::pi})) { + for (auto track : uniform_track_generator(ps, ts, p_mag)) { track.set_overstep_tolerance(overstep_tolerance); // Put it into vector of trajectories diff --git a/tests/common/include/tests/common/tools/intersectors/helix_cylinder_intersector.hpp b/tests/common/include/tests/common/tools/intersectors/helix_cylinder_intersector.hpp index c8245d8f6..43ab34f40 100644 --- a/tests/common/include/tests/common/tools/intersectors/helix_cylinder_intersector.hpp +++ b/tests/common/include/tests/common/tools/intersectors/helix_cylinder_intersector.hpp @@ -66,7 +66,7 @@ struct helix_cylinder_intersector // Guard against inifinite loops constexpr std::size_t max_n_tries{1000u}; // Tolerance for convergence - constexpr scalar_type tol{1e-4f}; + constexpr scalar_type tol{1e-3f}; // Get the surface placement const auto &sm = trf.matrix(); @@ -128,10 +128,12 @@ struct helix_cylinder_intersector const vector3 crp = vector::cross(h.pos(s) - sc, sz); const scalar_type denom{ 2.f * vector::dot(crp, vector::cross(h.dir(s), sz))}; + // No intersection can be found if dividing by zero if (denom == 0.f) { return ret; } + // x_n+1 = x_n - f(s) / f'(s) s_prev = s; s -= (vector::dot(crp, crp) - r * r) / denom; @@ -149,16 +151,6 @@ struct helix_cylinder_intersector is.local = mask.to_local_frame(trf, p3); is.status = mask.is_inside(is.local, mask_tolerance); - // Perform the r-check for Newton solution even if it is not - // required by the mask's shape - const bool r_check = - std::abs(r - is.local[2]) < - mask_tolerance + - 5.f * std::numeric_limits::epsilon(); - if (not r_check) { - is.status = intersection::status::e_outside; - } - // Compute some additional information if the intersection is valid if (is.status == intersection::status::e_inside) { is.sf_desc = sf; diff --git a/tests/common/include/tests/common/tools/particle_gun.hpp b/tests/common/include/tests/common/tools/particle_gun.hpp index b974b6d68..95ce87ba0 100644 --- a/tests/common/include/tests/common/tools/particle_gun.hpp +++ b/tests/common/include/tests/common/tools/particle_gun.hpp @@ -63,7 +63,8 @@ struct particle_gun { // Retrieve candidate(s) from the surface const auto sf = surface{detector, sf_desc}; sf.template visit_mask( - intersections, traj, sf_desc, tf_store, mask_tolerance); + intersections, traj, sf_desc, tf_store, + sf.is_portal() ? 0.f : mask_tolerance); // Candidate is invalid if it lies in the opposite direction for (auto &sfi : intersections) { diff --git a/tests/unit_tests/cpu/tools_particle_gun.cpp b/tests/unit_tests/cpu/tools_particle_gun.cpp index d2b014126..b8c30be22 100644 --- a/tests/unit_tests/cpu/tools_particle_gun.cpp +++ b/tests/unit_tests/cpu/tools_particle_gun.cpp @@ -43,7 +43,6 @@ GTEST_TEST(detray_tools, particle_gun) { unsigned int theta_steps{50u}; unsigned int phi_steps{50u}; - const point3 ori{0.f, 0.f, 0.f}; // Record ray tracing using detector_t = decltype(toy_det); @@ -52,8 +51,8 @@ GTEST_TEST(detray_tools, particle_gun) { std::vector>> expected; // Iterate through uniformly distributed momentum directions with ray for (const auto test_ray : - uniform_track_generator>( - theta_steps, phi_steps, ori)) { + uniform_track_generator>(phi_steps, + theta_steps)) { // Record all intersections and objects along the ray const auto intersection_record = @@ -66,7 +65,7 @@ GTEST_TEST(detray_tools, particle_gun) { std::size_t n_tracks{0u}; for (const auto track : uniform_track_generator>( - theta_steps, phi_steps, ori)) { + phi_steps, theta_steps)) { const detail::helix test_helix(track, &B); // Record all intersections and objects along the ray diff --git a/tests/unit_tests/cpu/tools_propagator.cpp b/tests/unit_tests/cpu/tools_propagator.cpp index ff3ed2038..3b0589007 100644 --- a/tests/unit_tests/cpu/tools_propagator.cpp +++ b/tests/unit_tests/cpu/tools_propagator.cpp @@ -160,7 +160,6 @@ class PropagatorWithRkStepper /// Track generator configuration unsigned int theta_steps{50u}; unsigned int phi_steps{50u}; - point3 ori{0.f, 0.f, 0.f}; scalar mom{10.f * unit::GeV}; /// Stepper configuration @@ -200,7 +199,7 @@ TEST_P(PropagatorWithRkStepper, rk4_propagator_const_bfield) { // Iterate through uniformly distributed momentum directions for (auto track : - uniform_track_generator(theta_steps, phi_steps, ori, mom)) { + uniform_track_generator(phi_steps, theta_steps, mom)) { // Generate second track state used for propagation with pathlimit track_t lim_track(track); @@ -241,7 +240,7 @@ TEST_P(PropagatorWithRkStepper, rk4_propagator_const_bfield) { // Propagate the entire detector ASSERT_TRUE(p.propagate(state, actor_states)) << print_insp_state.to_string() << std::endl; - // << state._navigation.inspector().to_string() << std::endl; + // << state._navigation.inspector().to_string() << std::endl; // Propagate with path limit ASSERT_NEAR(pathlimit_aborter_state.path_limit(), path_limit, tol); @@ -299,7 +298,7 @@ TEST_P(PropagatorWithRkStepper, rk4_propagator_inhom_bfield) { // Iterate through uniformly distributed momentum directions for (auto track : - uniform_track_generator(theta_steps, phi_steps, ori, mom)) { + uniform_track_generator(phi_steps, theta_steps, mom)) { // Genrate second track state used for propagation with pathlimit track_t lim_track(track); diff --git a/tests/unit_tests/cpu/tools_stepper.cpp b/tests/unit_tests/cpu/tools_stepper.cpp index 55b9acfbb..b05221c60 100644 --- a/tests/unit_tests/cpu/tools_stepper.cpp +++ b/tests/unit_tests/cpu/tools_stepper.cpp @@ -166,7 +166,6 @@ GTEST_TEST(detray_propagator, rk_stepper) { constexpr scalar stepsize_constr{0.5f * unit::mm}; // Track generator configuration - const point3 ori{0.f, 0.f, 0.f}; const scalar p_mag{10.f * unit::GeV}; constexpr unsigned int theta_steps = 100u; constexpr unsigned int phi_steps = 100u; @@ -174,7 +173,7 @@ GTEST_TEST(detray_propagator, rk_stepper) { // Iterate through uniformly distributed momentum directions for (auto track : uniform_track_generator>( - theta_steps, phi_steps, ori, p_mag)) { + phi_steps, theta_steps, p_mag)) { // Generate track state used for propagation with constrained step size free_track_parameters c_track(track); @@ -241,7 +240,7 @@ GTEST_TEST(detray_propagator, rk_stepper) { ASSERT_NEAR(crk_state.path_length(), 0.f, tol); const point3 backward_relative_error{1.f / (2.f * path_length) * - (rk_state().pos() - ori)}; + (rk_state().pos())}; // Make sure that relative error is smaller than the tolerance EXPECT_NEAR(getter::norm(backward_relative_error), 0.f, tol); @@ -271,7 +270,6 @@ TEST(detray_propagator, rk_stepper_inhomogeneous_bfield) { constexpr scalar stepsize_constr{0.5f * unit::mm}; // Track generator configuration - const point3 ori{0.f, 0.f, 0.f}; const scalar p_mag{10.f * unit::GeV}; constexpr unsigned int theta_steps = 100u; constexpr unsigned int phi_steps = 100u; @@ -279,7 +277,7 @@ TEST(detray_propagator, rk_stepper_inhomogeneous_bfield) { // Iterate through uniformly distributed momentum directions for (auto track : uniform_track_generator>( - theta_steps, phi_steps, ori, p_mag)) { + phi_steps, theta_steps, p_mag)) { // Generate track state used for propagation with constrained step size free_track_parameters c_track(track); @@ -328,7 +326,7 @@ TEST(detray_propagator, rk_stepper_inhomogeneous_bfield) { ASSERT_NEAR(crk_state.path_length(), 0.f, tol); const point3 backward_relative_error{1.f / (2.f * path_length) * - (rk_state().pos() - ori)}; + (rk_state().pos())}; // Make sure that relative error is smaller than the tolerance EXPECT_NEAR(getter::norm(backward_relative_error), 0.f, tol); diff --git a/tests/unit_tests/cpu/tools_track_generators.cpp b/tests/unit_tests/cpu/tools_track_generators.cpp index 1517181ab..c89780800 100644 --- a/tests/unit_tests/cpu/tools_track_generators.cpp +++ b/tests/unit_tests/cpu/tools_track_generators.cpp @@ -1,6 +1,6 @@ /** Detray library, part of the ACTS project (R&D line) * - * (c) 2022 CERN for the benefit of the ACTS project + * (c) 2022-2023 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -20,21 +20,25 @@ using point3 = test::point3; using transform3 = test::transform3; GTEST_TEST(detray_simulation, uniform_track_generator) { + using generator_t = + uniform_track_generator>; constexpr const scalar tol{1e-5f}; + constexpr const scalar epsilon{generator_t::configuration::epsilon}; constexpr std::size_t phi_steps{50u}; constexpr std::size_t theta_steps{50u}; std::array momenta{}; - // Loops of theta values ]0,pi[ + // Loop over theta values ]0,pi[ for (std::size_t itheta{0u}; itheta < theta_steps; ++itheta) { - const scalar theta{0.01f + static_cast(itheta) * - (constant::pi - 0.01f) / - static_cast(theta_steps)}; + const scalar theta{epsilon + + static_cast(itheta) * + (constant::pi - 2.f * epsilon) / + static_cast(theta_steps - 1u)}; - // Loops of phi values [-pi, pi] + // Loop over phi values [-pi, pi] for (std::size_t iphi{0u}; iphi < phi_steps; ++iphi) { // The direction const scalar phi{-constant::pi + @@ -55,10 +59,8 @@ GTEST_TEST(detray_simulation, uniform_track_generator) { // Now run the track generator and compare std::size_t n_tracks{0u}; - for (const auto track : - uniform_track_generator>( - theta_steps, phi_steps)) { - vector3 &expected = momenta[n_tracks]; + for (const auto track : generator_t(phi_steps, theta_steps)) { + vector3& expected = momenta[n_tracks]; vector3 result = track.mom(); // Compare with for loop @@ -74,9 +76,8 @@ GTEST_TEST(detray_simulation, uniform_track_generator) { // Generate rays n_tracks = 0u; - for (const auto r : uniform_track_generator>( - theta_steps, phi_steps)) { - vector3 &expected = momenta[n_tracks]; + for (const auto r : generator_t(phi_steps, theta_steps)) { + vector3& expected = momenta[n_tracks]; vector3 result = r.dir(); // Compare with for loop @@ -94,11 +95,9 @@ GTEST_TEST(detray_simulation, uniform_track_generator) { const vector3 B{0.f * unit::T, 0.f * unit::T, 2.f * unit::T}; n_tracks = 0u; - for (const auto track : - uniform_track_generator>( - theta_steps, phi_steps)) { + for (const auto track : generator_t(phi_steps, theta_steps)) { detail::helix helix_traj(track, &B); - vector3 &expected = momenta[n_tracks]; + vector3& expected = momenta[n_tracks]; vector3 result = helix_traj.dir(0.f); // Compare with for loop @@ -113,19 +112,77 @@ GTEST_TEST(detray_simulation, uniform_track_generator) { ASSERT_EQ(momenta.size(), n_tracks); } -GTEST_TEST(detray_simulation, uniform_track_generator_with_range) { +GTEST_TEST(detray_simulation, uniform_track_generator_eta) { + using generator_t = + uniform_track_generator>; constexpr const scalar tol{1e-5f}; - constexpr std::size_t theta_steps{2u}; - constexpr std::size_t phi_steps{4u}; + constexpr std::size_t phi_steps{50u}; + constexpr std::size_t eta_steps{50u}; + + std::array momenta{}; + + // Loop over eta values [-5, 5] + for (std::size_t ieta{0u}; ieta < eta_steps; ++ieta) { + const scalar eta{-5.f + static_cast(ieta) * (10.f) / + static_cast(eta_steps - 1u)}; + const scalar theta{2.f * std::atan(std::exp(-eta))}; + + // Loop over phi values [-pi, pi] + for (std::size_t iphi{0u}; iphi < phi_steps; ++iphi) { + // The direction + const scalar phi{-constant::pi + + static_cast(iphi) * + (2.f * constant::pi) / + static_cast(phi_steps)}; + + // intialize a track + vector3 mom{std::cos(phi) * std::sin(theta), + std::sin(phi) * std::sin(theta), std::cos(theta)}; + vector::normalize(mom); + free_track_parameters traj({0.f, 0.f, 0.f}, 0.f, mom, + -1.f); + + momenta[ieta * phi_steps + iphi] = traj.mom(); + } + } + + // Now run the track generator and compare + std::size_t n_tracks{0u}; + + auto trk_generator = generator_t{}; + trk_generator.config().phi_steps(phi_steps).eta_steps(eta_steps); + + for (const auto track : trk_generator) { + vector3& expected = momenta[n_tracks]; + vector3 result = track.mom(); + + // Compare with for loop + EXPECT_NEAR(getter::norm(expected - result), 0.f, tol) + << "Expected\tResult: \n" + << expected[0] << "\t" << result[0] << "\n" + << expected[1] << "\t" << result[1] << "\n" + << expected[2] << "\t" << result[2] << std::endl; + + ++n_tracks; + } + ASSERT_EQ(momenta.size(), n_tracks); +} + +GTEST_TEST(detray_simulation, uniform_track_generator_with_range) { + using generator_t = + uniform_track_generator>; + + constexpr const scalar tol{1e-5f}; std::vector> theta_phi; - for (const auto track : - uniform_track_generator>( - theta_steps, phi_steps, {0.f, 0.f, 0.f}, 1.f * unit::GeV, - {1.f, 2.f}, {-2.f, 2.f})) { + auto trk_gen_cfg = generator_t::configuration{}; + trk_gen_cfg.phi_range(-2.f, 2.f).phi_steps(4u); + trk_gen_cfg.theta_range(1.f, 2.f).theta_steps(2u); + + for (const auto track : generator_t{trk_gen_cfg}) { const auto dir = track.dir(); theta_phi.push_back({getter::theta(dir), getter::phi(dir)}); } @@ -139,13 +196,13 @@ GTEST_TEST(detray_simulation, uniform_track_generator_with_range) { EXPECT_NEAR(theta_phi[2][1], 0.f, tol); EXPECT_NEAR(theta_phi[3][0], 1.f, tol); EXPECT_NEAR(theta_phi[3][1], 1.f, tol); - EXPECT_NEAR(theta_phi[4][0], 1.5f, tol); + EXPECT_NEAR(theta_phi[4][0], 2.f, tol); EXPECT_NEAR(theta_phi[4][1], -2.f, tol); - EXPECT_NEAR(theta_phi[5][0], 1.5f, tol); + EXPECT_NEAR(theta_phi[5][0], 2.f, tol); EXPECT_NEAR(theta_phi[5][1], -1.f, tol); - EXPECT_NEAR(theta_phi[6][0], 1.5f, tol); + EXPECT_NEAR(theta_phi[6][0], 2.f, tol); EXPECT_NEAR(theta_phi[6][1], 0.f, tol); - EXPECT_NEAR(theta_phi[7][0], 1.5f, tol); + EXPECT_NEAR(theta_phi[7][0], 2.f, tol); EXPECT_NEAR(theta_phi[7][1], 1.f, tol); } @@ -156,6 +213,9 @@ GTEST_TEST(detray_simulation, random_track_generator_uniform) { using uniform_gen_t = random_numbers, std::seed_seq>; + using trk_generator_t = + random_track_generator, + uniform_gen_t>; // Tolerance depends on sample size constexpr scalar tol{0.05f}; @@ -165,14 +225,13 @@ GTEST_TEST(detray_simulation, random_track_generator_uniform) { constexpr std::size_t n_gen_tracks{10000u}; // Other params - point3 ori = {0.f, 0.f, 0.f}; - point3 ori_stddev = {0.1f * unit::mm, 0.f * unit::mm, - 0.2f * unit::mm}; - std::array mom_range = {1.f * unit::GeV, - 2.f * unit::GeV}; - std::array theta_range = {0.01f, constant::pi}; - std::array phi_range = {-0.9f * constant::pi, - 0.8f * constant::pi}; + trk_generator_t::configuration trk_gen_cfg{}; + trk_gen_cfg.n_tracks(n_gen_tracks); + trk_gen_cfg.phi_range(-0.9f * constant::pi, + 0.8f * constant::pi); + trk_gen_cfg.mom_range(1.f * unit::GeV, 2.f * unit::GeV); + trk_gen_cfg.origin_stddev({0.1f * unit::mm, 0.f * unit::mm, + 0.2f * unit::mm}); // Catch the results std::array x{}; @@ -182,11 +241,7 @@ GTEST_TEST(detray_simulation, random_track_generator_uniform) { std::array phi{}; std::array theta{}; - for (const auto track : - random_track_generator, - uniform_gen_t>(n_gen_tracks, ori, ori_stddev, - mom_range, theta_range, - phi_range)) { + for (const auto track : trk_generator_t{trk_gen_cfg}) { const auto pos = track.pos(); x[n_tracks] = pos[0]; y[n_tracks] = pos[1]; @@ -201,6 +256,11 @@ GTEST_TEST(detray_simulation, random_track_generator_uniform) { ASSERT_EQ(n_gen_tracks, n_tracks); // Check uniform distrubution + const auto& ori = trk_gen_cfg.origin(); + const auto& ori_stddev = trk_gen_cfg.origin_stddev(); + const auto& phi_range = trk_gen_cfg.phi_range(); + const auto& theta_range = trk_gen_cfg.theta_range(); + const auto& mom_range = trk_gen_cfg.mom_range(); // Mean EXPECT_NEAR(statistics::mean(x), ori[0], tol); @@ -232,6 +292,8 @@ GTEST_TEST(detray_simulation, random_track_generator_normal) { // Use deterministic random number generator for testing using normal_gen_t = random_numbers, std::seed_seq>; + using trk_generator_t = + random_track_generator, normal_gen_t>; // Tolerance depends on sample size constexpr scalar tol{0.05f}; @@ -241,14 +303,14 @@ GTEST_TEST(detray_simulation, random_track_generator_normal) { constexpr std::size_t n_gen_tracks{10000u}; // Other params - point3 ori = {0.f, 0.f, 0.f}; - point3 ori_stddev = {0.1f * unit::mm, 0.5f * unit::mm, - 0.3f * unit::mm}; - std::array mom_range = {1.f * unit::GeV, - 2.f * unit::GeV}; - std::array theta_range = {0.01f, constant::pi}; - std::array phi_range = {-0.9f * constant::pi, - 0.8f * constant::pi}; + trk_generator_t::configuration trk_gen_cfg{}; + trk_gen_cfg.n_tracks(n_gen_tracks); + trk_gen_cfg.phi_range(-0.9f * constant::pi, + 0.8f * constant::pi); + trk_gen_cfg.mom_range(1.f * unit::GeV, 2.f * unit::GeV); + trk_gen_cfg.origin({0.f, 0.f, 0.f}); + trk_gen_cfg.origin_stddev({0.1f * unit::mm, 0.5f * unit::mm, + 0.3f * unit::mm}); // Catch the results std::array x{}; @@ -258,11 +320,7 @@ GTEST_TEST(detray_simulation, random_track_generator_normal) { std::array phi{}; std::array theta{}; - for (const auto track : - random_track_generator, - normal_gen_t>(n_gen_tracks, ori, ori_stddev, - mom_range, theta_range, - phi_range)) { + for (const auto track : trk_generator_t{trk_gen_cfg}) { const auto pos = track.pos(); x[n_tracks] = pos[0]; y[n_tracks] = pos[1]; @@ -277,6 +335,11 @@ GTEST_TEST(detray_simulation, random_track_generator_normal) { ASSERT_EQ(n_gen_tracks, n_tracks); // check gaussian distribution - values are clamped to phi/theta range + const auto& ori = trk_gen_cfg.origin(); + const auto& ori_stddev = trk_gen_cfg.origin_stddev(); + const auto& phi_range = trk_gen_cfg.phi_range(); + const auto& theta_range = trk_gen_cfg.theta_range(); + const auto& mom_range = trk_gen_cfg.mom_range(); // Mean EXPECT_NEAR(statistics::mean(x), ori[0], tol); @@ -302,4 +365,4 @@ GTEST_TEST(detray_simulation, random_track_generator_normal) { EXPECT_NEAR(statistics::variance(theta), std::pow(0.5f / 3.0f * (theta_range[1] - theta_range[0]), 2.f), tol); -} \ No newline at end of file +} diff --git a/tests/unit_tests/device/cuda/navigator_cuda.cpp b/tests/unit_tests/device/cuda/navigator_cuda.cpp index 6597e0038..a622042b3 100644 --- a/tests/unit_tests/device/cuda/navigator_cuda.cpp +++ b/tests/unit_tests/device/cuda/navigator_cuda.cpp @@ -32,15 +32,13 @@ TEST(navigator_cuda, navigator) { vecmem::vector> tracks_host(&mng_mr); vecmem::vector> tracks_device(&mng_mr); - // Set origin position of tracks - const point3 ori{0.f, 0.f, 0.f}; + // Magnitude of total momentum of tracks const scalar p_mag{10.f * unit::GeV}; // Iterate through uniformly distributed momentum directions for (auto track : uniform_track_generator>( - theta_steps, phi_steps, ori, p_mag, {0.01f, constant::pi}, - {-constant::pi, constant::pi})) { + phi_steps, theta_steps, p_mag)) { track.set_overstep_tolerance(overstep_tolerance); tracks_host.push_back(track); diff --git a/tests/unit_tests/svgtools/intersections.cpp b/tests/unit_tests/svgtools/intersections.cpp index 3d086ac7c..00efe6c5a 100644 --- a/tests/unit_tests/svgtools/intersections.cpp +++ b/tests/unit_tests/svgtools/intersections.cpp @@ -48,15 +48,14 @@ int main(int, char**) { const auto svg_det = il.draw_detector("detector", context, view); // Creating the rays. - unsigned int theta_steps{10u}; - unsigned int phi_steps{10u}; - const typename detector_t::point3 ori{0.f, 0.f, 100.f}; + using generator_t = + detray::uniform_track_generator>; + auto trk_gen_cfg = generator_t::configuration{}; + trk_gen_cfg.origin({0.f, 0.f, 100.f}).phi_steps(10u).theta_steps(10u); std::size_t index = 0; // Iterate through uniformly distributed momentum directions with ray - for (const auto test_ray : - detray::uniform_track_generator>( - theta_steps, phi_steps, ori)) { + for (const auto test_ray : generator_t{trk_gen_cfg}) { // Record all intersections and objects along the ray const auto intersection_record = diff --git a/tests/validation/include/detray/validation/detail/navigation_check_helper.hpp b/tests/validation/include/detray/validation/detail/navigation_check_helper.hpp index 96bff3b21..127ec05c2 100644 --- a/tests/validation/include/detray/validation/detail/navigation_check_helper.hpp +++ b/tests/validation/include/detray/validation/detail/navigation_check_helper.hpp @@ -84,7 +84,7 @@ bool compare_traces(const inters_trace_t &intersection_trace, // Do a final check on the trace sizes const bool is_size{n_inters_nav == intersection_trace.size()}; EXPECT_TRUE(is_size) << "ERROR: Intersection traces found different number " - "of surfaces! Please check the last elements" + "of surfaces! Please check the last elements\n" << debug_stream.str(); if (not is_size) { return false; diff --git a/tests/validation/src/toy_detector_validation.cpp b/tests/validation/src/toy_detector_validation.cpp index 47b1f53b0..31c81a822 100644 --- a/tests/validation/src/toy_detector_validation.cpp +++ b/tests/validation/src/toy_detector_validation.cpp @@ -61,10 +61,7 @@ int main(int argc, char **argv) { cfg_hel_scan.name("toy_detector_helix_scan"); cfg_hel_scan.overstepping_tolerance(-100.f * unit::um); cfg_hel_scan.track_generator().p_mag(10.f * unit::GeV); - // Fails in single precision for more helices (unless is configured with - // only three edc layers) - // cfg_hel_scan.track_generator().theta_steps(100u).phi_steps(100u); - cfg_hel_scan.track_generator().theta_steps(30u).phi_steps(30u); + cfg_hel_scan.track_generator().theta_steps(100u).phi_steps(100u); detray::detail::register_checks(toy_det, toy_names, cfg_hel_scan); @@ -84,7 +81,7 @@ int main(int argc, char **argv) { // TODO: Fails due to mask tolerances for more helices, regardless of edc // configuration/precision // cfg_hel_nav.track_generator().theta_steps(100u).phi_steps(100u); - cfg_hel_nav.track_generator().theta_steps(30u).phi_steps(30u); + cfg_hel_nav.track_generator().theta_steps(50u).phi_steps(50u); detail::register_checks(toy_det, toy_names, cfg_hel_nav); diff --git a/tests/validation/src/wire_chamber_validation.cpp b/tests/validation/src/wire_chamber_validation.cpp index 96803d378..b7e8a20d6 100644 --- a/tests/validation/src/wire_chamber_validation.cpp +++ b/tests/validation/src/wire_chamber_validation.cpp @@ -60,9 +60,7 @@ int main(int argc, char **argv) { cfg_hel_scan.name("wire_chamber_helix_scan"); cfg_hel_scan.overstepping_tolerance(-100.f * unit::um); cfg_hel_scan.track_generator().p_mag(10.f * unit::GeV); - // Fails in single precision for more helices - // cfg_hel_scan.track_generator().theta_steps(100u).phi_steps(100u); - cfg_hel_scan.track_generator().theta_steps(3u).phi_steps(10u); + cfg_hel_scan.track_generator().theta_steps(100u).phi_steps(100u); detray::detail::register_checks(det, names, cfg_hel_scan); @@ -80,7 +78,7 @@ int main(int argc, char **argv) { cfg_hel_nav.track_generator() = cfg_hel_scan.track_generator(); // TODO: Fails for more helices // cfg_hel_nav.track_generator().theta_steps(100u).phi_steps(100u); - cfg_hel_nav.track_generator().theta_steps(3u).phi_steps(10u); + cfg_hel_nav.track_generator().theta_steps(25u).phi_steps(25u); detail::register_checks(det, names, cfg_hel_nav); diff --git a/tutorials/src/cpu/detector/detector_ray_scan.cpp b/tutorials/src/cpu/detector/detector_ray_scan.cpp index ffc9c665d..2249a18c1 100644 --- a/tutorials/src/cpu/detector/detector_ray_scan.cpp +++ b/tutorials/src/cpu/detector/detector_ray_scan.cpp @@ -56,13 +56,9 @@ int main() { // Index of the volume that the ray origin lies in detray::dindex start_index{0u}; - // Number of rays in theta and phi - unsigned int theta_steps{100u}; - unsigned int phi_steps{100u}; - // Origin of the rays - const detray::tutorial::point3 origin{0.f, 0.f, 0.f}; - auto ray_generator = - detray::uniform_track_generator(theta_steps, phi_steps, origin); + // Generate a number of rays + auto ray_generator = detray::uniform_track_generator{}; + ray_generator.config().theta_steps(100u).phi_steps(100u); // Run the check std::cout << "\nScanning detector (" << ray_generator.size() diff --git a/tutorials/src/cpu/propagation/navigation_inspection.cpp b/tutorials/src/cpu/propagation/navigation_inspection.cpp index 0420fe5eb..344dfde38 100644 --- a/tutorials/src/cpu/propagation/navigation_inspection.cpp +++ b/tutorials/src/cpu/propagation/navigation_inspection.cpp @@ -75,11 +75,10 @@ int main() { using ray_type = detray::detail::ray; constexpr std::size_t theta_steps{1}; constexpr std::size_t phi_steps{1}; - const detray::tutorial::point3 origin{0.f, 0.f, 0.f}; // Iterate through uniformly distributed momentum directions - for (const auto ray : detray::uniform_track_generator( - theta_steps, phi_steps, origin)) { + for (const auto ray : + detray::uniform_track_generator(phi_steps, theta_steps)) { // Shoot ray through the detector and record all surface intersections const auto intersection_trace = diff --git a/tutorials/src/device/cuda/propagation.cpp b/tutorials/src/device/cuda/propagation.cpp index b25f1ca94..0135e21e3 100644 --- a/tutorials/src/device/cuda/propagation.cpp +++ b/tutorials/src/device/cuda/propagation.cpp @@ -42,8 +42,7 @@ int main() { // Track directions to be generated constexpr unsigned int theta_steps{10u}; constexpr unsigned int phi_steps{10u}; - // Set origin and direction of tracks - const detray::tutorial::point3 origin{0.f, 0.f, 0.f}; + // Set momentum of tracks const detray::scalar p_mag{10.f * detray::unit::GeV}; // How much can the navigator overshoot on a given surface? constexpr detray::scalar overstep_tolerance{ @@ -52,7 +51,7 @@ int main() { // Genrate the tracks for (auto track : detray::uniform_track_generator< detray::free_track_parameters>( - theta_steps, phi_steps, origin, p_mag)) { + phi_steps, theta_steps, p_mag)) { // Set the oversetpping tolerance for every track track.set_overstep_tolerance(overstep_tolerance); // Put it into vector of tracks diff --git a/utils/include/detray/simulation/event_generator/random_track_generator.hpp b/utils/include/detray/simulation/event_generator/random_track_generator.hpp index e85d005b2..341fbdfee 100644 --- a/utils/include/detray/simulation/event_generator/random_track_generator.hpp +++ b/utils/include/detray/simulation/event_generator/random_track_generator.hpp @@ -82,6 +82,8 @@ class random_track_generator /// Configure how tracks are generated struct configuration { + /// Ensure sensible values at the theta bounds, even in single precision + static constexpr scalar epsilon{1e-2f}; /// Gaussian vertex smearing bool m_do_vtx_smearing = true; @@ -89,10 +91,12 @@ class random_track_generator /// How many tracks will be generated std::size_t m_n_tracks{10u}; - /// Range for theta and phi + /// Range for phi and theta std::array m_phi_range{-constant::pi, constant::pi}; - std::array m_theta_range{0.01f, constant::pi}; + std::array m_theta_range{epsilon, + constant::pi - epsilon}; + /// Momentum range std::array m_mom_range{1.f * unit::GeV, 1.f * unit::GeV}; @@ -104,35 +108,49 @@ class random_track_generator /// Setters /// @{ - configuration& do_vertex_smearing(bool b) { + DETRAY_HOST_DEVICE configuration& do_vertex_smearing(bool b) { m_do_vtx_smearing = b; return *this; } - configuration& n_tracks(std::size_t n) { + DETRAY_HOST_DEVICE configuration& n_tracks(std::size_t n) { + assert(n > 0); m_n_tracks = n; return *this; } - configuration& theta_range(scalar low, scalar high) { - m_theta_range = {low, high}; + DETRAY_HOST_DEVICE configuration& phi_range(scalar low, scalar high) { + assert(low <= high); + m_phi_range = {low, high}; return *this; } - configuration& phi_range(scalar low, scalar high) { - m_phi_range = {low, high}; + DETRAY_HOST_DEVICE configuration& theta_range(scalar low, scalar high) { + auto min_theta{ + std::clamp(low, epsilon, constant::pi - epsilon)}; + auto max_theta{ + std::clamp(high, epsilon, constant::pi - epsilon)}; + + assert(min_theta <= max_theta); + + m_theta_range = {min_theta, max_theta}; return *this; } - configuration& origin(point3 ori) { + DETRAY_HOST_DEVICE configuration& mom_range(scalar low, scalar high) { + assert(low <= high); + m_mom_range = {low, high}; + return *this; + } + DETRAY_HOST_DEVICE configuration& origin(point3 ori) { m_origin = ori; return *this; } - configuration& origin_stddev(point3 stddev) { + DETRAY_HOST_DEVICE configuration& origin_stddev(point3 stddev) { m_origin_stddev = stddev; return *this; } - configuration& time(scalar t) { + DETRAY_HOST_DEVICE configuration& time(scalar t) { m_time = t; return *this; } - configuration& charge(scalar q) { + DETRAY_HOST_DEVICE configuration& charge(scalar q) { m_charge = q; return *this; } @@ -140,23 +158,32 @@ class random_track_generator /// Getters /// @{ - constexpr bool do_vertex_smearing() const { return m_do_vtx_smearing; } - constexpr std::size_t n_tracks() const { return m_n_tracks; } - constexpr const std::array& theta_range() const { - return m_theta_range; + DETRAY_HOST_DEVICE constexpr bool do_vertex_smearing() const { + return m_do_vtx_smearing; + } + DETRAY_HOST_DEVICE constexpr std::size_t n_tracks() const { + return m_n_tracks; } - constexpr const std::array& phi_range() const { + DETRAY_HOST_DEVICE constexpr const std::array& phi_range() + const { return m_phi_range; } - constexpr const std::array& mom_range() const { + DETRAY_HOST_DEVICE constexpr const std::array& theta_range() + const { + return m_theta_range; + } + DETRAY_HOST_DEVICE constexpr const std::array& mom_range() + const { return m_mom_range; } - constexpr const point3& origin() const { return m_origin; } - constexpr const point3& origin_stddev() const { + DETRAY_HOST_DEVICE constexpr const point3& origin() const { + return m_origin; + } + DETRAY_HOST_DEVICE constexpr const point3& origin_stddev() const { return m_origin_stddev; } - constexpr scalar time() const { return m_time; } - constexpr scalar charge() const { return m_charge; } + DETRAY_HOST_DEVICE constexpr scalar time() const { return m_time; } + DETRAY_HOST_DEVICE constexpr scalar charge() const { return m_charge; } /// @} }; @@ -253,38 +280,35 @@ class random_track_generator constexpr random_track_generator() = default; /// Construct from external configuration + DETRAY_HOST_DEVICE constexpr random_track_generator(configuration cfg) : m_gen{}, m_cfg(cfg) {} - /// Paramtetrized constructor for fine-grained configurations + /// Paramtetrized constructor for quick construction of simple tasks + /// + /// @note For more complex tasks, use the @c configuration type /// /// @param n_tracks the number of steps in the theta space - /// @param trk_origin the starting point of the track - /// @param origin_stddev the standard deviation of origin /// @param mom_range the range of the track momentum (in GeV) - /// @param theta_range the range for theta values - /// @param phi_range the range for phi values - /// @param time time measurement (micro seconds) /// @param charge charge of particle (e) DETRAY_HOST_DEVICE random_track_generator( - std::size_t n_tracks, point3 trk_origin = {0.f, 0.f, 0.f}, - point3 origin_stddev = {0.f, 0.f, 0.f}, + std::size_t n_tracks, std::array mom_range = {1.f * unit::GeV, 1.f * unit::GeV}, - std::array theta_range = {0.01f, constant::pi}, - std::array phi_range = {-constant::pi, - constant::pi}, - scalar time = 0.f * unit::us, - scalar charge = -1.f * unit::e, bool do_vertex_smearing = true) - : m_gen{}, m_cfg{do_vertex_smearing, n_tracks, phi_range, - theta_range, mom_range, trk_origin, - origin_stddev, time, charge} {} + scalar charge = -1.f * unit::e) + : m_gen{}, m_cfg{} { + m_cfg.n_tracks(n_tracks); + m_cfg.mom_range(mom_range); + m_cfg.charge(charge); + } /// Move constructor + DETRAY_HOST_DEVICE random_track_generator(random_track_generator&& other) : m_gen(std::move(other.m_gen)), m_cfg(std::move(other.m_cfg)) {} /// Access the configuration + DETRAY_HOST_DEVICE constexpr configuration& config() { return m_cfg; } /// @returns the generator in initial state. @@ -304,4 +328,4 @@ class random_track_generator } }; -} // namespace detray \ No newline at end of file +} // namespace detray diff --git a/utils/include/detray/simulation/event_generator/uniform_track_generator.hpp b/utils/include/detray/simulation/event_generator/uniform_track_generator.hpp index 78327bc4b..4c1382366 100644 --- a/utils/include/detray/simulation/event_generator/uniform_track_generator.hpp +++ b/utils/include/detray/simulation/event_generator/uniform_track_generator.hpp @@ -39,14 +39,23 @@ class uniform_track_generator /// Configure how tracks are generated struct configuration { + /// Ensure sensible values at the theta bounds, even in single precision + static constexpr scalar epsilon{1e-2f}; /// Range for theta and phi - std::array m_theta_range{0.01f, constant::pi}; std::array m_phi_range{-constant::pi, constant::pi}; + std::array m_theta_range{epsilon, + constant::pi - epsilon}; + std::array m_eta_range{-5.f, 5.f}; + /// Angular step size - std::size_t m_theta_steps{50u}; std::size_t m_phi_steps{50u}; + std::size_t m_theta_steps{50u}; + + /// Do uniform eta steps instead of uniform theta steps + /// (use same number of steps and range) + bool m_uniform_eta{false}; /// Track origin point3 m_origin{0.f, 0.f, 0.f}; @@ -56,39 +65,74 @@ class uniform_track_generator scalar m_p_mag{1.f * unit::GeV}; /// Time parameter and charge of the track - scalar m_time{0.f * unit::us}, m_charge{-1.f * unit::e}; + scalar m_time{0.f * unit::us}; + scalar m_charge{-1.f * unit::e}; /// Setters /// @{ - configuration& theta_range(scalar low, scalar high) { - m_theta_range = {low, high}; + DETRAY_HOST_DEVICE configuration& phi_range(scalar low, scalar high) { + assert(low <= high); + m_phi_range = {low, high}; return *this; } - configuration& phi_range(scalar low, scalar high) { - m_phi_range = {low, high}; + DETRAY_HOST_DEVICE configuration& theta_range(scalar low, scalar high) { + auto min_theta{ + std::clamp(low, epsilon, constant::pi - epsilon)}; + auto max_theta{ + std::clamp(high, epsilon, constant::pi - epsilon)}; + + assert(min_theta <= max_theta); + + m_theta_range = {min_theta, max_theta}; + m_uniform_eta = false; return *this; } - configuration& theta_steps(std::size_t n) { - m_theta_steps = n; + DETRAY_HOST_DEVICE configuration& eta_range(scalar low, scalar high) { + // This value is more or less random + constexpr auto num_max{0.001f * std::numeric_limits::max()}; + auto min_eta{low > -num_max ? low : -num_max}; + auto max_eta{high < num_max ? high : num_max}; + + assert(min_eta <= max_eta); + + m_eta_range = {min_eta, max_eta}; + m_uniform_eta = true; return *this; } - configuration& phi_steps(std::size_t n) { + DETRAY_HOST_DEVICE configuration& phi_steps(std::size_t n) { + assert(n > 0); m_phi_steps = n; return *this; } - configuration& origin(point3 ori) { + DETRAY_HOST_DEVICE configuration& theta_steps(std::size_t n) { + assert(n > 0); + m_theta_steps = n; + m_uniform_eta = false; + return *this; + } + DETRAY_HOST_DEVICE configuration& eta_steps(std::size_t n) { + assert(n > 0); + m_theta_steps = n; + m_uniform_eta = true; + return *this; + } + DETRAY_HOST_DEVICE configuration& uniform_eta(bool b) { + m_uniform_eta = b; + return *this; + } + DETRAY_HOST_DEVICE configuration& origin(point3 ori) { m_origin = ori; return *this; } - configuration& p_mag(scalar p) { + DETRAY_HOST_DEVICE configuration& p_mag(scalar p) { m_p_mag = p; return *this; } - configuration& time(scalar t) { + DETRAY_HOST_DEVICE configuration& time(scalar t) { m_time = t; return *this; } - configuration& charge(scalar q) { + DETRAY_HOST_DEVICE configuration& charge(scalar q) { m_charge = q; return *this; } @@ -96,18 +140,33 @@ class uniform_track_generator /// Getters /// @{ - constexpr std::array theta_range() const { + DETRAY_HOST_DEVICE constexpr std::array phi_range() const { + return m_phi_range; + } + DETRAY_HOST_DEVICE constexpr std::array theta_range() const { return m_theta_range; } - constexpr std::array phi_range() const { - return m_phi_range; + DETRAY_HOST_DEVICE constexpr std::array eta_range() const { + return m_eta_range; + } + DETRAY_HOST_DEVICE constexpr std::size_t phi_steps() const { + return m_phi_steps; + } + DETRAY_HOST_DEVICE constexpr std::size_t theta_steps() const { + return m_theta_steps; + } + DETRAY_HOST_DEVICE constexpr std::size_t eta_steps() const { + return m_theta_steps; + } + DETRAY_HOST_DEVICE constexpr bool uniform_eta() const { + return m_uniform_eta; + } + DETRAY_HOST_DEVICE constexpr const point3& origin() const { + return m_origin; } - constexpr std::size_t theta_steps() const { return m_theta_steps; } - constexpr std::size_t phi_steps() const { return m_phi_steps; } - constexpr const point3& origin() const { return m_origin; } - constexpr scalar p_mag() const { return m_p_mag; } - constexpr scalar time() const { return m_time; } - constexpr scalar charge() const { return m_charge; } + DETRAY_HOST_DEVICE constexpr scalar p_mag() const { return m_p_mag; } + DETRAY_HOST_DEVICE constexpr scalar time() const { return m_time; } + DETRAY_HOST_DEVICE constexpr scalar charge() const { return m_charge; } /// @} }; @@ -127,12 +186,15 @@ class uniform_track_generator constexpr iterator(configuration cfg, std::size_t iph = 1u, std::size_t ith = 0u) : m_cfg{cfg}, - m_theta_step_size{(cfg.theta_range()[1] - cfg.theta_range()[0]) / - static_cast(cfg.theta_steps())}, m_phi_step_size{(cfg.phi_range()[1] - cfg.phi_range()[0]) / static_cast(cfg.phi_steps())}, + m_theta_step_size{(cfg.theta_range()[1] - cfg.theta_range()[0]) / + static_cast(cfg.theta_steps() - 1u)}, + m_eta_step_size{(cfg.eta_range()[1] - cfg.eta_range()[0]) / + static_cast(cfg.eta_steps() - 1u)}, m_phi{cfg.phi_range()[0]}, - m_theta{cfg.theta_range()[0]}, + m_theta{cfg.uniform_eta() ? get_theta(cfg.eta_range()[0]) + : cfg.theta_range()[0]}, i_phi{iph}, i_theta{ith} {} @@ -153,8 +215,10 @@ class uniform_track_generator /// @returns the generator at its next position. DETRAY_HOST_DEVICE constexpr auto operator++() -> iterator& { + // Check theta range according to step size if (i_theta < m_cfg.theta_steps()) { + // Check phi sub-range if (i_phi < m_cfg.phi_steps()) { // Calculate new phi in the given range @@ -169,8 +233,16 @@ class uniform_track_generator // Calculate new theta in the given range ++i_theta; - m_theta = m_cfg.theta_range()[0] + - static_cast(i_theta) * m_theta_step_size; + + if (m_cfg.uniform_eta()) { + const scalar eta = + m_cfg.eta_range()[0] + + static_cast(i_theta) * m_eta_step_size; + m_theta = get_theta(eta); + } else { + m_theta = m_cfg.theta_range()[0] + + static_cast(i_theta) * m_theta_step_size; + } } return *this; } @@ -179,29 +251,38 @@ class uniform_track_generator DETRAY_HOST_DEVICE track_t operator*() const { // Momentum direction from angles - vector3 mom{math_ns::cos(m_phi) * std::sin(m_theta), - std::sin(m_phi) * std::sin(m_theta), - math_ns::cos(m_theta)}; + vector3 p{math_ns::cos(m_phi) * std::sin(m_theta), + std::sin(m_phi) * std::sin(m_theta), + math_ns::cos(m_theta)}; // Magnitude of momentum - vector::normalize(mom); - mom = m_cfg.p_mag() * mom; + vector::normalize(p); + p = m_cfg.p_mag() * p; - return track_t{m_cfg.origin(), m_cfg.time(), mom, m_cfg.charge()}; + return track_t{m_cfg.origin(), m_cfg.time(), p, m_cfg.charge()}; } /// Current configuration configuration m_cfg{}; /// Angular step sizes - scalar m_theta_step_size{0.f}; scalar m_phi_step_size{0.f}; + scalar m_theta_step_size{0.f}; + scalar m_eta_step_size{0.f}; /// Phi and theta angles of momentum direction - scalar m_phi{-constant::pi}, m_theta{0.01f}; + scalar m_phi{-constant::pi}; + scalar m_theta{configuration::epsilon}; /// Iteration indices std::size_t i_phi{0u}; std::size_t i_theta{0u}; + + private: + /// @returns the theta angle for a given @param eta value + DETRAY_HOST_DEVICE + scalar get_theta(const scalar eta) { + return 2.f * std::atan(std::exp(-eta)); + } }; configuration m_cfg{}; @@ -216,30 +297,29 @@ class uniform_track_generator DETRAY_HOST_DEVICE constexpr uniform_track_generator(configuration cfg) : m_cfg{cfg} {} - /// Paramtetrized constructor for fine-grained configurations + /// Paramtetrized constructor for quick construction of simple tasks + /// + /// @note For more complex tasks, use the @c configuration type /// /// @param n_theta the number of steps in the theta space /// @param n_phi the number of steps in the phi space - /// @param trk_origin the starting point of the track - /// @param trk_mom magnitude of the track momentum (in GeV) - /// @param theta_range the range for theta values - /// @param phi_range the range for phi values - /// @param time time measurement (micro seconds) + /// @param p_mag magnitude of the track momentum (in GeV) + /// @param uniform_eta uniformly step through eta space instead of theta /// @param charge charge of particle (e) DETRAY_HOST_DEVICE - uniform_track_generator( - std::size_t n_theta, std::size_t n_phi, - point3 trk_origin = {0.f, 0.f, 0.f}, - scalar trk_mom = 1.f * unit::GeV, - std::array theta_range = {0.01f, constant::pi}, - std::array phi_range = {-constant::pi, - constant::pi}, - scalar time = 0.f * unit::us, - scalar charge = -1.f * unit::e) - : m_cfg{theta_range, phi_range, n_theta, n_phi, - trk_origin, trk_mom, time, charge} {} + uniform_track_generator(std::size_t n_phi, std::size_t n_theta, + scalar p_mag = 1.f * unit::GeV, + bool uniform_eta = false, + scalar charge = -1.f * unit::e) + : m_cfg{} { + m_cfg.phi_steps(n_phi).theta_steps(n_theta); + m_cfg.uniform_eta(uniform_eta); + m_cfg.p_mag(p_mag); + m_cfg.charge(charge); + } /// Move constructor + DETRAY_HOST_DEVICE uniform_track_generator(uniform_track_generator&& other) : m_cfg(std::move(other.m_cfg)) {} @@ -252,6 +332,7 @@ class uniform_track_generator } /// Access the configuration + DETRAY_HOST_DEVICE constexpr configuration& config() { return m_cfg; } /// @returns the generator in initial state: Default values reflect the @@ -270,7 +351,7 @@ class uniform_track_generator /// @returns the number of tracks that will be generated DETRAY_HOST_DEVICE constexpr auto size() const noexcept -> std::size_t { - return m_cfg.theta_steps() * m_cfg.phi_steps(); + return m_cfg.phi_steps() * m_cfg.theta_steps(); } };