cpu/propagator/backward_propagation.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/identifier.hpp"
0012 #include "detray/geometry/shapes/rectangle2D.hpp"
0013 #include "detray/navigation/caching_navigator.hpp"
0014 #include "detray/propagator/actors.hpp"
0015 #include "detray/propagator/propagator.hpp"
0016 #include "detray/propagator/rk_stepper.hpp"
0017 #include "detray/tracks/ray.hpp"
0018 #include "detray/tracks/tracks.hpp"
0019
0020 // Detray test include(s)
0021 #include "detray/test/common/bfield.hpp"
0022 #include "detray/test/common/build_telescope_detector.hpp"
0023 #include "detray/test/framework/types.hpp"
0024
0025 // Vecmem include(s)
0026 #include <vecmem/memory/host_memory_resource.hpp>
0027
0028 // google-test include(s).
0029 #include <gtest/gtest.h>
0030
0031 using namespace detray;
0032
0033 // Algebra types
0034 using test_algebra = test::algebra;
0035 using scalar = test::scalar;
0036 using point2 = test::point2;
0037 using vector3 = test::vector3;
0038
0039 // Test class for the backward propagation
0040 // Input tuple: < std::vector<plane positions>, tolerance >
0041 class BackwardPropagation
0042     : public ::testing::TestWithParam<std::tuple<std::vector<scalar>, scalar>> {
0043 };
0044
0045 TEST_P(BackwardPropagation, backward_propagation) {
0046   const scalar tol = std::get<1>(GetParam());
0047
0048   vecmem::host_memory_resource host_mr;
0049
0050   // Build in x-direction from given module positions
0051   detail::ray<test_algebra> traj{{0.f, 0.f, 0.f}, 0.f, {1.f, 0.f, 0.f}, -1.f};
0052   std::vector<scalar> positions = std::get<0>(GetParam());
0053
0054   tel_det_config<test_algebra, rectangle2D> tel_cfg{200.f * unit<scalar>::mm,
0055                                                     200.f * unit<scalar>::mm};
0056   tel_cfg.positions(positions).pilot_track(traj).mat_thickness(
0057       10.f * unit<scalar>::mm);
0058
0059   // Build telescope detector with rectangular planes
0060   const auto [det, names] = build_telescope_detector(host_mr, tel_cfg);
0061
0062   // Create b field
0063   using bfield_t = bfield::const_field_t<scalar>;
0064   vector3 B{0.f * unit<scalar>::T, 0.f * unit<scalar>::T,
0065             1.f * unit<scalar>::T};
0066   const bfield_t hom_bfield = create_const_field<scalar>(B);
0067   const bfield_t::view_t bfield_view(hom_bfield);
0068
0069   using navigator_t = caching_navigator<decltype(det)>;
0070   using rk_stepper_t = rk_stepper<bfield_t::view_t, test_algebra>;
0071   using actor_chain_t = actor_chain<actor::parameter_updater<
0072       test_algebra, actor::pointwise_material_interactor<test_algebra>>>;
0073   using propagator_t = propagator<rk_stepper_t, navigator_t, actor_chain_t>;
0074
0075   // Particle hypothesis
0076   pdg_particle<scalar> ptc = muon<scalar>();
0077
0078   // Bound vector
0079   bound_parameters_vector<test_algebra> bound_vector{};
0080   bound_vector.set_theta(constant<scalar>::pi_2);
0081   bound_vector.set_qop(ptc.charge() / (1.f * unit<scalar>::GeV));
0082
0083   // Bound covariance
0084   auto bound_cov = matrix::identity<
0085       typename bound_track_parameters<test_algebra>::covariance_type>();
0086
0087   // Bound track parameter
0088   const bound_track_parameters<test_algebra> bound_param0(
0089       det.surface(0u).identifier(), bound_vector, bound_cov);
0090
0091   propagation::config prop_cfg{};
0092   prop_cfg.stepping.rk_error_tol = 1e-7f * unit<float>::mm;
0093   prop_cfg.navigation.intersection.overstep_tolerance =
0094       -100.f * unit<float>::um;
0095   prop_cfg.navigation.estimate_scattering_noise = false;
0096   propagator_t p{prop_cfg};
0097
0098   // Actors
0099   actor::parameter_updater_state<test_algebra> updater_state{prop_cfg,
0100                                                              bound_param0};
0101   actor::pointwise_material_interactor<test_algebra>::state interactor{};
0102
0103   // Forward state
0104   propagator_t::state fw_state(bound_param0, bfield_view, det,
0105                                prop_cfg.context);
0106   fw_state.set_particle(ptc);
0107   fw_state.debug(true);
0108
0109   // Run propagator
0110   p.propagate(fw_state, detray::tie(interactor, updater_state));
0111
0112   // Bound state after propagation
0113   const auto bound_param1 = updater_state.bound_params();
0114
0115   // Check if the track reaches the final surface
0116   EXPECT_EQ(bound_param0.surface_link().volume(), 0u);
0117   EXPECT_EQ(bound_param0.surface_link().index(), 0u);
0118   EXPECT_EQ(bound_param1.surface_link().volume(), 0u);
0119   EXPECT_EQ(bound_param1.surface_link().id(), surface_id::e_sensitive);
0120   EXPECT_EQ(bound_param1.surface_link().index(), positions.size() - 1u);
0121
0122   // Backward state
0123   propagator_t::state bw_state(bound_param1, bfield_view, det,
0124                                prop_cfg.context);
0125   bw_state.set_particle(ptc);
0126   bw_state.debug(true);
0127   bw_state.navigation().set_direction(navigation::direction::e_backward);
0128   updater_state =
0129       actor::parameter_updater_state<test_algebra>{prop_cfg, bound_param1};
0130
0131   // Run propagator
0132   p.propagate(bw_state, detray::tie(interactor, updater_state));
0133
0134   // Bound state after propagation
0135   const auto& bound_param2 = updater_state.bound_params();
0136
0137   // Check if the track reaches the initial surface
0138   EXPECT_EQ(bound_param2.surface_link().volume(), 0u);
0139   EXPECT_EQ(bound_param2.surface_link().id(), surface_id::e_sensitive);
0140   EXPECT_EQ(bound_param2.surface_link().index(), 0u);
0141
0142   const auto bound_vec0 = bound_param0.vector();
0143   const auto bound_vec2 = bound_param2.vector();
0144
0145   // Check vector
0146   for (unsigned int i = 0u; i < e_bound_size; i++) {
0147     EXPECT_NEAR(getter::element(bound_vec0, i, 0),
0148                 getter::element(bound_vec2, i, 0), tol)
0149         << "index: " << i;
0150   }
0151
0152   const auto bound_cov0 = bound_param0.covariance();
0153   const auto bound_cov1 = bound_param1.covariance();
0154
0155   // Some sanity checks
0156   EXPECT_GT(bound_param0.p(ptc.charge()), bound_param1.p(ptc.charge()));
0157   EXPECT_GT(bound_param2.p(ptc.charge()), bound_param1.p(ptc.charge()));
0158
0159   EXPECT_GT(getter::element(bound_cov1, e_bound_qoverp, e_bound_qoverp),
0160             getter::element(bound_cov0, e_bound_qoverp, e_bound_qoverp));
0161   EXPECT_GT(getter::element(bound_cov1, e_bound_theta, e_bound_theta),
0162             getter::element(bound_cov0, e_bound_theta, e_bound_theta));
0163   EXPECT_GT(getter::element(bound_cov1, e_bound_phi, e_bound_phi),
0164             getter::element(bound_cov0, e_bound_phi, e_bound_phi));
0165 }
0166
0167 INSTANTIATE_TEST_SUITE_P(
0168     telescope, BackwardPropagation,
0169     ::testing::Values(
0170         std::make_tuple(std::vector<scalar>{0.f}, 1e-5f),
0171         std::make_tuple(std::vector<scalar>{0.f, 10.f}, 1e-3f),
0172         std::make_tuple(std::vector<scalar>{0.f, 10.f, 20.f, 30.f, 40.f, 50.f,
0173                                             60.f, 70.f, 80.f, 90.f, 100.f},
0174                         1e-2f)));