Merge pull request #560 from niermann999/feat-trackgen-eta

feat: regular eta steps for uniform track generator

core/include/detray/propagator/navigator.hpp

@@ -126,7 +126,7 @@ class navigator {
 
             sf.template visit_mask<intersection_initialize>(
                 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<intersection_update>(
-            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:

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<free_track_parameters<test::transform3>>;
+
+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<intersection2D<typename detector_t::surface_type,
+                               typename detector_t::transform3>>
+        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<intersection2D<typename detector_t::surface_type,
-                                   typename detector_t::transform3>>
-            intersections{};
-        std::size_t n_surfaces{0u};
-
-        // Iterate through uniformly distributed momentum directions
-        for (const auto track :
-             uniform_track_generator<free_track_parameters<test::transform3>>(
-                 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()) {

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<detail::ray<test::transform3>>;
+
 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<rectangle2D<>> 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<detail::ray<test::transform3>>(
-                 theta_steps, phi_steps, ori, 1.f)) {
+        for (const auto ray : ray_generator) {
 
             for (const auto &plane : planes) {
                 auto pi = rect.intersector<intersection2D<
@@ -115,16 +118,16 @@ void BM_INTERSECT_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<detail::ray<test::transform3>>(
-                 theta_steps, phi_steps, ori, 1.f)) {
+        for (const auto ray : ray_generator) {
 
             for (const auto &cylinder : cylinders) {
                 auto ci = cylinder.intersector<intersection_t>();
@@ -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<detail::ray<test::transform3>>(
-                 theta_steps, phi_steps, ori, 1.f)) {
+        for (const auto ray : ray_generator) {
 
             for (const auto &cylinder : cylinders) {
                 auto cpi = cylinder.intersector<intersection_t>();
@@ -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<detail::ray<test::transform3>>(
-                 theta_steps, phi_steps, ori, 1.f)) {
+        for (const auto ray : ray_generator) {
 
             for (const auto &cylinder : cylinders) {
                 auto cci = cylinder.intersector<intersection_t>();

tests/benchmarks/cuda/benchmark_propagator_cuda.cpp

@@ -31,13 +31,12 @@ toy_det_config toy_cfg{4u, 7u};
 
 void fill_tracks(vecmem::vector<free_track_parameters<transform3>> &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<scalar>::GeV};
 
     // Iterate through uniformly distributed momentum directions
     for (auto traj : uniform_track_generator<free_track_parameters<transform3>>(
-             theta_steps, phi_steps, ori, mom_mag)) {
+             phi_steps, theta_steps, mom_mag)) {
         tracks.push_back(traj);
     }
 }

tests/common/include/tests/common/test_base/propagator_test.hpp

@@ -139,14 +139,11 @@ inline vecmem::vector<track_t> generate_tracks(
     // Track collection
     vecmem::vector<track_t> 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<scalar>::GeV};
 
     // Iterate through uniformly distributed momentum directions
-    for (auto track : uniform_track_generator<track_t>(
-             ts, ps, ori, p_mag, {0.01f, constant<scalar>::pi},
-             {-constant<scalar>::pi, constant<scalar>::pi})) {
+    for (auto track : uniform_track_generator<track_t>(ps, ts, p_mag)) {
         track.set_overstep_tolerance(overstep_tolerance);
 
         // Put it into vector of trajectories

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<scalar_type>::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;

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<intersection_kernel_t>(
-                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) {

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<std::vector<std::pair<dindex, intersection_t>>> expected;
     //  Iterate through uniformly distributed momentum directions with ray
     for (const auto test_ray :
-         uniform_track_generator<detail::ray<transform3_type>>(
-             theta_steps, phi_steps, ori)) {
+         uniform_track_generator<detail::ray<transform3_type>>(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<free_track_parameters<transform3_type>>(
-             theta_steps, phi_steps, ori)) {
+             phi_steps, theta_steps)) {
         const detail::helix test_helix(track, &B);
 
         // Record all intersections and objects along the ray

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<scalar>::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<track_t>(theta_steps, phi_steps, ori, mom)) {
+         uniform_track_generator<track_t>(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<track_t>(theta_steps, phi_steps, ori, mom)) {
+         uniform_track_generator<track_t>(phi_steps, theta_steps, mom)) {
         // Genrate second track state used for propagation with pathlimit
         track_t lim_track(track);
 

tests/unit_tests/cpu/tools_stepper.cpp

@@ -166,15 +166,14 @@ GTEST_TEST(detray_propagator, rk_stepper) {
     constexpr scalar stepsize_constr{0.5f * unit<scalar>::mm};
 
     // Track generator configuration
-    const point3 ori{0.f, 0.f, 0.f};
     const scalar p_mag{10.f * unit<scalar>::GeV};
     constexpr unsigned int theta_steps = 100u;
     constexpr unsigned int phi_steps = 100u;
 
     // Iterate through uniformly distributed momentum directions
     for (auto track :
          uniform_track_generator<free_track_parameters<transform3>>(
-             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,15 +270,14 @@ TEST(detray_propagator, rk_stepper_inhomogeneous_bfield) {
     constexpr scalar stepsize_constr{0.5f * unit<scalar>::mm};
 
     // Track generator configuration
-    const point3 ori{0.f, 0.f, 0.f};
     const scalar p_mag{10.f * unit<scalar>::GeV};
     constexpr unsigned int theta_steps = 100u;
     constexpr unsigned int phi_steps = 100u;
 
     // Iterate through uniformly distributed momentum directions
     for (auto track :
          uniform_track_generator<free_track_parameters<transform3>>(
-             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);