Enable regular steps in eta for uniform_track_generator (this is need…

…ed 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

@@ -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);