Core/Surfaces/CylinderSurfaceTests.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/CylinderBounds.hpp"
0025 #include "Acts/Surfaces/CylinderSurface.hpp"
0026 #include "Acts/Surfaces/Surface.hpp"
0027 #include "Acts/Surfaces/SurfaceBounds.hpp"
0028 #include "Acts/Surfaces/SurfaceMergingException.hpp"
0029 #include "Acts/Utilities/AxisDefinitions.hpp"
0030 #include "Acts/Utilities/BinUtility.hpp"
0031 #include "Acts/Utilities/BinningType.hpp"
0032 #include "Acts/Utilities/Diagnostics.hpp"
0033 #include "Acts/Utilities/Intersection.hpp"
0034 #include "Acts/Utilities/Logger.hpp"
0035 #include "Acts/Utilities/Result.hpp"
0036 #include "Acts/Utilities/ThrowAssert.hpp"
0037 #include "Acts/Utilities/TransformHelpers.hpp"
0038 #include "Acts/Utilities/detail/periodic.hpp"
0039 #include "ActsTests/CommonHelpers/DetectorElementStub.hpp"
0040 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0041 #include "ActsTests/CommonHelpers/SurfaceDerivatives.hpp"
0042
0043 #include <cmath>
0044 #include <memory>
0045 #include <numbers>
0046 #include <string>
0047
0048 using namespace Acts;
0049 using namespace Acts::UnitLiterals;
0050
0051 namespace ActsTests {
0052
0053 auto logger = Acts::getDefaultLogger("UnitTests", Acts::Logging::VERBOSE);
0054
0055 // Create a test context
0056 GeometryContext testContext = GeometryContext::dangerouslyDefaultConstruct();
0057
0058 BOOST_AUTO_TEST_SUITE(SurfacesSuite)
0059 /// Unit test for creating compliant/non-compliant CylinderSurface object
0060 BOOST_AUTO_TEST_CASE(CylinderSurfaceConstruction) {
0061   /// Test default construction
0062   // default construction is deleted
0063
0064   /// Constructor with transform, radius and halfZ
0065   const double radius = 1.;
0066   const double halfZ = 10.;
0067   const double halfPhiSector = std::numbers::pi / 8.;
0068   const Translation3 translation{0., 1., 2.};
0069
0070   auto pTransform = Transform3(translation);
0071   BOOST_CHECK_EQUAL(
0072       Surface::makeShared<CylinderSurface>(pTransform, radius, halfZ)->type(),
0073       Surface::Cylinder);
0074
0075   /// Constructor with transform pointer, radius, halfZ and halfPhiSector
0076   BOOST_CHECK_EQUAL(Surface::makeShared<CylinderSurface>(pTransform, radius,
0077                                                          halfZ, halfPhiSector)
0078                         ->type(),
0079                     Surface::Cylinder);
0080
0081   /// Constructor with transform and CylinderBounds pointer
0082   auto pCylinderBounds = std::make_shared<const CylinderBounds>(radius, halfZ);
0083   BOOST_CHECK_EQUAL(
0084       Surface::makeShared<CylinderSurface>(pTransform, pCylinderBounds)->type(),
0085       Surface::Cylinder);
0086
0087   /// Copy constructor
0088   auto cylinderSurfaceObject =
0089       Surface::makeShared<CylinderSurface>(pTransform, radius, halfZ);
0090   auto copiedCylinderSurface =
0091       Surface::makeShared<CylinderSurface>(*cylinderSurfaceObject);
0092   BOOST_CHECK_EQUAL(copiedCylinderSurface->type(), Surface::Cylinder);
0093   BOOST_CHECK(*copiedCylinderSurface == *cylinderSurfaceObject);
0094
0095   /// Copied and transformed
0096   auto copiedTransformedCylinderSurface = Surface::makeShared<CylinderSurface>(
0097       testContext, *cylinderSurfaceObject, pTransform);
0098   BOOST_CHECK_EQUAL(copiedTransformedCylinderSurface->type(),
0099                     Surface::Cylinder);
0100
0101   /// Construct with nullptr bounds
0102   BOOST_CHECK_THROW(auto nullBounds = Surface::makeShared<CylinderSurface>(
0103                         Transform3::Identity(), nullptr),
0104                     AssertionFailureException);
0105 }
0106
0107 /// Unit test for testing CylinderSurface properties
0108 BOOST_AUTO_TEST_CASE(CylinderSurfaceProperties) {
0109   /// Test clone method
0110   const double radius = 1.;
0111   const double halfZ = 10.;
0112   const Translation3 translation{0., 1., 2.};
0113
0114   auto pTransform = Transform3(translation);
0115   auto cylinderSurfaceObject =
0116       Surface::makeShared<CylinderSurface>(pTransform, radius, halfZ);
0117
0118   /// Test type (redundant)
0119   BOOST_CHECK_EQUAL(cylinderSurfaceObject->type(), Surface::Cylinder);
0120
0121   /// Test referencePosition
0122   Vector3 referencePosition{0., 1., 2.};
0123   CHECK_CLOSE_ABS(cylinderSurfaceObject->referencePosition(
0124                       testContext, AxisDirection::AxisPhi),
0125                   referencePosition, 1e-9);
0126
0127   /// Test referenceFrame
0128   const double invSqrt2 = 1. / std::numbers::sqrt2;
0129   Vector3 globalPosition{invSqrt2, 1. - invSqrt2, 0.};
0130   Vector3 globalPositionZ{invSqrt2, 1. - invSqrt2, 2.};
0131   Vector3 momentum{15., 15., 15.};
0132   Vector3 momentum2{6.6, -3., 2.};
0133   RotationMatrix3 expectedFrame;
0134   expectedFrame << invSqrt2, 0., invSqrt2, invSqrt2, 0., -invSqrt2, 0., 1., 0.;
0135   // check without shift
0136   CHECK_CLOSE_OR_SMALL(cylinderSurfaceObject->referenceFrame(
0137                            testContext, globalPosition, momentum),
0138                        expectedFrame, 1e-6, 1e-9);
0139   // check with shift and different momentum
0140   CHECK_CLOSE_OR_SMALL(cylinderSurfaceObject->referenceFrame(
0141                            testContext, globalPositionZ, momentum2),
0142                        expectedFrame, 1e-6, 1e-9);
0143
0144   /// Test normal, given 3D position
0145   Vector3 origin{0., 0., 0.};
0146   Vector3 normal3D = {0., -1., 0.};
0147   CHECK_CLOSE_ABS(cylinderSurfaceObject->normal(testContext, origin), normal3D,
0148                   1e-9);
0149
0150   Vector3 pos45deg = {invSqrt2, 1 + invSqrt2, 0.};
0151   Vector3 pos45degZ = {invSqrt2, 1 + invSqrt2, 4.};
0152   Vector3 normal45deg = {invSqrt2, invSqrt2, 0.};
0153   // test the normal vector
0154   CHECK_CLOSE_ABS(cylinderSurfaceObject->normal(testContext, pos45deg),
0155                   normal45deg, 1e-6 * invSqrt2);
0156   // test that the normal vector is independent of z coordinate
0157   CHECK_CLOSE_ABS(cylinderSurfaceObject->normal(testContext, pos45degZ),
0158                   normal45deg, 1e-6 * invSqrt2);
0159
0160   /// Test normal given 2D rphi position
0161   Vector2 positionPiBy2(1., 0.);
0162   Vector3 normalAtPiBy2{std::cos(1.), std::sin(1.), 0.};
0163   CHECK_CLOSE_ABS(cylinderSurfaceObject->normal(testContext, positionPiBy2),
0164                   normalAtPiBy2, 1e-9);
0165
0166   /// Test rotational symmetry axis
0167   Vector3 symmetryAxis{0., 0., 1.};
0168   CHECK_CLOSE_ABS(cylinderSurfaceObject->rotSymmetryAxis(testContext),
0169                   symmetryAxis, 1e-9);
0170
0171   /// Test bounds
0172   BOOST_CHECK_EQUAL(cylinderSurfaceObject->bounds().type(),
0173                     SurfaceBounds::eCylinder);
0174
0175   /// Test localToGlobal
0176   Vector2 localPosition{0., 0.};
0177   globalPosition = cylinderSurfaceObject->localToGlobal(
0178       testContext, localPosition, momentum);
0179   Vector3 expectedPosition{1, 1, 2};
0180   BOOST_CHECK_EQUAL(globalPosition, expectedPosition);
0181
0182   /// Testing globalToLocal
0183   localPosition = cylinderSurfaceObject
0184                       ->globalToLocal(testContext, globalPosition, momentum)
0185                       .value();
0186   Vector2 expectedLocalPosition{0., 0.};
0187   BOOST_CHECK_EQUAL(localPosition, expectedLocalPosition);
0188
0189   /// Test isOnSurface
0190   Vector3 offSurface{100, 1, 2};
0191   BOOST_CHECK(cylinderSurfaceObject->isOnSurface(
0192       testContext, globalPosition, momentum, BoundaryTolerance::None()));
0193   BOOST_CHECK(cylinderSurfaceObject->isOnSurface(testContext, globalPosition,
0194                                                  BoundaryTolerance::None()));
0195   BOOST_CHECK(!cylinderSurfaceObject->isOnSurface(
0196       testContext, offSurface, momentum, BoundaryTolerance::None()));
0197   BOOST_CHECK(!cylinderSurfaceObject->isOnSurface(testContext, offSurface,
0198                                                   BoundaryTolerance::None()));
0199
0200   /// Intersection test
0201   Vector3 direction{-1., 0, 0};
0202   auto sfIntersection = cylinderSurfaceObject->intersect(
0203       testContext, offSurface, direction, BoundaryTolerance::Infinite());
0204   Intersection3D expectedIntersect{Vector3{1, 1, 2}, 99.,
0205                                    IntersectionStatus::reachable};
0206   BOOST_CHECK(sfIntersection[0].isValid());
0207   CHECK_CLOSE_ABS(sfIntersection[0].position(), expectedIntersect.position(),
0208                   1e-9);
0209   CHECK_CLOSE_ABS(sfIntersection[0].pathLength(),
0210                   expectedIntersect.pathLength(), 1e-9);
0211   // there is a second solution & and it should be valid
0212   BOOST_CHECK(sfIntersection[1].isValid());
0213   // And it's path should be further away then the primary solution
0214   double pn = sfIntersection[0].pathLength();
0215   double pa = sfIntersection[1].pathLength();
0216   BOOST_CHECK_LT(std::abs(pn), std::abs(pa));
0217
0218   /// Test pathCorrection
0219   CHECK_CLOSE_REL(cylinderSurfaceObject->pathCorrection(testContext, offSurface,
0220                                                         momentum.normalized()),
0221                   std::numbers::sqrt3, 0.01);
0222
0223   /// Test name
0224   BOOST_CHECK_EQUAL(cylinderSurfaceObject->name(),
0225                     std::string("Acts::CylinderSurface"));
0226
0227   /// Test dump
0228   boost::test_tools::output_test_stream dumpOutput;
0229   std::string expected =
0230       "Acts::CylinderSurface\n\
0231      Center position  (x, y, z) = (0.0000, 1.0000, 2.0000)\n\
0232      Rotation:             colX = (1.000000, 0.000000, 0.000000)\n\
0233                            colY = (0.000000, 1.000000, 0.000000)\n\
0234                            colZ = (0.000000, 0.000000, 1.000000)\n\
0235      Bounds  : Acts::CylinderBounds: (radius, halfLengthZ, halfPhiSector, averagePhi) = (1.0000000, 10.0000000, 3.1415927, 0.0000000)";
0236   dumpOutput << cylinderSurfaceObject->toStream(testContext);
0237   BOOST_CHECK(dumpOutput.is_equal(expected));
0238 }
0239
0240 BOOST_AUTO_TEST_CASE(CylinderSurfaceEqualityOperators) {
0241   const double radius = 1.;
0242   const double halfZ = 10.;
0243   const Translation3 translation{0., 1., 2.};
0244
0245   auto pTransform = Transform3(translation);
0246   auto cylinderSurfaceObject =
0247       Surface::makeShared<CylinderSurface>(pTransform, radius, halfZ);
0248
0249   auto cylinderSurfaceObject2 =
0250       Surface::makeShared<CylinderSurface>(pTransform, radius, halfZ);
0251
0252   /// Test equality operator
0253   BOOST_CHECK(*cylinderSurfaceObject == *cylinderSurfaceObject2);
0254
0255   BOOST_TEST_CHECKPOINT(
0256       "Create and then assign a CylinderSurface object to the existing one");
0257   /// Test assignment
0258   auto assignedCylinderSurface =
0259       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 6.6, 5.4);
0260   // TODO: surface copy-assignment is deprecated; switch to a makeShared copy
0261   // when it is removed (SonarCloud cpp:S3624).
0262   ACTS_PUSH_IGNORE_DEPRECATED()
0263   *assignedCylinderSurface = *cylinderSurfaceObject;
0264   ACTS_POP_IGNORE_DEPRECATED()
0265   /// Test equality of assigned to original
0266   BOOST_CHECK(*assignedCylinderSurface == *cylinderSurfaceObject);
0267 }
0268
0269 /// Unit test for testing CylinderSurface properties
0270 BOOST_AUTO_TEST_CASE(CylinderSurfaceExtent) {
0271   using enum AxisDirection;
0272
0273   // Some radius and half length
0274   const double radius = 1.;
0275   const double halfZ = 10.;
0276   const Translation3 translation{0., 0., 2.};  // != {0., 1., 2.}
0277
0278   auto pTransform = Transform3(translation);
0279   auto cylinderSurface =
0280       Surface::makeShared<CylinderSurface>(pTransform, radius, halfZ);
0281   // The Extent, let's measure it
0282   auto cylinderExtent =
0283       cylinderSurface->polyhedronRepresentation(testContext, 1).extent();
0284
0285   CHECK_CLOSE_ABS(-8, cylinderExtent.min(AxisZ), s_onSurfaceTolerance);
0286   CHECK_CLOSE_ABS(12, cylinderExtent.max(AxisZ), s_onSurfaceTolerance);
0287   CHECK_CLOSE_ABS(radius, cylinderExtent.min(AxisR), s_onSurfaceTolerance);
0288   CHECK_CLOSE_ABS(radius, cylinderExtent.max(AxisR), s_onSurfaceTolerance);
0289   CHECK_CLOSE_ABS(-radius, cylinderExtent.min(AxisX), s_onSurfaceTolerance);
0290   CHECK_CLOSE_ABS(radius, cylinderExtent.max(AxisX), s_onSurfaceTolerance);
0291   CHECK_CLOSE_ABS(-radius, cylinderExtent.min(AxisY), s_onSurfaceTolerance);
0292   CHECK_CLOSE_ABS(radius, cylinderExtent.max(AxisY), s_onSurfaceTolerance);
0293 }
0294
0295 /// Unit test for testing CylinderSurface alignment derivatives
0296 BOOST_AUTO_TEST_CASE(CylinderSurfaceAlignment) {
0297   const double radius = 1.;
0298   const double halfZ = 10.;
0299   const Translation3 translation{0., 1., 2.};
0300   // Rotate by 90 degrees about global z. The cylinder axis stays global z, so
0301   // the point below stays on the surface. The center columns are then the
0302   // local-frame components, not the global center derivative.
0303   const AngleAxis3 rotationAngle(std::numbers::pi / 2., Vector3::UnitZ());
0304   Transform3 transform =
0305       makeTransform3(rotationAngle.toRotationMatrix(), translation.vector());
0306   auto cylinderSurfaceObject =
0307       Surface::makeShared<CylinderSurface>(transform, radius, halfZ);
0308
0309   const auto& rotation = transform.rotation();
0310   // The local frame z axis
0311   const Vector3 localZAxis = rotation.col(2);
0312   // Check the local z axis is aligned to global z axis
0313   CHECK_CLOSE_ABS(localZAxis, Vector3(0., 0., 1.), 1e-15);
0314
0315   // On the cylinder, along the rotated local x axis.
0316   Vector3 globalPosition{0, 2, 2};
0317
0318   // Test the derivative of bound track parameters local position w.r.t.
0319   // position in local 3D Cartesian coordinates
0320   const auto& loc3DToLocBound =
0321       cylinderSurfaceObject->localCartesianToBoundLocalDerivative(
0322           testContext, globalPosition, Vector3::UnitZ());
0323   // Local position is (R, 0, 0), so d(phi)/d(local) is along local y.
0324   Matrix<2, 3> expLoc3DToLocBound = Matrix<2, 3>::Zero();
0325   expLoc3DToLocBound << 0, 1, 0, 0, 0, 1;
0326   CHECK_CLOSE_ABS(loc3DToLocBound, expLoc3DToLocBound, 1e-10);
0327
0328   // Local hit (1, 0, 2) and local direction (0.6, 0, 0.8). The denominator
0329   // is 0.6, so the local path derivative is (1, 0, 0, 0, 2, 0) / 0.6.
0330   const Vector3 hit = transform * Vector3{1., 0., 2.};
0331   const Vector3 direction = rotation * Vector3{0.6, 0., 0.8};
0332   const AlignmentToPathMatrix alignToPath =
0333       cylinderSurfaceObject->alignmentToPathDerivative(testContext, hit,
0334                                                        direction);
0335   AlignmentToPathMatrix expected = AlignmentToPathMatrix::Zero();
0336   expected << 1. / 0.6, 0., 0., 0., 2. / 0.6, 0.;
0337   CHECK_CLOSE_ABS(alignToPath, expected, 1e-10);
0338
0339   // A hit off the local x axis, checked against finite differences.
0340   const Vector3 skewedHit = transform * Vector3{0.6, 0.8, 1.5};
0341   const Vector3 skewedDir = rotation * Vector3{0.6, 0.3, 0.2};
0342   checkSurfaceDerivatives(testContext, *cylinderSurfaceObject, transform,
0343                           skewedHit, skewedDir, [&](const Transform3& tf) {
0344                             return Surface::makeShared<CylinderSurface>(
0345                                 tf, radius, halfZ);
0346                           });
0347 }
0348
0349 BOOST_AUTO_TEST_CASE(CylinderSurfaceBinningPosition) {
0350   using namespace Acts::UnitLiterals;
0351   Vector3 s{5_mm, 7_mm, 10_cm};
0352   Transform3 trf;
0353   trf = Translation3(s) * AngleAxis3{0.5, Vector3::UnitZ()};
0354
0355   double r = 300;
0356   double halfZ = 330;
0357   double averagePhi = 0.1;
0358
0359   auto bounds = std::make_shared<CylinderBounds>(r, halfZ, std::numbers::pi / 8,
0360                                                  averagePhi);
0361   auto cylinder = Acts::Surface::makeShared<CylinderSurface>(trf, bounds);
0362
0363   Vector3 exp = Vector3{r * std::cos(averagePhi), r * std::sin(averagePhi), 0};
0364   exp = trf * exp;
0365
0366   Vector3 bp = cylinder->referencePosition(testContext, AxisDirection::AxisR);
0367   CHECK_CLOSE_ABS(bp, exp, 1e-10);
0368   CHECK_CLOSE_ABS(
0369       cylinder->referencePositionValue(testContext, AxisDirection::AxisR),
0370       VectorHelpers::perp(exp), 1e-10);
0371
0372   bp = cylinder->referencePosition(testContext, AxisDirection::AxisRPhi);
0373   CHECK_CLOSE_ABS(bp, exp, 1e-10);
0374   CHECK_CLOSE_ABS(
0375       cylinder->referencePositionValue(testContext, AxisDirection::AxisRPhi),
0376       VectorHelpers::phi(exp) * VectorHelpers::perp(exp), 1e-10);
0377
0378   for (auto b : {AxisDirection::AxisX, AxisDirection::AxisY,
0379                  AxisDirection::AxisZ, AxisDirection::AxisEta,
0380                  AxisDirection::AxisTheta, AxisDirection::AxisMag}) {
0381     BOOST_TEST_CONTEXT("binValue: " << b) {
0382       BOOST_CHECK_EQUAL(cylinder->referencePosition(testContext, b),
0383                         cylinder->center(testContext));
0384     }
0385   }
0386 }
0387
0388 BOOST_AUTO_TEST_SUITE(CylinderSurfaceMerging)
0389
0390 BOOST_AUTO_TEST_CASE(InvalidDetectorElement) {
0391   DetectorElementStub detElem;
0392
0393   auto bounds = std::make_shared<CylinderBounds>(100_mm, 100_mm);
0394   auto cyl1 = Surface::makeShared<CylinderSurface>(bounds, detElem);
0395   auto cyl2 = Surface::makeShared<CylinderSurface>(bounds, detElem);
0396
0397   BOOST_CHECK_THROW(
0398       cyl1->mergedWith(*cyl2, Acts::AxisDirection::AxisR, false, *logger),
0399       SurfaceMergingException);
0400 }
0401
0402 BOOST_DATA_TEST_CASE(IncompatibleZDirection,
0403                      (boost::unit_test::data::xrange(-135, 180, 45) *
0404                       boost::unit_test::data::make(Vector3{0_mm, 0_mm, 0_mm},
0405                                                    Vector3{20_mm, 0_mm, 0_mm},
0406                                                    Vector3{0_mm, 20_mm, 0_mm},
0407                                                    Vector3{20_mm, 20_mm, 0_mm},
0408                                                    Vector3{0_mm, 0_mm, 20_mm})),
0409                      angle, offset) {
0410   Logging::ScopedFailureThreshold ft{Logging::FATAL};
0411
0412   Transform3 base =
0413       AngleAxis3(angle * 1_degree, Vector3::UnitX()) * Translation3(offset);
0414
0415   auto cyl = Surface::makeShared<CylinderSurface>(base, 30_mm, 100_mm);
0416   auto cyl2 = Surface::makeShared<CylinderSurface>(
0417       base * Translation3{Vector3::UnitZ() * 200_mm}, 30_mm, 100_mm);
0418
0419   BOOST_CHECK_THROW(
0420       cyl->mergedWith(*cyl2, Acts::AxisDirection::AxisPhi, false, *logger),
0421       SurfaceMergingException);
0422
0423   auto cylShiftedXy = Surface::makeShared<CylinderSurface>(
0424       base * Translation3{Vector3{1_mm, 2_mm, 200_mm}}, 30_mm, 100_mm);
0425   BOOST_CHECK_THROW(cyl->mergedWith(*cylShiftedXy, Acts::AxisDirection::AxisZ,
0426                                     false, *logger),
0427                     SurfaceMergingException);
0428
0429   auto cylRotatedX = Surface::makeShared<CylinderSurface>(
0430       base * AngleAxis3{10_degree, Vector3::UnitX()} *
0431           Translation3{Vector3::UnitZ() * 200_mm},
0432       30_mm, 100_mm);
0433   BOOST_CHECK_THROW(
0434       cyl->mergedWith(*cylRotatedX, Acts::AxisDirection::AxisZ, false, *logger),
0435       SurfaceMergingException);
0436
0437   // Cylinder with different radius
0438   auto cyl3 = Surface::makeShared<CylinderSurface>(
0439       base * Translation3{Vector3::UnitZ() * 200_mm}, 35_mm, 100_mm);
0440   BOOST_CHECK_THROW(
0441       cyl->mergedWith(*cyl3, Acts::AxisDirection::AxisZ, false, *logger),
0442       SurfaceMergingException);
0443
0444   // Cylinder with overlap in z
0445   auto cyl6 = Surface::makeShared<CylinderSurface>(
0446       base * Translation3{Vector3::UnitZ() * 150_mm}, 30_mm, 100_mm);
0447   BOOST_CHECK_THROW(
0448       cyl->mergedWith(*cyl6, Acts::AxisDirection::AxisZ, false, *logger),
0449       SurfaceMergingException);
0450
0451   // Cylinder with gap in z
0452   auto cyl7 = Surface::makeShared<CylinderSurface>(
0453       base * Translation3{Vector3::UnitZ() * 250_mm}, 30_mm, 100_mm);
0454   BOOST_CHECK_THROW(
0455       cyl->mergedWith(*cyl7, Acts::AxisDirection::AxisZ, false, *logger),
0456       SurfaceMergingException);
0457
0458   // Cylinder with phi sector and relative z rotation
0459   auto cyl8 = Surface::makeShared<CylinderSurface>(
0460       base * AngleAxis3(14_degree, Vector3::UnitZ()) *
0461           Translation3{Vector3::UnitZ() * 200_mm},
0462       30_mm, 100_mm, 10_degree, 40_degree);
0463   BOOST_CHECK_THROW(
0464       cyl->mergedWith(*cyl8, Acts::AxisDirection::AxisZ, false, *logger),
0465       SurfaceMergingException);
0466
0467   auto cylPhi1 = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
0468                                                       30_mm, 100_mm, 45_degree);
0469   auto cylPhi2 = Surface::makeShared<CylinderSurface>(
0470       Transform3{Translation3{Vector3::UnitZ() * 150_mm}}, 30_mm, 50_mm,
0471       55_degree);
0472   BOOST_CHECK_THROW(
0473       cylPhi1->mergedWith(*cylPhi2, Acts::AxisDirection::AxisZ, false, *logger),
0474       SurfaceMergingException);
0475 }
0476
0477 BOOST_DATA_TEST_CASE(ZDirection,
0478                      (boost::unit_test::data::xrange(-135, 180, 45) *
0479                       boost::unit_test::data::make(Vector3{0_mm, 0_mm, 0_mm},
0480                                                    Vector3{20_mm, 0_mm, 0_mm},
0481                                                    Vector3{0_mm, 20_mm, 0_mm},
0482                                                    Vector3{20_mm, 20_mm, 0_mm},
0483                                                    Vector3{0_mm, 0_mm, 20_mm})),
0484                      angle, offset) {
0485   Transform3 base =
0486       AngleAxis3(angle * 1_degree, Vector3::UnitX()) * Translation3(offset);
0487
0488   auto cyl = Surface::makeShared<CylinderSurface>(base, 30_mm, 100_mm);
0489
0490   auto cyl2 = Surface::makeShared<CylinderSurface>(
0491       base * AngleAxis3(14_degree, Vector3::UnitZ()) *
0492           Translation3{Vector3::UnitZ() * 200_mm},
0493       30_mm, 100_mm);
0494
0495   auto [cyl3, reversed] =
0496       cyl->mergedWith(*cyl2, Acts::AxisDirection::AxisZ, false, *logger);
0497   BOOST_REQUIRE_NE(cyl3, nullptr);
0498   BOOST_CHECK(!reversed);
0499
0500   auto [cyl3Reversed, reversed2] =
0501       cyl2->mergedWith(*cyl, Acts::AxisDirection::AxisZ, false, *logger);
0502   BOOST_REQUIRE_NE(cyl3Reversed, nullptr);
0503   BOOST_CHECK(cyl3->bounds() == cyl3Reversed->bounds());
0504   BOOST_CHECK(reversed2);
0505
0506   auto bounds = cyl3->bounds();
0507
0508   BOOST_CHECK_EQUAL(bounds.get(CylinderBounds::eR), 30_mm);
0509   BOOST_CHECK_EQUAL(bounds.get(CylinderBounds::eHalfLengthZ), 200_mm);
0510   BOOST_CHECK_EQUAL(bounds.get(CylinderBounds::eAveragePhi), 0_degree);
0511   BOOST_CHECK_EQUAL(bounds.get(CylinderBounds::eHalfPhiSector), 180_degree);
0512
0513   // Rotation in z depends on the ordering, the left side "wins"
0514   Transform3 expected12 = base * Translation3{Vector3::UnitZ() * 100_mm};
0515   BOOST_CHECK_EQUAL(expected12.matrix(),
0516                     cyl3->localToGlobalTransform(testContext).matrix());
0517
0518   Transform3 expected21 = base * AngleAxis3(14_degree, Vector3::UnitZ()) *
0519                           Translation3{Vector3::UnitZ() * 100_mm};
0520   CHECK_CLOSE_OR_SMALL(
0521       cyl3Reversed->localToGlobalTransform(testContext).matrix(),
0522       expected21.matrix(), 1e-6, 1e-10);
0523
0524   auto cylPhi1 = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
0525                                                       30_mm, 100_mm, 45_degree);
0526   auto cylPhi2 = Surface::makeShared<CylinderSurface>(
0527       Transform3{Translation3{Vector3::UnitZ() * 150_mm}}, 30_mm, 50_mm,
0528       45_degree);
0529
0530   auto [cylPhi12, reversedPhy12] =
0531       cylPhi1->mergedWith(*cylPhi2, Acts::AxisDirection::AxisZ, false, *logger);
0532
0533   BOOST_REQUIRE_NE(cylPhi12, nullptr);
0534   auto boundsPhi12 = cylPhi12->bounds();
0535   BOOST_CHECK_EQUAL(boundsPhi12.get(CylinderBounds::eR), 30_mm);
0536   BOOST_CHECK_EQUAL(boundsPhi12.get(CylinderBounds::eHalfLengthZ), 150_mm);
0537   BOOST_CHECK_EQUAL(boundsPhi12.get(CylinderBounds::eAveragePhi), 0_degree);
0538   BOOST_CHECK_EQUAL(boundsPhi12.get(CylinderBounds::eHalfPhiSector), 45_degree);
0539 }
0540
0541 BOOST_DATA_TEST_CASE(IncompatibleRPhiDirection,
0542                      (boost::unit_test::data::xrange(-135, 180, 45) *
0543                       boost::unit_test::data::make(Vector3{0_mm, 0_mm, 0_mm},
0544                                                    Vector3{20_mm, 0_mm, 0_mm},
0545                                                    Vector3{0_mm, 20_mm, 0_mm},
0546                                                    Vector3{20_mm, 20_mm, 0_mm},
0547                                                    Vector3{0_mm, 0_mm, 20_mm}) *
0548                       boost::unit_test::data::xrange(-1300, 1300, 104)),
0549                      angle, offset, phiShift) {
0550   Logging::ScopedFailureThreshold ft{Logging::FATAL};
0551   Transform3 base =
0552       AngleAxis3(angle * 1_degree, Vector3::UnitX()) * Translation3(offset);
0553
0554   auto a = [phiShift](double v) {
0555     return detail::radian_sym(v + phiShift * 1_degree);
0556   };
0557
0558   auto cylPhi = Surface::makeShared<CylinderSurface>(base, 30_mm, 100_mm,
0559                                                      10_degree, a(40_degree));
0560
0561   // Cylinder with overlap in phi
0562   auto cylPhi2 = Surface::makeShared<CylinderSurface>(base, 30_mm, 100_mm,
0563                                                       45_degree, a(85_degree));
0564   BOOST_CHECK_THROW(cylPhi->mergedWith(*cylPhi2, Acts::AxisDirection::AxisRPhi,
0565                                        false, *logger),
0566                     SurfaceMergingException);
0567
0568   // Cylinder with gap in phi
0569   auto cylPhi3 = Surface::makeShared<CylinderSurface>(base, 30_mm, 100_mm,
0570                                                       45_degree, a(105_degree));
0571   BOOST_CHECK_THROW(cylPhi->mergedWith(*cylPhi3, Acts::AxisDirection::AxisRPhi,
0572                                        false, *logger),
0573                     SurfaceMergingException);
0574
0575   // Cylinder with a z shift
0576   auto cylPhi4 = Surface::makeShared<CylinderSurface>(
0577       base * Translation3{Vector3::UnitZ() * 20_mm}, 30_mm, 100_mm, 45_degree,
0578       a(95_degree));
0579   BOOST_CHECK_THROW(cylPhi->mergedWith(*cylPhi4, Acts::AxisDirection::AxisRPhi,
0580                                        false, *logger),
0581                     SurfaceMergingException);
0582
0583   // Test phi sector with different z halflengths
0584   auto cylPhi5 = Surface::makeShared<CylinderSurface>(base, 30_mm, 110_mm,
0585                                                       45_degree, a(95_degree));
0586   BOOST_CHECK_THROW(cylPhi->mergedWith(*cylPhi5, Acts::AxisDirection::AxisRPhi,
0587                                        false, *logger),
0588                     SurfaceMergingException);
0589 }
0590
0591 BOOST_DATA_TEST_CASE(RPhiDirection,
0592                      (boost::unit_test::data::xrange(-135, 180, 45) *
0593                       boost::unit_test::data::make(Vector3{0_mm, 0_mm, 0_mm},
0594                                                    Vector3{20_mm, 0_mm, 0_mm},
0595                                                    Vector3{0_mm, 20_mm, 0_mm},
0596                                                    Vector3{20_mm, 20_mm, 0_mm},
0597                                                    Vector3{0_mm, 0_mm, 20_mm}) *
0598                       boost::unit_test::data::xrange(-1300, 1300, 104)),
0599                      angle, offset, phiShift) {
0600   Transform3 base =
0601       AngleAxis3(angle * 1_degree, Vector3::UnitX()) * Translation3(offset);
0602
0603   auto a = [phiShift](double v) {
0604     return detail::radian_sym(v + phiShift * 1_degree);
0605   };
0606
0607   BOOST_TEST_CONTEXT("Internal rotation") {
0608     auto cyl = Surface::makeShared<CylinderSurface>(base, 30_mm, 100_mm,
0609                                                     10_degree, a(40_degree));
0610     auto cyl2 = Surface::makeShared<CylinderSurface>(base, 30_mm, 100_mm,
0611                                                      45_degree, a(95_degree));
0612
0613     auto [cyl3, reversed] =
0614         cyl->mergedWith(*cyl2, Acts::AxisDirection::AxisRPhi, false, *logger);
0615     BOOST_REQUIRE_NE(cyl3, nullptr);
0616     BOOST_CHECK_EQUAL(base.matrix(),
0617                       cyl3->localToGlobalTransform(testContext).matrix());
0618     BOOST_CHECK(reversed);
0619
0620     auto [cyl3Reversed, reversed2] =
0621         cyl2->mergedWith(*cyl, Acts::AxisDirection::AxisRPhi, false, *logger);
0622     BOOST_REQUIRE_NE(cyl3Reversed, nullptr);
0623     BOOST_CHECK(*cyl3 == *cyl3Reversed);
0624     BOOST_CHECK(!reversed2);
0625
0626     const auto& bounds = cyl3->bounds();
0627
0628     BOOST_CHECK_SMALL(
0629         detail::difference_periodic(bounds.get(CylinderBounds::eAveragePhi),
0630                                     a(85_degree), 2 * std::numbers::pi),
0631         1e-6);
0632     BOOST_CHECK_CLOSE(bounds.get(CylinderBounds::eHalfPhiSector), 55_degree,
0633                       0.1);
0634
0635     auto cyl4 = Surface::makeShared<CylinderSurface>(base, 30_mm, 100_mm,
0636                                                      20_degree, a(170_degree));
0637     auto cyl5 = Surface::makeShared<CylinderSurface>(base, 30_mm, 100_mm,
0638                                                      10_degree, a(-160_degree));
0639     auto [cyl45, reversed45] =
0640         cyl4->mergedWith(*cyl5, Acts::AxisDirection::AxisRPhi, false, *logger);
0641     BOOST_REQUIRE_NE(cyl45, nullptr);
0642     BOOST_CHECK_EQUAL(base.matrix(),
0643                       cyl45->localToGlobalTransform(testContext).matrix());
0644     BOOST_CHECK(reversed45);
0645
0646     auto [cyl54, reversed54] =
0647         cyl5->mergedWith(*cyl4, Acts::AxisDirection::AxisRPhi, false, *logger);
0648     BOOST_REQUIRE_NE(cyl54, nullptr);
0649     BOOST_CHECK(!reversed54);
0650
0651     BOOST_CHECK(*cyl54 == *cyl45);
0652
0653     BOOST_CHECK_SMALL(detail::difference_periodic(
0654                           cyl45->bounds().get(CylinderBounds::eAveragePhi),
0655                           a(180_degree), 2 * std::numbers::pi),
0656                       1e-6);
0657     BOOST_CHECK_CLOSE(cyl45->bounds().get(CylinderBounds::eHalfPhiSector),
0658                       30_degree, 1e-6);
0659
0660     auto cyl6 = Surface::makeShared<CylinderSurface>(base, 30_mm, 100_mm,
0661                                                      90_degree, a(90_degree));
0662     auto cyl7 = Surface::makeShared<CylinderSurface>(base, 30_mm, 100_mm,
0663                                                      90_degree, a(-90_degree));
0664
0665     auto [cyl67, reversed67] =
0666         cyl6->mergedWith(*cyl7, Acts::AxisDirection::AxisRPhi, false, *logger);
0667     BOOST_REQUIRE_NE(cyl67, nullptr);
0668     BOOST_CHECK_EQUAL(base.matrix(),
0669                       cyl67->localToGlobalTransform(testContext).matrix());
0670
0671     auto [cyl76, reversed76] =
0672         cyl7->mergedWith(*cyl6, Acts::AxisDirection::AxisRPhi, false, *logger);
0673     BOOST_REQUIRE_NE(cyl76, nullptr);
0674     BOOST_CHECK_EQUAL(base.matrix(),
0675                       cyl76->localToGlobalTransform(testContext).matrix());
0676
0677     // The ordering in this case is effectively arbitrary, you get the ordering
0678     // you put in
0679     BOOST_CHECK(!reversed67);
0680     BOOST_CHECK(!reversed76);
0681
0682     BOOST_CHECK_SMALL(detail::difference_periodic(
0683                           cyl67->bounds().get(CylinderBounds::eAveragePhi),
0684                           a(180_degree), 2 * std::numbers::pi),
0685                       1e-6);
0686     BOOST_CHECK_CLOSE(cyl67->bounds().get(CylinderBounds::eHalfPhiSector),
0687                       180_degree, 1e-6);
0688   }
0689
0690   BOOST_TEST_CONTEXT("External rotation") {
0691     Transform3 trf1 = base * AngleAxis3(a(40_degree), Vector3::UnitZ());
0692     auto cyl1 = Surface::makeShared<CylinderSurface>(trf1, 30_mm, 100_mm,
0693                                                      10_degree, 0_degree);
0694
0695     Transform3 trf2 = base * AngleAxis3(a(95_degree), Vector3::UnitZ());
0696     auto cyl2 = Surface::makeShared<CylinderSurface>(trf2, 30_mm, 100_mm,
0697                                                      45_degree, 0_degree);
0698
0699     auto [cyl3, reversed] =
0700         cyl1->mergedWith(*cyl2, Acts::AxisDirection::AxisRPhi, true, *logger);
0701
0702     BOOST_REQUIRE_NE(cyl3, nullptr);
0703     Transform3 trfExpected12 =
0704         base * AngleAxis3(a(85_degree), Vector3::UnitZ());
0705     CHECK_CLOSE_OR_SMALL(cyl3->localToGlobalTransform(testContext).matrix(),
0706                          trfExpected12.matrix(), 1e-6, 1e-10);
0707     BOOST_CHECK(reversed);
0708
0709     BOOST_CHECK_EQUAL(cyl3->bounds().get(CylinderBounds::eAveragePhi), 0);
0710     BOOST_CHECK_CLOSE(cyl3->bounds().get(CylinderBounds::eHalfPhiSector),
0711                       55_degree, 1e-6);
0712
0713     Transform3 trf4 = base * AngleAxis3(a(170_degree), Vector3::UnitZ());
0714     auto cyl4 = Surface::makeShared<CylinderSurface>(trf4, 30_mm, 100_mm,
0715                                                      20_degree, 0_degree);
0716     Transform3 trf5 = base * AngleAxis3(a(-160_degree), Vector3::UnitZ());
0717     auto cyl5 = Surface::makeShared<CylinderSurface>(trf5, 30_mm, 100_mm,
0718                                                      10_degree, 0_degree);
0719     auto [cyl45, reversed45] =
0720         cyl4->mergedWith(*cyl5, Acts::AxisDirection::AxisRPhi, true, *logger);
0721     BOOST_REQUIRE_NE(cyl45, nullptr);
0722     Transform3 trfExpected45 =
0723         base * AngleAxis3(a(180_degree), Vector3::UnitZ());
0724     CHECK_CLOSE_OR_SMALL(cyl45->localToGlobalTransform(testContext).matrix(),
0725                          trfExpected45.matrix(), 1e-6, 1e-10);
0726     BOOST_CHECK(reversed45);
0727
0728     auto [cyl54, reversed54] =
0729         cyl5->mergedWith(*cyl4, Acts::AxisDirection::AxisRPhi, true, *logger);
0730     BOOST_REQUIRE_NE(cyl54, nullptr);
0731     BOOST_CHECK(!reversed54);
0732
0733     BOOST_CHECK(*cyl54 == *cyl45);
0734
0735     BOOST_CHECK_EQUAL(cyl45->bounds().get(CylinderBounds::eAveragePhi), 0);
0736     BOOST_CHECK_CLOSE(cyl45->bounds().get(CylinderBounds::eHalfPhiSector),
0737                       30_degree, 1e-6);
0738
0739     Transform3 trf6 = base * AngleAxis3(a(90_degree), Vector3::UnitZ());
0740     auto cyl6 = Surface::makeShared<CylinderSurface>(trf6, 30_mm, 100_mm,
0741                                                      90_degree, 0_degree);
0742     Transform3 trf7 = base * AngleAxis3(a(-90_degree), Vector3::UnitZ());
0743     auto cyl7 = Surface::makeShared<CylinderSurface>(trf7, 30_mm, 100_mm,
0744                                                      90_degree, 0_degree);
0745
0746     auto [cyl67, reversed67] =
0747         cyl6->mergedWith(*cyl7, Acts::AxisDirection::AxisRPhi, true, *logger);
0748     BOOST_REQUIRE_NE(cyl67, nullptr);
0749     Transform3 expected67 = trf6 * AngleAxis3(90_degree, Vector3::UnitZ());
0750     CHECK_CLOSE_OR_SMALL(cyl67->localToGlobalTransform(testContext).matrix(),
0751                          expected67.matrix(), 1e-6, 1e-10);
0752
0753     auto [cyl76, reversed76] =
0754         cyl7->mergedWith(*cyl6, Acts::AxisDirection::AxisRPhi, true, *logger);
0755     BOOST_REQUIRE_NE(cyl76, nullptr);
0756     Transform3 expected76 = trf7 * AngleAxis3(90_degree, Vector3::UnitZ());
0757     CHECK_CLOSE_OR_SMALL(cyl76->localToGlobalTransform(testContext).matrix(),
0758                          expected76.matrix(), 1e-6, 1e-10);
0759
0760     // The ordering in this case is effectively arbitrary, you get the ordering
0761     // you put in
0762     BOOST_CHECK(!reversed67);
0763     BOOST_CHECK(!reversed76);
0764
0765     BOOST_CHECK_EQUAL(cyl67->bounds().get(CylinderBounds::eAveragePhi), 0);
0766     BOOST_CHECK_CLOSE(cyl67->bounds().get(CylinderBounds::eHalfPhiSector),
0767                       180_degree, 0.1);
0768   }
0769 }
0770
0771 BOOST_AUTO_TEST_SUITE_END()
0772
0773 BOOST_AUTO_TEST_CASE(CylinderSurfaceMaterialAssignment) {
0774   auto surface =
0775       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 10., 100.);
0776   MaterialSlab slab(Material::fromMolarDensity(1., 2., 3., 4., 5.), 0.1);
0777
0778   // HomogeneousSurfaceMaterial is always valid
0779   auto homMat = std::make_shared<HomogeneousSurfaceMaterial>(slab);
0780   BOOST_CHECK_NO_THROW(surface->assignSurfaceMaterial(homMat));
0781   BOOST_CHECK_NE(surface->surfaceMaterial(), nullptr);
0782
0783   // nullptr clears the material
0784   BOOST_CHECK_NO_THROW(surface->assignSurfaceMaterial(nullptr));
0785   BOOST_CHECK_EQUAL(surface->surfaceMaterial(), nullptr);
0786
0787   // 1D AxisRPhi - valid for cylinder
0788   BinUtility buRPhi(10, 0.f, 62.8f, Acts::closed, AxisDirection::AxisRPhi);
0789   auto matRPhi = std::make_shared<BinnedSurfaceMaterial>(
0790       buRPhi, MaterialSlabVector(10, slab));
0791   BOOST_CHECK_NO_THROW(surface->assignSurfaceMaterial(matRPhi));
0792
0793   // 1D AxisPhi - valid for cylinder
0794   BinUtility buPhi(10, 0.f, 1.f, Acts::closed, AxisDirection::AxisPhi);
0795   auto matPhi = std::make_shared<BinnedSurfaceMaterial>(
0796       buPhi, MaterialSlabVector(10, slab));
0797   BOOST_CHECK_NO_THROW(surface->assignSurfaceMaterial(matPhi));
0798
0799   // 1D AxisZ - valid for cylinder
0800   BinUtility buZ(10, -100.f, 100.f, Acts::open, AxisDirection::AxisZ);
0801   auto matZ = std::make_shared<BinnedSurfaceMaterial>(
0802       buZ, MaterialSlabVector(10, slab));
0803   BOOST_CHECK_NO_THROW(surface->assignSurfaceMaterial(matZ));
0804
0805   // 2D {AxisRPhi, AxisZ} - valid for cylinder
0806   BinUtility bu2D(5, 0.f, 62.8f, Acts::closed, AxisDirection::AxisRPhi);
0807   bu2D += BinUtility(3, -100.f, 100.f, Acts::open, AxisDirection::AxisZ);
0808   auto mat2D = std::make_shared<BinnedSurfaceMaterial>(
0809       bu2D, MaterialSlabMatrix(3, MaterialSlabVector(5, slab)));
0810   BOOST_CHECK_NO_THROW(surface->assignSurfaceMaterial(mat2D));
0811
0812   // Wrong axis directions - should throw
0813   for (auto badDir :
0814        {AxisDirection::AxisR, AxisDirection::AxisX, AxisDirection::AxisY}) {
0815     BinUtility buBad(10, 0.f, 10.f, Acts::open, badDir);
0816     auto matBad = std::make_shared<BinnedSurfaceMaterial>(
0817         buBad, MaterialSlabVector(10, slab));
0818     BOOST_CHECK_THROW(surface->assignSurfaceMaterial(matBad),
0819                       std::invalid_argument);
0820   }
0821 }
0822
0823 BOOST_AUTO_TEST_SUITE_END()
0824
0825 }  // namespace ActsTests