cpu/propagator/jacobian_polar.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/propagator/detail/jacobian_polar.hpp"
0011
0012 #include "detray/geometry/mask.hpp"
0013 #include "detray/geometry/shapes/ring2D.hpp"
0014 #include "detray/propagator/detail/jacobian_engine.hpp"
0015 #include "detray/tracks/tracks.hpp"
0016
0017 // Detray test include(s)
0018 #include "detray/test/framework/types.hpp"
0019
0020 // GTest include(s).
0021 #include <gtest/gtest.h>
0022
0023 using namespace detray;
0024
0025 using test_algebra = test::algebra;
0026 using scalar = test::scalar;
0027 using point3 = test::point3;
0028 using vector3 = test::vector3;
0029 using transform3 = test::transform3;
0030
0031 const scalar isclose{1e-5f};
0032
0033 // This test polar2D coordinate
0034 GTEST_TEST(detray_propagator, jacobian_polar2D) {
0035   using jac_engine = detail::jacobian_engine<test_algebra>;
0036   using frame_type = polar2D<test_algebra>;
0037
0038   // Preparation work
0039   const vector3 z = {0.f, 0.f, 1.f};
0040   const vector3 x = {1.f, 0.f, 0.f};
0041   const point3 t = {2.f, 3.f, 4.f};
0042   const transform3 trf(t, z, x);
0043   const point3 global1 = {4.f, 7.f, 4.f};
0044   const vector3 mom = {1.f, 2.f, 3.f};
0045   const scalar time{0.1f};
0046   const scalar charge{static_cast<scalar>(-1.)};
0047
0048   const scalar r{2.f};
0049   mask<ring2D, test_algebra> rng{0u, 0.f, r};
0050
0051   // Free track parameter
0052   const free_track_parameters<test_algebra> free_params(global1, time, mom,
0053                                                         charge);
0054   const auto bound_vec =
0055       detail::free_to_bound_vector<polar2D<test_algebra>>(trf, free_params);
0056   const auto free_params2 = detail::bound_to_free_vector(trf, rng, bound_vec);
0057
0058   // Check if the bound vector is correct
0059   ASSERT_NEAR(getter::element(bound_vec.bound_local(), e_bound_loc0),
0060               std::sqrt(20.f), isclose);
0061   ASSERT_NEAR(getter::element(bound_vec.bound_local(), e_bound_loc1),
0062               std::atan2(4.f, 2.f), isclose);
0063   ASSERT_NEAR(bound_vec.phi(), 1.1071487f, isclose);     // atan(2)
0064   ASSERT_NEAR(bound_vec.theta(), 0.64052231f, isclose);  // atan(sqrt(5)/3)
0065   ASSERT_NEAR(bound_vec.qop(), -1.f / 3.7416574f, isclose);
0066   ASSERT_NEAR(bound_vec.time(), 0.1f, isclose);
0067
0068   // Check if the same free vector is obtained
0069   for (unsigned int i = 0u; i < 8u; i++) {
0070     ASSERT_NEAR(free_params[i], free_params2[i], isclose);
0071   }
0072
0073   // Test Jacobian transformation
0074   const bound_matrix<test_algebra> J =
0075       jac_engine::free_to_bound_jacobian<frame_type>(trf, free_params) *
0076       jac_engine::bound_to_free_jacobian<frame_type>(trf, rng, bound_vec);
0077
0078   for (unsigned int i = 0u; i < 6u; i++) {
0079     for (unsigned int j = 0u; j < 6u; j++) {
0080       if (i == j) {
0081         EXPECT_NEAR(getter::element(J, i, j), 1.f, isclose)
0082             << "at: " << i << ", " << j;
0083       } else {
0084         EXPECT_NEAR(getter::element(J, i, j), 0.f, isclose)
0085             << "at: " << i << ", " << j;
0086       }
0087     }
0088   }
0089 }