Core/Surfaces/PointSurfaceTests.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 #include <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Alignment.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/Geometry/GeometryContext.hpp"
0015 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0016 #include "Acts/Surfaces/PointSurface.hpp"
0017 #include "Acts/Surfaces/Surface.hpp"
0018 #include "Acts/Utilities/Intersection.hpp"
0019 #include "Acts/Utilities/Result.hpp"
0020 #include "Acts/Utilities/UnitVectors.hpp"
0021 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0022
0023 #include <memory>
0024 #include <numbers>
0025 #include <tuple>
0026 #include <vector>
0027
0028 using namespace Acts;
0029
0030 namespace ActsTests {
0031
0032 // Create a test context
0033 GeometryContext tgContext = GeometryContext::dangerouslyDefaultConstruct();
0034
0035 BOOST_AUTO_TEST_SUITE(SurfacesSuite)
0036
0037 /// Construction, type, name, bounds
0038 BOOST_AUTO_TEST_CASE(PointSurfaceConstruction) {
0039   Vector3 center{1., 2., 3.};
0040
0041   // unbounded
0042   auto unbounded = Surface::makeShared<PointSurface>(center);
0043   BOOST_CHECK_EQUAL(unbounded->type(), Surface::Point);
0044   BOOST_CHECK_EQUAL(unbounded->name(), "Acts::PointSurface");
0045   CHECK_CLOSE_ABS(unbounded->center(tgContext), center, 1e-12);
0046   // unbounded surface reports the shared boundless bounds
0047   BOOST_CHECK_EQUAL(unbounded->bounds().type(), SurfaceBounds::eBoundless);
0048   BOOST_CHECK(unbounded->boundsPtr() == nullptr);
0049
0050   // bounded
0051   auto bounded = Surface::makeShared<PointSurface>(center, 5.);
0052   BOOST_CHECK_EQUAL(bounded->type(), Surface::Point);
0053   BOOST_CHECK_EQUAL(bounded->bounds().type(), SurfaceBounds::ePoint);
0054
0055   // copy
0056   auto copied = Surface::makeShared<PointSurface>(*bounded);
0057   BOOST_CHECK_EQUAL(copied->type(), Surface::Point);
0058   BOOST_CHECK(*copied == *bounded);
0059 }
0060
0061 /// The measurement frame normal equals the (normalized) direction
0062 BOOST_AUTO_TEST_CASE(PointSurfaceNormalIsDirection) {
0063   auto surface = Surface::makeShared<PointSurface>(Vector3(0., 0., 0.));
0064   Vector3 dir = Vector3(1., 2., 3.).normalized();
0065   CHECK_CLOSE_ABS(surface->normal(tgContext, Vector3::Zero(), dir), dir, 1e-12);
0066
0067   auto rframe = surface->referenceFrame(tgContext, Vector3::Zero(), dir);
0068   CHECK_CLOSE_ABS(rframe.col(2), dir, 1e-12);
0069   // orthonormal frame
0070   CHECK_CLOSE_ABS(rframe.col(0).dot(rframe.col(1)), 0., 1e-12);
0071   CHECK_CLOSE_ABS(rframe.col(0).norm(), 1., 1e-12);
0072   CHECK_CLOSE_ABS(rframe.col(1).norm(), 1., 1e-12);
0073 }
0074
0075 /// Round-trip: intersect -> globalToLocal -> localToGlobal
0076 BOOST_AUTO_TEST_CASE(PointSurfaceRoundTrip) {
0077   Vector3 center{2., -1., 0.5};
0078   auto surface = Surface::makeShared<PointSurface>(center);
0079
0080   auto roundTrip = [&](const Vector3& pos, const Vector3& dir) {
0081     Intersection3D intersection =
0082         surface->intersect(tgContext, pos, dir).closest();
0083     Vector3 global = intersection.position();
0084     Vector2 local = *surface->globalToLocal(tgContext, global, dir);
0085     Vector3 global2 = surface->localToGlobal(tgContext, local, dir);
0086     return std::make_tuple(global, local, global2);
0087   };
0088
0089   {
0090     Vector3 pos = {-0.5, 0.3, 1.2};
0091     Vector3 dir = Vector3(-0.2, -0.7, 0.4).normalized();
0092     auto [global, local, global2] = roundTrip(pos, dir);
0093     CHECK_CLOSE_ABS(global, global2, 1e-10);
0094   }
0095   {
0096     Vector3 pos = {5., 5., 5.};
0097     Vector3 dir = Vector3(1., -1., 0.2).normalized();
0098     auto [global, local, global2] = roundTrip(pos, dir);
0099     CHECK_CLOSE_ABS(global, global2, 1e-10);
0100   }
0101 }
0102
0103 /// Round-trip stability including large eta and direction parallel to Z
0104 BOOST_AUTO_TEST_CASE(PointSurfaceRoundTripEtaStability) {
0105   Vector3 center = {3., 0., 0.};
0106   auto surface = Surface::makeShared<PointSurface>(center);
0107
0108   const std::vector<double> etas = {0, 1, 2, 3, 5, 10};
0109
0110   for (double eta : etas) {
0111     Vector3 dir = makeDirectionFromPhiEta(std::numbers::pi / 4., eta);
0112     Vector3 pos = center + Vector3(0.7, -0.4, 0.9);
0113
0114     Intersection3D intersection =
0115         surface->intersect(tgContext, pos, dir).closest();
0116     Vector3 global = intersection.position();
0117     Vector2 local = *surface->globalToLocal(tgContext, global, dir);
0118     Vector3 global2 = surface->localToGlobal(tgContext, local, dir);
0119
0120     CHECK_CLOSE_ABS(global, global2, 1e-9);
0121     // the PCA residual must be perpendicular to the direction
0122     CHECK_CLOSE_ABS((global - center).dot(dir), 0., 1e-9);
0123   }
0124
0125   // direction (anti-)parallel to Z (non-standard frame branch)
0126   for (Vector3 dir : {Vector3(0., 0., 1.), Vector3(0., 0., -1.)}) {
0127     Vector3 pos = center + Vector3(1., 2., 4.);
0128     Intersection3D intersection =
0129         surface->intersect(tgContext, pos, dir).closest();
0130     Vector3 global = intersection.position();
0131     Vector2 local = *surface->globalToLocal(tgContext, global, dir);
0132     Vector3 global2 = surface->localToGlobal(tgContext, local, dir);
0133     CHECK_CLOSE_ABS(global, global2, 1e-9);
0134     CHECK_CLOSE_ABS((global - center).dot(dir), 0., 1e-9);
0135   }
0136 }
0137
0138 /// The intersection is the point of closest approach; no degeneracy
0139 BOOST_AUTO_TEST_CASE(PointSurfaceIntersection) {
0140   Vector3 center{1., 1., 1.};
0141   auto surface = Surface::makeShared<PointSurface>(center);
0142
0143   Vector3 pos{0., 0., 0.};
0144   // Try several directions, including axis-aligned ones (would be degenerate
0145   // for a LineSurface but are fine here).
0146   const std::vector<Vector3> dirs = {
0147       Vector3(1., 0., 0.), Vector3(0., 1., 0.), Vector3(0., 0., 1.),
0148       Vector3(1., 1., 1.).normalized(), Vector3(-1., 2., -0.5).normalized()};
0149
0150   for (const Vector3& dir : dirs) {
0151     Intersection3D intersection =
0152         surface->intersect(tgContext, pos, dir).closest();
0153     BOOST_CHECK(intersection.isValid());
0154
0155     double u = (center - pos).dot(dir);
0156     Vector3 expected = pos + u * dir;
0157     CHECK_CLOSE_ABS(intersection.position(), expected, 1e-10);
0158     CHECK_CLOSE_ABS(intersection.pathLength(), u, 1e-10);
0159     // residual perpendicular to the direction
0160     CHECK_CLOSE_ABS((intersection.position() - center).dot(dir), 0., 1e-10);
0161   }
0162
0163   // starting exactly at the PCA reports onSurface
0164   Vector3 dir = Vector3(1., 0., 0.);
0165   Vector3 pca = pos + (center - pos).dot(dir) * dir;
0166   Intersection3D onSurf = surface->intersect(tgContext, pca, dir).closest();
0167   BOOST_CHECK_EQUAL(onSurf.status(), IntersectionStatus::onSurface);
0168 }
0169
0170 /// Bounds are respected in the intersection boundary check
0171 BOOST_AUTO_TEST_CASE(PointSurfaceIntersectionBounds) {
0172   Vector3 center{0., 0., 0.};
0173   // maximum distance 1
0174   auto bounded = Surface::makeShared<PointSurface>(center, 1.);
0175   auto unbounded = Surface::makeShared<PointSurface>(center);
0176
0177   // A track passing at transverse distance ~2 from the point
0178   Vector3 pos{2., 0., 0.};
0179   Vector3 dir = Vector3(0., 0., 1.);
0180
0181   auto within = BoundaryTolerance::None();
0182
0183   Intersection3D boundedHit =
0184       bounded->intersect(tgContext, pos, dir, within).closest();
0185   BOOST_CHECK_EQUAL(boundedHit.status(), IntersectionStatus::unreachable);
0186
0187   Intersection3D unboundedHit =
0188       unbounded->intersect(tgContext, pos, dir, within).closest();
0189   BOOST_CHECK(unboundedHit.isValid());
0190
0191   // A track passing close to the point is inside the bounds
0192   Vector3 posClose{0.3, 0., 0.};
0193   Intersection3D closeHit =
0194       bounded->intersect(tgContext, posClose, dir, within).closest();
0195   BOOST_CHECK(closeHit.isValid());
0196 }
0197
0198 // Build a free vector from a bound vector consistent with the surface frame
0199 FreeVector boundToFree(const PointSurface& surface, const BoundVector& b) {
0200   Vector3 dir = makeDirectionFromPhiTheta(b[eBoundPhi], b[eBoundTheta]);
0201   Vector3 pos = surface.localToGlobal(
0202       tgContext, Vector2(b[eBoundLoc0], b[eBoundLoc1]), dir);
0203   FreeVector f = FreeVector::Zero();
0204   f.segment<3>(eFreePos0) = pos;
0205   f[eFreeTime] = b[eBoundTime];
0206   f.segment<3>(eFreeDir0) = dir;
0207   f[eFreeQOverP] = b[eBoundQOverP];
0208   return f;
0209 }
0210
0211 /// boundToFreeJacobian validated against finite differences at nonzero local
0212 BOOST_AUTO_TEST_CASE(PointSurfaceBoundToFreeJacobian) {
0213   auto surface = Surface::makeShared<PointSurface>(Vector3(1., -2., 0.5));
0214
0215   // pick a non-forward direction (standard frame branch) and nonzero local so
0216   // the frame-rotation correction is exercised
0217   double phi = 0.6;
0218   double theta = std::numbers::pi / 3.;
0219   BoundVector b;
0220   b << 0.7, -0.4, phi, theta, 0.5, 1.5;
0221
0222   FreeVector f0 = boundToFree(*surface, b);
0223   Vector3 pos0 = f0.segment<3>(eFreePos0);
0224   Vector3 dir0 = f0.segment<3>(eFreeDir0);
0225
0226   BoundToFreeMatrix jac = surface->boundToFreeJacobian(tgContext, pos0, dir0);
0227
0228   double h = 1e-6;
0229   for (int i = 0; i < static_cast<int>(eBoundSize); ++i) {
0230     BoundVector bp = b;
0231     BoundVector bm = b;
0232     bp[i] += h;
0233     bm[i] -= h;
0234     FreeVector numeric =
0235         (boundToFree(*surface, bp) - boundToFree(*surface, bm)) / (2 * h);
0236     for (int j = 0; j < static_cast<int>(eFreeSize); ++j) {
0237       BOOST_CHECK_SMALL(numeric[j] - jac(j, i), 1e-6);
0238     }
0239   }
0240 }
0241
0242 /// freeToPathDerivative validated against finite differences
0243 BOOST_AUTO_TEST_CASE(PointSurfaceFreeToPathDerivative) {
0244   Vector3 center{1., -2., 0.5};
0245   auto surface = Surface::makeShared<PointSurface>(center);
0246
0247   Vector3 dir = Vector3(0.3, -0.5, 0.8).normalized();
0248   // start on the surface (PCA) so the assert holds
0249   Vector3 pos = center + Vector3(0.4, 0.1, -0.2);
0250   // project pos onto the measurement plane through the center
0251   pos -= (pos - center).dot(dir) * dir;
0252
0253   FreeToPathMatrix analytic =
0254       surface->freeToPathDerivative(tgContext, pos, dir);
0255
0256   auto pathU = [&](const Vector3& p, const Vector3& d) {
0257     return (center - p).dot(d);
0258   };
0259
0260   double h = 1e-7;
0261   for (int k = 0; k < 3; ++k) {
0262     Vector3 pp = pos, pm = pos;
0263     pp[k] += h;
0264     pm[k] -= h;
0265     double numeric = (pathU(pp, dir) - pathU(pm, dir)) / (2 * h);
0266     BOOST_CHECK_SMALL(numeric - analytic(0, eFreePos0 + k), 1e-6);
0267   }
0268   for (int k = 0; k < 3; ++k) {
0269     Vector3 dp = dir, dm = dir;
0270     dp[k] += h;
0271     dm[k] -= h;
0272     double numeric = (pathU(pos, dp) - pathU(pos, dm)) / (2 * h);
0273     BOOST_CHECK_SMALL(numeric - analytic(0, eFreeDir0 + k), 1e-6);
0274   }
0275 }
0276
0277 /// alignmentToPathDerivative: only the center contributes (= direction)
0278 BOOST_AUTO_TEST_CASE(PointSurfaceAlignmentToPathDerivative) {
0279   Vector3 center{1., -2., 0.5};
0280   auto surface = Surface::makeShared<PointSurface>(center);
0281   Vector3 dir = Vector3(0.3, -0.5, 0.8).normalized();
0282   Vector3 pos = center + Vector3(0.4, 0.1, -0.2);
0283   pos -= (pos - center).dot(dir) * dir;
0284
0285   AlignmentToPathMatrix analytic =
0286       surface->alignmentToPathDerivative(tgContext, pos, dir);
0287
0288   // center block equals the direction
0289   CHECK_CLOSE_ABS(analytic.segment<3>(eAlignmentCenter0).transpose(), dir,
0290                   1e-10);
0291   // rotation block is zero
0292   CHECK_CLOSE_ABS(analytic.segment<3>(eAlignmentRotation0).norm(), 0., 1e-12);
0293 }
0294
0295 BOOST_AUTO_TEST_SUITE_END()
0296
0297 }  // namespace ActsTests