navigation/intersection/helix_intersector.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/navigation/intersection/helix_intersector.hpp"
0011
0012 #include "detray/definitions/units.hpp"
0013 #include "detray/geometry/mask.hpp"
0014 #include "detray/geometry/shapes/concentric_cylinder2D.hpp"
0015 #include "detray/geometry/shapes/line.hpp"
0016 #include "detray/geometry/shapes/rectangle2D.hpp"
0017 #include "detray/geometry/shapes/unmasked.hpp"
0018 #include "detray/geometry/surface_descriptor.hpp"
0019 #include "detray/tracks/helix.hpp"
0020 #include "detray/tracks/tracks.hpp"
0021
0022 // Detray test include(s)
0023 #include "detray/test/framework/types.hpp"
0024
0025 // Google Test include(s).
0026 #include <gtest/gtest.h>
0027
0028 // System include(s)
0029 #include <limits>
0030
0031 using namespace detray;
0032
0033 namespace {
0034
0035 // Three-dimensional definitions
0036 using test_algebra = test::algebra;
0037 using transform3_t = test::transform3;
0038 using vector3 = test::vector3;
0039 using point3 = test::point3;
0040 using scalar = test::scalar;
0041 using helix_t = detray::detail::helix<test_algebra>;
0042 using intersection_t = intersection2D<surface_descriptor<>, test_algebra,
0043                                       intersection::contains_pos>;
0044
0045 constexpr auto not_defined{detail::invalid_value<scalar>()};
0046 constexpr scalar tol{1e-4f};
0047
0048 // z axis
0049 const vector3 z_axis{0.f, 0.f, 1.f};
0050
0051 // Track defined on origin point
0052 const free_track_parameters<test_algebra> free_trk({0.f, 0.f, 0.f}, 0.f,
0053                                                    {0.1f * unit<scalar>::GeV,
0054                                                     static_cast<scalar>(0.f),
0055                                                     static_cast<scalar>(0.f)},
0056                                                    -1.f);
0057
0058 // Magnetic field
0059 const vector3 B{static_cast<scalar>(0.f), static_cast<scalar>(0.f),
0060                 1.f * unit<scalar>::T};
0061 const vector3 B_0{0.f * unit<scalar>::T, tol * unit<scalar>::T,
0062                   tol * unit<scalar>::T};
0063
0064 // Test helix
0065 const helix_t hlx(free_trk, B);
0066
0067 // Path along the helix
0068 const scalar path = 10.f * unit<scalar>::cm;
0069
0070 // Transform translation vector
0071 const vector3 trl = hlx(path);
0072
0073 // Surface normal vector
0074 const vector3 w = hlx.dir(path);
0075
0076 }  // anonymous namespace
0077
0078 /// This defines the local frame test suite
0079 GTEST_TEST(detray_intersection, helix_plane_intersector_no_bfield) {
0080   // Create a shifted plane
0081   const transform3_t shifted(vector3{3.f, 2.f, 10.f});
0082
0083   // Test helix
0084   const point3 pos{2.f, 1.f, 0.f};
0085   const vector3 mom{0.f, 0.f, 1.f};
0086   const detail::helix<test_algebra> h({pos, 0.f, mom, -1.f}, B_0);
0087
0088   // The same test but bound to local frame
0089   helix_intersector<unmasked<2>, test_algebra> pi;
0090   mask<unmasked<2>, test_algebra> unmasked_bound{};
0091   const auto hit_bound =
0092       pi(h, surface_descriptor<>{}, unmasked_bound, shifted, tol);
0093
0094   ASSERT_TRUE(hit_bound.is_inside());
0095   // Global intersection information - unchanged
0096   const auto global0 =
0097       unmasked_bound.to_global_frame(shifted, hit_bound.local());
0098   ASSERT_NEAR(global0[0], 2.f, tol);
0099   ASSERT_NEAR(global0[1], 1.f, tol);
0100   ASSERT_NEAR(global0[2], 10.f, tol);
0101   // Local intersection information
0102   ASSERT_NEAR(hit_bound.local()[0], -1.f, tol);
0103   ASSERT_NEAR(hit_bound.local()[1], -1.f, tol);
0104
0105   // The same test but bound to local frame & masked - inside
0106   mask<rectangle2D, test_algebra> rect_for_inside{0u, 3.f, 3.f};
0107   const auto hit_bound_inside =
0108       pi(h, surface_descriptor<>{}, rect_for_inside, shifted, tol);
0109   ASSERT_TRUE(hit_bound_inside.is_inside());
0110   // Global intersection information - unchanged
0111   const auto global1 =
0112       rect_for_inside.to_global_frame(shifted, hit_bound_inside.local());
0113   ASSERT_NEAR(global1[0], 2.f, tol);
0114   ASSERT_NEAR(global1[1], 1.f, tol);
0115   ASSERT_NEAR(global1[2], 10.f, tol);
0116   // Local intersection infoimation - unchanged
0117   ASSERT_NEAR(hit_bound_inside.local()[0], -1.f, tol);
0118   ASSERT_NEAR(hit_bound_inside.local()[1], -1.f, tol);
0119
0120   // The same test but bound to local frame & masked - outside
0121   mask<rectangle2D, test_algebra> rect_for_outside{0u, 0.5f, 3.5f};
0122   const auto hit_bound_outside =
0123       pi(h, surface_descriptor<>{}, rect_for_outside, shifted, tol);
0124   ASSERT_FALSE(hit_bound_outside.is_inside());
0125   const auto global2 =
0126       rect_for_outside.to_global_frame(shifted, hit_bound_outside.local());
0127   // Global intersection information - unchanged
0128   ASSERT_NEAR(global2[0], 2.f, tol);
0129   ASSERT_NEAR(global2[1], 1.f, tol);
0130   ASSERT_NEAR(global2[2], 10.f, tol);
0131   // Local intersection infoimation - unchanged
0132   ASSERT_NEAR(hit_bound_outside.local()[0], -1.f, tol);
0133   ASSERT_NEAR(hit_bound_outside.local()[1], -1.f, tol);
0134 }
0135
0136 /// Test the intersection between a helical trajectory and a plane
0137 GTEST_TEST(detray_intersection, helix_plane_intersector) {
0138   // Vector on the surface
0139   const vector3 v = vector::cross(z_axis, w);
0140
0141   // Transform matrix
0142   const transform3_t trf(trl, w, v);
0143
0144   // Rectangle surface
0145   const mask<rectangle2D, test_algebra> rectangle{0u, 10.f * unit<scalar>::cm,
0146                                                   10.f * unit<scalar>::cm};
0147
0148   const helix_intersector<rectangle2D, test_algebra> hpi;
0149
0150   // Get the intersection on the next surface
0151   const auto is = hpi(hlx, surface_descriptor<>{}, rectangle, trf, tol);
0152
0153   // Check the values
0154   EXPECT_TRUE(is.is_inside());
0155   EXPECT_NEAR(is.path(), path, tol);
0156   EXPECT_NEAR(is.local()[0], 0.f, tol);
0157   EXPECT_NEAR(is.local()[1], 0.f, tol);
0158
0159   const auto global = rectangle.to_global_frame(trf, is.local());
0160   EXPECT_NEAR(global[0], trl[0], tol);
0161   EXPECT_NEAR(global[1], trl[1], tol);
0162   EXPECT_NEAR(global[2], trl[2], tol);
0163 }
0164
0165 /// This checks the closest solution of a helix-cylinder intersection
0166 GTEST_TEST(detray_intersection, helix_cylinder_intersector_no_bfield) {
0167   const scalar r{4.f * unit<scalar>::mm};
0168   const scalar hz{10.f * unit<scalar>::mm};
0169
0170   // Create a translated cylinder and test untersection
0171   const transform3_t shifted(vector3{3.f, 2.f, 10.f});
0172   helix_intersector<cylinder2D, test_algebra> hi;
0173
0174   // Test helix
0175   const point3 pos{3.f, 2.f, 5.f};
0176   const vector3 mom{1.f, 0.f, 0.f};
0177   const helix_t h({pos, 0.f * unit<scalar>::s, mom, -1 * unit<scalar>::e}, B_0);
0178
0179   // Intersect
0180   mask<cylinder2D, test_algebra, std::uint_least16_t> cylinder{0u, r, -hz, hz};
0181   const auto hits_bound = hi(h, surface_descriptor<>{}, cylinder, shifted, tol);
0182
0183   // No magnetic field, so the solutions must be the same as for a ray
0184
0185   // second intersection lies in front of the track
0186   EXPECT_TRUE(hits_bound[0].is_inside());
0187   EXPECT_FALSE(hits_bound[0].is_along());
0188
0189   const auto global0 = cylinder.to_global_frame(shifted, hits_bound[0].local());
0190
0191   EXPECT_NEAR(global0[0], -1.f, tol);
0192   EXPECT_NEAR(global0[1], 2.f, tol);
0193   EXPECT_NEAR(global0[2], 5.f, tol);
0194   ASSERT_TRUE(hits_bound[0].local()[0] != not_defined &&
0195               hits_bound[0].local()[1] != not_defined);
0196   // p2[0] = r * phi : 180deg in the opposite direction with r = 4
0197   EXPECT_NEAR(hits_bound[0].local()[0], 4.f * constant<scalar>::pi, tol);
0198   EXPECT_NEAR(hits_bound[0].local()[1], -5.f, tol);
0199
0200   // first intersection lies behind the track
0201   const auto global1 = cylinder.to_global_frame(shifted, hits_bound[1].local());
0202   EXPECT_TRUE(hits_bound[1].is_inside());
0203   EXPECT_TRUE(hits_bound[1].is_along());
0204   EXPECT_NEAR(global1[0], 7.f, tol);
0205   EXPECT_NEAR(global1[1], 2.f, tol);
0206   EXPECT_NEAR(global1[2], 5.f, tol);
0207   ASSERT_TRUE(hits_bound[1].local()[0] != not_defined &&
0208               hits_bound[1].local()[1] != not_defined);
0209   EXPECT_NEAR(hits_bound[1].local()[0], 0.f, tol);
0210   EXPECT_NEAR(hits_bound[1].local()[1], -5., tol);
0211 }
0212
0213 /// Test the intersection between a helical trajectory and a cylinder
0214 GTEST_TEST(detray_intersection, helix_cylinder_intersector) {
0215   // Transform matrix
0216   const transform3_t trf(trl, z_axis, w);
0217
0218   // Cylinder surface (5 cm radius)
0219   const scalar r{4.f * unit<scalar>::cm};
0220   const scalar hz{10.f * unit<scalar>::cm};
0221   const mask<cylinder2D, test_algebra> cylinder{0u, r, -hz, hz};
0222
0223   const helix_intersector<cylinder2D, test_algebra> hci;
0224
0225   // Get the intersection on the next surface
0226   const auto is = hci(hlx, surface_descriptor<>{}, cylinder, trf, tol);
0227
0228   // First solution
0229   const vector3 pos_near = hlx.pos(is[0].path());
0230   const vector3 loc_near = pos_near - trl;
0231   const scalar phi_near = std::acos(vector::dot(w, loc_near) /
0232                                     (vector::norm(w) * vector::norm(loc_near)));
0233
0234   EXPECT_TRUE(is[0].is_inside());
0235   // Not precise due to helix curvature
0236   EXPECT_NEAR(is[0].path(), path - r, 5000.f * tol);
0237   EXPECT_NEAR(is[0].local()[0], r * phi_near, tol);
0238   EXPECT_NEAR(is[0].local()[1], 0.f, tol);
0239   const auto global0 = cylinder.to_global_frame(trf, is[0].local());
0240   EXPECT_NEAR(global0[0], pos_near[0], tol);
0241   EXPECT_NEAR(global0[1], pos_near[1], tol);
0242   EXPECT_NEAR(global0[2], pos_near[2], tol);
0243
0244   // Second solution
0245   const vector3 pos_far = hlx.pos(is[1].path());
0246   const vector3 loc_far = pos_far - trl;
0247   const scalar phi_far = std::acos(vector::dot(w, loc_far) /
0248                                    (vector::norm(w) * vector::norm(loc_far)));
0249
0250   EXPECT_TRUE(is[1].is_inside());
0251   // Not precise due to helix curvature
0252   EXPECT_NEAR(is[1].path(), path + r, 5000.f * tol);
0253   EXPECT_NEAR(is[1].local()[0], r * phi_far, tol);
0254   EXPECT_NEAR(is[1].local()[1], 0.f, tol);
0255   const auto global1 = cylinder.to_global_frame(trf, is[1].local());
0256   EXPECT_NEAR(global1[0], pos_far[0], tol);
0257   EXPECT_NEAR(global1[1], pos_far[1], tol);
0258   EXPECT_NEAR(global1[2], pos_far[2], tol);
0259 }
0260
0261 /// This checks the closest solution of a helix-concentric cylinder intersection
0262 GTEST_TEST(detray_intersection,
0263            helix_concentric_cylinder_intersector_no_bfield) {
0264   const scalar r{4.f * unit<scalar>::mm};
0265   const scalar hz{10.f * unit<scalar>::mm};
0266
0267   // Create a translated cylinder and test untersection
0268   const transform3_t identity{};
0269   helix_intersector<concentric_cylinder2D, test_algebra> c_hi;
0270
0271   // Test helix
0272   const point3 pos{0.f, 0.f, -5.f};
0273   const vector3 mom{1.f, 0.f, 0.f};
0274   const helix_t h({pos, 0.f * unit<scalar>::s, mom, -1 * unit<scalar>::e}, B_0);
0275
0276   // Intersect
0277   mask<concentric_cylinder2D, test_algebra, std::uint_least16_t> cylinder{
0278       0u, r, -hz, hz};
0279   const auto hits_bound =
0280       c_hi(h, surface_descriptor<>{}, cylinder, identity, tol);
0281
0282   // No magnetic field, so the solutions must be the same as for a ray
0283
0284   // second intersection lies in front of the track
0285   EXPECT_TRUE(hits_bound.is_inside());
0286   EXPECT_TRUE(hits_bound.is_along());
0287
0288   const auto global0 = cylinder.to_global_frame(identity, hits_bound.local());
0289
0290   EXPECT_NEAR(global0[0], 4.f, tol);
0291   EXPECT_NEAR(global0[1], 0.f, tol);
0292   EXPECT_NEAR(global0[2], -5.f, tol);
0293   EXPECT_NEAR(hits_bound.local()[0], 0.f, tol);
0294   EXPECT_NEAR(hits_bound.local()[1], -5.f, tol);
0295 }
0296
0297 /// Test the intersection between a helical trajectory and a line
0298 GTEST_TEST(detray_intersection, helix_line_intersector) {
0299   // Intersector object
0300   const helix_intersector<line_circular, test_algebra> hli;
0301
0302   // Get radius of track
0303   const scalar R{hlx.radius()};
0304
0305   // Path length for pi/4
0306   const scalar s0 = constant<scalar>::pi_2 * R;
0307
0308   // Wire properties
0309   const scalar scope = 2.f * unit<scalar>::cm;
0310   const scalar half_z = std::numeric_limits<scalar>::max();
0311
0312   // Straw wire
0313   const mask<line_circular, test_algebra> straw_tube{0u, scope, half_z};
0314
0315   // Cell wire
0316   const mask<line_square, test_algebra> drift_cell{0u, scope, half_z};
0317
0318   // Offset to shift the translation of transform matrix
0319   const scalar offset = 1.f * unit<scalar>::cm;
0320
0321   //---------------------
0322   // Forward direction
0323   //---------------------
0324
0325   // Reference point for transform matrix
0326   const point3 r0_fw = hlx.pos(s0);
0327
0328   // Translation is shifted from reference point
0329   const point3 trl_fw = r0_fw + vector3{offset, static_cast<scalar>(0.f),
0330                                         static_cast<scalar>(0.f)};
0331
0332   // Transform matrix
0333   const transform3_t trf_fw(trl_fw, z_axis, hlx.dir(s0));
0334
0335   // Get the intersection on the next surface
0336   auto is = hli(hlx, surface_descriptor<>{}, straw_tube, trf_fw, tol);
0337
0338   EXPECT_NEAR(is.path(), s0, tol);
0339   // track (helix) is at the left side w.r.t wire
0340   EXPECT_NEAR(is.local()[0], offset, tol);
0341   EXPECT_NEAR(is.local()[1], 0.f, tol);
0342   EXPECT_TRUE(is.is_inside());
0343   EXPECT_TRUE(is.is_along());
0344
0345   // Get the intersection on the next surface
0346   is = hli(hlx, surface_descriptor<>{}, drift_cell, trf_fw, tol);
0347
0348   EXPECT_NEAR(is.path(), s0, tol);
0349   // track (helix) is at the left side w.r.t wire
0350   EXPECT_NEAR(is.local()[0], offset, tol);
0351   EXPECT_NEAR(is.local()[1], 0.f, tol);
0352   EXPECT_TRUE(is.is_inside());
0353   EXPECT_TRUE(is.is_along());
0354
0355   //---------------------
0356   // Backward direction
0357   //---------------------
0358
0359   // Reference point for transform matrix
0360   const point3 r0_bw = hlx.pos(-s0);
0361
0362   // Translation is shifted from reference point
0363   const point3 trl_bw = r0_bw + vector3{offset, static_cast<scalar>(0.f),
0364                                         static_cast<scalar>(0.f)};
0365
0366   // Transform matrix
0367   const transform3_t trf_bw(trl_bw, z_axis, hlx.dir(-s0));
0368
0369   // Get the intersection on the next surface
0370   is = hli(hlx, surface_descriptor<>{}, straw_tube, trf_bw, tol);
0371
0372   EXPECT_NEAR(is.path(), -s0, tol);
0373   // track (helix) is at the right side w.r.t wire
0374   EXPECT_NEAR(is.local()[0], -offset, tol);
0375   EXPECT_NEAR(is.local()[1], 0.f, tol);
0376   EXPECT_TRUE(is.is_inside());
0377   EXPECT_FALSE(is.is_along());
0378
0379   // Get the intersection on the next surface
0380   is = hli(hlx, surface_descriptor<>{}, drift_cell, trf_bw, tol);
0381
0382   EXPECT_NEAR(is.path(), -s0, tol);
0383   // track (helix) is at the right side w.r.t wire
0384   EXPECT_NEAR(is.local()[0], -offset, tol);
0385   EXPECT_NEAR(is.local()[1], 0.f, tol);
0386   EXPECT_TRUE(is.is_inside());
0387   EXPECT_FALSE(is.is_along());
0388 }