ref: Reuse ray/helix scan results and switch to random track generator (#713)

This PR adds a simple whiteboard to the validation so that the ray and helix scan traces can be reused by the navigation tests, which reduces execution time of the tests a lot. Also allows to add new scans in the future by moving the surface loop into a distinct functor.

Switches the navigation validation to use the random_track_generator to avoid effects from detector symmetry and cleans that generator ups a bit (mostly setting a default seed that is not zero).

With the random tracks, the number of mismatches between the Newton intersector and the navigator increased, so I adapted the mask tolerances in both of them dynamically. In the Newton solver, the an estimation of the approximation error gets projected onto the tangential surface of the mask at the track position and in the navigator, the mask tolerance gets scaled by the distance of the track to the surface (per mille of the path) and then clamped to a configurable min-max range. That has the advantage, that at long distances to the surface, when the ray is a poor approximation to the curved track, the mask tolerance is very large, while it goes to zero as we approach the candidate, ensuring a precise decision. The minimal mask tolerance is the machine epsilon by default, but be set in a propagation executable e.g. to the mean position uncertainty of the stepper or to what the CKF needs.

Fixes a few smaller bugs and tightens the portals around the sensitive surfaces in the toy detector for better navigation. Clamps the phi values in the track generators to -pi, pi.

The detector integration tests now run 10000 helices at p_T=1GeV, except for the wire chamber which starts failing for p_T<3Gev. That needs some more investigation.

I also discovered, that the CUDA propagation test was essentially not testing anything, since the host propagation silently failed. I set that test back up, but found discrepancies in the number of track position in single precision (double precision ran fine), so I configured fewer tracks and added step size constraints until it worked. That test should probably get an overhaul at some point.

Author: niermann999

core/include/detray/navigation/detail/helix.hpp

@@ -67,7 +67,7 @@ class helix {
         _h0 = vector::normalize(*_mag_field);
 
         // Normalized tangent vector
-        _t0 = dir;
+        _t0 = vector::normalize(dir);
 
         // Momentum
         const vector3_type mom =

core/include/detray/navigation/intersection/helix_cylinder_intersector.hpp

@@ -62,9 +62,14 @@ struct helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t>
     DETRAY_HOST_DEVICE inline std::array<intersection_type<surface_descr_t>, 2>
     operator()(const helix_type &h, const surface_descr_t &sf_desc,
                const mask_t &mask, const transform3_type &trf,
-               const scalar_type mask_tolerance = 0.f,
+               const std::array<scalar_type, 2u> mask_tolerance =
+                   {detail::invalid_value<scalar_type>(),
+                    detail::invalid_value<scalar_type>()},
                const scalar_type = 0.f) const {
 
+        assert((mask_tolerance[0] == mask_tolerance[1]) &&
+               "Helix intersectors use only one mask tolerance value");
+
         std::array<intersection_type<surface_descr_t>, 2> ret;
 
         // Guard against inifinite loops
@@ -117,7 +122,7 @@ struct helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t>
         for (unsigned int i = 0u; i < n_runs; ++i) {
 
             scalar_type &s = paths[i];
-            intersection_type<surface_descr_t> &is = ret[i];
+            intersection_type<surface_descr_t> &sfi = ret[i];
 
             // Path length in the previous iteration step
             scalar_type s_prev{0.f};
@@ -150,26 +155,48 @@ struct helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t>
             }
 
             // Build intersection struct from helix parameters
-            is.path = s;
+            sfi.path = s;
             const auto p3 = h.pos(s);
-            is.local = mask.to_local_frame(trf, p3);
-            is.status = mask.is_inside(is.local, mask_tolerance);
+            sfi.local = mask.to_local_frame(trf, p3);
+            sfi.cos_incidence_angle = vector::dot(
+                mask.local_frame().normal(trf, sfi.local), h.dir(s));
+
+            scalar_type tol{mask_tolerance[1]};
+            if (detail::is_invalid_value(tol)) {
+                // Due to floating point errors this can be negative if cos ~ 1
+                const scalar_type sin_inc2{math::abs(
+                    1.f - sfi.cos_incidence_angle * sfi.cos_incidence_angle)};
+
+                tol = math::abs((s - s_prev) * math::sqrt(sin_inc2));
+            }
+            sfi.status = mask.is_inside(sfi.local, tol);
 
             // Compute some additional information if the intersection is valid
-            if (is.status == intersection::status::e_inside) {
-                is.sf_desc = sf_desc;
-                is.direction = math::signbit(s)
-                                   ? intersection::direction::e_opposite
-                                   : intersection::direction::e_along;
-                is.volume_link = mask.volume_link();
+            if (sfi.status == intersection::status::e_inside) {
+                sfi.sf_desc = sf_desc;
+                sfi.direction = math::signbit(s)
+                                    ? intersection::direction::e_opposite
+                                    : intersection::direction::e_along;
+                sfi.volume_link = mask.volume_link();
             }
         }
 
         return ret;
     }
 
+    /// Interface to use fixed mask tolerance
+    template <typename surface_descr_t, typename mask_t>
+    DETRAY_HOST_DEVICE inline std::array<intersection_type<surface_descr_t>, 2>
+    operator()(const helix_type &h, const surface_descr_t &sf_desc,
+               const mask_t &mask, const transform3_type &trf,
+               const scalar_type mask_tolerance,
+               const scalar_type = 0.f) const {
+        return this->operator()(h, sf_desc, mask, trf, {mask_tolerance, 0.f},
+                                0.f);
+    }
+
     /// Tolerance for convergence
-    scalar_type convergence_tolerance{1e-3f};
+    scalar_type convergence_tolerance{1.f * unit<scalar_type>::um};
 };
 
 template <typename algebra_t>

core/include/detray/navigation/intersection/helix_line_intersector.hpp

@@ -54,9 +54,15 @@ struct helix_intersector_impl<line2D<algebra_t>, algebra_t> {
     template <typename surface_descr_t, typename mask_t>
     DETRAY_HOST_DEVICE inline intersection_type<surface_descr_t> operator()(
         const helix_type &h, const surface_descr_t &sf_desc, const mask_t &mask,
-        const transform3_type &trf, const scalar_type mask_tolerance = 0.f,
+        const transform3_type &trf,
+        const std::array<scalar_type, 2u> mask_tolerance =
+            {detail::invalid_value<scalar_type>(),
+             detail::invalid_value<scalar_type>()},
         const scalar_type = 0.f) const {
 
+        assert((mask_tolerance[0] == mask_tolerance[1]) &&
+               "Helix intersectors use only one mask tolerance value");
+
         intersection_type<surface_descr_t> sfi;
 
         // Guard against inifinite loops
@@ -147,9 +153,17 @@ struct helix_intersector_impl<line2D<algebra_t>, algebra_t> {
         // Build intersection struct from helix parameters
         sfi.path = s;
         sfi.local = mask.to_local_frame(trf, h.pos(s), h.dir(s));
-
-        const point3_type local = mask.to_local_frame(trf, h.pos(s), h.dir(s));
-        sfi.status = mask.is_inside(local, mask_tolerance);
+        sfi.cos_incidence_angle =
+            vector::dot(mask.local_frame().normal(trf, sfi.local), h.dir(s));
+        scalar_type tol{mask_tolerance[1]};
+        if (detail::is_invalid_value(tol)) {
+            // Due to floating point errors this can be negative if cos ~ 1
+            const scalar_type sin_inc2{math::abs(
+                1.f - sfi.cos_incidence_angle * sfi.cos_incidence_angle)};
+
+            tol = math::abs((s - s_prev) * math::sqrt(sin_inc2));
+        }
+        sfi.status = mask.is_inside(sfi.local, tol);
 
         // Compute some additional information if the intersection is valid
         if (sfi.status == intersection::status::e_inside) {
@@ -163,8 +177,18 @@ struct helix_intersector_impl<line2D<algebra_t>, algebra_t> {
         return sfi;
     }
 
+    /// Interface to use fixed mask tolerance
+    template <typename surface_descr_t, typename mask_t>
+    DETRAY_HOST_DEVICE inline intersection_type<surface_descr_t> operator()(
+        const helix_type &h, const surface_descr_t &sf_desc, const mask_t &mask,
+        const transform3_type &trf, const scalar_type mask_tolerance,
+        const scalar_type = 0.f) const {
+        return this->operator()(h, sf_desc, mask, trf, {mask_tolerance, 0.f},
+                                0.f);
+    }
+
     /// Tolerance for convergence
-    scalar_type convergence_tolerance{1e-3f};
+    scalar_type convergence_tolerance{1.f * unit<scalar_type>::um};
 };
 
 }  // namespace detray

core/include/detray/navigation/intersection/helix_plane_intersector.hpp

@@ -56,9 +56,15 @@ struct helix_intersector_impl<cartesian2D<algebra_t>, algebra_t> {
     template <typename surface_descr_t, typename mask_t>
     DETRAY_HOST_DEVICE inline intersection_type<surface_descr_t> operator()(
         const helix_type &h, const surface_descr_t &sf_desc, const mask_t &mask,
-        const transform3_type &trf, const scalar_type mask_tolerance = 0.f,
+        const transform3_type &trf,
+        const std::array<scalar_type, 2u> mask_tolerance =
+            {detail::invalid_value<scalar_type>(),
+             detail::invalid_value<scalar_type>()},
         const scalar_type = 0.f) const {
 
+        assert((mask_tolerance[0] == mask_tolerance[1]) &&
+               "Helix intersectors use only one mask tolerance value");
+
         intersection_type<surface_descr_t> sfi;
 
         // Guard against inifinite loops
@@ -100,7 +106,18 @@ struct helix_intersector_impl<cartesian2D<algebra_t>, algebra_t> {
         // Build intersection struct from helix parameters
         sfi.path = s;
         sfi.local = mask.to_local_frame(trf, h.pos(s), h.dir(s));
-        sfi.status = mask.is_inside(sfi.local, mask_tolerance);
+        sfi.cos_incidence_angle =
+            vector::dot(mask.local_frame().normal(trf, sfi.local), h.dir(s));
+
+        scalar_type tol{mask_tolerance[1]};
+        if (detail::is_invalid_value(tol)) {
+            // Due to floating point errors this can be negative if cos ~ 1
+            const scalar_type sin_inc2{math::abs(
+                1.f - sfi.cos_incidence_angle * sfi.cos_incidence_angle)};
+
+            tol = math::abs((s - s_prev) * math::sqrt(sin_inc2));
+        }
+        sfi.status = mask.is_inside(sfi.local, tol);
 
         // Compute some additional information if the intersection is valid
         if (sfi.status == intersection::status::e_inside) {
@@ -114,8 +131,18 @@ struct helix_intersector_impl<cartesian2D<algebra_t>, algebra_t> {
         return sfi;
     }
 
+    /// Interface to use fixed mask tolerance
+    template <typename surface_descr_t, typename mask_t>
+    DETRAY_HOST_DEVICE inline intersection_type<surface_descr_t> operator()(
+        const helix_type &h, const surface_descr_t &sf_desc, const mask_t &mask,
+        const transform3_type &trf, const scalar_type mask_tolerance,
+        const scalar_type = 0.f) const {
+        return this->operator()(h, sf_desc, mask, trf, {mask_tolerance, 0.f},
+                                0.f);
+    }
+
     /// Tolerance for convergence
-    scalar_type convergence_tolerance{1e-3f};
+    scalar_type convergence_tolerance{1.f * unit<scalar_type>::um};
 };
 
 template <typename algebra_t>

core/include/detray/navigation/intersection/ray_concentric_cylinder_intersector.hpp

@@ -55,7 +55,9 @@ struct ray_concentric_cylinder_intersector {
     template <typename surface_descr_t, typename mask_t>
     DETRAY_HOST_DEVICE inline intersection_type<surface_descr_t> operator()(
         const ray_type &ray, const surface_descr_t &sf, const mask_t &mask,
-        const transform3_type & /*trf*/, const scalar_type mask_tolerance = 0.f,
+        const transform3_type & /*trf*/,
+        const std::array<scalar_type, 2u> mask_tolerance =
+            {0.f, 1.f * unit<scalar_type>::mm},
         const scalar_type overstep_tol = 0.f) const {
 
         intersection_type<surface_descr_t> is;
@@ -110,7 +112,11 @@ struct ray_concentric_cylinder_intersector {
                 is.path = t01[cindex];
                 // In this case, the point has to be in cylinder3 coordinates
                 // for the r-check
-                is.status = mask.is_inside(is.local, mask_tolerance);
+                // Tolerance: per mille of the distance
+                is.status = mask.is_inside(
+                    is.local, math::max(mask_tolerance[0],
+                                        math::min(mask_tolerance[1],
+                                                  1e-3f * math::abs(is.path))));
 
                 // prepare some additional information in case the intersection
                 // is valid
@@ -131,6 +137,16 @@ struct ray_concentric_cylinder_intersector {
         return is;
     }
 
+    /// Interface to use fixed mask tolerance
+    template <typename surface_descr_t, typename mask_t>
+    DETRAY_HOST_DEVICE inline intersection_type<surface_descr_t> operator()(
+        const ray_type &ray, const surface_descr_t &sf, const mask_t &mask,
+        const transform3_type &trf, const scalar_type mask_tolerance,
+        const scalar_type overstep_tol = 0.f) const {
+        return this->operator()(ray, sf, mask, trf, {mask_tolerance, 0.f},
+                                overstep_tol);
+    }
+
     /// Operator function to find intersections between a ray and a 2D cylinder
     ///
     /// @tparam mask_t is the input mask type
@@ -146,7 +162,8 @@ struct ray_concentric_cylinder_intersector {
     DETRAY_HOST_DEVICE inline void update(
         const ray_type &ray, intersection_type<surface_descr_t> &sfi,
         const mask_t &mask, const transform3_type &trf,
-        const scalar_type mask_tolerance = 0.f,
+        const std::array<scalar_type, 2u> &mask_tolerance =
+            {0.f, 1.f * unit<scalar_type>::mm},
         const scalar_type overstep_tol = 0.f) const {
         sfi = this->operator()(ray, sfi.sf_desc, mask, trf, mask_tolerance,
                                overstep_tol)[0];

core/include/detray/navigation/intersection/ray_cylinder_intersector.hpp

@@ -57,7 +57,8 @@ struct ray_intersector_impl<cylindrical2D<algebra_t>, algebra_t> {
     DETRAY_HOST_DEVICE inline std::array<intersection_type<surface_descr_t>, 2>
     operator()(const ray_type &ray, const surface_descr_t &sf,
                const mask_t &mask, const transform3_type &trf,
-               const scalar_type mask_tolerance = 0.f,
+               const std::array<scalar_type, 2u> mask_tolerance =
+                   {0.f, 100.f * unit<scalar_type>::um},
                const scalar_type overstep_tol = 0.f) const {
 
         // One or both of these solutions might be invalid
@@ -89,6 +90,17 @@ struct ray_intersector_impl<cylindrical2D<algebra_t>, algebra_t> {
         return ret;
     }
 
+    /// Interface to use fixed mask tolerance
+    template <typename surface_descr_t, typename mask_t>
+    DETRAY_HOST_DEVICE inline std::array<intersection_type<surface_descr_t>, 2>
+    operator()(const ray_type &ray, const surface_descr_t &sf,
+               const mask_t &mask, const transform3_type &trf,
+               const scalar_type mask_tolerance,
+               const scalar_type overstep_tol = 0.f) const {
+        return this->operator()(ray, sf, mask, trf, {mask_tolerance, 0.f},
+                                overstep_tol);
+    }
+
     /// Operator function to find intersections between a ray and a 2D cylinder
     ///
     /// @tparam mask_t is the input mask type
@@ -103,7 +115,8 @@ struct ray_intersector_impl<cylindrical2D<algebra_t>, algebra_t> {
     DETRAY_HOST_DEVICE inline void update(
         const ray_type &ray, intersection_type<surface_descr_t> &sfi,
         const mask_t &mask, const transform3_type &trf,
-        const scalar_type mask_tolerance = 0.f,
+        const std::array<scalar_type, 2u> mask_tolerance =
+            {0.f, 1.f * unit<scalar_type>::mm},
         const scalar_type overstep_tol = 0.f) const {
 
         // One or both of these solutions might be invalid
@@ -155,8 +168,12 @@ struct ray_intersector_impl<cylindrical2D<algebra_t>, algebra_t> {
     DETRAY_HOST_DEVICE inline intersection_type<surface_descr_t>
     build_candidate(const ray_type &ray, mask_t &mask,
                     const transform3_type &trf, const scalar_type path,
-                    const scalar_type mask_tolerance = 0.f,
-                    const scalar_type overstep_tol = 0.f) const {
+                    const std::array<scalar_type, 2u> mask_tolerance,
+                    const scalar_type overstep_tol) const {
+
+        assert((mask_tolerance[0] <= mask_tolerance[1]) &&
+               "Minimal mask tolerance needs to be smaller or equal maximal "
+               "mask tolerance");
 
         intersection_type<surface_descr_t> is;
 
@@ -170,7 +187,11 @@ struct ray_intersector_impl<cylindrical2D<algebra_t>, algebra_t> {
             const point3_type p3 = ro + is.path * rd;
 
             is.local = mask.to_local_frame(trf, p3);
-            is.status = mask.is_inside(is.local, mask_tolerance);
+            // Tolerance: per mille of the distance
+            is.status = mask.is_inside(
+                is.local, math::max(mask_tolerance[0],
+                                    math::min(mask_tolerance[1],
+                                              1e-3f * math::abs(is.path))));
 
             // prepare some additional information in case the intersection
             // is valid

core/include/detray/navigation/intersection/ray_cylinder_portal_intersector.hpp

@@ -62,7 +62,9 @@ struct ray_intersector_impl<concentric_cylindrical2D<algebra_t>, algebra_t>
     template <typename surface_descr_t, typename mask_t>
     DETRAY_HOST_DEVICE inline intersection_type<surface_descr_t> operator()(
         const ray_type &ray, const surface_descr_t &sf, const mask_t &mask,
-        const transform3_type &trf, const scalar_type mask_tolerance = 0.f,
+        const transform3_type &trf,
+        const std::array<scalar_type, 2u> mask_tolerance =
+            {0.f, 1.f * unit<scalar_type>::mm},
         const scalar_type overstep_tol = 0.f) const {
 
         intersection_type<surface_descr_t> is;
@@ -87,6 +89,16 @@ struct ray_intersector_impl<concentric_cylindrical2D<algebra_t>, algebra_t>
         return is;
     }
 
+    /// Interface to use fixed mask tolerance
+    template <typename surface_descr_t, typename mask_t>
+    DETRAY_HOST_DEVICE inline intersection_type<surface_descr_t> operator()(
+        const ray_type &ray, const surface_descr_t &sf, const mask_t &mask,
+        const transform3_type &trf, const scalar_type mask_tolerance,
+        const scalar_type overstep_tol = 0.f) const {
+        return this->operator()(ray, sf, mask, trf, {mask_tolerance, 0.f},
+                                overstep_tol);
+    }
+
     /// Operator function to find intersections between a ray and a 2D cylinder
     ///
     /// @tparam mask_t is the input mask type
@@ -101,7 +113,8 @@ struct ray_intersector_impl<concentric_cylindrical2D<algebra_t>, algebra_t>
     DETRAY_HOST_DEVICE inline void update(
         const ray_type &ray, intersection_type<surface_descr_t> &sfi,
         const mask_t &mask, const transform3_type &trf,
-        const scalar_type mask_tolerance = 0.f,
+        const std::array<scalar_type, 2u> &mask_tolerance =
+            {0.f, 1.f * unit<scalar_type>::mm},
         const scalar_type overstep_tol = 0.f) const {
         sfi = this->operator()(ray, sfi.sf_desc, mask, trf, mask_tolerance,
                                overstep_tol);