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(); } };