navigation/intersection/helix_trajectory.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/tracks/helix.hpp"
0012 #include "detray/tracks/tracks.hpp"
0013
0014 // Detray test include(s)
0015 #include "detray/test/framework/types.hpp"
0016
0017 // GTest include(s)
0018 #include <gtest/gtest.h>
0019
0020 using namespace detray;
0021
0022 using scalar = test::scalar;
0023
0024 constexpr const scalar tol{1e-5f};
0025
0026 // This tests the base functionality of the Helix Gun
0027 GTEST_TEST(detray_intersection, helix_trajectory) {
0028   using test_algebra = test::algebra;
0029   using vector3 = test::vector3;
0030   using point3 = test::point3;
0031
0032   const point3 pos{0.f, 0.f, 0.f};
0033   const scalar time{0.f};
0034   const vector3 mom{static_cast<scalar>(1.f), static_cast<scalar>(0.f),
0035                     1.f * unit<scalar>::GeV};
0036   const scalar q{static_cast<scalar>(-1.) * unit<scalar>::e};
0037
0038   // vertex
0039   free_track_parameters<test_algebra> vertex(pos, time, mom, q);
0040
0041   // magnetic field
0042   const vector3 B{static_cast<scalar>(0.f), static_cast<scalar>(0.f),
0043                   1.f * unit<scalar>::T};
0044
0045   const scalar p_mag{vector::norm(mom)};
0046   const scalar B_mag{vector::norm(B)};
0047   const scalar pz_along{vector::dot(mom, vector::normalize(B))};
0048   const scalar pt{std::sqrt(p_mag * p_mag - pz_along * pz_along)};
0049
0050   // helix trajectory
0051   detail::helix helix_traj(vertex, B);
0052   EXPECT_NEAR(helix_traj.time(), 0.f, tol);
0053   EXPECT_NEAR(helix_traj.qop(), -constant<scalar>::inv_sqrt2, tol);
0054
0055   // radius of helix
0056   scalar R{helix_traj.radius()};
0057   EXPECT_NEAR(R, pt / B_mag, tol);
0058
0059   // (1 T of b_z, 1 GeV/c of p_T) ==> 3.33564095 m of radius
0060   EXPECT_NEAR(R, 3.33564095f * unit<scalar>::m, 10.f * tol);
0061
0062   // Path length for one loop
0063   scalar S = 2.f * p_mag / B_mag * constant<scalar>::pi;
0064
0065   // After half turn
0066   point3 half_loop_pos = helix_traj(S / 2.f);
0067   EXPECT_NEAR(half_loop_pos[0], 0.f, R * tol);
0068   EXPECT_NEAR(half_loop_pos[1], 2.f * R, R * tol);
0069   EXPECT_NEAR(half_loop_pos[2], pz_along / B_mag * constant<scalar>::pi,
0070               R * tol);
0071
0072   point3 half_loop_dir = helix_traj.dir(S / 2.f);
0073   EXPECT_NEAR(half_loop_dir[0], -vertex.dir()[0], R * tol);
0074   EXPECT_NEAR(half_loop_dir[1], -vertex.dir()[1], R * tol);
0075   EXPECT_NEAR(half_loop_dir[2], vertex.dir()[2], R * tol);
0076
0077   // After half turn in the opposite direction
0078   half_loop_pos = helix_traj(-S / 2.f);
0079   EXPECT_NEAR(half_loop_pos[0], 0.f, R * tol);
0080   EXPECT_NEAR(half_loop_pos[1], 2.f * R, R * tol);
0081   EXPECT_NEAR(half_loop_pos[2], -pz_along / B_mag * constant<scalar>::pi,
0082               R * tol);
0083
0084   half_loop_dir = helix_traj.dir(-S / 2.f);
0085   EXPECT_NEAR(half_loop_dir[0], -vertex.dir()[0], R * tol);
0086   EXPECT_NEAR(half_loop_dir[1], -vertex.dir()[1], R * tol);
0087   EXPECT_NEAR(half_loop_dir[2], vertex.dir()[2], R * tol);
0088
0089   // After one full turn
0090   point3 one_loop_pos = helix_traj(S);
0091   EXPECT_NEAR(one_loop_pos[0], 0.f, R * tol);
0092   EXPECT_NEAR(one_loop_pos[1], 0.f, R * tol);
0093   EXPECT_NEAR(one_loop_pos[2], 2.f * pz_along / B_mag * constant<scalar>::pi,
0094               R * tol);
0095
0096   point3 one_loop_dir = helix_traj.dir(S);
0097   EXPECT_NEAR(one_loop_dir[0], vertex.dir()[0], R * tol);
0098   EXPECT_NEAR(one_loop_dir[1], vertex.dir()[1], R * tol);
0099   EXPECT_NEAR(one_loop_dir[2], vertex.dir()[2], R * tol);
0100
0101   // After one full turn in the opposite direction
0102   one_loop_pos = helix_traj(-S);
0103   EXPECT_NEAR(one_loop_pos[0], 0.f, R * tol);
0104   EXPECT_NEAR(one_loop_pos[1], 0.f, R * tol);
0105   EXPECT_NEAR(one_loop_pos[2], -2.f * pz_along / B_mag * constant<scalar>::pi,
0106               R * tol);
0107
0108   one_loop_dir = helix_traj.dir(-S);
0109   EXPECT_NEAR(one_loop_dir[0], vertex.dir()[0], R * tol);
0110   EXPECT_NEAR(one_loop_dir[1], vertex.dir()[1], R * tol);
0111   EXPECT_NEAR(one_loop_dir[2], vertex.dir()[2], R * tol);
0112
0113   /*********************************
0114    * Same test with oppsite charge
0115    *********************************/
0116
0117   free_track_parameters<test_algebra> vertex2(pos, time, mom, -q);
0118
0119   // helix trajectory
0120   detail::helix helix_traj2(vertex2, B);
0121
0122   EXPECT_NEAR(R, helix_traj2.radius(), tol);
0123
0124   // After half turn
0125   half_loop_pos = helix_traj2(S / 2.f);
0126   EXPECT_NEAR(half_loop_pos[0], 0.f, R * tol);
0127   EXPECT_NEAR(half_loop_pos[1], -2.f * R, R * tol);
0128   EXPECT_NEAR(half_loop_pos[2], pz_along / B_mag * constant<scalar>::pi,
0129               R * tol);
0130
0131   half_loop_dir = helix_traj2.dir(S / 2.f);
0132   EXPECT_NEAR(half_loop_dir[0], -vertex2.dir()[0], R * tol);
0133   EXPECT_NEAR(half_loop_dir[1], -vertex2.dir()[1], R * tol);
0134   EXPECT_NEAR(half_loop_dir[2], vertex2.dir()[2], R * tol);
0135
0136   // After half turn in the opposite direction
0137   half_loop_pos = helix_traj2(-S / 2.f);
0138   EXPECT_NEAR(half_loop_pos[0], 0.f, R * tol);
0139   EXPECT_NEAR(half_loop_pos[1], -2.f * R, R * tol);
0140   EXPECT_NEAR(half_loop_pos[2], -pz_along / B_mag * constant<scalar>::pi,
0141               R * tol);
0142
0143   half_loop_dir = helix_traj.dir(-S / 2.f);
0144   EXPECT_NEAR(half_loop_dir[0], -vertex.dir()[0], R * tol);
0145   EXPECT_NEAR(half_loop_dir[1], -vertex.dir()[1], R * tol);
0146   EXPECT_NEAR(half_loop_dir[2], vertex.dir()[2], R * tol);
0147
0148   // After one full turn
0149   one_loop_pos = helix_traj2(S);
0150   EXPECT_NEAR(one_loop_pos[0], 0.f, R * tol);
0151   EXPECT_NEAR(one_loop_pos[1], 0.f, R * tol);
0152   EXPECT_NEAR(one_loop_pos[2], 2.f * pz_along / B_mag * constant<scalar>::pi,
0153               R * tol);
0154
0155   one_loop_dir = helix_traj2.dir(S);
0156   EXPECT_NEAR(one_loop_dir[0], vertex2.dir()[0], R * tol);
0157   EXPECT_NEAR(one_loop_dir[1], vertex2.dir()[1], R * tol);
0158   EXPECT_NEAR(one_loop_dir[2], vertex2.dir()[2], R * tol);
0159
0160   // After one full turn in the opposite direction
0161   one_loop_pos = helix_traj2(-S);
0162   EXPECT_NEAR(one_loop_pos[0], 0.f, R * tol);
0163   EXPECT_NEAR(one_loop_pos[1], 0.f, R * tol);
0164   EXPECT_NEAR(one_loop_pos[2], -2.f * pz_along / B_mag * constant<scalar>::pi,
0165               R * tol);
0166
0167   one_loop_dir = helix_traj2.dir(-S);
0168   EXPECT_NEAR(one_loop_dir[0], vertex2.dir()[0], R * tol);
0169   EXPECT_NEAR(one_loop_dir[1], vertex2.dir()[1], R * tol);
0170   EXPECT_NEAR(one_loop_dir[2], vertex2.dir()[2], R * tol);
0171 }
0172
0173 GTEST_TEST(detray_intersection, helix_trajectory_small_pT) {
0174   using test_algebra = test::algebra;
0175   using vector3 = test::vector3;
0176   using point3 = test::point3;
0177
0178   const point3 pos{0.f, 0.f, 0.f};
0179   const scalar time{0.f};
0180   const vector3 mom{static_cast<scalar>(0.f), tol, 1.f * unit<scalar>::GeV};
0181   const scalar q{static_cast<scalar>(-1.) * unit<scalar>::e};
0182
0183   // vertex
0184   free_track_parameters<test_algebra> vertex(pos, time, mom, q);
0185
0186   // magnetic field
0187   const vector3 B{static_cast<scalar>(0.f), static_cast<scalar>(0.f),
0188                   1.f * unit<scalar>::T};
0189
0190   // helix trajectory
0191   detail::helix helix_traj(vertex, B);
0192   EXPECT_NEAR(helix_traj.time(), 0.f, tol);
0193   EXPECT_NEAR(helix_traj.qop(), -1.f, tol);
0194
0195   // After 10 mm
0196   const scalar path_length{10.f * unit<scalar>::mm};
0197   const point3 helix_pos = helix_traj(path_length);
0198   const point3 true_pos = pos + path_length * vector::normalize(mom);
0199
0200   EXPECT_NEAR(true_pos[0], helix_pos[0], tol);
0201   EXPECT_NEAR(true_pos[1], helix_pos[1], tol);
0202   EXPECT_NEAR(true_pos[2], helix_pos[2], tol);
0203 }
0204
0205 GTEST_TEST(detray_intersection, helix_direction_stability) {
0206   using test_algebra = test::algebra;
0207   using vector3 = test::vector3;
0208   using point3 = test::point3;
0209
0210   // magnetic field
0211   const vector3 B{static_cast<scalar>(0.f), static_cast<scalar>(0.f),
0212                   1.f * unit<scalar>::T};
0213
0214   const point3 pos{0.f, 0.f, 0.f};
0215   const scalar time{0.f};
0216   const vector3 mom{1.f * unit<scalar>::GeV, 1.f * unit<scalar>::GeV,
0217                     1.f * unit<scalar>::GeV};
0218   const scalar q{static_cast<scalar>(-1.) * unit<scalar>::e};
0219
0220   // vertex
0221   free_track_parameters<test_algebra> vertex(pos, time, mom, q);
0222
0223   // helix trajectory
0224   detail::helix hlx(vertex, B);
0225
0226   for (int i = 0; i < 100; i++) {
0227     const auto d = hlx.dir(static_cast<scalar>(i) * 10.f);
0228     ASSERT_FLOAT_EQ(static_cast<float>(vector::theta(d)),
0229                     static_cast<float>(vector::theta(mom)));
0230   }
0231 }