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