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