cpu/propagator/covariance_transport.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/mask.hpp"
0012 #include "detray/geometry/shapes.hpp"
0013 #include "detray/geometry/shapes/unbounded.hpp"
0014 #include "detray/geometry/surface_descriptor.hpp"
0015 #include "detray/navigation/caching_navigator.hpp"
0016 #include "detray/propagator/actor_chain.hpp"
0017 #include "detray/propagator/actors/parameter_updater.hpp"
0018 #include "detray/propagator/concepts.hpp"
0019 #include "detray/propagator/line_stepper.hpp"
0020 #include "detray/propagator/propagator.hpp"
0021 #include "detray/tracks/ray.hpp"
0022 #include "detray/tracks/tracks.hpp"
0023 #include "detray/utils/axis_rotation.hpp"
0024
0025 // Detray test include(s)
0026 #include "detray/test/common/build_telescope_detector.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 // google-test include(s).
0034 #include <gtest/gtest.h>
0035
0036 using namespace detray;
0037
0038 // Algebra types
0039 using test_algebra = test::algebra;
0040 using scalar = test::scalar;
0041 using point2 = test::point2;
0042 using vector3 = test::vector3;
0043
0044 // Mask types to be tested
0045 // @TODO: Remove unbounded tag
0046 using annulus_type =
0047     detray::mask<detray::unbounded<detray::annulus2D>, test_algebra>;
0048 using rectangle_type = detray::mask<detray::rectangle2D, test_algebra>;
0049 using trapezoid_type = detray::mask<detray::trapezoid2D, test_algebra>;
0050 using ring_type = detray::mask<detray::ring2D, test_algebra>;
0051 using cylinder_type = detray::mask<detray::cylinder2D, test_algebra>;
0052 using straw_tube_type = detray::mask<detray::line_circular, test_algebra>;
0053 using drift_cell_type = detray::mask<detray::line_square, test_algebra>;
0054
0055 constexpr scalar tol{1e-6f};
0056 constexpr std::size_t cache_size{navigation::default_cache_size};
0057
0058 GTEST_TEST(detray_propagator, covariance_transport) {
0059   static_assert(
0060       detray::concepts::actor<actor::parameter_transporter<test_algebra>>);
0061   static_assert(detray::concepts::actor<actor::parameter_setter<test_algebra>>);
0062
0063   vecmem::host_memory_resource host_mr;
0064
0065   // Build in x-direction from given module positions
0066   detail::ray<test_algebra> traj{{0.f, 0.f, 0.f}, 0.f, {1.f, 0.f, 0.f}, -1.f};
0067   std::vector<scalar> positions = {0.f, 10.f, 20.f, 30.f, 40.f, 50.f, 60.f};
0068
0069   tel_det_config<test_algebra, rectangle2D> tel_cfg{200.f * unit<scalar>::mm,
0070                                                     200.f * unit<scalar>::mm};
0071   tel_cfg.positions(positions).pilot_track(traj);
0072
0073   // Build telescope detector with unbounded planes
0074   const auto [det, names] =
0075       build_telescope_detector<test_algebra>(host_mr, tel_cfg);
0076
0077   using navigator_t =
0078       caching_navigator<decltype(det), cache_size, navigation::print_inspector>;
0079   using cline_stepper_t = line_stepper<test_algebra>;
0080   using actor_chain_t = actor_chain<actor::parameter_updater<test_algebra>>;
0081   using propagator_t = propagator<cline_stepper_t, navigator_t, actor_chain_t>;
0082
0083   // Bound vector
0084   bound_parameters_vector<test_algebra> bound_vector{};
0085   bound_vector.set_theta(constant<scalar>::pi_4);
0086   bound_vector.set_qop(-0.1f);
0087
0088   // Bound covariance
0089   auto bound_cov = matrix::identity<
0090       typename bound_track_parameters<test_algebra>::covariance_type>();
0091
0092   // Note: Set angle error as ZERO, to constrain the loc0 and loc1 divergence
0093   getter::element(bound_cov, e_bound_phi, e_bound_phi) = 0.f;
0094   getter::element(bound_cov, e_bound_theta, e_bound_theta) = 0.f;
0095
0096   // Bound track parameter
0097   const bound_track_parameters<test_algebra> bound_param0(
0098       det.surface(0u).identifier(), bound_vector, bound_cov);
0099
0100   propagation::config prop_cfg{};
0101   prop_cfg.navigation.intersection.overstep_tolerance =
0102       -100.f * unit<float>::um;
0103   propagator_t p{prop_cfg};
0104   propagator_t::state propagation(bound_param0, det, prop_cfg.context);
0105
0106   // Run propagator
0107   actor::parameter_updater_state<test_algebra> updater_state{prop_cfg,
0108                                                              bound_param0};
0109   updater_state.always_update();
0110
0111   EXPECT_TRUE(p.propagate(propagation, detray::tie(updater_state)))
0112       << propagation.navigation().inspector().to_string();
0113
0114   // Bound state after one turn propagation
0115   const auto& bound_param1 = updater_state.bound_params();
0116
0117   // Check if the track reaches the final surface
0118   EXPECT_EQ(bound_param0.surface_link().volume(), 0u);
0119   EXPECT_EQ(bound_param0.surface_link().index(), 0u);
0120   EXPECT_EQ(bound_param1.surface_link().volume(), 0u);
0121   EXPECT_EQ(bound_param1.surface_link().id(), surface_id::e_sensitive);
0122   EXPECT_EQ(bound_param1.surface_link().index(), 6u);
0123
0124   const auto bound_cov0 = bound_param0.covariance();
0125   const auto bound_cov1 = bound_param1.covariance();
0126
0127   // Check covaraince
0128   for (unsigned int i = 0u; i < e_bound_size; i++) {
0129     for (unsigned int j = 0u; j < e_bound_size; j++) {
0130       EXPECT_NEAR(getter::element(bound_cov0, i, j),
0131                   getter::element(bound_cov1, i, j), tol);
0132     }
0133   }
0134 }