fix: To make traccc update work (#725)

niermann999 · web-flow · commit 9dd741988c35 · 2024-04-29T21:49:57.000Z
Some fixes to make the current traccc update work: widening the angle space for the track generators and fixing a resulting bug in the quadratic equation. Additionally fixed a bug in the navigation that could lead to a navigation abort. Also harmonizing the track generator config interface with traccc.
diff --git a/core/include/detray/grids/grid2.hpp b/core/include/detray/grids/grid2.hpp
@@ -209,7 +209,8 @@ class grid2 {
      *
      * @return the const reference to the value in this bin
      **/
-    template <typename point2_t>
+    template <typename point2_t,
+              std::enable_if_t<!std::is_scalar_v<point2_t>, bool> = true>
     DETRAY_HOST_DEVICE typename serialized_storage::const_reference bin(
         const point2_t &p2) const {
         return _data_serialized[_serializer.template serialize<axis_p0_type,
@@ -223,7 +224,8 @@ class grid2 {
      *
      * @return the const reference to the value in this bin
      **/
-    template <typename point2_t>
+    template <typename point2_t,
+              std::enable_if_t<!std::is_scalar_v<point2_t>, bool> = true>
     DETRAY_HOST_DEVICE typename serialized_storage::reference bin(
         const point2_t &p2) {
         return _data_serialized[_serializer.template serialize<axis_p0_type,
diff --git a/core/include/detray/navigation/detail/helix.hpp b/core/include/detray/navigation/detail/helix.hpp
@@ -54,7 +54,7 @@ class helix {
     ///
     /// @param pos the the origin of the helix
     /// @param time the time parameter
-    /// @param dir the initial direction for the helix
+    /// @param dir the initial normalized direction for the helix
     /// @param q the charge of the particle
     /// @param mag_field the magnetic field vector
     DETRAY_HOST_DEVICE
@@ -67,11 +67,14 @@ class helix {
         _h0 = vector::normalize(*_mag_field);
 
         // Normalized tangent vector
-        _t0 = vector::normalize(dir);
+        _t0 = dir;
+
+        assert((math::abs(getter::norm(_t0) - 1.f) < 1e-5f) &&
+               "The helix direction must be normalized");
 
         // Momentum
         const vector3_type mom =
-            1.f / static_cast<scalar_type>(math::abs(qop)) * dir;
+            1.f / static_cast<scalar_type>(math::abs(qop)) * _t0;
 
         // Normalized _h0 X _t0
         _n0 = vector::normalize(vector::cross(_h0, _t0));
diff --git a/core/include/detray/navigation/intersection/helix_cylinder_intersector.hpp b/core/include/detray/navigation/intersection/helix_cylinder_intersector.hpp
@@ -117,8 +117,15 @@ struct helix_intersector_impl<cylindrical2D<algebra_t>, algebra_t>
 
         // Obtain both possible solutions by looping over the (different)
         // starting positions
-        unsigned int n_runs =
-            math::abs(paths[0] - paths[1]) < convergence_tolerance ? 1u : 2u;
+        unsigned int n_runs = static_cast<unsigned int>(qe.solutions());
+
+        // Even if the ray is parallel to the cylinder, the helix might still
+        // hit it
+        if (qe.solutions() == 0) {
+            n_runs = 2u;
+            paths[0] = r;
+            paths[1] = -r;
+        }
         for (unsigned int i = 0u; i < n_runs; ++i) {
 
             scalar_type &s = paths[i];
diff --git a/core/include/detray/navigation/intersection/intersection.hpp b/core/include/detray/navigation/intersection/intersection.hpp
@@ -101,6 +101,7 @@ struct intersection2D {
                                     const intersection2D &is) {
         out_stream << "dist:" << is.path
                    << "\tsurface: " << is.sf_desc.barcode()
+                   << ", type: " << static_cast<int>(is.sf_desc.mask().id())
                    << ", links to vol:" << is.volume_link << ")"
                    << ", loc [" << is.local[0] << ", " << is.local[1] << ", "
                    << is.local[2] << "], ";
diff --git a/core/include/detray/navigation/navigator.hpp b/core/include/detray/navigation/navigator.hpp
@@ -610,33 +610,34 @@ class navigator {
             // Update next candidate: If not reachable, 'high trust' is broken
             if (not update_candidate(*navigation.next(), track, det, cfg)) {
                 navigation.m_status = navigation::status::e_unknown;
-                navigation.set_no_trust();
-                return;
-            }
+                navigation.set_fair_trust();
+            } else {
 
-            // Update navigation flow on the new candidate information
-            update_navigation_state(cfg, propagation);
+                // Update navigation flow on the new candidate information
+                update_navigation_state(cfg, propagation);
 
-            navigation.run_inspector(cfg, "Update complete: high trust: ");
+                navigation.run_inspector(cfg, "Update complete: high trust: ");
 
-            // The work is done if: the track has not reached a surface yet or
-            // trust is gone (portal was reached or the cache is broken).
-            if (navigation.status() == navigation::status::e_towards_object or
-                navigation.trust_level() ==
-                    navigation::trust_level::e_no_trust) {
-                return;
-            }
+                // The work is done if: the track has not reached a surface yet
+                // or trust is gone (portal was reached or the cache is broken).
+                if (navigation.status() ==
+                        navigation::status::e_towards_object or
+                    navigation.trust_level() ==
+                        navigation::trust_level::e_no_trust) {
+                    return;
+                }
 
-            // Else: Track is on module.
-            // Ready the next candidate after the current module
-            if (update_candidate(*navigation.next(), track, det, cfg)) {
-                return;
-            }
+                // Else: Track is on module.
+                // Ready the next candidate after the current module
+                if (update_candidate(*navigation.next(), track, det, cfg)) {
+                    return;
+                }
 
-            // If next candidate is not reachable, don't 'return', but
-            // escalate the trust level.
-            // This will run into the fair trust case below.
-            navigation.set_fair_trust();
+                // If next candidate is not reachable, don't 'return', but
+                // escalate the trust level.
+                // This will run into the fair trust case below.
+                navigation.set_fair_trust();
+            }
         }
 
         // Re-evaluate all currently available candidates and sort again
@@ -660,7 +661,9 @@ class navigator {
 
             navigation.run_inspector(cfg, "Update complete: fair trust: ");
 
-            return;
+            if (!navigation.is_exhausted()) {
+                return;
+            }
         }
 
         // Actor flagged cache as broken (other cases of 'no trust' are
diff --git a/core/include/detray/utils/quadratic_equation.hpp b/core/include/detray/utils/quadratic_equation.hpp
@@ -39,8 +39,12 @@ class quadratic_equation {
         const scalar_t tolerance = std::numeric_limits<scalar_t>::epsilon()) {
         // linear case
         if (math::abs(a) <= tolerance) {
-            m_solutions = 1;
-            m_values[0] = -c / b;
+            if (math::abs(b) <= tolerance) {
+                m_solutions = 0;
+            } else {
+                m_solutions = 1;
+                m_values[0] = -c / b;
+            }
         } else {
             const scalar_t discriminant{b * b - 4.f * a * c};
             // If there is more than one solution, then a != 0 and q != 0
diff --git a/tests/integration_tests/cpu/detectors/toy_detector_navigation.cpp b/tests/integration_tests/cpu/detectors/toy_detector_navigation.cpp
@@ -62,6 +62,7 @@ int main(int argc, char **argv) {
     cfg_hel_scan.name("toy_detector_helix_scan");
     cfg_hel_scan.whiteboard(white_board);
     cfg_hel_scan.track_generator().n_tracks(10000u);
+    cfg_hel_scan.track_generator().eta_range(-4.f, 4.f);
     cfg_hel_scan.track_generator().p_T(1.f * unit<scalar_t>::GeV);
 
     detail::register_checks<test::helix_scan>(toy_det, toy_names, cfg_hel_scan);
@@ -83,8 +84,6 @@ int main(int argc, char **argv) {
     cfg_hel_nav.name("toy_detector_helix_navigation");
     cfg_hel_nav.whiteboard(white_board);
     cfg_hel_nav.propagation().navigation.search_window = {3u, 3u};
-    // For one surface the toy detector seems to need a stricter tolerance
-    cfg_hel_nav.propagation().navigation.min_mask_tolerance = 1e-5f;
 
     detail::register_checks<test::helix_navigation>(toy_det, toy_names,
                                                     cfg_hel_nav);
diff --git a/tests/integration_tests/cpu/detectors/wire_chamber_navigation.cpp b/tests/integration_tests/cpu/detectors/wire_chamber_navigation.cpp
@@ -62,16 +62,17 @@ int main(int argc, char **argv) {
     cfg_hel_scan.name("wire_chamber_helix_scan");
     cfg_hel_scan.whiteboard(white_board);
     cfg_hel_scan.track_generator().n_tracks(10000u);
+    cfg_hel_scan.track_generator().eta_range(-1.f, 1.f);
     // TODO: Fails for smaller momenta
-    cfg_hel_scan.track_generator().p_T(3.f * unit<scalar_t>::GeV);
+    cfg_hel_scan.track_generator().p_T(5.f * unit<scalar_t>::GeV);
 
     detail::register_checks<test::helix_scan>(det, names, cfg_hel_scan);
 
     // Comparison of straight line navigation with ray scan
     test::straight_line_navigation<wire_chamber_t>::config cfg_str_nav{};
     cfg_str_nav.name("wire_chamber_straight_line_navigation");
     cfg_str_nav.whiteboard(white_board);
-    cfg_str_nav.propagation().navigation.search_window = {2u, 2u};
+    cfg_str_nav.propagation().navigation.search_window = {3u, 3u};
     auto mask_tolerance = cfg_ray_scan.mask_tolerance();
     cfg_str_nav.propagation().navigation.min_mask_tolerance = mask_tolerance[0];
     cfg_str_nav.propagation().navigation.max_mask_tolerance = mask_tolerance[1];
diff --git a/tests/integration_tests/device/cuda/propagator_cuda.cpp b/tests/integration_tests/device/cuda/propagator_cuda.cpp
@@ -35,6 +35,7 @@ TEST_P(CudaPropConstBFieldMng, propagator) {
     propagator_test_config cfg{};
     cfg.track_generator.phi_steps(20).theta_steps(20);
     cfg.track_generator.p_tot(10.f * unit<scalar_t>::GeV);
+    cfg.track_generator.eta_range(-3.f, 3.f);
     cfg.propagation.navigation.search_window = {3u, 3u};
     // Configuration for non-z-aligned B-fields
     cfg.propagation.navigation.overstep_tolerance = std::get<0>(GetParam());
@@ -69,6 +70,7 @@ TEST_P(CudaPropConstBFieldCpy, propagator) {
     propagator_test_config cfg{};
     cfg.track_generator.phi_steps(20u).theta_steps(20u);
     cfg.track_generator.p_tot(10.f * unit<scalar_t>::GeV);
+    cfg.track_generator.eta_range(-3.f, 3.f);
     cfg.propagation.navigation.search_window = {3u, 3u};
     // Configuration for non-z-aligned B-fields
     cfg.propagation.navigation.overstep_tolerance = std::get<0>(GetParam());
diff --git a/tests/unit_tests/cpu/simulation/detector_scanner.cpp b/tests/unit_tests/cpu/simulation/detector_scanner.cpp
@@ -26,7 +26,7 @@ using namespace detray;
 using algebra_t = test::algebra;
 using vector3 = test::vector3;
 
-constexpr const scalar tol{1e-3f};
+constexpr const scalar tol{1e-7f};
 
 /// Brute force test: Intersect toy geometry and compare between ray and helix
 /// without B-field
@@ -74,10 +74,13 @@ GTEST_TEST(detray_simulation, detector_scanner) {
 
         // Should have encountered the same number of tracks (vulnerable to
         // floating point errors)
-        EXPECT_EQ(expected[n_tracks].size(), intersection_trace.size());
+        EXPECT_EQ(expected[n_tracks].size(), intersection_trace.size())
+            << test_helix;
 
         // Check every single recorded intersection
-        for (std::size_t i = 0u; i < intersection_trace.size(); ++i) {
+        for (std::size_t i = 0u;
+             i < std::min(expected[n_tracks].size(), intersection_trace.size());
+             ++i) {
             if (expected[n_tracks][i].vol_idx !=
                 intersection_trace[i].vol_idx) {
                 // Intersection record at portal bound might be flipped
diff --git a/tests/unit_tests/cpu/simulation/track_generators.cpp b/tests/unit_tests/cpu/simulation/track_generators.cpp
@@ -29,26 +29,24 @@ GTEST_TEST(detray_simulation, uniform_track_generator) {
         uniform_track_generator<free_track_parameters<algebra_t>>;
 
     constexpr const scalar_t tol{1e-5f};
-    constexpr const scalar_t epsilon{generator_t::configuration::epsilon};
+    constexpr const scalar_t max_pi{generator_t::configuration::k_max_pi};
 
     constexpr std::size_t phi_steps{50u};
     constexpr std::size_t theta_steps{50u};
 
     std::array<vector3, phi_steps * theta_steps> momenta{};
 
-    // Loop over theta values ]0,pi[
+    // Loop over theta values [0,pi)
     for (std::size_t itheta{0u}; itheta < theta_steps; ++itheta) {
-        const scalar_t theta{epsilon +
-                             static_cast<scalar_t>(itheta) *
-                                 (constant<scalar_t>::pi - 2.f * epsilon) /
-                                 static_cast<scalar_t>(theta_steps - 1u)};
+        const scalar_t theta{static_cast<scalar_t>(itheta) * max_pi /
+                             static_cast<scalar_t>(theta_steps - 1u)};
 
-        // Loop over phi values [-pi, pi]
+        // Loop over phi values [-pi, pi)
         for (std::size_t iphi{0u}; iphi < phi_steps; ++iphi) {
             // The direction
             const scalar_t phi{-constant<scalar_t>::pi +
                                static_cast<scalar_t>(iphi) *
-                                   (2.f * constant<scalar_t>::pi) /
+                                   (constant<scalar_t>::pi + max_pi) /
                                    static_cast<scalar_t>(phi_steps)};
 
             // intialize a track
@@ -125,6 +123,7 @@ GTEST_TEST(detray_simulation, uniform_track_generator_eta) {
         uniform_track_generator<free_track_parameters<algebra_t>>;
 
     constexpr const scalar_t tol{1e-5f};
+    constexpr const scalar_t max_pi{generator_t::configuration::k_max_pi};
 
     constexpr std::size_t phi_steps{50u};
     constexpr std::size_t eta_steps{50u};
@@ -137,12 +136,12 @@ GTEST_TEST(detray_simulation, uniform_track_generator_eta) {
                                       static_cast<scalar_t>(eta_steps - 1u)};
         const scalar_t theta{2.f * std::atan(std::exp(-eta))};
 
-        // Loop over phi values [-pi, pi]
+        // Loop over phi values [-pi, pi)
         for (std::size_t iphi{0u}; iphi < phi_steps; ++iphi) {
             // The direction
             const scalar_t phi{-constant<scalar_t>::pi +
                                static_cast<scalar_t>(iphi) *
-                                   (2.f * constant<scalar_t>::pi) /
+                                   (constant<scalar_t>::pi + max_pi) /
                                    static_cast<scalar_t>(phi_steps)};
 
             // intialize a track
@@ -224,7 +223,7 @@ GTEST_TEST(detray_simulation, random_track_generator_uniform) {
         random_track_generator<free_track_parameters<algebra_t>, uniform_gen_t>;
 
     // Tolerance depends on sample size
-    constexpr scalar_t tol{0.05f};
+    constexpr scalar_t tol{0.02f};
 
     // Track counter
     std::size_t n_tracks{0u};
@@ -236,6 +235,7 @@ GTEST_TEST(detray_simulation, random_track_generator_uniform) {
     trk_gen_cfg.seed(42u);
     trk_gen_cfg.phi_range(-0.9f * constant<scalar_t>::pi,
                           0.8f * constant<scalar_t>::pi);
+    trk_gen_cfg.eta_range(-4.f, 4.f);
     trk_gen_cfg.mom_range(1.f * unit<scalar_t>::GeV, 2.f * unit<scalar_t>::GeV);
     trk_gen_cfg.origin_stddev({0.1f * unit<scalar_t>::mm,
                                0.f * unit<scalar_t>::mm,
@@ -268,6 +268,8 @@ GTEST_TEST(detray_simulation, random_track_generator_uniform) {
     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();
+    ASSERT_NEAR(theta_range[0], 0.0366f, tol);
+    ASSERT_NEAR(theta_range[1], 3.105f, tol);
     const auto& mom_range = trk_gen_cfg.mom_range();
 
     // Mean
@@ -304,7 +306,7 @@ GTEST_TEST(detray_simulation, random_track_generator_normal) {
         random_track_generator<free_track_parameters<algebra_t>, normal_gen_t>;
 
     // Tolerance depends on sample size
-    constexpr scalar_t tol{0.05f};
+    constexpr scalar_t tol{0.02f};
 
     // Track counter
     std::size_t n_tracks{0u};
@@ -316,6 +318,7 @@ GTEST_TEST(detray_simulation, random_track_generator_normal) {
     trk_gen_cfg.seed(42u);
     trk_gen_cfg.phi_range(-0.9f * constant<scalar_t>::pi,
                           0.8f * constant<scalar_t>::pi);
+    trk_gen_cfg.eta_range(-4.f, 4.f);
     trk_gen_cfg.mom_range(1.f * unit<scalar_t>::GeV, 2.f * unit<scalar_t>::GeV);
     trk_gen_cfg.origin({0.f, 0.f, 0.f});
     trk_gen_cfg.origin_stddev({0.1f * unit<scalar_t>::mm,
@@ -349,6 +352,8 @@ GTEST_TEST(detray_simulation, random_track_generator_normal) {
     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();
+    ASSERT_NEAR(theta_range[0], 0.0366f, tol);
+    ASSERT_NEAR(theta_range[1], 3.105f, tol);
     const auto& mom_range = trk_gen_cfg.mom_range();
 
     // Mean
diff --git a/utils/include/detray/simulation/event_generator/random_track_generator.hpp b/utils/include/detray/simulation/event_generator/random_track_generator.hpp
diff --git a/utils/include/detray/simulation/event_generator/uniform_track_generator.hpp b/utils/include/detray/simulation/event_generator/uniform_track_generator.hpp
