Core/Surfaces/LineSurfaceTests.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/tools/output_test_stream.hpp>
0010 #include <boost/test/unit_test.hpp>
0011
0012 #include "Acts/Definitions/Algebra.hpp"
0013 #include "Acts/Definitions/Alignment.hpp"
0014 #include "Acts/Definitions/TrackParametrization.hpp"
0015 #include "Acts/EventData/ParticleHypothesis.hpp"
0016 #include "Acts/EventData/detail/GenerateParameters.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/Material/BinnedSurfaceMaterial.hpp"
0019 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
0020 #include "Acts/Material/Material.hpp"
0021 #include "Acts/Material/MaterialSlab.hpp"
0022 #include "Acts/Propagator/Propagator.hpp"
0023 #include "Acts/Propagator/StraightLineStepper.hpp"
0024 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0025 #include "Acts/Surfaces/LineBounds.hpp"
0026 #include "Acts/Surfaces/Surface.hpp"
0027 #include "Acts/Utilities/AxisDefinitions.hpp"
0028 #include "Acts/Utilities/BinUtility.hpp"
0029 #include "Acts/Utilities/BinningType.hpp"
0030 #include "Acts/Utilities/Intersection.hpp"
0031 #include "Acts/Utilities/Result.hpp"
0032 #include "Acts/Utilities/ThrowAssert.hpp"
0033 #include "Acts/Utilities/UnitVectors.hpp"
0034 #include "Acts/Utilities/VectorHelpers.hpp"
0035 #include "ActsTests/CommonHelpers/DetectorElementStub.hpp"
0036 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0037 #include "ActsTests/CommonHelpers/LineSurfaceStub.hpp"
0038 #include "ActsTests/CommonHelpers/PredefinedMaterials.hpp"
0039
0040 #include <cmath>
0041 #include <memory>
0042 #include <numbers>
0043 #include <optional>
0044 #include <vector>
0045
0046 using namespace Acts;
0047
0048 namespace ActsTests {
0049
0050 // Create a test context
0051 GeometryContext tgContext = GeometryContext::dangerouslyDefaultConstruct();
0052
0053 BOOST_AUTO_TEST_SUITE(SurfacesSuite)
0054
0055 /// Unit test for creating compliant/non-compliant LineSurface object
0056 BOOST_AUTO_TEST_CASE(LineSurface_Constructors_test) {
0057   /// Test default construction
0058   // default construction is deleted
0059
0060   Translation3 translation{0., 1., 2.};
0061   Transform3 transform(translation);
0062   auto pTransform = Transform3(translation);
0063   const double radius = 2.;
0064   const double halfZ = 20.;
0065   BOOST_CHECK(LineSurfaceStub(pTransform, radius, halfZ).constructedOk());
0066
0067   /// ctor with nullptr for LineBounds
0068   BOOST_CHECK(LineSurfaceStub(pTransform).constructedOk());
0069
0070   /// ctor with LineBounds
0071   auto pLineBounds = std::make_shared<const LineBounds>(2., 10.);
0072   BOOST_CHECK(LineSurfaceStub(pTransform, pLineBounds).constructedOk());
0073
0074   /// ctor with LineBounds, detector element, Identifier
0075   auto pMaterial =
0076       std::make_shared<const HomogeneousSurfaceMaterial>(makePercentSlab());
0077   DetectorElementStub detElement{pTransform, pLineBounds, 0.2, pMaterial};
0078   BOOST_CHECK(LineSurfaceStub(pLineBounds, detElement).constructedOk());
0079   LineSurfaceStub lineToCopy(pTransform, 2., 20.);
0080
0081   /// Copy ctor
0082   BOOST_CHECK(LineSurfaceStub(lineToCopy).constructedOk());
0083
0084   /// Copied and transformed ctor
0085   BOOST_CHECK(
0086       LineSurfaceStub(tgContext, lineToCopy, transform).constructedOk());
0087
0088   /// Construct with nullptr bounds
0089   DetectorElementStub detElem;
0090   BOOST_CHECK_THROW(LineSurfaceStub nullBounds(nullptr, detElem),
0091                     AssertionFailureException);
0092
0093   BOOST_TEST_MESSAGE(
0094       "All LineSurface constructors are callable without problem");
0095 }
0096
0097 /// Unit tests of all named methods
0098 BOOST_AUTO_TEST_CASE(LineSurface_allNamedMethods_test) {
0099   // referencePosition()
0100   Translation3 translation{0., 1., 2.};
0101   Transform3 transform(translation);
0102   LineSurfaceStub line(transform, 2., 20.);
0103   Vector3 referencePosition{0., 1., 2.};
0104   CHECK_CLOSE_ABS(referencePosition,
0105                   line.referencePosition(tgContext, AxisDirection::AxisX),
0106                   1e-6);
0107
0108   // bounds()
0109   auto pLineBounds = std::make_shared<const LineBounds>(2., 10.);
0110   LineSurfaceStub boundedLine(transform, pLineBounds);
0111   const LineBounds& bounds =
0112       dynamic_cast<const LineBounds&>(boundedLine.bounds());
0113   BOOST_CHECK_EQUAL(bounds, LineBounds(2., 10.));
0114
0115   // globalToLocal()
0116   Vector3 gpos{0., 1., 0.};
0117   const Vector3 mom{20., 0., 0.};  // needs more realistic parameters
0118   Vector2 localPosition =
0119       line.globalToLocal(tgContext, gpos, mom.normalized()).value();
0120   const Vector2 expectedResult{0, -2};
0121   CHECK_CLOSE_ABS(expectedResult, localPosition, 1e-6);
0122
0123   // intersection
0124   {
0125     const Vector3 direction{0., 1., 2.};
0126     Intersection3D sfIntersection =
0127         line.intersect(tgContext, {0., 0., 0.}, direction.normalized(),
0128                        BoundaryTolerance::Infinite())
0129             .closest();
0130     BOOST_CHECK(sfIntersection.isValid());
0131     Vector3 expectedIntersection(0, 1., 2.);
0132     CHECK_CLOSE_ABS(sfIntersection.position(), expectedIntersection,
0133                     1e-6);  // need more tests..
0134   }
0135
0136   // isOnSurface
0137   const Vector3 insidePosition{0., 2.5, 0.};
0138   BOOST_CHECK(line.isOnSurface(
0139       tgContext, insidePosition, mom,
0140       BoundaryTolerance::Infinite()));  // need better test here
0141   const Vector3 outsidePosition{100., 100., 200.};
0142   BOOST_CHECK(!line.isOnSurface(tgContext, outsidePosition, mom,
0143                                 BoundaryTolerance::None()));
0144
0145   // localToGlobal
0146   Vector3 returnedGlobalPosition{0., 0., 0.};
0147   // Vector2 localPosition{0., 0.};
0148   const Vector3 momentum{300., 200., 0.};  // find better values!
0149   returnedGlobalPosition =
0150       line.localToGlobal(tgContext, localPosition, momentum.normalized());
0151   const Vector3 expectedGlobalPosition{0, 1, 0};
0152   CHECK_CLOSE_ABS(returnedGlobalPosition, expectedGlobalPosition, 1e-6);
0153
0154   // referenceFrame
0155   Vector3 globalPosition{0., 0., 0.};
0156   auto returnedRotationMatrix =
0157       line.referenceFrame(tgContext, globalPosition, momentum.normalized());
0158   double v0 = std::cos(std::atan(2. / 3.));
0159   double v1 = std::sin(std::atan(2. / 3.));
0160   RotationMatrix3 expectedRotationMatrix;
0161   expectedRotationMatrix << -v1, 0., v0, v0, 0., v1, 0., 1., -0.;
0162   CHECK_CLOSE_OR_SMALL(returnedRotationMatrix, expectedRotationMatrix, 1e-6,
0163                        1e-9);
0164
0165   // name()
0166   boost::test_tools::output_test_stream output;
0167   output << line.name();
0168   BOOST_CHECK(output.is_equal("Acts::LineSurface"));
0169
0170   // normal
0171   // arbitrary position, because should be irrelevant
0172   Vector3 position{5, 5, 5};  // should be irrelevant
0173   {
0174     Vector3 direction{1, 0, 0};
0175     CHECK_CLOSE_ABS(line.normal(tgContext, position, direction), direction,
0176                     1e-6);
0177   }
0178   {
0179     Vector3 direction = Vector3{1, 0, 0.1}.normalized();
0180     CHECK_CLOSE_ABS(line.normal(tgContext, position, direction),
0181                     Vector3::UnitX(), 1e-6);
0182   }
0183   {
0184     Vector3 direction{-1, 0, 0};
0185     CHECK_CLOSE_ABS(line.normal(tgContext, position, direction), direction,
0186                     1e-6);
0187   }
0188   {
0189     Vector3 direction{0, 1, 0};
0190     CHECK_CLOSE_ABS(line.normal(tgContext, position, direction), direction,
0191                     1e-6);
0192   }
0193
0194   // pathCorrection
0195   Vector3 any3DVector = Vector3::Random();
0196   CHECK_CLOSE_REL(line.pathCorrection(tgContext, any3DVector, any3DVector), 1.,
0197                   1e-6);
0198 }
0199
0200 /// Unit test for testing LineSurface assignment
0201 BOOST_AUTO_TEST_CASE(LineSurface_assignment_test) {
0202   Translation3 translation{0., 1., 2.};
0203   Transform3 transform(translation);
0204   LineSurfaceStub originalLine(transform, 2., 20.);
0205   LineSurfaceStub assignedLine(transform, 1., 1.);
0206   BOOST_CHECK(assignedLine != originalLine);  // operator != from base
0207   assignedLine = originalLine;
0208   BOOST_CHECK(assignedLine == originalLine);  // operator == from base
0209 }
0210
0211 /// Unit test for testing LineSurface alignment derivatives
0212 BOOST_AUTO_TEST_CASE(LineSurfaceAlignment) {
0213   Translation3 translation{0., 1., 2.};
0214   Transform3 transform(translation);
0215   LineSurfaceStub line(transform, 2., 20.);
0216
0217   const auto& rotation = transform.rotation();
0218   // The local frame z axis
0219   const Vector3 localZAxis = rotation.col(2);
0220   // Check the local z axis is aligned to global z axis
0221   CHECK_CLOSE_ABS(localZAxis, Vector3(0., 0., 1.), 1e-15);
0222
0223   // Define the track (global) position and direction
0224   Vector3 globalPosition{1, 2, 4};
0225   Vector3 momentum{-1, 1, 1};
0226   Vector3 direction = momentum.normalized();
0227
0228   // (a) Test the derivative of path length w.r.t. alignment parameters
0229   const AlignmentToPathMatrix& alignToPath =
0230       line.alignmentToPathDerivative(tgContext, globalPosition, direction);
0231   // The expected results
0232   AlignmentToPathMatrix expAlignToPath = AlignmentToPathMatrix::Zero();
0233   const double value = std::numbers::sqrt3 / 2;
0234   expAlignToPath << -value, value, 0, -3 * value, -value, 0;
0235   // Check if the calculated derivative is as expected
0236   CHECK_CLOSE_ABS(alignToPath, expAlignToPath, 1e-10);
0237
0238   // (b) Test the derivative of bound track parameters local position w.r.t.
0239   // position in local 3D Cartesian coordinates
0240   const auto& loc3DToLocBound =
0241       line.localCartesianToBoundLocalDerivative(tgContext, globalPosition);
0242   // Check if the result is as expected
0243   Matrix<2, 3> expLoc3DToLocBound = Matrix<2, 3>::Zero();
0244   expLoc3DToLocBound << 1 / std::numbers::sqrt2, 1 / std::numbers::sqrt2, 0, 0,
0245       0, 1;
0246   CHECK_CLOSE_ABS(loc3DToLocBound, expLoc3DToLocBound, 1e-10);
0247 }
0248
0249 BOOST_AUTO_TEST_CASE(LineSurfaceTransformRoundTrip) {
0250   LineSurfaceStub surface(Transform3::Identity());
0251
0252   auto roundTrip = [&surface](const Vector3& pos, const Vector3& dir) {
0253     Intersection3D intersection =
0254         surface.intersect(tgContext, pos, dir).closest();
0255     Vector3 global = intersection.position();
0256     Vector2 local = *surface.globalToLocal(tgContext, global, dir);
0257     Vector3 global2 = surface.localToGlobal(tgContext, local, dir);
0258     return std::make_tuple(global, local, global2);
0259   };
0260
0261   {
0262     Vector3 pos = {-0.02801, 0.00475611, 0.285106};
0263     Vector3 dir = Vector3(-0.03951, -0.221457, -0.564298).normalized();
0264
0265     auto [global, local, global2] = roundTrip(pos, dir);
0266
0267     CHECK_CLOSE_ABS(global, global2, 1e-10);
0268   }
0269
0270   {
0271     Vector3 pos = {-64.2892, 65.2697, -0.839014};
0272     Vector3 dir = Vector3(-0.236602, -0.157616, 0.956786).normalized();
0273
0274     auto [global, local, global2] = roundTrip(pos, dir);
0275
0276     CHECK_CLOSE_ABS(global, global2, 1e-10);
0277   }
0278 }
0279
0280 BOOST_AUTO_TEST_CASE(LineSurfaceTransformRoundTripEtaStability) {
0281   LineSurfaceStub surface(Transform3::Identity());
0282
0283   // eta=6 is already crashing
0284   const std::vector<double> etas = {0, 1, 2, 3, 4, 5};
0285
0286   for (double eta : etas) {
0287     Vector3 pca = {5, 0, 0};
0288     Vector3 dir = makeDirectionFromPhiEta(std::numbers::pi / 2., eta);
0289     Vector3 pos = pca + dir;
0290
0291     Intersection3D intersection =
0292         surface.intersect(tgContext, pos, dir).closest();
0293
0294     Vector3 global = intersection.position();
0295     Vector2 local = *surface.globalToLocal(tgContext, global, dir);
0296     Vector3 global2 = surface.localToGlobal(tgContext, local, dir);
0297
0298     CHECK_CLOSE_ABS(global, global2, 1e-10);
0299     CHECK_CLOSE_ABS(pca, global2, 1e-10);
0300   }
0301 }
0302
0303 BOOST_AUTO_TEST_CASE(LineSurfaceIntersection) {
0304   using namespace Acts::UnitLiterals;
0305
0306   double eps = 1e-10;
0307
0308   Vector3 direction = Vector3(1, 1, 100).normalized();
0309   BoundVector boundVector;
0310   boundVector << 1_cm, 1_cm, VectorHelpers::phi(direction),
0311       VectorHelpers::theta(direction), 1, 0;
0312   double pathLimit = 1_cm;
0313
0314   auto surface = std::make_shared<LineSurfaceStub>(Transform3::Identity());
0315
0316   BoundTrackParameters initialParams{surface, boundVector, std::nullopt,
0317                                      ParticleHypothesis::pion()};
0318
0319   using Propagator = Propagator<StraightLineStepper>;
0320   using PropagatorOptions = Propagator::Options<>;
0321
0322   Propagator propagator({});
0323
0324   BoundTrackParameters displacedParameters =
0325       BoundTrackParameters::createCurvilinear(Vector4::Zero(), Vector3::Zero(),
0326                                               1, std::nullopt,
0327                                               ParticleHypothesis::pion());
0328   {
0329     PropagatorOptions options(tgContext, {});
0330     options.direction = Acts::Direction::Backward();
0331     options.pathLimit = pathLimit;
0332
0333     auto result = propagator.propagate(initialParams, options);
0334     BOOST_CHECK(result.ok());
0335     BOOST_CHECK(result.value().endParameters);
0336
0337     displacedParameters = result.value().endParameters.value();
0338   }
0339
0340   Intersection3D intersection =
0341       surface
0342           ->intersect(tgContext, displacedParameters.position(tgContext),
0343                       displacedParameters.direction())
0344           .closest();
0345   CHECK_CLOSE_ABS(intersection.pathLength(), pathLimit, eps);
0346
0347   BoundTrackParameters endParameters{surface,
0348                                      detail::Test::someBoundParametersA(),
0349                                      std::nullopt, ParticleHypothesis::pion()};
0350   {
0351     PropagatorOptions options(tgContext, {});
0352     options.direction = Acts::Direction::Forward();
0353     options.stepping.maxStepSize = 1_mm;
0354
0355     auto result = propagator.propagate(displacedParameters, *surface, options);
0356     BOOST_CHECK(result.ok());
0357     BOOST_CHECK(result.value().endParameters);
0358     CHECK_CLOSE_ABS(result.value().pathLength, pathLimit, eps);
0359     endParameters = result.value().endParameters.value();
0360   }
0361
0362   CHECK_CLOSE_ABS(initialParams.parameters(), endParameters.parameters(), eps);
0363 }
0364
0365 BOOST_AUTO_TEST_CASE(LineSurfaceMaterialAssignment) {
0366   LineSurfaceStub surface(Transform3::Identity(), 1., 100.);
0367   MaterialSlab slab(Material::fromMolarDensity(1., 2., 3., 4., 5.), 0.1);
0368
0369   // HomogeneousSurfaceMaterial is always valid
0370   auto homMat = std::make_shared<HomogeneousSurfaceMaterial>(slab);
0371   BOOST_CHECK_NO_THROW(surface.assignSurfaceMaterial(homMat));
0372   BOOST_CHECK_NE(surface.surfaceMaterial(), nullptr);
0373
0374   // nullptr clears the material
0375   BOOST_CHECK_NO_THROW(surface.assignSurfaceMaterial(nullptr));
0376   BOOST_CHECK_EQUAL(surface.surfaceMaterial(), nullptr);
0377
0378   // 1D AxisR - valid for line surface
0379   BinUtility buR(10, 0.f, 1.f, Acts::open, AxisDirection::AxisR);
0380   auto matR = std::make_shared<BinnedSurfaceMaterial>(
0381       buR, MaterialSlabVector(10, slab));
0382   BOOST_CHECK_NO_THROW(surface.assignSurfaceMaterial(matR));
0383
0384   // 1D AxisZ - valid for line surface
0385   BinUtility buZ(10, -100.f, 100.f, Acts::open, AxisDirection::AxisZ);
0386   auto matZ = std::make_shared<BinnedSurfaceMaterial>(
0387       buZ, MaterialSlabVector(10, slab));
0388   BOOST_CHECK_NO_THROW(surface.assignSurfaceMaterial(matZ));
0389
0390   // 2D {AxisR, AxisZ} - valid for line surface
0391   BinUtility bu2D(5, 0.f, 1.f, Acts::open, AxisDirection::AxisR);
0392   bu2D += BinUtility(3, -100.f, 100.f, Acts::open, AxisDirection::AxisZ);
0393   auto mat2D = std::make_shared<BinnedSurfaceMaterial>(
0394       bu2D, MaterialSlabMatrix(3, MaterialSlabVector(5, slab)));
0395   BOOST_CHECK_NO_THROW(surface.assignSurfaceMaterial(mat2D));
0396
0397   // Wrong axis directions - should throw
0398   for (auto badDir : {AxisDirection::AxisRPhi, AxisDirection::AxisPhi,
0399                       AxisDirection::AxisX, AxisDirection::AxisY}) {
0400     BinUtility buBad(10, 0.f, 10.f, Acts::open, badDir);
0401     auto matBad = std::make_shared<BinnedSurfaceMaterial>(
0402         buBad, MaterialSlabVector(10, slab));
0403     BOOST_CHECK_THROW(surface.assignSurfaceMaterial(matBad),
0404                       std::invalid_argument);
0405   }
0406 }
0407
0408 BOOST_AUTO_TEST_SUITE_END()
0409
0410 }  // namespace ActsTests