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/intersection/helix_intersector.hpp"
0016 #include "detray/propagator/detail/jacobian_engine.hpp"
0017 #include "detray/test/framework/types.hpp"
0018 #include "detray/tracks/tracks.hpp"
0019 #include "detray/utils/axis_rotation.hpp"
0020
0021 // Vecmem include(s)
0022 #include <vecmem/memory/host_memory_resource.hpp>
0023
0024 // google-test include(s).
0025 #include <gtest/gtest.h>
0026
0027 using namespace detray;
0028
0029 // Algebra types
0030 using test_algebra = test::algebra;
0031 using scalar = test::scalar;
0032 using transform3 = test::transform3;
0033 using vector3 = test::vector3;
0034 using intersection_t = intersection2D<surface_descriptor<>, test_algebra,
0035                                       intersection::contains_pos>;
0036
0037 // Mask types to be tested
0038 // @TODO: Remove unbounded tag
0039 using annulus_type =
0040     detray::mask<detray::unbounded<detray::annulus2D>, test_algebra>;
0041 using rectangle_type = detray::mask<detray::rectangle2D, test_algebra>;
0042 using trapezoid_type = detray::mask<detray::trapezoid2D, test_algebra>;
0043 using ring_type = detray::mask<detray::ring2D, test_algebra>;
0044 using cylinder_type = detray::mask<detray::cylinder2D, test_algebra>;
0045 using straw_tube_type = detray::mask<detray::line_circular, test_algebra>;
0046 using drift_cell_type = detray::mask<detray::line_square, test_algebra>;
0047
0048 // Test class for covariance transport
0049 template <typename T>
0050 class detray_propagation_HelixCovarianceTransportValidation
0051     : public ::testing::Test {
0052  public:
0053   // Environment Setup
0054   const vector3 B{static_cast<scalar>(0.f), static_cast<scalar>(0.f),
0055                   1.f * unit<scalar>::T};
0056   const vector3 y_axis{0.f, 1.f, 0.f};
0057   const vector3 z_axis{0.f, 0.f, 1.f};
0058   static constexpr const scalar mask_tolerance{1e-3f};
0059   static constexpr const scalar tolerance{3e-2f};
0060
0061   // Test types
0062   using mask_type = T;
0063   using local_frame_type = typename mask_type::local_frame;
0064
0065   // First mask at the origin is always rectangle
0066   using first_mask_type = rectangle_type;
0067   using first_local_frame_type = typename first_mask_type::local_frame;
0068
0069   // Algebra type
0070   using algebra_type = typename local_frame_type::algebra_type;
0071
0072   // Vector and matrix types
0073   using bound_param_vector_t =
0074       typename bound_track_parameters<algebra_type>::parameter_vector_type;
0075   using bound_matrix_t = bound_matrix<algebra_type>;
0076   using bound_to_free_matrix_t = bound_to_free_matrix<algebra_type>;
0077
0078   using free_matrix_t = free_matrix<algebra_type>;
0079   using free_to_bound_matrix_t = free_to_bound_matrix<algebra_type>;
0080
0081   std::tuple<annulus_type, rectangle_type, trapezoid_type, ring_type,
0082              cylinder_type, straw_tube_type, drift_cell_type>
0083       masks = std::make_tuple(
0084           annulus_type{0u, 7.2f * unit<scalar>::mm, 12.0f * unit<scalar>::mm,
0085                        0.74195f, 1.33970f, 0.f, -2.f, 2.f},
0086           rectangle_type{0u, 50 * unit<scalar>::mm, 50 * unit<scalar>::mm},
0087           trapezoid_type{0u, 50 * unit<scalar>::mm, 100 * unit<scalar>::mm,
0088                          50 * unit<scalar>::mm,
0089                          1.f / (2.f * 50 * unit<scalar>::mm)},
0090           ring_type{0u, 0.f, 50.f * unit<scalar>::mm},
0091           cylinder_type{0u, 0.3f * unit<scalar>::mm, -100.f * unit<scalar>::mm,
0092                         100.f * unit<scalar>::mm},
0093           straw_tube_type{0u, 50.f * unit<scalar>::mm,
0094                           100.f * unit<scalar>::mm},
0095           drift_cell_type{0u, 50.f * unit<scalar>::mm,
0096                           100.f * unit<scalar>::mm});
0097
0098   // Create transform matrices
0099   std::vector<transform3> create_transforms(
0100       detail::helix<algebra_type>& reference_helix,
0101       const std::size_t n_planes) {
0102     std::vector<transform3> trfs;
0103
0104     // Step size between two neighbor planes
0105     const scalar S{2.f * constant<scalar>::pi /
0106                    std::abs(reference_helix.qop() * reference_helix.B())};
0107     const scalar step_size = S / static_cast<scalar>(n_planes);
0108
0109     for (std::size_t i = 0u; i < n_planes; i++) {
0110       const scalar s = step_size * static_cast<scalar>(i);
0111
0112       // Translation of the new surface
0113       vector3 trl = reference_helix(s);
0114
0115       // Normal vector of the new surface
0116       vector3 w = reference_helix.dir(s);
0117       vector3 v = vector::cross(z_axis, w);
0118
0119       if (i > 0u &&
0120           (std::is_same_v<local_frame_type, cylindrical2D<algebra_type>> ||
0121            std::is_same_v<local_frame_type, line2D<algebra_type>>)) {
0122         const vector3 r_axis = vector::cross(w, z_axis);
0123
0124         // Rotate for cylinder and wire
0125         axis_rotation<algebra_type> axis_rot(r_axis,
0126                                              constant<scalar>::pi / 2.f);
0127
0128         // Test masks are rotated
0129         w = axis_rot(w);
0130         EXPECT_NEAR(vector::norm(r_axis), 1.f, tolerance);
0131
0132         v = vector::cross(r_axis, w);
0133
0134         /*
0135         const vector3 offset{0.f, 0.f, 10.f * unit<scalar>::mm};
0136         trl = trl + offset;
0137         */
0138       } else {
0139         // @TODO why does this offset (in y-direction) fail the test???
0140         // const vector3 offset{0.f, 10.f * unit<scalar>::mm, 10.f *
0141         // unit<scalar>::mm};
0142         const vector3 offset{static_cast<scalar>(0.f), static_cast<scalar>(0.f),
0143                              10.f * unit<scalar>::mm};
0144         trl = trl + offset;
0145       }
0146
0147       // Add transform matrix
0148       trfs.emplace_back(trl, w, v);
0149     }
0150
0151     return trfs;
0152   }
0153
0154   // Error propagation
0155   template <typename departure_mask_type, typename destination_mask_type>
0156   bound_track_parameters<algebra_type> propagate(
0157       const bound_track_parameters<algebra_type>& bound_params,
0158       const transform3& trf_0, const transform3& trf_1,
0159       const departure_mask_type& mask_0, const destination_mask_type& mask_1,
0160       scalar& total_path_length, std::vector<intersection_t>& sfis) {
0161     using departure_frame = typename departure_mask_type::local_frame;
0162     using destination_frame = typename destination_mask_type::local_frame;
0163
0164     using jacobian_engine = detail::jacobian_engine<algebra_type>;
0165
0166     const bound_param_vector_t& bound_vec_0 = bound_params;
0167     const bound_matrix_t& bound_cov_0 = bound_params.covariance();
0168
0169     // Free vector at the departure surface
0170     const auto free_trk_0 =
0171         detail::bound_to_free_vector(trf_0, mask_0, bound_vec_0);
0172
0173     // Helix from the departure surface
0174     detail::helix<algebra_type> hlx(free_trk_0, B);
0175
0176     // Bound-to-free jacobian at the departure surface
0177     const bound_to_free_matrix_t bound_to_free_jacobi =
0178         jacobian_engine::template bound_to_free_jacobian<departure_frame>(
0179             trf_0, mask_0, bound_vec_0);
0180
0181     // Get the intersection on the next surface
0182     helix_intersector<typename destination_mask_type::shape, algebra_type>
0183         helix_inters{};
0184     helix_inters.run_rtsafe = false;
0185     const intersection_t is = get_intersection(helix_inters(
0186         hlx, surface_descriptor<>{}, mask_1, trf_1, this->mask_tolerance));
0187     // Check for successful intersection
0188     EXPECT_TRUE(is.is_inside()) << is;
0189
0190     sfis.push_back(is);
0191
0192     // Helical path length between two surfaces
0193     const auto path_length = is.path();
0194
0195     // Add the path length to the total path length
0196     total_path_length += path_length;
0197
0198     // Free transport jacobian between two surfaces
0199     const free_matrix_t transport_jacobi = hlx.jacobian(path_length);
0200
0201     // r at the destination surface
0202     const vector3 r = hlx.pos(path_length);
0203
0204     // dr/ds, or t at the destination surface
0205     const vector3 t = hlx.dir(path_length);
0206
0207     // d^2r/ds^2, or dt/ds at the destination surface
0208     const vector3 dtds = hlx.qop() * vector::cross(t, B);
0209
0210     // d(qop)/ds, which is zero in this test
0211     const scalar dqopds = 0.f;
0212
0213     // Free track at the destination surface
0214     free_track_parameters<algebra_type> free_trk_1;
0215     free_trk_1.set_pos(r);
0216     free_trk_1.set_dir(t);
0217     free_trk_1.set_qop(free_trk_0.qop());
0218
0219     // Path correction
0220     const free_matrix_t path_correction =
0221         jacobian_engine::template path_correction<destination_frame>(
0222             r, t, dtds, dqopds, trf_1);
0223
0224     // Correction term for the path variation
0225     const free_matrix_t correction_term =
0226         matrix::identity<free_matrix_t>() + path_correction;
0227
0228     // Free-to-bound jacobian at the destination surface
0229     const free_to_bound_matrix_t free_to_bound_jacobi =
0230         jacobian_engine::template free_to_bound_jacobian<destination_frame>(
0231             trf_1, free_trk_1);
0232
0233     // Bound vector at the destination surface
0234     const bound_param_vector_t bound_vec_1 =
0235         detail::free_to_bound_vector<destination_frame>(trf_1, free_trk_1);
0236
0237     // Full jacobian
0238     const bound_matrix_t full_jacobi = free_to_bound_jacobi * correction_term *
0239                                        transport_jacobi * bound_to_free_jacobi;
0240
0241     // Update the covariance at the destination surface
0242     const bound_matrix_t bound_cov_1 =
0243         full_jacobi * bound_cov_0 * matrix::transpose(full_jacobi);
0244
0245     bound_track_parameters<algebra_type> ret;
0246     ret.set_parameter_vector(bound_vec_1);
0247     ret.set_covariance(bound_cov_1);
0248
0249     return ret;
0250   }
0251
0252   intersection_t get_intersection(
0253       const std::array<intersection_t, 2>& inters) const {
0254     return inters[0];
0255   }
0256
0257   intersection_t get_intersection(const intersection_t& inters) const {
0258     return inters;
0259   }
0260 };
0261
0262 using TestTypes =
0263     ::testing::Types<annulus_type, rectangle_type, trapezoid_type, ring_type,
0264                      cylinder_type, straw_tube_type, drift_cell_type>;
0265 TYPED_TEST_SUITE(detray_propagation_HelixCovarianceTransportValidation,
0266                  TestTypes, );
0267
0268 TYPED_TEST(detray_propagation_HelixCovarianceTransportValidation,
0269            one_loop_test) {
0270   // @NOTE: The test with high energy (>1 GeV) might fail due
0271   // to the numerical instability
0272   free_track_parameters<test_algebra> free_trk(
0273       {static_cast<scalar>(0.f), static_cast<scalar>(0.f),
0274        static_cast<scalar>(0.f)},
0275       static_cast<scalar>(0.f),
0276       {0.1f * unit<scalar>::GeV, static_cast<scalar>(0.f),
0277        static_cast<scalar>(0.f)},
0278       static_cast<scalar>(-1.f));
0279
0280   detail::helix<test_algebra> reference_helix(free_trk, this->B);
0281
0282   const std::size_t n_planes = 10u;
0283   std::vector<transform3> trfs =
0284       this->create_transforms(reference_helix, n_planes);
0285   ASSERT_EQ(trfs.size(), 10u);
0286
0287   // Set the initial bound vector
0288   bound_parameters_vector<test_algebra> bound_vec_0 =
0289       detail::free_to_bound_vector<
0290           typename TestFixture::first_local_frame_type>(trfs[0], free_trk);
0291
0292   // Set the initial bound covariance
0293   auto bound_cov_0 = matrix::zero<
0294       typename bound_track_parameters<test_algebra>::covariance_type>();
0295   getter::element(bound_cov_0, e_bound_loc0, e_bound_loc0) = 1.f;
0296   getter::element(bound_cov_0, e_bound_loc1, e_bound_loc1) = 1.f;
0297   getter::element(bound_cov_0, e_bound_phi, e_bound_phi) = 1.f;
0298   // Set theta error to zero, to suppress the loc1 (z) divergence
0299   getter::element(bound_cov_0, e_bound_theta, e_bound_theta) = 0.f;
0300   getter::element(bound_cov_0, e_bound_qoverp, e_bound_qoverp) = 1.f;
0301   getter::element(bound_cov_0, e_bound_time, e_bound_time) = 0.f;
0302
0303   // Set bound track parameters
0304   bound_track_parameters<test_algebra> bound_params;
0305   bound_params.set_parameter_vector(bound_vec_0);
0306   bound_params.set_covariance(bound_cov_0);
0307
0308   // Create masks
0309   const auto first_mask = std::get<rectangle_type>(this->masks);
0310   const auto test_mask = std::get<typename TestFixture::mask_type>(this->masks);
0311
0312   // Total path length, just for testing purpose
0313   scalar total_path_length = 0.f;
0314
0315   // Intersections for testing purporse
0316   std::vector<intersection_t> sfis;
0317
0318   // Iterate over the planes until we reach the first plane (one loop)
0319   for (std::size_t i_p = 0u; i_p < n_planes; i_p++) {
0320     if (i_p == 0) {
0321       bound_params =
0322           this->propagate(bound_params, trfs[i_p], trfs[i_p + 1], first_mask,
0323                           test_mask, total_path_length, sfis);
0324     } else if (i_p > 0 && i_p < n_planes - 1) {
0325       bound_params =
0326           this->propagate(bound_params, trfs[i_p], trfs[i_p + 1], test_mask,
0327                           test_mask, total_path_length, sfis);
0328
0329     } else if (i_p == n_planes - 1) {
0330       bound_params =
0331           this->propagate(bound_params, trfs[i_p], trfs[0], test_mask,
0332                           first_mask, total_path_length, sfis);
0333     }
0334   }
0335
0336   // Check if the total path length is the expected value
0337   ASSERT_TRUE(total_path_length > 1e-3);
0338   ASSERT_NEAR(total_path_length,
0339               2.f * constant<scalar>::pi /
0340                   std::abs(reference_helix.qop() * reference_helix.B()),
0341               this->tolerance);
0342
0343   ASSERT_EQ(sfis.size(), n_planes);
0344   for (std::size_t i = 0u; i < n_planes; i++) {
0345     EXPECT_TRUE(sfis[i].is_inside());
0346     EXPECT_TRUE(sfis[i].is_along());
0347   }
0348
0349   // Check if the same vector is obtained after one loop
0350   for (unsigned int i = 0u; i < e_bound_size; i++) {
0351     EXPECT_NEAR(bound_vec_0[i], bound_params[i], this->tolerance)
0352         << "i: " << i << "\n"
0353         << std::setprecision(8) << bound_params;
0354   }
0355
0356   // Check if the same covariance is obtained after one loop
0357   for (unsigned int i = 0u; i < e_bound_size; i++) {
0358     for (unsigned int j = 0u; j < e_bound_size; j++) {
0359       EXPECT_NEAR(getter::element(bound_cov_0, i, j),
0360                   getter::element(bound_params.covariance(), i, j),
0361                   this->tolerance)
0362           << "i: " << i << "\nj: " << j << "\n"
0363           << std::setprecision(8) << bound_params;
0364     }
0365   }
0366 }