Core/Surfaces/DiscSurfaceTests.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/data/test_case.hpp>
0010 #include <boost/test/tools/output_test_stream.hpp>
0011 #include <boost/test/unit_test.hpp>
0012
0013 #include "Acts/Definitions/Algebra.hpp"
0014 #include "Acts/Definitions/Alignment.hpp"
0015 #include "Acts/Definitions/Tolerance.hpp"
0016 #include "Acts/Definitions/Units.hpp"
0017 #include "Acts/Geometry/Extent.hpp"
0018 #include "Acts/Geometry/GeometryContext.hpp"
0019 #include "Acts/Geometry/Polyhedron.hpp"
0020 #include "Acts/Surfaces/AnnulusBounds.hpp"
0021 #include "Acts/Surfaces/DiscSurface.hpp"
0022 #include "Acts/Surfaces/RadialBounds.hpp"
0023 #include "Acts/Surfaces/Surface.hpp"
0024 #include "Acts/Surfaces/SurfaceBounds.hpp"
0025 #include "Acts/Surfaces/SurfaceMergingException.hpp"
0026 #include "Acts/Tests/CommonHelpers/DetectorElementStub.hpp"
0027 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0028 #include "Acts/Utilities/Intersection.hpp"
0029 #include "Acts/Utilities/Result.hpp"
0030 #include "Acts/Utilities/ThrowAssert.hpp"
0031 #include "Acts/Utilities/detail/periodic.hpp"
0032
0033 #include <cmath>
0034 #include <memory>
0035 #include <numbers>
0036 #include <string>
0037
0038 using namespace Acts::UnitLiterals;
0039
0040 namespace Acts::Test {
0041 // Create a test context
0042 GeometryContext tgContext = GeometryContext();
0043 auto logger = Acts::getDefaultLogger("UnitTests", Acts::Logging::VERBOSE);
0044
0045 BOOST_AUTO_TEST_SUITE(Surfaces)
0046 /// Unit tests for creating DiscSurface object
0047 BOOST_AUTO_TEST_CASE(DiscSurfaceConstruction) {
0048   /// Test default construction
0049   // default construction is deleted
0050
0051   const double rMin = 1.;
0052   const double rMax = 5.;
0053   const double halfPhiSector = std::numbers::pi / 8.;
0054
0055   /// Test DiscSurface constructor with default halfPhiSector
0056   BOOST_CHECK_NO_THROW(
0057       Surface::makeShared<DiscSurface>(Transform3::Identity(), rMin, rMax));
0058
0059   /// Test DiscSurface constructor with a transform specified
0060   Translation3 translation{0., 1., 2.};
0061   auto pTransform = Transform3(translation);
0062   BOOST_CHECK_NO_THROW(
0063       Surface::makeShared<DiscSurface>(pTransform, rMin, rMax, halfPhiSector));
0064
0065   /// Copy constructed DiscSurface
0066   auto anotherDiscSurface =
0067       Surface::makeShared<DiscSurface>(pTransform, rMin, rMax, halfPhiSector);
0068   // N.B. Just using
0069   // BOOST_CHECK_NO_THROW(Surface::makeShared<DiscSurface>(anotherDiscSurface))
0070   // tries to call the (deleted) default constructor.
0071   auto copiedSurface = Surface::makeShared<DiscSurface>(*anotherDiscSurface);
0072   BOOST_TEST_MESSAGE("Copy constructed DiscSurface ok");
0073
0074   /// Copied and transformed DiscSurface
0075   BOOST_CHECK_NO_THROW(Surface::makeShared<DiscSurface>(
0076       tgContext, *anotherDiscSurface, pTransform));
0077
0078   /// Construct with nullptr bounds
0079   DetectorElementStub detElem;
0080   BOOST_CHECK_THROW(
0081       auto nullBounds = Surface::makeShared<DiscSurface>(nullptr, detElem),
0082       AssertionFailureException);
0083 }
0084
0085 /// Unit tests of all named methods
0086 BOOST_AUTO_TEST_CASE(DiscSurfaceProperties) {
0087   const double rMin = 1.;
0088   const double rMax = 5.;
0089   const double halfPhiSector = std::numbers::pi / 8.;
0090
0091   const Vector3 origin3D{0, 0, 0};
0092
0093   auto discSurfaceObject = Surface::makeShared<DiscSurface>(
0094       Transform3::Identity(), rMin, rMax, halfPhiSector);
0095
0096   /// Test type
0097   BOOST_CHECK_EQUAL(discSurfaceObject->type(), Surface::Disc);
0098
0099   /// Test normal, no local position specified
0100   Vector3 zAxis{0, 0, 1};
0101   BOOST_CHECK_EQUAL(discSurfaceObject->normal(tgContext), zAxis);
0102
0103   /// Test normal, local position specified
0104   Vector2 lpos(2., 0.05);
0105   BOOST_CHECK_EQUAL(discSurfaceObject->normal(tgContext, lpos), zAxis);
0106
0107   /// Test referencePosition
0108   BOOST_CHECK_EQUAL(
0109       discSurfaceObject->referencePosition(tgContext, AxisDirection::AxisRPhi),
0110       origin3D);
0111
0112   /// Test bounds
0113   BOOST_CHECK_EQUAL(discSurfaceObject->bounds().type(), SurfaceBounds::eDisc);
0114
0115   Vector3 ignoredMomentum{0., 0., 0.};
0116   /// Test isOnSurface()
0117   Vector3 point3DNotInSector{0., 1.2, 0};
0118   Vector3 point3DOnSurface{1.2, 0., 0};
0119   BOOST_CHECK(!discSurfaceObject->isOnSurface(tgContext, point3DNotInSector,
0120                                               ignoredMomentum,
0121                                               BoundaryTolerance::None()));
0122   BOOST_CHECK(!discSurfaceObject->isOnSurface(tgContext, point3DNotInSector,
0123                                               BoundaryTolerance::None()));
0124   BOOST_CHECK(discSurfaceObject->isOnSurface(
0125       tgContext, point3DOnSurface, ignoredMomentum, BoundaryTolerance::None()));
0126   BOOST_CHECK(discSurfaceObject->isOnSurface(tgContext, point3DOnSurface,
0127                                              BoundaryTolerance::None()));
0128
0129   /// Test localToGlobal
0130   Vector3 returnedPosition{10.9, 8.7, 6.5};
0131   Vector3 expectedPosition{1.2, 0, 0};
0132   Vector2 rPhiOnDisc{1.2, 0.};
0133   Vector2 rPhiNotInSector{
0134       1.2, std::numbers::pi};  // outside sector at Phi=0, +/- pi/8
0135   returnedPosition =
0136       discSurfaceObject->localToGlobal(tgContext, rPhiOnDisc, ignoredMomentum);
0137   CHECK_CLOSE_ABS(returnedPosition, expectedPosition, 1e-6);
0138
0139   returnedPosition = discSurfaceObject->localToGlobal(
0140       tgContext, rPhiNotInSector, ignoredMomentum);
0141   Vector3 expectedNonPosition{-1.2, 0, 0};
0142   CHECK_CLOSE_ABS(returnedPosition, expectedNonPosition, 1e-6);
0143
0144   /// Test globalToLocal
0145   Vector2 returnedLocalPosition{33., 44.};
0146   Vector2 expectedLocalPosition{1.2, 0.};
0147   returnedLocalPosition =
0148       discSurfaceObject
0149           ->globalToLocal(tgContext, point3DOnSurface, ignoredMomentum)
0150           .value();
0151   CHECK_CLOSE_ABS(returnedLocalPosition, expectedLocalPosition, 1e-6);
0152
0153   // Global to local does not check inside bounds
0154   returnedLocalPosition =
0155       discSurfaceObject
0156           ->globalToLocal(tgContext, point3DNotInSector, ignoredMomentum)
0157           .value();
0158
0159   Vector3 pointOutsideR{0., 100., 0};
0160   returnedLocalPosition =
0161       discSurfaceObject
0162           ->globalToLocal(tgContext, pointOutsideR, ignoredMomentum)
0163           .value();
0164
0165   /// Test localPolarToCartesian
0166   Vector2 rPhi1_1{std::numbers::sqrt2, std::numbers::pi / 4.};
0167   Vector2 cartesian1_1{1., 1.};
0168   CHECK_CLOSE_REL(discSurfaceObject->localPolarToCartesian(rPhi1_1),
0169                   cartesian1_1, 1e-6);
0170
0171   /// Test localCartesianToPolar
0172   CHECK_CLOSE_REL(discSurfaceObject->localCartesianToPolar(cartesian1_1),
0173                   rPhi1_1, 1e-6);
0174
0175   /// Test localPolarToLocalCartesian
0176   CHECK_CLOSE_REL(discSurfaceObject->localPolarToLocalCartesian(rPhi1_1),
0177                   cartesian1_1, 1e-6);
0178
0179   /// Test localCartesianToGlobal
0180   Vector3 cartesian3D1_1{1., 1., 0.};
0181   CHECK_CLOSE_ABS(
0182       discSurfaceObject->localCartesianToGlobal(tgContext, cartesian1_1),
0183       cartesian3D1_1, 1e-6);
0184
0185   /// Test globalToLocalCartesian
0186   CHECK_CLOSE_REL(
0187       discSurfaceObject->globalToLocalCartesian(tgContext, cartesian3D1_1),
0188       cartesian1_1, 1e-6);
0189
0190   /// Test pathCorrection
0191   double projected3DMomentum = std::numbers::sqrt3 * 1.e6;
0192   Vector3 momentum{projected3DMomentum, projected3DMomentum,
0193                    projected3DMomentum};
0194   Vector3 ignoredPosition = discSurfaceObject->center(tgContext);
0195   CHECK_CLOSE_REL(discSurfaceObject->pathCorrection(tgContext, ignoredPosition,
0196                                                     momentum.normalized()),
0197                   std::numbers::sqrt3, 0.01);
0198
0199   /// intersection test
0200   Vector3 globalPosition{1.2, 0., -10.};
0201   Vector3 direction{0., 0., 1.};  // must be normalised
0202   Vector3 expected{1.2, 0., 0.};
0203
0204   // intersect is a struct of (Vector3) position, pathLength, distance and
0205   // (bool) valid, it's contained in a Surface intersection
0206   auto sfIntersection = discSurfaceObject
0207                             ->intersect(tgContext, globalPosition, direction,
0208                                         BoundaryTolerance::Infinite())
0209                             .closest();
0210   Intersection3D expectedIntersect{Vector3{1.2, 0., 0.}, 10.,
0211                                    IntersectionStatus::reachable};
0212   BOOST_CHECK(sfIntersection.isValid());
0213   CHECK_CLOSE_ABS(sfIntersection.position(), expectedIntersect.position(),
0214                   1e-9);
0215   CHECK_CLOSE_ABS(sfIntersection.pathLength(), expectedIntersect.pathLength(),
0216                   1e-9);
0217   BOOST_CHECK_EQUAL(&sfIntersection.surface(), discSurfaceObject.get());
0218
0219   /// Test name
0220   boost::test_tools::output_test_stream nameOuput;
0221   nameOuput << discSurfaceObject->name();
0222   BOOST_CHECK(nameOuput.is_equal("Acts::DiscSurface"));
0223 }
0224
0225 /// Unit test for testing DiscSurface assignment and equality
0226 BOOST_AUTO_TEST_CASE(DiscSurfaceAssignment) {
0227   const double rMin = 1.;
0228   const double rMax = 5.;
0229   const double halfPhiSector = std::numbers::pi / 8.;
0230
0231   auto discSurfaceObject = Surface::makeShared<DiscSurface>(
0232       Transform3::Identity(), rMin, rMax, halfPhiSector);
0233   auto assignedDisc =
0234       Surface::makeShared<DiscSurface>(Transform3::Identity(), 2.2, 4.4, 0.07);
0235
0236   BOOST_CHECK_NO_THROW(*assignedDisc = *discSurfaceObject);
0237   BOOST_CHECK((*assignedDisc) == (*discSurfaceObject));
0238 }
0239
0240 /// Unit test for testing DiscSurface assignment and equality
0241 BOOST_AUTO_TEST_CASE(DiscSurfaceExtent) {
0242   const double rMin = 1.;
0243   const double rMax = 5.;
0244
0245   auto pDisc =
0246       Surface::makeShared<DiscSurface>(Transform3::Identity(), 0., rMax);
0247   auto pDiscExtent = pDisc->polyhedronRepresentation(tgContext, 1).extent();
0248
0249   CHECK_CLOSE_ABS(0., pDiscExtent.min(AxisDirection::AxisZ),
0250                   s_onSurfaceTolerance);
0251   CHECK_CLOSE_ABS(0., pDiscExtent.max(AxisDirection::AxisZ),
0252                   s_onSurfaceTolerance);
0253   CHECK_CLOSE_ABS(0., pDiscExtent.min(AxisDirection::AxisR),
0254                   s_onSurfaceTolerance);
0255   CHECK_CLOSE_ABS(rMax, pDiscExtent.max(AxisDirection::AxisR),
0256                   s_onSurfaceTolerance);
0257   CHECK_CLOSE_ABS(-rMax, pDiscExtent.min(AxisDirection::AxisX),
0258                   s_onSurfaceTolerance);
0259   CHECK_CLOSE_ABS(rMax, pDiscExtent.max(AxisDirection::AxisX),
0260                   s_onSurfaceTolerance);
0261   CHECK_CLOSE_ABS(-rMax, pDiscExtent.min(AxisDirection::AxisY),
0262                   s_onSurfaceTolerance);
0263   CHECK_CLOSE_ABS(rMax, pDiscExtent.max(AxisDirection::AxisY),
0264                   s_onSurfaceTolerance);
0265   CHECK_CLOSE_ABS(-std::numbers::pi, pDiscExtent.min(AxisDirection::AxisPhi),
0266                   s_onSurfaceTolerance);
0267   CHECK_CLOSE_ABS(std::numbers::pi, pDiscExtent.max(AxisDirection::AxisPhi),
0268                   s_onSurfaceTolerance);
0269
0270   auto pRing =
0271       Surface::makeShared<DiscSurface>(Transform3::Identity(), rMin, rMax);
0272   auto pRingExtent = pRing->polyhedronRepresentation(tgContext, 1).extent();
0273
0274   CHECK_CLOSE_ABS(0., pRingExtent.min(AxisDirection::AxisZ),
0275                   s_onSurfaceTolerance);
0276   CHECK_CLOSE_ABS(0., pRingExtent.max(AxisDirection::AxisZ),
0277                   s_onSurfaceTolerance);
0278   CHECK_CLOSE_ABS(rMin, pRingExtent.min(AxisDirection::AxisR),
0279                   s_onSurfaceTolerance);
0280   CHECK_CLOSE_ABS(rMax, pRingExtent.max(AxisDirection::AxisR),
0281                   s_onSurfaceTolerance);
0282   CHECK_CLOSE_ABS(-rMax, pRingExtent.min(AxisDirection::AxisX),
0283                   s_onSurfaceTolerance);
0284   CHECK_CLOSE_ABS(rMax, pRingExtent.max(AxisDirection::AxisX),
0285                   s_onSurfaceTolerance);
0286   CHECK_CLOSE_ABS(-rMax, pRingExtent.min(AxisDirection::AxisY),
0287                   s_onSurfaceTolerance);
0288   CHECK_CLOSE_ABS(rMax, pRingExtent.max(AxisDirection::AxisY),
0289                   s_onSurfaceTolerance);
0290 }
0291
0292 /// Unit test for testing DiscSurface alignment derivatives
0293 BOOST_AUTO_TEST_CASE(DiscSurfaceAlignment) {
0294   Translation3 translation{0., 1., 2.};
0295   Transform3 transform(translation);
0296   const double rMin = 1.;
0297   const double rMax = 5.;
0298   const double halfPhiSector = std::numbers::pi / 8.;
0299
0300   auto discSurfaceObject =
0301       Surface::makeShared<DiscSurface>(transform, rMin, rMax, halfPhiSector);
0302
0303   const auto& rotation = transform.rotation();
0304   // The local frame z axis
0305   const Vector3 localZAxis = rotation.col(2);
0306   // Check the local z axis is aligned to global z axis
0307   CHECK_CLOSE_ABS(localZAxis, Vector3(0., 0., 1.), 1e-15);
0308
0309   // Define the track (global) position and direction
0310   Vector3 globalPosition{0, 4, 2};
0311   Vector3 momentum{0, 0, 1};
0312   Vector3 direction = momentum.normalized();
0313
0314   // (a) Test the derivative of path length w.r.t. alignment parameters
0315   const AlignmentToPathMatrix& alignToPath =
0316       discSurfaceObject->alignmentToPathDerivative(tgContext, globalPosition,
0317                                                    direction);
0318   // The expected results
0319   AlignmentToPathMatrix expAlignToPath = AlignmentToPathMatrix::Zero();
0320   expAlignToPath << 0, 0, 1, 3, 0, 0;
0321   // Check if the calculated derivative is as expected
0322   CHECK_CLOSE_ABS(alignToPath, expAlignToPath, 1e-10);
0323
0324   // (b) Test the derivative of bound track parameters local position w.r.t.
0325   // position in local 3D Cartesian coordinates
0326   const auto& loc3DToLocBound =
0327       discSurfaceObject->localCartesianToBoundLocalDerivative(tgContext,
0328                                                               globalPosition);
0329   // Check if the result is as expected
0330   ActsMatrix<2, 3> expLoc3DToLocBound = ActsMatrix<2, 3>::Zero();
0331   expLoc3DToLocBound << 0, 1, 0, -1. / 3, 0, 0;
0332   CHECK_CLOSE_ABS(loc3DToLocBound, expLoc3DToLocBound, 1e-10);
0333 }
0334
0335 BOOST_AUTO_TEST_CASE(DiscSurfaceBinningPosition) {
0336   using namespace Acts::UnitLiterals;
0337   Vector3 s{5_mm, 7_mm, 10_cm};
0338   Transform3 trf;
0339   trf = Translation3(s) * AngleAxis3{0.5, Vector3::UnitZ()};
0340
0341   double minR = 300;
0342   double maxR = 330;
0343
0344   {
0345     // Radial Bounds
0346     auto bounds =
0347         std::make_shared<RadialBounds>(minR, maxR, std::numbers::pi / 8, 0.1);
0348     auto disc = Acts::Surface::makeShared<Acts::DiscSurface>(trf, bounds);
0349
0350     Vector3 bp = disc->referencePosition(tgContext, AxisDirection::AxisR);
0351     double r = (bounds->rMax() + bounds->rMin()) / 2.0;
0352     double phi = bounds->get(RadialBounds::eAveragePhi);
0353     Vector3 exp = Vector3{r * std::cos(phi), r * std::sin(phi), 0};
0354     exp = trf * exp;
0355
0356     BOOST_CHECK_EQUAL(bp, exp);
0357     BOOST_CHECK_EQUAL(
0358         disc->referencePositionValue(tgContext, AxisDirection::AxisR),
0359         VectorHelpers::perp(exp));
0360
0361     bp = disc->referencePosition(tgContext, AxisDirection::AxisPhi);
0362     BOOST_CHECK_EQUAL(bp, exp);
0363     BOOST_CHECK_EQUAL(
0364         disc->referencePositionValue(tgContext, AxisDirection::AxisPhi),
0365         VectorHelpers::phi(exp));
0366
0367     for (auto b :
0368          {AxisDirection::AxisX, AxisDirection::AxisY, AxisDirection::AxisZ,
0369           AxisDirection::AxisEta, AxisDirection::AxisRPhi,
0370           AxisDirection::AxisTheta, AxisDirection::AxisMag}) {
0371       BOOST_TEST_CONTEXT("binValue: " << b) {
0372         BOOST_CHECK_EQUAL(disc->referencePosition(tgContext, b),
0373                           disc->center(tgContext));
0374       }
0375     }
0376   }
0377
0378   {
0379     // Annulus Bounds
0380     double minPhiRel = -0.3;
0381     double maxPhiRel = 0.2;
0382     Vector2 origin{5_mm, 5_mm};
0383     auto bounds = std::make_shared<AnnulusBounds>(minR, maxR, minPhiRel,
0384                                                   maxPhiRel, origin);
0385
0386     auto disc = Acts::Surface::makeShared<Acts::DiscSurface>(trf, bounds);
0387
0388     Vector3 bp = disc->referencePosition(tgContext, AxisDirection::AxisR);
0389     double r = (bounds->rMax() + bounds->rMin()) / 2.0;
0390     double phi = bounds->get(AnnulusBounds::eAveragePhi);
0391     Vector3 exp = Vector3{r * std::cos(phi), r * std::sin(phi), 0};
0392     exp = trf * exp;
0393
0394     BOOST_CHECK_EQUAL(bp, exp);
0395
0396     bp = disc->referencePosition(tgContext, AxisDirection::AxisPhi);
0397     BOOST_CHECK_EQUAL(bp, exp);
0398
0399     for (auto b :
0400          {AxisDirection::AxisX, AxisDirection::AxisY, AxisDirection::AxisZ,
0401           AxisDirection::AxisEta, AxisDirection::AxisRPhi,
0402           AxisDirection::AxisTheta, AxisDirection::AxisMag}) {
0403       BOOST_TEST_CONTEXT("binValue: " << b) {
0404         BOOST_CHECK_EQUAL(disc->referencePosition(tgContext, b),
0405                           disc->center(tgContext));
0406       }
0407     }
0408   }
0409 }
0410
0411 BOOST_AUTO_TEST_SUITE(DiscSurfaceMerging)
0412
0413 namespace {
0414 std::shared_ptr<DiscSurface> makeDisc(const Transform3& transform, double rMin,
0415                                       double rMax,
0416                                       double halfPhi = std::numbers::pi,
0417                                       double avgPhi = 0) {
0418   return Surface::makeShared<DiscSurface>(
0419       transform, std::make_shared<RadialBounds>(rMin, rMax, halfPhi, avgPhi));
0420 }
0421
0422 }  // namespace
0423
0424 BOOST_AUTO_TEST_CASE(IncompatibleBounds) {
0425   Logging::ScopedFailureThreshold ft{Logging::FATAL};
0426   Transform3 base = Transform3::Identity();
0427   auto discRadial = makeDisc(base, 30_mm, 100_mm);
0428   auto discTrap =
0429       Surface::makeShared<DiscSurface>(base, 20_mm, 40_mm, 100_mm, 150_mm);
0430   auto discTrap2 =
0431       Surface::makeShared<DiscSurface>(base, 20_mm, 40_mm, 30_mm, 100_mm);
0432
0433   BOOST_CHECK_THROW(
0434       discRadial->mergedWith(*discTrap, AxisDirection::AxisR, false, *logger),
0435
0436       SurfaceMergingException);
0437
0438   BOOST_CHECK_THROW(
0439       discTrap2->mergedWith(*discTrap, AxisDirection::AxisR, false, *logger),
0440       SurfaceMergingException);
0441 }
0442
0443 BOOST_AUTO_TEST_CASE(InvalidDetectorElement) {
0444   DetectorElementStub detElem;
0445
0446   auto bounds1 = std::make_shared<RadialBounds>(30_mm, 100_mm);
0447   auto disc1 = Surface::makeShared<DiscSurface>(bounds1, detElem);
0448
0449   auto bounds2 = std::make_shared<RadialBounds>(100_mm, 150_mm);
0450   auto disc2 = Surface::makeShared<DiscSurface>(bounds2, detElem);
0451
0452   BOOST_CHECK_THROW(
0453       disc1->mergedWith(*disc2, AxisDirection::AxisR, false, *logger),
0454       SurfaceMergingException);
0455 }
0456
0457 BOOST_DATA_TEST_CASE(IncompatibleRDirection,
0458                      (boost::unit_test::data::xrange(-135, 180, 45) *
0459                       boost::unit_test::data::make(Vector3{0_mm, 0_mm, 0_mm},
0460                                                    Vector3{20_mm, 0_mm, 0_mm},
0461                                                    Vector3{0_mm, 20_mm, 0_mm},
0462                                                    Vector3{20_mm, 20_mm, 0_mm},
0463                                                    Vector3{0_mm, 0_mm, 20_mm})),
0464                      angle, offset) {
0465   Logging::ScopedFailureThreshold ft{Logging::FATAL};
0466
0467   Transform3 base =
0468       AngleAxis3(angle * 1_degree, Vector3::UnitX()) * Translation3(offset);
0469
0470   auto disc = makeDisc(base, 30_mm, 100_mm);
0471
0472   // Disc with overlap in r
0473   auto discOverlap = makeDisc(base, 90_mm, 150_mm);
0474   BOOST_CHECK_THROW(disc->mergedWith(*discOverlap, Acts::AxisDirection::AxisR,
0475                                      false, *logger),
0476                     SurfaceMergingException);
0477
0478   // Disc with gap in r
0479   auto discGap = makeDisc(base, 110_mm, 150_mm);
0480   BOOST_CHECK_THROW(
0481       disc->mergedWith(*discGap, Acts::AxisDirection::AxisR, false, *logger),
0482       SurfaceMergingException);
0483
0484   auto discShiftedZ = Surface::makeShared<DiscSurface>(
0485       base * Translation3{Vector3::UnitZ() * 10_mm}, 100_mm, 150_mm);
0486   BOOST_CHECK_THROW(disc->mergedWith(*discShiftedZ, Acts::AxisDirection::AxisR,
0487                                      false, *logger),
0488                     SurfaceMergingException);
0489
0490   auto discShiftedXy = makeDisc(
0491       base * Translation3{Vector3{1_mm, 2_mm, 200_mm}}, 100_mm, 150_mm);
0492   BOOST_CHECK_THROW(disc->mergedWith(*discShiftedXy, Acts::AxisDirection::AxisZ,
0493                                      false, *logger),
0494                     SurfaceMergingException);
0495
0496   auto discRotatedZ = makeDisc(base * AngleAxis3{10_degree, Vector3::UnitZ()},
0497                                100_mm, 150_mm, 60_degree, 0_degree);
0498   BOOST_CHECK_THROW(disc->mergedWith(*discRotatedZ, Acts::AxisDirection::AxisR,
0499                                      false, *logger),
0500                     SurfaceMergingException);
0501
0502   auto discRotatedX =
0503       makeDisc(base * AngleAxis3{10_degree, Vector3::UnitX()}, 100_mm, 150_mm);
0504   BOOST_CHECK_THROW(disc->mergedWith(*discRotatedX, Acts::AxisDirection::AxisR,
0505                                      false, *logger),
0506                     SurfaceMergingException);
0507
0508   // Test not same phi sector
0509   auto discPhi1 = makeDisc(base, 30_mm, 100_mm, 10_degree, 40_degree);
0510   auto discPhi2 = makeDisc(base, 100_mm, 160_mm, 20_degree, 40_degree);
0511   auto discPhi3 = makeDisc(base, 100_mm, 160_mm, 10_degree, 50_degree);
0512   BOOST_CHECK_THROW(
0513       discPhi1->mergedWith(*discPhi2, AxisDirection::AxisR, false, *logger),
0514       SurfaceMergingException);
0515
0516   BOOST_CHECK_THROW(
0517       discPhi1->mergedWith(*discPhi3, AxisDirection::AxisR, false, *logger),
0518       SurfaceMergingException);
0519 }
0520
0521 BOOST_DATA_TEST_CASE(RDirection,
0522                      (boost::unit_test::data::xrange(-135, 180, 45) *
0523                       boost::unit_test::data::make(Vector3{0_mm, 0_mm, 0_mm},
0524                                                    Vector3{20_mm, 0_mm, 0_mm},
0525                                                    Vector3{0_mm, 20_mm, 0_mm},
0526                                                    Vector3{20_mm, 20_mm, 0_mm},
0527                                                    Vector3{0_mm, 0_mm, 20_mm})),
0528                      angle, offset) {
0529   Transform3 base =
0530       AngleAxis3(angle * 1_degree, Vector3::UnitX()) * Translation3(offset);
0531
0532   auto disc = makeDisc(base, 30_mm, 100_mm);
0533
0534   auto disc2 =
0535       makeDisc(base * AngleAxis3(14_degree, Vector3::UnitZ()), 100_mm, 150_mm);
0536
0537   auto [disc3, reversed] =
0538       disc->mergedWith(*disc2, Acts::AxisDirection::AxisR, false, *logger);
0539   BOOST_REQUIRE_NE(disc3, nullptr);
0540   BOOST_CHECK(!reversed);
0541
0542   auto [disc3Reversed, reversed2] =
0543       disc2->mergedWith(*disc, Acts::AxisDirection::AxisR, false, *logger);
0544   BOOST_REQUIRE_NE(disc3Reversed, nullptr);
0545   BOOST_CHECK(disc3->bounds() == disc3Reversed->bounds());
0546   BOOST_CHECK(reversed2);
0547
0548   const auto* bounds = dynamic_cast<const RadialBounds*>(&disc3->bounds());
0549   BOOST_REQUIRE_NE(bounds, nullptr);
0550
0551   BOOST_CHECK_EQUAL(bounds->get(RadialBounds::eMinR), 30_mm);
0552   BOOST_CHECK_EQUAL(bounds->get(RadialBounds::eMaxR), 150_mm);
0553
0554   // Disc did not move
0555   BOOST_CHECK_EQUAL(base.matrix(), disc3->transform(tgContext).matrix());
0556
0557   // Rotation in z depends on the ordering, the left side "wins"
0558   Transform3 expected12 = base;
0559   BOOST_CHECK_EQUAL(expected12.matrix(), disc3->transform(tgContext).matrix());
0560
0561   Transform3 expected21 = base * AngleAxis3(14_degree, Vector3::UnitZ());
0562   CHECK_CLOSE_OR_SMALL(disc3Reversed->transform(tgContext).matrix(),
0563                        expected21.matrix(), 1e-6, 1e-10);
0564
0565   // Test r merging with phi sectors (matching)
0566   auto discPhi1 = makeDisc(base, 30_mm, 100_mm, 10_degree, 40_degree);
0567   auto discPhi2 = makeDisc(base, 100_mm, 160_mm, 10_degree, 40_degree);
0568   auto [discPhi12, reversedPhi12] =
0569       discPhi1->mergedWith(*discPhi2, AxisDirection::AxisR, false, *logger);
0570   BOOST_REQUIRE_NE(discPhi12, nullptr);
0571
0572   const auto* boundsPhi12 =
0573       dynamic_cast<const RadialBounds*>(&discPhi12->bounds());
0574   BOOST_REQUIRE_NE(boundsPhi12, nullptr);
0575
0576   BOOST_CHECK_EQUAL(boundsPhi12->get(RadialBounds::eMinR), 30_mm);
0577   BOOST_CHECK_EQUAL(boundsPhi12->get(RadialBounds::eMaxR), 160_mm);
0578   BOOST_CHECK_EQUAL(boundsPhi12->get(RadialBounds::eHalfPhiSector), 10_degree);
0579 }
0580
0581 BOOST_DATA_TEST_CASE(IncompatiblePhiDirection,
0582                      (boost::unit_test::data::xrange(-135, 180, 45) *
0583                       boost::unit_test::data::make(Vector3{0_mm, 0_mm, 0_mm},
0584                                                    Vector3{20_mm, 0_mm, 0_mm},
0585                                                    Vector3{0_mm, 20_mm, 0_mm},
0586                                                    Vector3{20_mm, 20_mm, 0_mm},
0587                                                    Vector3{0_mm, 0_mm, 20_mm}) *
0588                       boost::unit_test::data::xrange(-1300, 1300, 104)),
0589                      angle, offset, phiShift) {
0590   Logging::ScopedFailureThreshold ft{Logging::FATAL};
0591   Transform3 base =
0592       AngleAxis3(angle * 1_degree, Vector3::UnitX()) * Translation3(offset);
0593
0594   auto a = [phiShift](double v) {
0595     return detail::radian_sym(v + phiShift * 1_degree);
0596   };
0597
0598   auto discPhi = makeDisc(base, 30_mm, 100_mm, 10_degree, a(40_degree));
0599
0600   // Disc with overlap in phi
0601   auto discPhi2 = makeDisc(base, 30_mm, 100_mm, 45_degree, a(85_degree));
0602   BOOST_CHECK_THROW(discPhi->mergedWith(*discPhi2, Acts::AxisDirection::AxisPhi,
0603                                         false, *logger),
0604                     SurfaceMergingException);
0605
0606   // Disc with gap in phi
0607   auto discPhi3 = makeDisc(base, 30_mm, 100_mm, 45_degree, a(105_degree));
0608   BOOST_CHECK_THROW(discPhi->mergedWith(*discPhi3, Acts::AxisDirection::AxisPhi,
0609                                         false, *logger),
0610                     SurfaceMergingException);
0611
0612   // Disc with a z shift
0613   auto discPhi4 = makeDisc(base * Translation3{Vector3::UnitZ() * 20_mm}, 30_mm,
0614                            100_mm, 45_degree, a(95_degree));
0615   BOOST_CHECK_THROW(discPhi->mergedWith(*discPhi4, Acts::AxisDirection::AxisPhi,
0616                                         false, *logger),
0617                     SurfaceMergingException);
0618
0619   // Disc with different r bounds: could be merged in r but not in phi
0620   auto discPhi5 = makeDisc(base, 100_mm, 150_mm, 45_degree, a(95_degree));
0621   BOOST_CHECK_THROW(discPhi->mergedWith(*discPhi5, Acts::AxisDirection::AxisPhi,
0622                                         false, *logger),
0623                     SurfaceMergingException);
0624 }
0625
0626 BOOST_DATA_TEST_CASE(PhiDirection,
0627                      (boost::unit_test::data::xrange(-135, 180, 45) *
0628                       boost::unit_test::data::make(Vector3{0_mm, 0_mm, 0_mm},
0629                                                    Vector3{20_mm, 0_mm, 0_mm},
0630                                                    Vector3{0_mm, 20_mm, 0_mm},
0631                                                    Vector3{20_mm, 20_mm, 0_mm},
0632                                                    Vector3{0_mm, 0_mm, 20_mm}) *
0633                       boost::unit_test::data::xrange(-1300, 1300, 104)),
0634                      angle, offset, phiShift) {
0635   Transform3 base =
0636       AngleAxis3(angle * 1_degree, Vector3::UnitX()) * Translation3(offset);
0637
0638   auto a = [phiShift](double v) {
0639     return detail::radian_sym(v + phiShift * 1_degree);
0640   };
0641
0642   BOOST_TEST_CONTEXT("Internal rotation") {
0643     auto disc = makeDisc(base, 30_mm, 100_mm, 10_degree, a(40_degree));
0644     auto disc2 = makeDisc(base, 30_mm, 100_mm, 45_degree, a(95_degree));
0645
0646     auto [disc3, reversed] =
0647         disc->mergedWith(*disc2, Acts::AxisDirection::AxisPhi, false, *logger);
0648     BOOST_REQUIRE_NE(disc3, nullptr);
0649     BOOST_CHECK_EQUAL(base.matrix(), disc3->transform(tgContext).matrix());
0650     BOOST_CHECK(reversed);
0651
0652     auto [disc3Reversed, reversed2] =
0653         disc2->mergedWith(*disc, Acts::AxisDirection::AxisPhi, false, *logger);
0654     BOOST_REQUIRE_NE(disc3Reversed, nullptr);
0655     BOOST_CHECK(*disc3 == *disc3Reversed);
0656     BOOST_CHECK(!reversed2);
0657
0658     const auto* bounds = dynamic_cast<const RadialBounds*>(&disc3->bounds());
0659     BOOST_REQUIRE_NE(bounds, nullptr);
0660
0661     BOOST_CHECK_SMALL(
0662         detail::difference_periodic(bounds->get(RadialBounds::eAveragePhi),
0663                                     a(85_degree), 2 * std::numbers::pi),
0664         1e-6);
0665     BOOST_CHECK_CLOSE(bounds->get(RadialBounds::eHalfPhiSector), 55_degree,
0666                       1e-6);
0667
0668     auto disc4 = makeDisc(base, 30_mm, 100_mm, 20_degree, a(170_degree));
0669     auto disc5 = makeDisc(base, 30_mm, 100_mm, 10_degree, a(-160_degree));
0670     auto [disc45, reversed45] =
0671         disc4->mergedWith(*disc5, Acts::AxisDirection::AxisPhi, false, *logger);
0672     BOOST_REQUIRE_NE(disc45, nullptr);
0673     BOOST_CHECK_EQUAL(base.matrix(), disc45->transform(tgContext).matrix());
0674     BOOST_CHECK(reversed45);
0675
0676     auto [disc54, reversed54] =
0677         disc5->mergedWith(*disc4, Acts::AxisDirection::AxisPhi, false, *logger);
0678     BOOST_REQUIRE_NE(disc54, nullptr);
0679     BOOST_CHECK(!reversed54);
0680
0681     BOOST_CHECK(*disc54 == *disc45);
0682
0683     const auto* bounds45 = dynamic_cast<const RadialBounds*>(&disc45->bounds());
0684     BOOST_REQUIRE_NE(bounds, nullptr);
0685
0686     BOOST_CHECK_SMALL(
0687         detail::difference_periodic(bounds45->get(RadialBounds::eAveragePhi),
0688                                     a(180_degree), 2 * std::numbers::pi),
0689         1e-6);
0690     BOOST_CHECK_CLOSE(bounds45->get(RadialBounds::eHalfPhiSector), 30_degree,
0691                       1e-6);
0692
0693     auto disc6 = makeDisc(base, 30_mm, 100_mm, 90_degree, a(0_degree));
0694     auto disc7 = makeDisc(base, 30_mm, 100_mm, 90_degree, a(180_degree));
0695
0696     auto [disc67, reversed67] =
0697         disc6->mergedWith(*disc7, Acts::AxisDirection::AxisPhi, false, *logger);
0698     BOOST_REQUIRE_NE(disc67, nullptr);
0699     CHECK_CLOSE_OR_SMALL(disc67->transform(tgContext).matrix(), base.matrix(),
0700                          1e-6, 1e-10);
0701     BOOST_CHECK(!reversed67);
0702
0703     auto [disc76, reversed76] =
0704         disc7->mergedWith(*disc6, Acts::AxisDirection::AxisPhi, false, *logger);
0705     BOOST_REQUIRE_NE(disc76, nullptr);
0706     // surfaces are not equal because bounds are not equal
0707     BOOST_CHECK(*disc76 != *disc67);
0708     // bounds are different because of avg phi
0709     BOOST_CHECK_NE(disc76->bounds(), disc67->bounds());
0710     // transforms should be the same
0711     BOOST_CHECK_EQUAL(disc76->transform(tgContext).matrix(),
0712                       disc67->transform(tgContext).matrix());
0713     // not reversed either because you get the ordering you put in
0714     BOOST_CHECK(!reversed76);
0715
0716     const auto* bounds67 = dynamic_cast<const RadialBounds*>(&disc67->bounds());
0717     BOOST_REQUIRE_NE(bounds67, nullptr);
0718     BOOST_CHECK_SMALL(
0719         detail::difference_periodic(bounds67->get(RadialBounds::eAveragePhi),
0720                                     a(90_degree), 2 * std::numbers::pi),
0721         1e-6);
0722     BOOST_CHECK_CLOSE(bounds67->get(RadialBounds::eHalfPhiSector), 180_degree,
0723                       1e-6);
0724   }
0725
0726   BOOST_TEST_CONTEXT("External rotation") {
0727     Transform3 trf1 = base * AngleAxis3(a(40_degree), Vector3::UnitZ());
0728     auto disc = makeDisc(trf1, 30_mm, 100_mm, 10_degree, 0_degree);
0729     Transform3 trf2 = base * AngleAxis3(a(95_degree), Vector3::UnitZ());
0730     auto disc2 = makeDisc(trf2, 30_mm, 100_mm, 45_degree, 0_degree);
0731
0732     auto [disc3, reversed] =
0733         disc->mergedWith(*disc2, Acts::AxisDirection::AxisPhi, true, *logger);
0734     BOOST_REQUIRE_NE(disc3, nullptr);
0735     Transform3 trfExpected12 =
0736         base * AngleAxis3(a(85_degree), Vector3::UnitZ());
0737     CHECK_CLOSE_OR_SMALL(disc3->transform(tgContext).matrix(),
0738                          trfExpected12.matrix(), 1e-6, 1e-10);
0739     BOOST_CHECK(reversed);
0740
0741     auto [disc3Reversed, reversed2] =
0742         disc2->mergedWith(*disc, Acts::AxisDirection::AxisPhi, true, *logger);
0743     BOOST_REQUIRE_NE(disc3Reversed, nullptr);
0744     BOOST_CHECK(*disc3 == *disc3Reversed);
0745     BOOST_CHECK(!reversed2);
0746
0747     const auto* bounds = dynamic_cast<const RadialBounds*>(&disc3->bounds());
0748     BOOST_REQUIRE_NE(bounds, nullptr);
0749
0750     BOOST_CHECK_EQUAL(bounds->get(RadialBounds::eAveragePhi), 0);
0751     BOOST_CHECK_CLOSE(bounds->get(RadialBounds::eHalfPhiSector), 55_degree,
0752                       1e-6);
0753
0754     Transform3 trf4 = base * AngleAxis3(a(170_degree), Vector3::UnitZ());
0755     auto disc4 = makeDisc(trf4, 30_mm, 100_mm, 20_degree, 0_degree);
0756     Transform3 trf5 = base * AngleAxis3(a(-160_degree), Vector3::UnitZ());
0757     auto disc5 = makeDisc(trf5, 30_mm, 100_mm, 10_degree, 0_degree);
0758     auto [disc45, reversed45] =
0759         disc4->mergedWith(*disc5, Acts::AxisDirection::AxisPhi, true, *logger);
0760     BOOST_REQUIRE_NE(disc45, nullptr);
0761     Transform3 trfExpected45 =
0762         base * AngleAxis3(a(180_degree), Vector3::UnitZ());
0763     CHECK_CLOSE_OR_SMALL(disc45->transform(tgContext).matrix(),
0764                          trfExpected45.matrix(), 1e-6, 1e-10);
0765     BOOST_CHECK(reversed45);
0766
0767     auto [disc54, reversed54] =
0768         disc5->mergedWith(*disc4, Acts::AxisDirection::AxisPhi, true, *logger);
0769     BOOST_REQUIRE_NE(disc54, nullptr);
0770     BOOST_CHECK(!reversed54);
0771
0772     BOOST_CHECK(*disc54 == *disc45);
0773
0774     const auto* bounds45 = dynamic_cast<const RadialBounds*>(&disc45->bounds());
0775     BOOST_REQUIRE_NE(bounds, nullptr);
0776
0777     BOOST_CHECK_EQUAL(bounds45->get(RadialBounds::eAveragePhi), 0);
0778     BOOST_CHECK_CLOSE(bounds45->get(RadialBounds::eHalfPhiSector), 30_degree,
0779                       1e-6);
0780
0781     Transform3 trf6 = base * AngleAxis3(a(0_degree), Vector3::UnitZ());
0782     auto disc6 = makeDisc(trf6, 30_mm, 100_mm, 90_degree, 0_degree);
0783     Transform3 trf7 = base * AngleAxis3(a(180_degree), Vector3::UnitZ());
0784     auto disc7 = makeDisc(trf7, 30_mm, 100_mm, 90_degree, 0_degree);
0785     auto [disc67, reversed67] =
0786         disc6->mergedWith(*disc7, Acts::AxisDirection::AxisPhi, true, *logger);
0787     BOOST_REQUIRE_NE(disc67, nullptr);
0788     Transform3 trfExpected67 =
0789         base * AngleAxis3(a(90_degree), Vector3::UnitZ());
0790     CHECK_CLOSE_OR_SMALL(disc67->transform(tgContext).matrix(),
0791                          trfExpected67.matrix(), 1e-6, 1e-10);
0792     BOOST_CHECK(!reversed67);
0793
0794     auto [disc76, reversed76] =
0795         disc7->mergedWith(*disc6, Acts::AxisDirection::AxisPhi, true, *logger);
0796     BOOST_REQUIRE_NE(disc76, nullptr);
0797     // surfaces are not equal due to different transforms
0798     BOOST_CHECK(*disc76 != *disc67);
0799     BOOST_CHECK_NE(disc76->transform(tgContext).matrix(),
0800                    disc67->transform(tgContext).matrix());
0801     // bounds should be equal however
0802     BOOST_CHECK_EQUAL(disc76->bounds(), disc67->bounds());
0803
0804     // not reversed either because you get the ordering you put in
0805     BOOST_CHECK(!reversed76);
0806
0807     const auto* bounds67 = dynamic_cast<const RadialBounds*>(&disc67->bounds());
0808     BOOST_REQUIRE_NE(bounds67, nullptr);
0809     BOOST_CHECK_EQUAL(bounds67->get(RadialBounds::eAveragePhi), 0);
0810     BOOST_CHECK_CLOSE(bounds67->get(RadialBounds::eHalfPhiSector), 180_degree,
0811                       1e-6);
0812   }
0813 }
0814
0815 BOOST_AUTO_TEST_SUITE_END()
0816
0817 BOOST_AUTO_TEST_SUITE_END()
0818
0819 }  // namespace Acts::Test