cpu/propagator/propagation_direct_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/navigation/caching_navigator.hpp"
0012 #include "detray/navigation/direct_navigator.hpp"
0013 #include "detray/propagator/actors.hpp"
0014 #include "detray/propagator/actors/parameter_updater.hpp"
0015 #include "detray/propagator/propagator.hpp"
0016 #include "detray/propagator/rk_stepper.hpp"
0017 #include "detray/tracks/tracks.hpp"
0018 #include "detray/tracks/trajectories.hpp"
0019
0020 // Detray test include(s)
0021 #include "detray/test/common/bfield.hpp"
0022 #include "detray/test/common/build_toy_detector.hpp"
0023 #include "detray/test/common/build_wire_chamber.hpp"
0024 #include "detray/test/common/track_generators.hpp"
0025 #include "detray/test/framework/types.hpp"
0026 #include "detray/test/utils/inspectors.hpp"
0027
0028 // Vecmem include(s)
0029 #include <vecmem/memory/host_memory_resource.hpp>
0030
0031 // GTest include(s)
0032 #include <gtest/gtest.h>
0033
0034 using namespace detray;
0035
0036 using test_algebra = test::algebra;
0037 using scalar = test::scalar;
0038 using vector3 = test::vector3;
0039
0040 constexpr std::size_t cache_size{navigation::default_cache_size};
0041
0042 /// Fixture for Runge-Kutta Propagation with direct navigator and toy detector
0043 class PropagatorWithRkStepperDirectNavigatorToyDetector
0044     : public ::testing::TestWithParam<std::tuple<scalar, vector3>> {};
0045
0046 /// Test propagation in a constant magnetic field using a Runge-Kutta stepper
0047 TEST_P(PropagatorWithRkStepperDirectNavigatorToyDetector, direct_navigator) {
0048   constexpr scalar tol{2e-3f};
0049
0050   // Track generator configuration
0051   using generator_t =
0052       uniform_track_generator<free_track_parameters<test_algebra>>;
0053
0054   // Toy detector
0055   using detector_t = detector<test::toy_metadata>;
0056   using surface_t = typename detector_t::surface_type;
0057
0058   // Runge-Kutta propagation
0059   using navigator_t =
0060       caching_navigator<detector_t, cache_size, navigation::print_inspector>;
0061   using direct_navigator_t = direct_navigator<detector_t>;
0062
0063   // Constant magnetic field and stepper type
0064   using bfield_t = bfield::const_field_t<scalar>;
0065   using stepper_t = rk_stepper<bfield_t::view_t, test_algebra>;
0066
0067   // Include helix actor to check track position/covariance
0068   using actor_chain_t = actor_chain<
0069       actor::parameter_updater<
0070           test_algebra, actor::pointwise_material_interactor<test_algebra>>,
0071       actor::surface_sequencer<surface_t>>;
0072
0073   using propagator_t = propagator<stepper_t, navigator_t, actor_chain_t>;
0074   using direct_propagator_t =
0075       propagator<stepper_t, direct_navigator_t, actor_chain_t>;
0076
0077   // Memory resource
0078   vecmem::host_memory_resource host_mr;
0079
0080   // Build toy detector
0081   toy_det_config<scalar> toy_cfg =
0082       toy_det_config<scalar>{}.n_brl_layers(4u).n_edc_layers(7u);
0083
0084   // TODO: Test the collection of portal material
0085   toy_cfg.use_material_maps(false);
0086   const auto [det, names] = build_toy_detector<test_algebra>(host_mr, toy_cfg);
0087
0088   // Build mangetic field
0089   const bfield_t bfield = create_const_field<scalar>(std::get<1>(GetParam()));
0090   const bfield_t::view_t bfield_view(bfield);
0091
0092   // Track generator config
0093   generator_t::configuration trk_gen_cfg{};
0094   trk_gen_cfg.phi_steps(50u).theta_steps(50u);
0095   trk_gen_cfg.p_tot(10.f * unit<scalar>::GeV);
0096
0097   // Propagation configuration
0098   propagation::config cfg{};
0099   cfg.navigation.intersection.overstep_tolerance =
0100       static_cast<float>(std::get<0>(GetParam()));
0101   cfg.navigation.search_window = {3u, 3u};
0102   cfg.navigation.estimate_scattering_noise = false;
0103   propagator_t p{cfg};
0104   direct_propagator_t direct_p{cfg};
0105
0106   // Iterate through uniformly distributed momentum directions
0107   for (auto track : generator_t{trk_gen_cfg}) {
0108     // Build actor states: the helix inspector can be shared
0109     actor::parameter_updater_state<test_algebra> updater_state{cfg};
0110     actor::parameter_updater_state<test_algebra> fw_updater_state{cfg};
0111     actor::parameter_updater_state<test_algebra> bw_updater_state{cfg};
0112
0113     actor::pointwise_material_interactor<test_algebra>::state
0114         interactor_state{};
0115     vecmem::data::vector_buffer<surface_t> seqs_buffer{
0116         100u, host_mr, vecmem::data::buffer_type::resizable};
0117     vecmem::data::vector_buffer<surface_t> seqs_forward_buffer{
0118         100u, host_mr, vecmem::data::buffer_type::resizable};
0119     vecmem::data::vector_buffer<surface_t> seqs_backward_buffer{
0120         100u, host_mr, vecmem::data::buffer_type::resizable};
0121     vecmem::copy m_copy;
0122     m_copy.setup(seqs_buffer)->wait();
0123     m_copy.setup(seqs_forward_buffer)->wait();
0124     m_copy.setup(seqs_backward_buffer)->wait();
0125
0126     vecmem::device_vector<surface_t> seqs_device(seqs_buffer);
0127     vecmem::device_vector<surface_t> seqs_forward_device(seqs_forward_buffer);
0128     vecmem::device_vector<surface_t> seqs_backward_device(seqs_backward_buffer);
0129
0130     actor::surface_sequencer<surface_t>::state sequencer_state(seqs_device);
0131     actor::surface_sequencer<surface_t>::state sequencer_forward_state(
0132         seqs_forward_device);
0133     actor::surface_sequencer<surface_t>::state sequencer_backward_state(
0134         seqs_backward_device);
0135
0136     auto actor_states =
0137         detray::tie(updater_state, interactor_state, sequencer_state);
0138
0139     propagator_t::state state(track, bfield_view, det);
0140     navigator_t::state& navigation = state.navigation();
0141
0142     // Propagate the entire detector
0143     ASSERT_TRUE(p.propagate(state, actor_states))
0144         << navigation.inspector().to_string();
0145
0146     if (!seqs_device.empty()) {
0147       auto direct_forward_actor_states = detray::tie(
0148           fw_updater_state, interactor_state, sequencer_forward_state);
0149       auto direct_backward_actor_states = detray::tie(
0150           bw_updater_state, interactor_state, sequencer_backward_state);
0151
0152       direct_propagator_t::state direct_forward_state(track, bfield_view, det,
0153                                                       seqs_buffer);
0154
0155       ASSERT_TRUE(direct_p.propagate(direct_forward_state,
0156                                      direct_forward_actor_states));
0157
0158       // Check if all surfaces in the sequence are encountered
0159       ASSERT_TRUE(direct_forward_state.navigation().finished());
0160       ASSERT_EQ(sequencer_state.sequence().size(),
0161                 sequencer_forward_state.sequence().size());
0162       for (unsigned int i = 0; i < sequencer_state.sequence().size(); i++) {
0163         ASSERT_EQ(sequencer_state.sequence().at(i),
0164                   sequencer_forward_state.sequence().at(i));
0165       }
0166
0167       const auto ptc = state.stepping().particle_hypothesis();
0168       ASSERT_EQ(
0169           ptc.pdg_num(),
0170           direct_forward_state.stepping().particle_hypothesis().pdg_num());
0171       const auto q = ptc.charge();
0172
0173       // The initial momentum should be higher than the momentum at the
0174       // last surface
0175       ASSERT_TRUE(track.p(q) > updater_state.bound_params().p(q));
0176       ASSERT_NEAR(static_cast<float>(updater_state.bound_params().p(q)),
0177                   static_cast<float>(fw_updater_state.bound_params().p(q)),
0178                   static_cast<float>(updater_state.bound_params().p(q)) * tol);
0179
0180       direct_propagator_t::state direct_backward_state(
0181           fw_updater_state.bound_params(), bfield_view, det, seqs_buffer);
0182       direct_backward_state.navigation().set_direction(
0183           detray::navigation::direction::e_backward);
0184       direct_backward_state.navigation().reset();
0185
0186       ASSERT_TRUE(direct_p.propagate(direct_backward_state,
0187                                      direct_backward_actor_states));
0188       // Check if all surfaces in the sequence are encountered
0189       ASSERT_TRUE(direct_backward_state.navigation().finished());
0190       ASSERT_EQ(sequencer_state.sequence().size(),
0191                 sequencer_backward_state.sequence().size());
0192       for (unsigned int i = 0u; i < sequencer_state.sequence().size(); i++) {
0193         unsigned int j = sequencer_state.sequence().size() - 1u - i;
0194         ASSERT_EQ(sequencer_state.sequence().at(i),
0195                   sequencer_backward_state.sequence().at(j));
0196       }
0197
0198       ASSERT_NEAR(static_cast<float>(track.p(q)),
0199                   static_cast<float>(bw_updater_state.bound_params().p(q)),
0200                   static_cast<float>(track.p(q)) * tol);
0201       ASSERT_TRUE(bw_updater_state.bound_params().p(q) >
0202                   fw_updater_state.bound_params().p(q));
0203     }
0204   }
0205 }
0206
0207 INSTANTIATE_TEST_SUITE_P(
0208     detray_propagator_direct_navigator_toy_detector,
0209     PropagatorWithRkStepperDirectNavigatorToyDetector,
0210     ::testing::Values(
0211         std::make_tuple(-100.f * unit<scalar>::um,
0212                         vector3{0.f * unit<scalar>::T, 0.f * unit<scalar>::T,
0213                                 2.f * unit<scalar>::T}),
0214         std::make_tuple(-400.f * unit<scalar>::um,
0215                         vector3{0.f * unit<scalar>::T, 1.f * unit<scalar>::T,
0216                                 1.f * unit<scalar>::T}),
0217         std::make_tuple(-400.f * unit<scalar>::um,
0218                         vector3{1.f * unit<scalar>::T, 0.f * unit<scalar>::T,
0219                                 1.f * unit<scalar>::T}),
0220         std::make_tuple(-600.f * unit<scalar>::um,
0221                         vector3{1.f * unit<scalar>::T, 1.f * unit<scalar>::T,
0222                                 1.f * unit<scalar>::T})));
0223
0224 /// Fixture for Runge-Kutta Propagation with direct navigator and wire chamber
0225 class PropagatorWithRkStepperDirectNavigatorWireChamber
0226     : public ::testing::TestWithParam<std::tuple<scalar, vector3>> {};
0227
0228 /// Test propagation in a constant magnetic field using a Runge-Kutta stepper
0229 TEST_P(PropagatorWithRkStepperDirectNavigatorWireChamber, direct_navigator) {
0230   constexpr scalar tol{2e-5f};
0231
0232   // Track generator configuration
0233   using generator_t =
0234       uniform_track_generator<free_track_parameters<test_algebra>>;
0235
0236   // Toy detector
0237   using detector_t = detector<test::wire_chamber_metadata>;
0238   using surface_t = typename detector_t::surface_type;
0239
0240   // Default navigator for comparison
0241   using navigator_t =
0242       caching_navigator<detector_t, cache_size, navigation::print_inspector>;
0243   using direct_navigator_t = direct_navigator<detector_t>;
0244
0245   // Constant magnetic field and stepper type
0246   using bfield_t = bfield::const_field_t<scalar>;
0247   using stepper_t = rk_stepper<bfield_t::view_t, test_algebra>;
0248
0249   // Include helix actor to check track position/covariance
0250   using actor_chain_t = actor_chain<
0251       actor::parameter_updater<
0252           test_algebra, actor::pointwise_material_interactor<test_algebra>>,
0253       actor::surface_sequencer<surface_t>>;
0254
0255   using propagator_t = propagator<stepper_t, navigator_t, actor_chain_t>;
0256   using direct_propagator_t =
0257       propagator<stepper_t, direct_navigator_t, actor_chain_t>;
0258
0259   // Memory resource
0260   vecmem::host_memory_resource host_mr;
0261
0262   // Build wire chamber
0263   wire_chamber_config<scalar> wire_cfg{};
0264   auto [det, names] = build_wire_chamber<test::algebra>(host_mr, wire_cfg);
0265
0266   // Build mangetic field
0267   const bfield_t bfield = create_const_field<scalar>(std::get<1>(GetParam()));
0268   const bfield_t::view_t bfield_view(bfield);
0269
0270   // Truth track generation
0271   generator_t::configuration trk_gen_cfg{};
0272   trk_gen_cfg.phi_steps(50u).theta_steps(50u);
0273   trk_gen_cfg.p_tot(10.f * unit<scalar>::GeV);
0274
0275   // Propagation configuration
0276   propagation::config cfg{};
0277   cfg.navigation.intersection.overstep_tolerance =
0278       static_cast<float>(std::get<0>(GetParam()));
0279   cfg.navigation.search_window = {5u, 5u};
0280   cfg.navigation.estimate_scattering_noise = false;
0281   propagator_t p{cfg};
0282   direct_propagator_t direct_p{cfg};
0283
0284   // Iterate through uniformly distributed momentum directions
0285   for (auto track : generator_t{trk_gen_cfg}) {
0286     // Build actor states: the helix inspector can be shared
0287     actor::parameter_updater_state<test_algebra> updater_state{cfg};
0288     actor::parameter_updater_state<test_algebra> fw_updater_state{cfg};
0289     actor::parameter_updater_state<test_algebra> bw_updater_state{cfg};
0290     actor::pointwise_material_interactor<test_algebra>::state
0291         interactor_state{};
0292
0293     vecmem::data::vector_buffer<surface_t> seqs_buffer{
0294         100u, host_mr, vecmem::data::buffer_type::resizable};
0295     vecmem::data::vector_buffer<surface_t> seqs_forward_buffer{
0296         100u, host_mr, vecmem::data::buffer_type::resizable};
0297     vecmem::data::vector_buffer<surface_t> seqs_backward_buffer{
0298         100u, host_mr, vecmem::data::buffer_type::resizable};
0299     vecmem::copy m_copy;
0300     m_copy.setup(seqs_buffer)->wait();
0301     m_copy.setup(seqs_forward_buffer)->wait();
0302     m_copy.setup(seqs_backward_buffer)->wait();
0303
0304     vecmem::device_vector<surface_t> seqs_device(seqs_buffer);
0305     vecmem::device_vector<surface_t> seqs_forward_device(seqs_forward_buffer);
0306     vecmem::device_vector<surface_t> seqs_backward_device(seqs_backward_buffer);
0307
0308     actor::surface_sequencer<surface_t>::state sequencer_state(seqs_device);
0309     actor::surface_sequencer<surface_t>::state sequencer_forward_state(
0310         seqs_forward_device);
0311     actor::surface_sequencer<surface_t>::state sequencer_backward_state(
0312         seqs_backward_device);
0313
0314     auto actor_states =
0315         detray::tie(updater_state, interactor_state, sequencer_state);
0316
0317     propagator_t::state state(track, bfield_view, det);
0318     navigator_t::state& navigation = state.navigation();
0319
0320     // Propagate the entire detector
0321     ASSERT_TRUE(p.propagate(state, actor_states))
0322         << navigation.inspector().to_string();
0323
0324     if (!seqs_device.empty()) {
0325       auto direct_forward_actor_states = detray::tie(
0326           fw_updater_state, interactor_state, sequencer_forward_state);
0327       auto direct_backward_actor_states = detray::tie(
0328           bw_updater_state, interactor_state, sequencer_backward_state);
0329
0330       direct_propagator_t::state direct_forward_state(track, bfield_view, det,
0331                                                       seqs_buffer);
0332
0333       ASSERT_TRUE(direct_p.propagate(direct_forward_state,
0334                                      direct_forward_actor_states));
0335
0336       // Check if all surfaces in the sequence are encountered
0337       ASSERT_TRUE(direct_forward_state.navigation().finished());
0338       ASSERT_EQ(sequencer_state.sequence().size(),
0339                 sequencer_forward_state.sequence().size());
0340       for (unsigned int i = 0; i < sequencer_state.sequence().size(); i++) {
0341         ASSERT_EQ(sequencer_state.sequence().at(i),
0342                   sequencer_forward_state.sequence().at(i));
0343       }
0344
0345       const auto ptc = state.stepping().particle_hypothesis();
0346       ASSERT_EQ(
0347           ptc.pdg_num(),
0348           direct_forward_state.stepping().particle_hypothesis().pdg_num());
0349       const auto q = ptc.charge();
0350
0351       // The initial momentum should be higher than or equal to the
0352       // momentum at the last surface
0353       ASSERT_GE(track.p(q), updater_state.bound_params().p(q));
0354       ASSERT_NEAR(static_cast<float>(updater_state.bound_params().p(q)),
0355                   static_cast<float>(fw_updater_state.bound_params().p(q)),
0356                   static_cast<float>(updater_state.bound_params().p(q)) * tol);
0357
0358       direct_propagator_t::state direct_backward_state(
0359           fw_updater_state.bound_params(), bfield_view, det, seqs_buffer);
0360       direct_backward_state.navigation().set_direction(
0361           detray::navigation::direction::e_backward);
0362       direct_backward_state.navigation().reset();
0363
0364       ASSERT_TRUE(direct_p.propagate(direct_backward_state,
0365                                      direct_backward_actor_states));
0366
0367       // Check if all surfaces in the sequence are encountered
0368       ASSERT_TRUE(direct_backward_state.navigation().finished());
0369
0370       ASSERT_EQ(sequencer_state.sequence().size(),
0371                 sequencer_backward_state.sequence().size());
0372       for (unsigned int i = 0u; i < sequencer_state.sequence().size(); i++) {
0373         unsigned int j = sequencer_state.sequence().size() - 1u - i;
0374         EXPECT_EQ(sequencer_state.sequence().at(i),
0375                   sequencer_backward_state.sequence().at(j));
0376       }
0377
0378       ASSERT_NEAR(static_cast<float>(track.p(q)),
0379                   static_cast<float>(bw_updater_state.bound_params().p(q)),
0380                   static_cast<float>(track.p(q)) * tol);
0381       ASSERT_GE(bw_updater_state.bound_params().p(q),
0382                 fw_updater_state.bound_params().p(q));
0383     }
0384   }
0385 }
0386
0387 INSTANTIATE_TEST_SUITE_P(
0388     detray_propagator_direct_navigator_wire_chamber,
0389     PropagatorWithRkStepperDirectNavigatorWireChamber,
0390     ::testing::Values(
0391         std::make_tuple(-100.f * unit<scalar>::um,
0392                         vector3{0.f * unit<scalar>::T, 0.f * unit<scalar>::T,
0393                                 2.f * unit<scalar>::T}),
0394         std::make_tuple(-400.f * unit<scalar>::um,
0395                         vector3{0.f * unit<scalar>::T, 1.f * unit<scalar>::T,
0396                                 1.f * unit<scalar>::T}),
0397         std::make_tuple(-400.f * unit<scalar>::um,
0398                         vector3{1.f * unit<scalar>::T, 0.f * unit<scalar>::T,
0399                                 1.f * unit<scalar>::T}),
0400         std::make_tuple(-600.f * unit<scalar>::um,
0401                         vector3{1.f * unit<scalar>::T, 1.f * unit<scalar>::T,
0402                                 1.f * unit<scalar>::T})));