navigation/intersection/intersection_kernel.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 // Project include(s)
0010 #include "detray/navigation/detail/intersection_kernel.hpp"
0011
0012 #include "detray/core/detail/multi_store.hpp"
0013 #include "detray/core/detail/single_store.hpp"
0014 #include "detray/geometry/mask.hpp"
0015 #include "detray/geometry/shapes.hpp"
0016 #include "detray/geometry/surface_descriptor.hpp"
0017 #include "detray/navigation/intersection/helix_intersector.hpp"
0018 #include "detray/navigation/intersection/ray_intersector.hpp"
0019 #include "detray/tracks/tracks.hpp"
0020 #include "detray/tracks/trajectories.hpp"
0021 #include "detray/utils/ranges.hpp"
0022
0023 // Detray test include(s)
0024 #include "detray/test/framework/types.hpp"
0025
0026 // Vecmem include(s)
0027 #include <vecmem/memory/host_memory_resource.hpp>
0028
0029 // Google Test include(s)
0030 #include <gtest/gtest.h>
0031
0032 using namespace detray;
0033
0034 using test_algebra = test::algebra;
0035 using scalar = test::scalar;
0036 using vector3 = test::vector3;
0037 using point3 = test::point3;
0038 using transform_link_t = dindex;
0039
0040 namespace {
0041
0042 constexpr const scalar is_close{1e-5f};
0043 constexpr const scalar external_tol{1.f * unit<scalar>::mm};
0044
0045 const intersection::config intr_cfg{
0046     .min_mask_tolerance = 1e-3f * unit<float>::mm,
0047     .max_mask_tolerance = 1e-3f * unit<float>::mm,
0048     .mask_tolerance_scalor = 0.f,
0049     .overstep_tolerance = 0.f};
0050
0051 enum class mask_id : unsigned int {
0052   e_rectangle2D = 0u,
0053   e_trapezoid2D = 1u,
0054   e_annulus2D = 2u,
0055   e_cylinder2D = 3u,
0056   e_concentric_cylinder2D = 4u,
0057 };
0058
0059 enum class material_id : unsigned int {
0060   e_material_slab = 0u,
0061 };
0062
0063 }  // anonymous namespace
0064
0065 /// Surface components:
0066 using volume_link_t = std::uint_least16_t;
0067
0068 /// - masks, with mask identifiers 0,1,2
0069 using rectangle_t = mask<rectangle2D, test_algebra, volume_link_t>;
0070 using trapezoid_t = mask<trapezoid2D, test_algebra, volume_link_t>;
0071 using annulus_t = mask<annulus2D, test_algebra, volume_link_t>;
0072 using cylinder_t = mask<cylinder2D, test_algebra, volume_link_t>;
0073 using cylinder_portal_t =
0074     mask<concentric_cylinder2D, test_algebra, volume_link_t>;
0075
0076 using mask_container_t =
0077     regular_multi_store<mask_id, empty_context, dtuple, dvector, rectangle_t,
0078                         trapezoid_t, annulus_t, cylinder_t, cylinder_portal_t>;
0079 using mask_link_t = typename mask_container_t::single_link;
0080 using material_link_t = dtyped_index<material_id, dindex>;
0081
0082 using transform_container_t =
0083     single_store<test::transform3, dvector, geometry_context>;
0084
0085 /// The Surface definition:
0086 using surface_t =
0087     surface_descriptor<mask_link_t, material_link_t, transform_link_t>;
0088 using surface_container_t = dvector<surface_t>;
0089
0090 // TODO: How about merging ray and helix tests into one to remove the code
0091 // repetition?
0092
0093 // This tests the construction of a surface
0094 GTEST_TEST(detray_intersection, intersection_kernel_ray) {
0095   vecmem::host_memory_resource host_mr;
0096
0097   // The transforms & their store
0098   typename transform_container_t::context_type static_context{};
0099   transform_container_t transform_store;
0100   // Transforms of the rectangle, trapezoid, annulus and the two cylinders
0101   transform_store.emplace_back(static_context, vector3{0.f, 0.f, 10.f});
0102   transform_store.emplace_back(static_context, vector3{0.f, 0.f, 20.f});
0103   transform_store.emplace_back(static_context, vector3{0.f, -20.f, 30.f});
0104   // 90deg rotation around y-axis
0105   transform_store.emplace_back(static_context, vector3{0.f, 0.f, 50.f},
0106                                vector3{1.f, 0.f, 0.f}, vector3{0.f, 0.f, -1.f});
0107   // Identity transform for concentric cylinder
0108   transform_store.emplace_back(static_context, vector3{0.f, 0.f, 0.f});
0109   // Shifted rectangle
0110   transform_store.emplace_back(static_context, vector3{12.f, 0.f, 1000.f});
0111
0112   // The masks & their store
0113   mask_container_t mask_store(host_mr);
0114
0115   const rectangle_t rect{0u, 10.f, 10.f};
0116   const trapezoid_t trap{0u, 10.f, 20.f, 30.f, 1.f / 60.f};
0117   const annulus_t annl{0u, 15.f, 55.f, 0.75f, 1.95f, 0.f, 2.f, -2.f};
0118   const cylinder_t cyl{0u, 5.f, -10.f, 10.f};
0119   const cylinder_portal_t cyl_portal{0u, 1.f, 0.f, 1000.f};
0120
0121   mask_store.template push_back<mask_id::e_rectangle2D>(rect, empty_context{});
0122   mask_store.template push_back<mask_id::e_rectangle2D>(rect, empty_context{});
0123   mask_store.template push_back<mask_id::e_trapezoid2D>(trap, empty_context{});
0124   mask_store.template push_back<mask_id::e_annulus2D>(annl, empty_context{});
0125   mask_store.template push_back<mask_id::e_cylinder2D>(cyl, empty_context{});
0126   mask_store.template push_back<mask_id::e_concentric_cylinder2D>(
0127       cyl_portal, empty_context{});
0128
0129   // The surfaces and their store
0130   surface_t rectangle_surface(0u, {mask_id::e_rectangle2D, 0u},
0131                               {material_id::e_material_slab, 0u}, 0u,
0132                               surface_id::e_sensitive);
0133   surface_t trapezoid_surface(1u, {mask_id::e_trapezoid2D, 0u},
0134                               {material_id::e_material_slab, 1u}, 0u,
0135                               surface_id::e_sensitive);
0136   surface_t annulus_surface(2u, {mask_id::e_annulus2D, 0u},
0137                             {material_id::e_material_slab, 2u}, 0u,
0138                             surface_id::e_sensitive);
0139   surface_t cyl_surface(3u, {mask_id::e_cylinder2D, 0u},
0140                         {material_id::e_material_slab, 2u}, 0u,
0141                         surface_id::e_passive);
0142   surface_t cyl_portal_surface(4u, {mask_id::e_concentric_cylinder2D, 0u},
0143                                {material_id::e_material_slab, 2u}, 0u,
0144                                surface_id::e_portal);
0145   surface_t rectangle_shifted_surface(5u, {mask_id::e_rectangle2D, 1u},
0146                                       {material_id::e_material_slab, 0u}, 0u,
0147                                       surface_id::e_sensitive);
0148   // Easier test debugging
0149   rectangle_surface.set_index(0u);
0150   trapezoid_surface.set_index(1u);
0151   annulus_surface.set_index(2u);
0152   cyl_surface.set_index(3u);
0153   cyl_portal_surface.set_index(4u);
0154   rectangle_shifted_surface.set_index(5u);
0155
0156   surface_container_t surfaces = {
0157       rectangle_surface, trapezoid_surface,  annulus_surface,
0158       cyl_surface,       cyl_portal_surface, rectangle_shifted_surface};
0159
0160   const point3 pos{0.f, 0.f, 0.f};
0161   const vector3 mom{0.01f, 0.01f, 10.f};
0162   const free_track_parameters<test_algebra> track(pos, 0.f, mom, -1.f);
0163
0164   // Validation data
0165   const point3 expected_rectangle{0.01f, 0.01f, 10.f};
0166   const point3 expected_trapezoid{0.02f, 0.02f, 20.f};
0167   const point3 expected_annulus{0.03f, 0.03f, 30.f};
0168   const point3 expected_cylinder1{0.045f, 0.045f, 45.0f};
0169   const point3 expected_cylinder2{0.055f, 0.055f, 55.0f};
0170   const point3 expected_cylinder_pt{0.7071068f, 0.7071068f, 707.1068f};
0171   const point3 expected_rectangle_shifted{1.f, 1.f, 1000.f};
0172
0173   const std::vector<point3> expected_points = {
0174       expected_rectangle,        expected_trapezoid, expected_annulus,
0175       expected_cylinder1,        expected_cylinder2, expected_cylinder_pt,
0176       expected_rectangle_shifted};
0177
0178   // Initialize kernel
0179   std::vector<
0180       intersection2D<surface_t, test_algebra, intersection::contains_pos>>
0181       sfi_init;
0182   sfi_init.reserve(expected_points.size());
0183
0184   for (const auto &surface : surfaces) {
0185     mask_store.visit<detail::intersection_initialize<ray_intersector>>(
0186         surface.mask(), sfi_init, detail::ray(track), surface, transform_store,
0187         static_context, intr_cfg, external_tol);
0188   }
0189
0190   EXPECT_EQ(expected_points.size(), sfi_init.size());
0191
0192   // Also check intersections
0193   for (std::size_t i = 0u;
0194        i < math::min(expected_points.size(), sfi_init.size()); ++i) {
0195     EXPECT_TRUE(sfi_init[i].is_along());
0196     EXPECT_EQ(sfi_init[i].volume_link(), 0u);
0197
0198     vector3 global{0.f, 0.f, 0.f};
0199
0200     const surface_t sf_desc = sfi_init[i].surface();
0201     if (sf_desc.mask().id() == mask_id::e_rectangle2D) {
0202       global = rect.to_global_frame(transform_store.at(sf_desc.transform()),
0203                                     sfi_init[i].local());
0204     } else if (sf_desc.mask().id() == mask_id::e_trapezoid2D) {
0205       global = trap.to_global_frame(transform_store.at(sf_desc.transform()),
0206                                     sfi_init[i].local());
0207     } else if (sf_desc.mask().id() == mask_id::e_annulus2D) {
0208       global = annl.to_global_frame(transform_store.at(sf_desc.transform()),
0209                                     sfi_init[i].local());
0210     } else if (sf_desc.mask().id() == mask_id::e_cylinder2D) {
0211       global = cyl.to_global_frame(transform_store.at(sf_desc.transform()),
0212                                    sfi_init[i].local());
0213     } else if (sf_desc.mask().id() == mask_id::e_concentric_cylinder2D) {
0214       global = cyl_portal.to_global_frame(
0215           transform_store.at(sf_desc.transform()), sfi_init[i].local());
0216     }
0217
0218     EXPECT_NEAR(global[0], expected_points[i][0], 1e-3f)
0219         << " at point " << i << ": surface " << sf_desc.identifier();
0220     EXPECT_NEAR(global[1], expected_points[i][1], 1e-3f)
0221         << " at point " << i << ": surface " << sf_desc.identifier();
0222     EXPECT_NEAR(global[2], expected_points[i][2], 1e-3f)
0223         << " at point " << i << ": surface " << sf_desc.identifier();
0224   }
0225
0226   // Update kernel
0227   // @fixme: The intersection update kernel does not work for non-portal
0228   // cylinders, since it assigns the closest intersection to both
0229   // solutions
0230   /*std::vector<intersection2D_point<surface_t, test_algebra>> sfi_update;
0231   sfi_update.resize(5);
0232
0233   for (const auto [idx, surface] : detray::views::enumerate(surfaces)) {
0234       sfi_update[idx].set_surface(surface);
0235       mask_store.visit<detail::intersection_update>(
0236           surface.mask(), detail::ray(track), sfi_update[idx],
0237           transform_store);
0238
0239       if(!sfi_update[idx].is_inside()) {
0240   continue; } ASSERT_TRUE(sfi_update[idx].is_along()) << " at surface " <<
0241   sfi_update[idx]
0242   << ", " << sfi_init[idx]; ASSERT_EQ(sfi_update[idx].volume_link(), 0u);
0243       ASSERT_NEAR(sfi_update[idx].p3[0], expected_points[idx][0],
0244   is_close)
0245       << " at surface " << sfi_update[idx] << ", " << sfi_init[idx];
0246       ASSERT_NEAR(sfi_update[idx].p3[1], expected_points[idx][1],
0247   is_close) << " at surface " << sfi_update[idx] << ", " << sfi_init[idx];
0248       ASSERT_NEAR(sfi_update[idx].p3[2], expected_points[idx][2],
0249   is_close)
0250       << " at surface " << sfi_update[idx] << ", " << sfi_init[idx];
0251   }
0252
0253   // Compare
0254   ASSERT_EQ(sfi_init.size(), 5u);
0255   ASSERT_EQ(sfi_update.size(), 5u);
0256   for (unsigned int i = 0u; i < 5u; i++) {
0257       ASSERT_EQ(sfi_init[i].p3, sfi_update[i].p3);
0258       ASSERT_EQ(sfi_init[i].local(), sfi_update[i].local());
0259       ASSERT_EQ(sfi_init[i].path(), sfi_update[i].path());
0260   }*/
0261 }
0262
0263 /// Reuse the intersection kernel test for particle gun
0264 GTEST_TEST(detray_intersection, intersection_kernel_helix) {
0265   vecmem::host_memory_resource host_mr;
0266
0267   // The transforms & their store
0268   typename transform_container_t::context_type static_context{};
0269   transform_container_t transform_store;
0270   // Transforms of the rectangle, trapezoid and annulus
0271   transform_store.emplace_back(static_context, vector3{0.f, 0.f, 10.f});
0272   transform_store.emplace_back(static_context, vector3{0.f, 0.f, 20.f});
0273   transform_store.emplace_back(static_context, vector3{0.f, -20.f, 30.f});
0274   // The masks & their store
0275   mask_container_t mask_store(host_mr);
0276
0277   const rectangle_t rect{0u, 10.f, 10.f};
0278   const trapezoid_t trap{0u, 10.f, 20.f, 30.f, 1.f / 60.f};
0279   const annulus_t annl{0u, 15.f, 55.f, 0.75f, 1.95f, 0.f, 2.f, -2.f};
0280   mask_store.template push_back<mask_id::e_rectangle2D>(rect, empty_context{});
0281   mask_store.template push_back<mask_id::e_trapezoid2D>(trap, empty_context{});
0282   mask_store.template push_back<mask_id::e_annulus2D>(annl, empty_context{});
0283
0284   // The surfaces and their store
0285   const surface_t rectangle_surface(0u, {mask_id::e_rectangle2D, 0u},
0286                                     {material_id::e_material_slab, 0u}, 0u,
0287                                     surface_id::e_sensitive);
0288   const surface_t trapezoid_surface(1u, {mask_id::e_trapezoid2D, 0u},
0289                                     {material_id::e_material_slab, 1u}, 0u,
0290                                     surface_id::e_sensitive);
0291   const surface_t annulus_surface(2u, {mask_id::e_annulus2D, 0u},
0292                                   {material_id::e_material_slab, 2u}, 0u,
0293                                   surface_id::e_sensitive);
0294   surface_container_t surfaces = {rectangle_surface, trapezoid_surface,
0295                                   annulus_surface};
0296   const point3 pos{0.f, 0.f, 0.f};
0297   const vector3 mom{0.01f, 0.01f, 10.f};
0298   const vector3 B{0.f * unit<scalar>::T, 0.f * unit<scalar>::T,
0299                   is_close * unit<scalar>::T};
0300   const detail::helix<test_algebra> h({pos, 0.f, mom, -1.f}, B);
0301
0302   // Validation data
0303   const point3 expected_rectangle{0.01f, 0.01f, 10.f};
0304   const point3 expected_trapezoid{0.02f, 0.02f, 20.f};
0305   const point3 expected_annulus{0.03f, 0.03f, 30.f};
0306   const std::vector<point3> expected_points = {
0307       expected_rectangle, expected_trapezoid, expected_annulus};
0308   std::vector<
0309       intersection2D<surface_t, test_algebra, intersection::contains_pos>>
0310       sfi_helix{};
0311   sfi_helix.reserve(expected_points.size());
0312
0313   // Try the intersections - with automated dispatching via the kernel
0314   for (const auto [sf_idx, surface] : detray::views::enumerate(surfaces)) {
0315     mask_store.visit<detail::intersection_initialize<helix_intersector>>(
0316         surface.mask(), sfi_helix, h, surface, transform_store, static_context,
0317         intr_cfg, scalar{0.f});
0318
0319     vector3 global{0.f, 0.f, 0.f};
0320
0321     if (surface.mask().id() == mask_id::e_rectangle2D) {
0322       global =
0323           rect.to_global_frame(transform_store.at(0), sfi_helix[0].local());
0324     } else if (surface.mask().id() == mask_id::e_trapezoid2D) {
0325       global =
0326           trap.to_global_frame(transform_store.at(1), sfi_helix[0].local());
0327     } else if (surface.mask().id() == mask_id::e_annulus2D) {
0328       global =
0329           annl.to_global_frame(transform_store.at(2), sfi_helix[0].local());
0330     }
0331
0332     ASSERT_NEAR(global[0], expected_points[sf_idx][0], is_close);
0333     ASSERT_NEAR(global[1], expected_points[sf_idx][1], is_close);
0334     ASSERT_NEAR(global[2], expected_points[sf_idx][2], is_close);
0335
0336     sfi_helix.clear();
0337   }
0338 }