cpu/detectors/telescope_detector.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/definitions/detail/qualifiers.hpp"
0011 #include "detray/definitions/units.hpp"
0012 #include "detray/geometry/mask.hpp"
0013 #include "detray/geometry/shapes/unbounded.hpp"
0014 #include "detray/navigation/caching_navigator.hpp"
0015 #include "detray/propagator/line_stepper.hpp"
0016 #include "detray/propagator/propagation_config.hpp"
0017 #include "detray/propagator/rk_stepper.hpp"
0018 #include "detray/tracks/tracks.hpp"
0019 #include "detray/utils/consistency_checker.hpp"
0020
0021 // Detray test include(s)
0022 #include "detray/test/common/build_telescope_detector.hpp"
0023 #include "detray/test/framework/types.hpp"
0024 #include "detray/test/utils/inspectors.hpp"
0025
0026 // Vecmem include(s)
0027 #include <vecmem/memory/host_memory_resource.hpp>
0028
0029 // Covfie include(s)
0030 #include <covfie/core/backend/primitive/constant.hpp>
0031 #include <covfie/core/field.hpp>
0032 #include <covfie/core/vector.hpp>
0033
0034 // GTest include
0035 #include <gtest/gtest.h>
0036
0037 // System include(s)
0038 #include <utility>
0039
0040 namespace detray {
0041
0042 namespace {
0043
0044 using test_algebra = test::algebra;
0045 using scalar = test::scalar;
0046 using vector3 = test::vector3;
0047
0048 // dummy propagator state
0049 template <typename stepping_t, typename navigation_t>
0050 struct prop_state {
0051   using context_t = typename navigation_t::detector_type::geometry_context;
0052
0053   template <typename track_t, typename field_type>
0054   prop_state(const track_t &t_in, const field_type &field,
0055              const typename navigation_t::detector_type &det,
0056              const context_t &ctx = {})
0057       : m_stepping(t_in, field), m_navigation(det), m_context(ctx) {}
0058
0059   constexpr navigation_t &navigation() { return m_navigation; }
0060   constexpr stepping_t &stepping() { return m_stepping; }
0061
0062   stepping_t m_stepping;
0063   navigation_t m_navigation;
0064   context_t m_context;
0065 };
0066
0067 inline constexpr bool verbose_check = true;
0068
0069 const propagation::config prop_cfg{};
0070
0071 }  // anonymous namespace
0072
0073 }  // namespace detray
0074
0075 // This tests the construction and general methods of the navigator
0076 GTEST_TEST(detray_detectors, telescope_detector) {
0077   using namespace detray;
0078
0079   // Use rectangle surfaces
0080   mask<rectangle2D, test_algebra> rectangle{0u, 20.f * unit<scalar>::mm,
0081                                             20.f * unit<scalar>::mm};
0082   tel_det_config<test_algebra> tel_cfg{rectangle};
0083
0084   using const_bfield_bknd_t =
0085       covfie::backend::constant<covfie::vector::vector_d<scalar, 3>,
0086                                 covfie::vector::vector_d<scalar, 3>>;
0087   using b_field_t = covfie::field<const_bfield_bknd_t>;
0088
0089   using rk_stepper_t = rk_stepper<b_field_t::view_t, test_algebra>;
0090   using inspector_t = navigation::print_inspector;
0091   constexpr std::size_t cache_size{navigation::default_cache_size};
0092
0093   // Test tolerance
0094   constexpr scalar tol{1e-4f};
0095
0096   vecmem::host_memory_resource host_mr;
0097
0098   // B-fields
0099   vector3 B_z{static_cast<scalar>(0.f), static_cast<scalar>(0.f),
0100               1.f * unit<scalar>::T};
0101   vector3 B_x{1.f * unit<scalar>::T, static_cast<scalar>(0.f),
0102               static_cast<scalar>(0.f)};
0103   b_field_t b_field_z{covfie::make_parameter_pack(
0104       b_field_t::backend_t::configuration_t{B_z[0], B_z[1], B_z[2]})};
0105   b_field_t b_field_x{covfie::make_parameter_pack(
0106       b_field_t::backend_t::configuration_t{B_x[0], B_x[1], B_x[2]})};
0107   b_field_t::view_t b_field_z_view(b_field_z);
0108   b_field_t::view_t b_field_x_view(b_field_x);
0109
0110   // steppers
0111   rk_stepper_t rk_stepper_z;
0112   rk_stepper_t rk_stepper_x;
0113
0114   //
0115   // telescope along z
0116   //
0117
0118   // Build from given module positions
0119   std::vector<scalar> positions = {0.f,   50.f,  100.f, 150.f, 200.f, 250.f,
0120                                    300.f, 350.f, 400.f, 450.f, 500.f};
0121   // Build telescope detector with unbounded planes
0122   const auto [z_tel_det1, z_tel_names1] =
0123       build_telescope_detector<test_algebra>(host_mr,
0124                                              tel_cfg.positions(positions));
0125
0126   // Some general checks
0127   const auto vol0 = tracking_volume{z_tel_det1, 0u};
0128   ASSERT_EQ(vol0.portals().size(), 6u);
0129   ASSERT_EQ(vol0.surfaces().size(), positions.size() + 6u);
0130   ASSERT_EQ(vol0.template surfaces<surface_id::e_sensitive>().size(),
0131             positions.size());
0132
0133   // Test this only once, it is the same for all telescope detectors
0134   EXPECT_EQ(z_tel_names1.get_detector_name(), "telescope_detector");
0135   EXPECT_EQ(z_tel_names1.at(0u), "telescope_world_0");
0136
0137   // Check general consistency of the detector
0138   detail::check_consistency(z_tel_det1, verbose_check, z_tel_names1);
0139
0140   // Build the same telescope detector with rectangular planes and given
0141   // length/number of surfaces
0142   tel_cfg.positions({}).n_surfaces(11u).length(500.f * unit<scalar>::mm);
0143   const auto [z_tel_det2, z_tel_names2] =
0144       build_telescope_detector<test_algebra>(host_mr, tel_cfg);
0145
0146   // Check general consistency of the detector
0147   detail::check_consistency(z_tel_det2, verbose_check, z_tel_names2);
0148
0149   // Compare
0150   for (std::size_t i{0u}; i < z_tel_det1.surfaces().size(); ++i) {
0151     geometry::identifier geo_id{};
0152     geo_id.set_volume(0u).set_index(i);
0153     geo_id.set_id((i == z_tel_det1.surfaces().size() - 1u)
0154                       ? surface_id::e_portal
0155                       : surface_id::e_sensitive);
0156     EXPECT_TRUE(z_tel_det1.surface(geo_id) == z_tel_det2.surface(geo_id));
0157   }
0158
0159   //
0160   // telescope along x
0161   //
0162
0163   // Same telescope, but in x direction and created from custom stepper
0164   detail::ray<test_algebra> x_track({0.f, 0.f, 0.f}, 0.f, {1.f, 0.f, 0.f},
0165                                     -1.f);
0166
0167   const auto [x_tel_det, x_tel_names] = build_telescope_detector<test_algebra>(
0168       host_mr, tel_cfg.pilot_track(x_track));
0169
0170   // Check general consistency of the detector
0171   detail::check_consistency(x_tel_det, verbose_check, x_tel_names);
0172
0173   //
0174   // test propagation in all telescope detector instances
0175   //
0176
0177   // Telescope navigation should be symmetric in x and z
0178   vector3 pos = {0.f, 0.f, 0.f};
0179   vector3 mom = {0.f, 0.f, 1.f};
0180   free_track_parameters<test_algebra> test_track_z1(pos, 0.f, mom, -1.f);
0181   free_track_parameters<test_algebra> test_track_z2(pos, 0.f, mom, -1.f);
0182   mom = {1.f, 0.f, 0.f};
0183   free_track_parameters<test_algebra> test_track_x(pos, 0.f, mom, -1.f);
0184
0185   // navigators
0186   caching_navigator<decltype(z_tel_det1), cache_size, inspector_t> navigator_z1;
0187   caching_navigator<decltype(z_tel_det2), cache_size, inspector_t> navigator_z2;
0188   caching_navigator<decltype(x_tel_det), cache_size, inspector_t> navigator_x;
0189   using navigation_state_t = decltype(navigator_z1)::state;
0190   using stepping_state_t = rk_stepper_t::state;
0191
0192   // propagation states
0193   prop_state<stepping_state_t, navigation_state_t> propgation_z1(
0194       test_track_z1, b_field_z_view, z_tel_det1);
0195   prop_state<stepping_state_t, navigation_state_t> propgation_z2(
0196       test_track_z2, b_field_z_view, z_tel_det2);
0197   prop_state<stepping_state_t, navigation_state_t> propgation_x(
0198       test_track_x, b_field_x_view, x_tel_det);
0199
0200   stepping_state_t &stepping_z1 = propgation_z1.stepping();
0201   stepping_state_t &stepping_z2 = propgation_z2.stepping();
0202   stepping_state_t &stepping_x = propgation_x.stepping();
0203
0204   navigation_state_t &navigation_z1 = propgation_z1.navigation();
0205   navigation_state_t &navigation_z2 = propgation_z2.navigation();
0206   navigation_state_t &navigation_x = propgation_x.navigation();
0207
0208   // propagate all telescopes
0209   navigator_z1.init(stepping_z1(), navigation_z1, prop_cfg.navigation,
0210                     prop_cfg.context);
0211   navigator_z2.init(stepping_z2(), navigation_z2, prop_cfg.navigation,
0212                     prop_cfg.context);
0213   navigator_x.init(stepping_x(), navigation_x, prop_cfg.navigation,
0214                    prop_cfg.context);
0215
0216   bool heartbeat_z1 = navigation_z1.is_alive();
0217   bool heartbeat_z2 = navigation_z2.is_alive();
0218   bool heartbeat_x = navigation_x.is_alive();
0219
0220   bool do_reset_z1{true};
0221   bool do_reset_z2{true};
0222   bool do_reset_x{true};
0223
0224   while (heartbeat_z1 && heartbeat_z2 && heartbeat_x) {
0225     // check that all propagation flows are still running
0226     EXPECT_TRUE(heartbeat_z1);
0227     EXPECT_TRUE(heartbeat_z2);
0228     EXPECT_TRUE(heartbeat_x);
0229
0230     heartbeat_z1 =
0231         heartbeat_z1 && rk_stepper_z.step(navigation_z1(), stepping_z1,
0232                                           prop_cfg.stepping, do_reset_z1);
0233     heartbeat_z2 =
0234         heartbeat_z2 && rk_stepper_z.step(navigation_z2(), stepping_z2,
0235                                           prop_cfg.stepping, do_reset_z2);
0236     heartbeat_x =
0237         heartbeat_x && rk_stepper_x.step(navigation_x(), stepping_x,
0238                                          prop_cfg.stepping, do_reset_x);
0239
0240     navigation_z1.set_high_trust();
0241     navigation_z2.set_high_trust();
0242     navigation_x.set_high_trust();
0243
0244     do_reset_z1 = navigator_z1.update(stepping_z1(), navigation_z1,
0245                                       prop_cfg.navigation, prop_cfg.context);
0246     do_reset_z2 = navigator_z2.update(stepping_z2(), navigation_z2,
0247                                       prop_cfg.navigation, prop_cfg.context);
0248     do_reset_x = navigator_x.update(stepping_x(), navigation_x,
0249                                     prop_cfg.navigation, prop_cfg.context);
0250
0251     // Also reset when reached a surface
0252     do_reset_z1 = do_reset_z1 || navigation_z1.is_on_surface();
0253     do_reset_z2 = do_reset_z2 || navigation_z2.is_on_surface();
0254     do_reset_x = do_reset_x || navigation_x.is_on_surface();
0255
0256     heartbeat_z1 = heartbeat_z1 && navigation_z1.is_alive();
0257     heartbeat_z2 = heartbeat_z2 && navigation_z2.is_alive();
0258     heartbeat_x = heartbeat_x && navigation_x.is_alive();
0259
0260     // The track path lengths should match between all propagations
0261     EXPECT_NEAR(
0262         std::abs(stepping_z1.path_length() - stepping_z2.path_length()) /
0263             stepping_z1.path_length(),
0264         0.f, tol);
0265     EXPECT_NEAR(std::abs(stepping_z1.path_length() - stepping_x.path_length()) /
0266                     stepping_x.path_length(),
0267                 0.f, tol);
0268     // The track positions in z should match exactly
0269     EXPECT_NEAR(vector::norm(stepping_z1().pos() - stepping_z2().pos()) /
0270                     vector::norm(stepping_z1().pos()),
0271                 0.f, tol);
0272     EXPECT_NEAR(vector::norm(stepping_z1().dir() - stepping_z2().dir()) /
0273                     vector::norm(stepping_z1().dir()),
0274                 0.f, tol);
0275   }
0276
0277   // check that all propagation flows exited successfully
0278   ASSERT_TRUE(navigation_z1.finished())
0279       << navigation_z1.inspector().to_string();
0280   ASSERT_TRUE(navigation_z2.finished())
0281       << navigation_z2.inspector().to_string();
0282   ASSERT_TRUE(navigation_x.finished()) << navigation_x.inspector().to_string();
0283
0284   //
0285   // Build a telescope along a bent track
0286   //
0287   pos = {0.f, 0.f, 0.f};
0288   mom = {0.f, 1.f, 0.f};
0289
0290   auto pilot_track = free_track_parameters<test_algebra>(pos, 0.f, mom, -1.f);
0291
0292   detail::helix<test_algebra> helix_bz(pilot_track, B_z);
0293
0294   tel_det_config htel_cfg{rectangle, helix_bz};
0295   htel_cfg.n_surfaces(11u).length(500.f * unit<scalar>::mm);
0296   const auto [tel_detector, tel_names] =
0297       build_telescope_detector<test_algebra>(host_mr, htel_cfg);
0298
0299   // Check general consistency of the detector
0300   detail::check_consistency(tel_detector, verbose_check, tel_names);
0301
0302   // make at least sure it is navigatable
0303   caching_navigator<decltype(tel_detector), cache_size, inspector_t>
0304       tel_navigator;
0305
0306   prop_state<stepping_state_t, navigation_state_t> tel_propagation(
0307       pilot_track, b_field_z_view, tel_detector);
0308   navigation_state_t &tel_navigation = tel_propagation.navigation();
0309   stepping_state_t &tel_stepping = tel_propagation.stepping();
0310
0311   // run propagation
0312   tel_navigator.init(tel_stepping(), tel_navigation, prop_cfg.navigation,
0313                      prop_cfg.context);
0314   bool heartbeat_tel = tel_navigation.is_alive();
0315
0316   bool do_reset_tel{true};
0317
0318   while (heartbeat_tel) {
0319     heartbeat_tel =
0320         heartbeat_tel && rk_stepper_z.step(tel_navigation(), tel_stepping,
0321                                            prop_cfg.stepping, do_reset_tel);
0322
0323     tel_navigation.set_high_trust();
0324
0325     do_reset_tel = tel_navigator.update(tel_stepping(), tel_navigation,
0326                                         prop_cfg.navigation, prop_cfg.context);
0327
0328     do_reset_tel = do_reset_tel || navigation_z1.is_on_surface();
0329     heartbeat_tel = heartbeat_tel && tel_navigation.is_alive();
0330   }
0331   // check that propagation was successful
0332   ASSERT_TRUE(tel_navigation.finished())
0333       << tel_navigation.inspector().to_string();
0334 }