cpu/propagator/propagation_caching_navigator.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/units.hpp"
0011 #include "detray/geometry/tracking_surface.hpp"
0012 #include "detray/navigation/caching_navigator.hpp"
0013 #include "detray/propagator/actors.hpp"
0014 #include "detray/propagator/base_actor.hpp"
0015 #include "detray/propagator/line_stepper.hpp"
0016 #include "detray/propagator/propagator.hpp"
0017 #include "detray/propagator/rk_stepper.hpp"
0018 #include "detray/tracks/tracks.hpp"
0019 #include "detray/tracks/trajectories.hpp"
0020 #include "detray/utils/logging.hpp"
0021
0022 // Detray test include(s)
0023 #include "detray/test/common/bfield.hpp"
0024 #include "detray/test/common/build_toy_detector.hpp"
0025 #include "detray/test/common/build_wire_chamber.hpp"
0026 #include "detray/test/common/track_generators.hpp"
0027 #include "detray/test/framework/types.hpp"
0028 #include "detray/test/utils/inspectors.hpp"
0029
0030 // Vecmem include(s)
0031 #include <vecmem/memory/host_memory_resource.hpp>
0032
0033 // GTest include(s)
0034 #include <gtest/gtest.h>
0035
0036 using namespace detray;
0037
0038 using test_algebra = test::algebra;
0039 using scalar = test::scalar;
0040 using point3 = test::point3;
0041 using vector3 = test::vector3;
0042
0043 constexpr scalar tol{1e-3f};
0044 constexpr scalar path_limit{5.f * unit<scalar>::cm};
0045 constexpr std::size_t cache_size{navigation::default_cache_size};
0046
0047 namespace detray::actor {
0048
0049 /// Compare helical track positions for stepper
0050 struct helix_inspector : public base_actor {
0051   /// Keeps the state of a helix gun to calculate track positions
0052   struct state {
0053     scalar path_from_surface{0.f};
0054   };
0055
0056   /// Check that the stepper remains on the right helical track for its pos.
0057   template <typename propagator_state_t>
0058   DETRAY_HOST_DEVICE void operator()(
0059       state& inspector_state, const propagator_state_t& prop_state,
0060       parameter_transporter_result<test_algebra>& res) const {
0061     const auto& navigation = prop_state.navigation();
0062     const auto& stepping = prop_state.stepping();
0063     const auto& bound_params = res.destination_params();
0064
0065     typename propagator_state_t::detector_type::geometry_context ctx{};
0066
0067     // Nothing has happened yet (first call of actor chain)
0068     if (stepping.path_length() < tol ||
0069         inspector_state.path_from_surface < tol) {
0070       return;
0071     }
0072
0073     if (bound_params.surface_link().is_invalid()) {
0074       return;
0075     }
0076
0077     // Surface
0078     const auto sf =
0079         tracking_surface{navigation.detector(), bound_params.surface_link()};
0080
0081     const free_track_parameters<test_algebra> free_params =
0082         sf.bound_to_free_vector(ctx, bound_params);
0083
0084     const auto last_pos = free_params.pos();
0085
0086     const auto bvec =
0087         stepping.field().at(last_pos[0], last_pos[1], last_pos[2]);
0088     const vector3 b{bvec[0], bvec[1], bvec[2]};
0089
0090     detail::helix<test_algebra> hlx(free_params, b);
0091
0092     const auto true_pos = hlx(inspector_state.path_from_surface);
0093
0094     const point3 relative_error{1.f / inspector_state.path_from_surface *
0095                                 (stepping().pos() - true_pos)};
0096
0097     ASSERT_NEAR(vector::norm(relative_error), 0.f, tol);
0098
0099     auto true_J = hlx.jacobian(inspector_state.path_from_surface);
0100
0101     for (unsigned int i = 0u; i < e_free_size; i++) {
0102       for (unsigned int j = 0u; j < e_free_size; j++) {
0103         ASSERT_NEAR(getter::element(stepping.transport_jacobian(), i, j),
0104                     getter::element(true_J, i, j),
0105                     inspector_state.path_from_surface * tol * 10.f);
0106       }
0107     }
0108     // The propagation does not start on a surface, skip the initial path
0109     if (!bound_params.surface_link().is_invalid()) {
0110       inspector_state.path_from_surface += stepping.step_size();
0111     }
0112     // Reset path from surface
0113     if (navigation.is_on_sensitive() || navigation.encountered_sf_material()) {
0114       inspector_state.path_from_surface = 0.f;
0115     }
0116   }
0117 };
0118
0119 }  // namespace detray::actor
0120
0121 /// Test basic functionality of the propagator using a straight line stepper
0122 GTEST_TEST(detray_propagator, propagator_line_stepper) {
0123   vecmem::host_memory_resource host_mr;
0124   toy_det_config<scalar> toy_cfg{};
0125   toy_cfg.use_material_maps(false);
0126   const auto [d, names] = build_toy_detector<test_algebra>(host_mr, toy_cfg);
0127
0128   using navigator_t =
0129       caching_navigator<decltype(d), cache_size, navigation::print_inspector>;
0130   using stepper_t = line_stepper<test_algebra>;
0131   using propagator_t = propagator<stepper_t, navigator_t, actor_chain<>>;
0132
0133   const point3 pos{0.f, 0.f, 0.f};
0134   const vector3 mom{1.f, 1.f, 0.f};
0135   free_track_parameters<test_algebra> track(pos, 0.f, mom, -1.f);
0136
0137   propagation::config prop_cfg{};
0138   propagator_t p{prop_cfg};
0139
0140   propagator_t::state state(track, d, prop_cfg.context);
0141
0142   EXPECT_TRUE(p.propagate(state))
0143       << state.navigation().inspector().to_string() << std::endl;
0144 }
0145
0146 /// Fixture for Runge-Kutta Propagation
0147 class PropagatorWithRkStepper : public ::testing::TestWithParam<
0148                                     std::tuple<scalar, scalar, test::vector3>> {
0149  public:
0150   using generator_t =
0151       uniform_track_generator<free_track_parameters<test_algebra>>;
0152
0153   /// Set the test environment up
0154   void SetUp() override {
0155     overstep_tol = std::get<0>(GetParam());
0156     step_constr = std::get<1>(GetParam());
0157
0158     trk_gen_cfg.phi_steps(50u).theta_steps(50u);
0159     trk_gen_cfg.p_tot(10.f * unit<scalar>::GeV);
0160   }
0161
0162   /// Clean up
0163   void TearDown() override { /* Do nothing */ }
0164
0165  protected:
0166   /// Detector configuration
0167   vecmem::host_memory_resource host_mr;
0168
0169   /// Toy detector configuration
0170   toy_det_config<scalar> toy_cfg =
0171       toy_det_config<scalar>{}.n_brl_layers(4u).n_edc_layers(7u);
0172
0173   /// Track generator configuration
0174   generator_t::configuration trk_gen_cfg{};
0175
0176   /// Stepper configuration
0177   scalar overstep_tol{std::numeric_limits<scalar>::max()};
0178   scalar step_constr{std::numeric_limits<scalar>::max()};
0179 };
0180
0181 /// Test propagation in a constant magnetic field using a Runge-Kutta stepper
0182 TEST_P(PropagatorWithRkStepper, rk4_propagator_const_bfield) {
0183   // Constant magnetic field type
0184   using bfield_t = bfield::const_field_t<scalar>;
0185
0186   // Toy detector
0187   using detector_t = detector<test::toy_metadata>;
0188
0189   // Runge-Kutta propagation
0190   using navigator_t =
0191       caching_navigator<detector_t, cache_size, navigation::print_inspector>;
0192   using track_t = free_track_parameters<test_algebra>;
0193   using constraints_t = constrained_step<scalar>;
0194   using policy_t = stepper_rk_policy<scalar>;
0195   using stepper_t =
0196       rk_stepper<bfield_t::view_t, test_algebra, constraints_t, policy_t>;
0197   // Include helix actor to check track position/covariance
0198   using actor_chain_t = actor_chain<
0199       actor::pathlimit_aborter<scalar>,
0200       actor::parameter_updater<
0201           test_algebra, actor::pointwise_material_interactor<test_algebra>,
0202           actor::helix_inspector>>;
0203   using propagator_t = propagator<stepper_t, navigator_t, actor_chain_t>;
0204
0205   // Build detector
0206   toy_cfg.use_material_maps(false);
0207   toy_cfg.mapped_material(detray::vacuum<scalar>());
0208   const auto [det, names] = build_toy_detector<test_algebra>(host_mr, toy_cfg);
0209
0210   const bfield_t bfield = create_const_field<scalar>(std::get<2>(GetParam()));
0211
0212   // Propagator is built from the stepper and navigator
0213   propagation::config cfg{};
0214   cfg.navigation.intersection.overstep_tolerance =
0215       static_cast<float>(overstep_tol);
0216   cfg.navigation.search_window = {3u, 3u};
0217   cfg.navigation.estimate_scattering_noise = false;
0218   propagator_t p{cfg};
0219
0220   // Iterate through uniformly distributed momentum directions
0221   for (auto track : generator_t{trk_gen_cfg}) {
0222     assert(track.qop() != 0.f);
0223
0224     // Generate second track state used for propagation with pathlimit
0225     track_t lim_track(track);
0226
0227     // Build actor states: the helix inspector can be shared
0228     auto actor_states = actor_chain_t::make_default_actor_states();
0229     auto actor_states_lim = actor_chain_t::make_default_actor_states();
0230
0231     // Make sure the lim state is being terminated
0232     auto& pathlimit_aborter_state =
0233         detail::get<actor::pathlimit_aborter<scalar>::state>(actor_states_lim);
0234     pathlimit_aborter_state.set_path_limit(path_limit);
0235
0236     // Init propagator states
0237     typename propagator_t::stepper_type::magnetic_field_type bfield_view(
0238         bfield);
0239     propagator_t::state state(track, bfield_view, det);
0240     propagator_t::state sync_state(track, bfield_view, det);
0241     propagator_t::state lim_state(lim_track, bfield_view, det);
0242
0243     state.debug(true);
0244     sync_state.debug(true);
0245     lim_state.debug(true);
0246
0247     // Set step constraints
0248     state.stepping().template set_constraint<step::constraint::e_accuracy>(
0249         step_constr);
0250     sync_state.stepping().template set_constraint<step::constraint::e_accuracy>(
0251         step_constr);
0252     lim_state.stepping().template set_constraint<step::constraint::e_accuracy>(
0253         step_constr);
0254
0255     // No multiple scattering simulated in this test
0256     using updater_state_t = actor::parameter_updater_state<test_algebra>;
0257     detail::get<updater_state_t>(actor_states)
0258         .noise_estimation_cfg()
0259         .estimate_scattering_noise = false;
0260     detail::get<updater_state_t>(actor_states_lim)
0261         .noise_estimation_cfg()
0262         .estimate_scattering_noise = false;
0263
0264     // Propagate the entire detector
0265     ASSERT_TRUE(
0266         p.propagate(state, actor_chain_t::setup_actor_states(actor_states)))
0267         << state.navigation().inspector().to_string() << std::endl;
0268
0269     // Propagate with path limit
0270     ASSERT_FALSE(p.propagate(
0271         lim_state, actor_chain_t::setup_actor_states(actor_states_lim)))
0272         << lim_state.navigation().inspector().to_string() << std::endl;
0273
0274     ASSERT_GE(std::abs(path_limit), lim_state.stepping().abs_path_length())
0275         << "Absolute path length: " << lim_state.stepping().abs_path_length()
0276         << ", path limit: " << path_limit << std::endl;
0277     //<< state.navigation().inspector().to_string() << std::endl;
0278   }
0279 }
0280
0281 /// Test propagation in an inhomogeneous magnetic field using a Runge-Kutta
0282 /// stepper
0283 TEST_P(PropagatorWithRkStepper, rk4_propagator_inhom_bfield) {
0284   // Magnetic field map using nearest neighbor interpolation
0285   using bfield_t = bfield::inhom_field_t<scalar>;
0286
0287   // Toy detector
0288   using detector_t = detector<test::toy_metadata>;
0289
0290   // Runge-Kutta propagation
0291   using navigator_t =
0292       caching_navigator<detector_t, cache_size, navigation::print_inspector>;
0293   using track_t = free_track_parameters<test_algebra>;
0294   using constraints_t = constrained_step<scalar>;
0295   using policy_t = stepper_rk_policy<scalar>;
0296   using stepper_t =
0297       rk_stepper<bfield_t::view_t, test_algebra, constraints_t, policy_t>;
0298   // Include helix actor to check track position/covariance
0299   using actor_chain_t = actor_chain<
0300       actor::pathlimit_aborter<scalar>,
0301       actor::parameter_updater<
0302           test_algebra, actor::pointwise_material_interactor<test_algebra>>>;
0303   using propagator_t = propagator<stepper_t, navigator_t, actor_chain_t>;
0304
0305   // Build detector and magnetic field
0306   toy_cfg.use_material_maps(false);
0307   const auto [det, names] = build_toy_detector<test_algebra>(host_mr, toy_cfg);
0308   const bfield_t bfield = create_inhom_field<scalar>();
0309
0310   // Propagator is built from the stepper and navigator
0311   propagation::config cfg{};
0312   cfg.navigation.intersection.overstep_tolerance =
0313       static_cast<float>(overstep_tol);
0314   cfg.navigation.search_window = {3u, 3u};
0315   cfg.navigation.estimate_scattering_noise = false;
0316   propagator_t p{cfg};
0317
0318   // Iterate through uniformly distributed momentum directions
0319   for (auto track : generator_t{trk_gen_cfg}) {
0320     // Generate second track state used for propagation with pathlimit
0321     track_t lim_track(track);
0322
0323     // Build actor states: the helix inspector can be shared
0324     actor::pathlimit_aborter<scalar>::state unlimted_aborter_state{};
0325     actor::pathlimit_aborter<scalar>::state pathlimit_aborter_state{path_limit};
0326     actor::parameter_updater_state<test_algebra> updater_state{cfg};
0327     actor::pointwise_material_interactor<test_algebra>::state
0328         interactor_state{};
0329
0330     // Create actor states tuples
0331     auto actor_states =
0332         detray::tie(unlimted_aborter_state, updater_state, interactor_state);
0333     auto lim_actor_states =
0334         detray::tie(pathlimit_aborter_state, updater_state, interactor_state);
0335
0336     // Init propagator states
0337     typename propagator_t::stepper_type::magnetic_field_type bfield_view(
0338         bfield);
0339     propagator_t::state state(track, bfield_view, det);
0340     propagator_t::state lim_state(lim_track, bfield_view, det);
0341
0342     // Set step constraints
0343     state.stepping().template set_constraint<step::constraint::e_accuracy>(
0344         step_constr);
0345     lim_state.stepping().template set_constraint<step::constraint::e_accuracy>(
0346         step_constr);
0347
0348     // Propagate the entire detector
0349     state.debug(true);
0350     ASSERT_TRUE(p.propagate(state, actor_states))
0351         << state.navigation().inspector().to_string() << std::endl;
0352
0353     // Propagate with path limit
0354     ASSERT_NEAR(pathlimit_aborter_state.path_limit(), path_limit, tol);
0355     lim_state.debug(true);
0356     ASSERT_FALSE(p.propagate(lim_state, lim_actor_states))
0357         << lim_state.navigation().inspector().to_string() << std::endl;
0358
0359     ASSERT_TRUE(lim_state.stepping().path_length() < std::abs(path_limit) + tol)
0360         << "path length: " << lim_state.stepping().path_length()
0361         << ", path limit: " << path_limit << std::endl;
0362     //<< state.navigation().inspector().to_string() << std::endl;
0363   }
0364 }
0365
0366 // No step size constraint
0367 INSTANTIATE_TEST_SUITE_P(
0368     detray_propagator_validation1, PropagatorWithRkStepper,
0369     ::testing::Values(std::make_tuple(-100.f * unit<scalar>::um,
0370                                       std::numeric_limits<scalar>::max(),
0371                                       vector3{0.f * unit<scalar>::T,
0372                                               0.f * unit<scalar>::T,
0373                                               2.f * unit<scalar>::T})));
0374
0375 // Add some restrictions for more frequent navigation updates in the cases of
0376 // non-z-aligned B-fields
0377 INSTANTIATE_TEST_SUITE_P(
0378     detray_propagator_validation2, PropagatorWithRkStepper,
0379     ::testing::Values(std::make_tuple(-400.f * unit<scalar>::um,
0380                                       std::numeric_limits<scalar>::max(),
0381                                       vector3{0.f * unit<scalar>::T,
0382                                               1.f * unit<scalar>::T,
0383                                               1.f * unit<scalar>::T})));
0384
0385 INSTANTIATE_TEST_SUITE_P(
0386     detray_propagator_validation3, PropagatorWithRkStepper,
0387     ::testing::Values(std::make_tuple(-400.f * unit<scalar>::um,
0388                                       std::numeric_limits<scalar>::max(),
0389                                       vector3{1.f * unit<scalar>::T,
0390                                               0.f * unit<scalar>::T,
0391                                               1.f * unit<scalar>::T})));
0392
0393 INSTANTIATE_TEST_SUITE_P(
0394     detray_propagator_validation4, PropagatorWithRkStepper,
0395     ::testing::Values(std::make_tuple(-600.f * unit<scalar>::um,
0396                                       std::numeric_limits<scalar>::max(),
0397                                       vector3{1.f * unit<scalar>::T,
0398                                               1.f * unit<scalar>::T,
0399                                               1.f * unit<scalar>::T})));