Core/Geometry/PortalLinkTests.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/context.hpp>
0011 #include <boost/test/tools/old/interface.hpp>
0012 #include <boost/test/unit_test.hpp>
0013 #include <boost/test/unit_test_suite.hpp>
0014
0015 #include "Acts/Definitions/Algebra.hpp"
0016 #include "Acts/Definitions/Units.hpp"
0017 #include "Acts/Geometry/CompositePortalLink.hpp"
0018 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
0019 #include "Acts/Geometry/GridPortalLink.hpp"
0020 #include "Acts/Geometry/TrackingVolume.hpp"
0021 #include "Acts/Geometry/TrivialPortalLink.hpp"
0022 #include "Acts/Surfaces/CylinderSurface.hpp"
0023 #include "Acts/Surfaces/RadialBounds.hpp"
0024 #include "Acts/Surfaces/SurfaceMergingException.hpp"
0025 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0026 #include "Acts/Utilities/AxisDefinitions.hpp"
0027 #include "Acts/Utilities/ThrowAssert.hpp"
0028
0029 #include <cstdio>
0030 #include <iostream>
0031 #include <memory>
0032 #include <numbers>
0033 #include <stdexcept>
0034
0035 using namespace Acts::UnitLiterals;
0036
0037 namespace Acts::Test {
0038
0039 auto logger = Acts::getDefaultLogger("UnitTests", Acts::Logging::VERBOSE);
0040
0041 struct Fixture {
0042   Logging::Level m_level;
0043   Fixture() {
0044     m_level = Acts::Logging::getFailureThreshold();
0045     Acts::Logging::setFailureThreshold(Acts::Logging::FATAL);
0046   }
0047
0048   ~Fixture() { Acts::Logging::setFailureThreshold(m_level); }
0049 };
0050
0051 std::shared_ptr<TrackingVolume> makeDummyVolume() {
0052   return std::make_shared<TrackingVolume>(
0053       Transform3::Identity(),
0054       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
0055 }
0056
0057 GeometryContext gctx;
0058
0059 template <typename T>
0060 std::unique_ptr<T> copy(const std::unique_ptr<T>& p) {
0061   return std::make_unique<T>(*p);
0062 }
0063
0064 template <typename link_t>
0065 void visitBins(const link_t& link,
0066                const std::function<void(const TrackingVolume*)>& func) {
0067   auto& grid = link.grid();
0068   auto loc = grid.numLocalBins();
0069   if constexpr (std::decay_t<decltype(grid)>::DIM == 1) {
0070     for (std::size_t i = 1; i <= loc[0]; i++) {
0071       func(grid.atLocalBins({i}));
0072     }
0073   } else {
0074     for (std::size_t i = 1; i <= loc[0]; i++) {
0075       for (std::size_t j = 1; j <= loc[1]; j++) {
0076         func(grid.atLocalBins({i, j}));
0077       }
0078     }
0079   }
0080 }
0081
0082 BOOST_FIXTURE_TEST_SUITE(Geometry, Fixture)
0083
0084 BOOST_AUTO_TEST_SUITE(GridConstruction)
0085
0086 BOOST_AUTO_TEST_CASE(Cylinder) {
0087   BOOST_TEST_CONTEXT("1D") {
0088     auto cyl = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
0089                                                     30_mm, 100_mm);
0090
0091     // Volume for bin testing
0092     auto vol = std::make_shared<TrackingVolume>(
0093         Transform3::Identity(),
0094         std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
0095
0096     // Incompatible binning
0097     BOOST_CHECK_THROW(GridPortalLink::make(cyl, AxisDirection::AxisZ,
0098                                            Axis{AxisBound, 0, 5, 5}),
0099                       std::invalid_argument);
0100
0101     auto grid1dCyl = GridPortalLink::make(cyl, AxisDirection::AxisZ,
0102                                           Axis{AxisBound, -100_mm, 100_mm, 10});
0103     BOOST_REQUIRE(grid1dCyl);
0104     grid1dCyl->setVolume(vol.get());
0105
0106     // Throws because non-closed axis
0107     BOOST_CHECK_THROW(GridPortalLink::make(cyl, AxisDirection::AxisRPhi,
0108                                            Axis{AxisBound, -180_degree * 30_mm,
0109                                                 180_degree * 30_mm, 10}),
0110                       std::invalid_argument);
0111
0112     auto grid1dCylRPhi = GridPortalLink::make(
0113         cyl, AxisDirection::AxisRPhi,
0114         Axis{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 10});
0115     BOOST_REQUIRE_NE(grid1dCylRPhi, nullptr);
0116     grid1dCylRPhi->setVolume(vol.get());
0117
0118     Axis axisExpected{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 10};
0119     BOOST_CHECK_EQUAL(grid1dCylRPhi->grid().axes().size(), 1);
0120     const auto& axis = *grid1dCylRPhi->grid().axes().front();
0121     BOOST_CHECK_EQUAL(axis, axisExpected);
0122
0123     // Another cylinder, shifted in z
0124     auto cyl2 = Surface::makeShared<CylinderSurface>(
0125         Transform3{Translation3{Vector3::UnitZ() * 150_mm}}, 30_mm, 50_mm);
0126
0127     auto grid1dCyl2 = GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0128                                            Axis{AxisBound, -50_mm, 50_mm, 5});
0129
0130     // Test exception on cylinder with non-zero average phi
0131     auto cylNonZeroAverage = Surface::makeShared<CylinderSurface>(
0132         Transform3::Identity(), 30_mm, 100_mm, 20_degree, 45_degree);
0133     BOOST_CHECK_THROW(
0134         GridPortalLink::make(cylNonZeroAverage, AxisDirection::AxisZ,
0135                              Axis{AxisBound, -100_mm, 100_mm, 10}),
0136         std::invalid_argument);
0137
0138     auto checkAllBins = [&](const auto& link) {
0139       visitBins(link, [&](const TrackingVolume* content) {
0140         BOOST_CHECK_EQUAL(content, vol.get());
0141       });
0142     };
0143
0144     checkAllBins(*grid1dCyl);
0145     checkAllBins(*grid1dCylRPhi);
0146
0147     // Extend to a 2D grid with auto phi binning
0148
0149     auto grid2dCyl1 = grid1dCyl->extendTo2d(nullptr);
0150     BOOST_REQUIRE(grid2dCyl1);
0151     BOOST_CHECK_EQUAL(grid2dCyl1->grid().axes().size(), 2);
0152     BOOST_CHECK_EQUAL(grid2dCyl1->surface().bounds(), cyl->bounds());
0153     const auto* axis1 = grid2dCyl1->grid().axes().front();
0154     const auto* axis2 = grid2dCyl1->grid().axes().back();
0155
0156     Axis axis1Expected{AxisClosed, -std::numbers::pi * 30_mm,
0157                        std::numbers::pi * 30_mm, 1};
0158     BOOST_CHECK_EQUAL(*axis1, axis1Expected);
0159     Axis axis2Expected{AxisBound, -100_mm, 100_mm, 10};
0160     BOOST_CHECK_EQUAL(*axis2, axis2Expected);
0161
0162     auto& concrete = dynamic_cast<
0163         GridPortalLinkT<decltype(axis1Expected), decltype(axis2Expected)>&>(
0164         *grid2dCyl1);
0165
0166     checkAllBins(concrete);
0167
0168     Axis axis1Explicit{AxisClosed, -std::numbers::pi * 30_mm,
0169                        std::numbers::pi * 30_mm, 13};
0170     auto grid2dCyl1Explicit = grid1dCyl->extendTo2d(&axis1Explicit);
0171     BOOST_REQUIRE(grid2dCyl1Explicit);
0172     BOOST_CHECK_EQUAL(grid2dCyl1Explicit->grid().axes().size(), 2);
0173     axis1 = grid2dCyl1Explicit->grid().axes().front();
0174     axis2 = grid2dCyl1Explicit->grid().axes().back();
0175
0176     BOOST_CHECK_EQUAL(*axis1, axis1Explicit);
0177     BOOST_CHECK_EQUAL(*axis2, axis2Expected);
0178
0179     auto& concrete2 = dynamic_cast<
0180         GridPortalLinkT<decltype(axis1Explicit), decltype(axis2Expected)>&>(
0181         *grid2dCyl1Explicit);
0182
0183     checkAllBins(concrete2);
0184
0185     auto cylPhi = Surface::makeShared<CylinderSurface>(
0186         Transform3::Identity(), 30_mm, 100_mm, 45_degree);
0187     std::unique_ptr<GridPortalLink> grid1dCylPhi = GridPortalLink::make(
0188         cylPhi, AxisDirection::AxisZ, Axis{AxisBound, -100_mm, 100_mm, 10});
0189
0190     grid1dCylPhi->setVolume(vol.get());
0191
0192     // Check that phi sector portal does not accept closed axis
0193     BOOST_CHECK_THROW(GridPortalLink::make(cylPhi, AxisDirection::AxisRPhi,
0194                                            Axis{AxisClosed, -45_degree * 30_mm,
0195                                                 45_degree * 30_mm, 10}),
0196                       std::invalid_argument);
0197
0198     auto grid2dCylPhi = grid1dCylPhi->extendTo2d(nullptr);
0199     BOOST_CHECK_EQUAL(grid2dCylPhi->grid().axes().size(), 2);
0200     BOOST_CHECK_EQUAL(grid2dCylPhi->surface().bounds(), cylPhi->bounds());
0201     const auto* axis1Phi = grid2dCylPhi->grid().axes().front();
0202     const auto* axis2Phi = grid2dCylPhi->grid().axes().back();
0203
0204     Axis axis1PhiExpected{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 1};
0205     BOOST_CHECK_EQUAL(*axis1Phi, axis1PhiExpected);
0206     Axis axis2PhiExpected{AxisBound, -100_mm, 100_mm, 10};
0207     BOOST_CHECK_EQUAL(*axis2Phi, axis2PhiExpected);
0208
0209     auto& concrete3 =
0210         dynamic_cast<GridPortalLinkT<decltype(axis1PhiExpected),
0211                                      decltype(axis2PhiExpected)>&>(
0212             *grid2dCylPhi);
0213
0214     checkAllBins(concrete3);
0215
0216     Axis axis1PhiExplicit{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 13};
0217     auto grid2dCylPhiExplicit = grid1dCylPhi->extendTo2d(&axis1PhiExplicit);
0218     BOOST_REQUIRE(grid2dCylPhiExplicit);
0219     BOOST_CHECK_EQUAL(grid2dCylPhiExplicit->grid().axes().size(), 2);
0220     axis1Phi = grid2dCylPhiExplicit->grid().axes().front();
0221     axis2Phi = grid2dCylPhiExplicit->grid().axes().back();
0222     BOOST_CHECK_EQUAL(*axis1Phi, axis1PhiExplicit);
0223     BOOST_CHECK_EQUAL(*axis2Phi, axis2PhiExpected);
0224
0225     auto& concrete4 =
0226         dynamic_cast<GridPortalLinkT<decltype(axis1PhiExplicit),
0227                                      decltype(axis2PhiExpected)>&>(
0228             *grid2dCylPhiExplicit);
0229
0230     checkAllBins(concrete4);
0231   }
0232
0233   BOOST_TEST_CONTEXT("2D") {
0234     auto cyl = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
0235                                                     30_mm, 100_mm, 45_degree);
0236
0237     // z bad, rphi bad
0238     BOOST_CHECK_THROW(GridPortalLink::make(cyl, Axis{AxisBound, 1, 2, 5},
0239                                            Axis{AxisBound, 3_mm, 4_mm, 5}),
0240                       std::invalid_argument);
0241
0242     // z good, rphi bad
0243     BOOST_CHECK_THROW(GridPortalLink::make(cyl, Axis{AxisBound, 3_mm, 4_mm, 5},
0244                                            Axis{AxisBound, -100_mm, 100_m, 5}),
0245                       std::invalid_argument);
0246
0247     // z bad, rphi good
0248     BOOST_CHECK_THROW(
0249         GridPortalLink::make(
0250             cyl, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5},
0251             Axis{AxisBound, -80_mm, 100_mm, 5}),
0252         std::invalid_argument);
0253
0254     auto grid1 = GridPortalLink::make(
0255         cyl, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5},
0256         Axis{AxisBound, -100_mm, 100_mm, 5});
0257
0258     auto cylFull = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
0259                                                         30_mm, 100_mm);
0260
0261     // Throws because non-closed axis
0262     BOOST_CHECK_THROW(GridPortalLink::make(cylFull,
0263                                            Axis{AxisBound, -180_degree * 30_mm,
0264                                                 180_degree * 30_mm, 5},
0265                                            Axis{AxisBound, -100_mm, 100_mm, 5}),
0266                       std::invalid_argument);
0267
0268     auto gridFull = GridPortalLink::make(
0269         cylFull, Axis{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 5},
0270         Axis{AxisBound, -100_mm, 100_mm, 5});
0271
0272     BOOST_CHECK_EQUAL(gridFull->grid().axes().size(), 2);
0273     BOOST_CHECK_EQUAL(gridFull->grid().axes().size(), 2);
0274     Axis axisFullExpected{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm,
0275                           5};
0276     const auto& axisFull = *gridFull->grid().axes().front();
0277     BOOST_CHECK_EQUAL(axisFull, axisFullExpected);
0278   }
0279 }
0280
0281 BOOST_AUTO_TEST_CASE(Disc) {
0282   using enum AxisType;
0283   using enum AxisBoundaryType;
0284   BOOST_TEST_CONTEXT("1D") {
0285     auto disc1 =
0286         Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 100_mm);
0287
0288     // Volume for bin testing
0289     auto vol = std::make_shared<TrackingVolume>(
0290         Transform3::Identity(),
0291         std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
0292
0293     BOOST_CHECK_THROW(GridPortalLink::make(disc1, AxisDirection::AxisZ,
0294                                            Axis{AxisBound, 30_mm, 100_mm, 3}),
0295                       std::invalid_argument);
0296
0297     // Check exception for full disc and non-closed phi axis
0298     BOOST_CHECK_THROW(
0299         GridPortalLink::make(disc1, AxisDirection::AxisPhi,
0300                              Axis{AxisBound, -180_degree, 180_degree, 3}),
0301         std::invalid_argument);
0302
0303     auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisR,
0304                                       Axis{AxisBound, 30_mm, 100_mm, 3});
0305     BOOST_REQUIRE_NE(grid1, nullptr);
0306     BOOST_CHECK_EQUAL(grid1->grid().axes().size(), 1);
0307     const auto& axis = *grid1->grid().axes().front();
0308     Axis axis1Expected{AxisBound, 30_mm, 100_mm, 3};
0309     BOOST_CHECK_EQUAL(axis, axis1Expected);
0310
0311     grid1->setVolume(vol.get());
0312
0313     auto discPhi = Surface::makeShared<DiscSurface>(Transform3::Identity(),
0314                                                     30_mm, 100_mm, 45_degree);
0315
0316     // Check thet disc with phi sector does not accept closed axis
0317     BOOST_CHECK_THROW(
0318         GridPortalLink::make(discPhi, AxisDirection::AxisPhi,
0319                              Axis{AxisClosed, -45_degree, 45_degree, 3}),
0320         std::invalid_argument);
0321
0322     auto gridPhi =
0323         GridPortalLink::make(discPhi, AxisDirection::AxisPhi,
0324                              Axis{AxisBound, -45_degree, 45_degree, 3});
0325     BOOST_REQUIRE_NE(gridPhi, nullptr);
0326     gridPhi->setVolume(vol.get());
0327
0328     // Test exception on disc with non-zero average phi
0329     auto discNonZeroAverage = Surface::makeShared<DiscSurface>(
0330         Transform3::Identity(),
0331         std::make_shared<RadialBounds>(30_mm, 100_mm, 45_degree, 75_degree));
0332     BOOST_CHECK_THROW(
0333         GridPortalLink::make(discNonZeroAverage, AxisDirection::AxisR,
0334                              Axis{AxisBound, 30_mm, 100_mm, 3}),
0335         std::invalid_argument);
0336
0337     BOOST_CHECK_EQUAL(gridPhi->grid().axes().size(), 1);
0338     const auto& axisPhi = *gridPhi->grid().axes().front();
0339     Axis axisPhi1Expected{AxisBound, -45_degree, 45_degree, 3};
0340     BOOST_CHECK_EQUAL(axisPhi, axisPhi1Expected);
0341
0342     auto checkAllBins = [&](const auto& grid) {
0343       visitBins(grid, [&](const TrackingVolume* content) {
0344         BOOST_CHECK_EQUAL(content, vol.get());
0345       });
0346     };
0347
0348     checkAllBins(*grid1);
0349     checkAllBins(*gridPhi);
0350
0351     // Test making 2D grids from the 1D ones
0352     auto grid2d = grid1->extendTo2d(nullptr);
0353     BOOST_REQUIRE(grid2d);
0354     BOOST_CHECK_EQUAL(grid2d->grid().axes().size(), 2);
0355     const auto* axis1 = grid2d->grid().axes().front();
0356     BOOST_CHECK_EQUAL(*axis1, axis1Expected);
0357     const auto* axis2 = grid2d->grid().axes().back();
0358     BOOST_CHECK_CLOSE(axis2->getMin(), -180_degree, 1e-6);
0359     BOOST_CHECK_CLOSE(axis2->getMax(), 180_degree, 1e-6);
0360     BOOST_CHECK_EQUAL(axis2->getNBins(), 1);
0361     BOOST_CHECK_EQUAL(axis2->getType(), AxisType::Equidistant);
0362     BOOST_CHECK_EQUAL(axis2->getBoundaryType(), AxisBoundaryType::Closed);
0363
0364     checkAllBins(
0365         dynamic_cast<GridPortalLinkT<decltype(axisPhi1Expected),
0366                                      Axis<Equidistant, Closed>>&>(*grid2d));
0367
0368     Axis axis2Explicit{AxisClosed, -180_degree, 180_degree, 3};
0369     auto grid2dExplicit = grid1->extendTo2d(&axis2Explicit);
0370     BOOST_REQUIRE(grid2dExplicit);
0371     BOOST_CHECK_EQUAL(grid2dExplicit->grid().axes().size(), 2);
0372     axis1 = grid2dExplicit->grid().axes().front();
0373     axis2 = grid2dExplicit->grid().axes().back();
0374     BOOST_CHECK_EQUAL(*axis1, axis1Expected);
0375     BOOST_CHECK_EQUAL(*axis2, axis2Explicit);
0376
0377     checkAllBins(dynamic_cast<GridPortalLinkT<decltype(axisPhi1Expected),
0378                                               decltype(axis2Explicit)>&>(
0379         *grid2dExplicit));
0380
0381     auto gridPhiBinnedInR = GridPortalLink::make(
0382         discPhi, AxisDirection::AxisR, Axis{AxisBound, 30_mm, 100_mm, 3});
0383     gridPhiBinnedInR->setVolume(vol.get());
0384     auto grid2dPhiNonClosed = gridPhiBinnedInR->extendTo2d(nullptr);
0385     BOOST_REQUIRE(grid2dPhiNonClosed);
0386     BOOST_CHECK_EQUAL(grid2dPhiNonClosed->grid().axes().size(), 2);
0387     Axis gridPhiBinnedInRExpected{AxisBound, 30_mm, 100_mm, 3};
0388     BOOST_CHECK_EQUAL(*grid2dPhiNonClosed->grid().axes().front(),
0389                       gridPhiBinnedInRExpected);
0390     const auto* axisPhiNonClosed = grid2dPhiNonClosed->grid().axes().back();
0391     BOOST_CHECK_CLOSE(axisPhiNonClosed->getMin(), -45_degree, 1e-6);
0392     BOOST_CHECK_CLOSE(axisPhiNonClosed->getMax(), 45_degree, 1e-6);
0393     BOOST_CHECK_EQUAL(axisPhiNonClosed->getNBins(), 1);
0394     BOOST_CHECK_EQUAL(axisPhiNonClosed->getType(), AxisType::Equidistant);
0395     BOOST_CHECK_EQUAL(axisPhiNonClosed->getBoundaryType(),
0396                       AxisBoundaryType::Bound);
0397
0398     checkAllBins(
0399         dynamic_cast<GridPortalLinkT<decltype(gridPhiBinnedInRExpected),
0400                                      Axis<Equidistant, Bound>>&>(
0401             *grid2dPhiNonClosed));
0402
0403     Axis axisPhiNonClosedExplicit{AxisBound, -45_degree, 45_degree, 3};
0404     auto grid2dPhiNonClosedExplicit =
0405         gridPhiBinnedInR->extendTo2d(&axisPhiNonClosedExplicit);
0406     BOOST_REQUIRE(grid2dPhiNonClosedExplicit);
0407     BOOST_CHECK_EQUAL(grid2dPhiNonClosedExplicit->grid().axes().size(), 2);
0408     axisPhiNonClosed = grid2dPhiNonClosedExplicit->grid().axes().back();
0409     BOOST_CHECK_EQUAL(*axisPhiNonClosed, axisPhiNonClosedExplicit);
0410     BOOST_CHECK_EQUAL(*grid2dPhiNonClosedExplicit->grid().axes().front(),
0411                       gridPhiBinnedInRExpected);
0412
0413     checkAllBins(
0414         dynamic_cast<GridPortalLinkT<decltype(gridPhiBinnedInRExpected),
0415                                      decltype(axisPhiNonClosedExplicit)>&>(
0416             *grid2dPhiNonClosedExplicit));
0417
0418     auto grid2dPhi = gridPhi->extendTo2d(nullptr);
0419     BOOST_REQUIRE(grid2dPhi);
0420     BOOST_CHECK_EQUAL(grid2dPhi->grid().axes().size(), 2);
0421     Axis axis2dPhiExpected{AxisBound, 30_mm, 100_mm, 1};
0422     BOOST_CHECK_EQUAL(*grid2dPhi->grid().axes().front(), axis2dPhiExpected);
0423     BOOST_CHECK_EQUAL(*grid2dPhi->grid().axes().back(), axisPhi1Expected);
0424
0425     checkAllBins(
0426         dynamic_cast<GridPortalLinkT<decltype(axis2dPhiExpected),
0427                                      decltype(axisPhi1Expected)>&>(*grid2dPhi));
0428
0429     Axis axis2dPhiExplicit{AxisBound, 30_mm, 100_mm, 3};
0430     auto grid2dPhiExplicit = gridPhi->extendTo2d(&axis2dPhiExplicit);
0431     BOOST_REQUIRE(grid2dPhiExplicit);
0432     BOOST_CHECK_EQUAL(grid2dPhiExplicit->grid().axes().size(), 2);
0433     BOOST_CHECK_EQUAL(*grid2dPhiExplicit->grid().axes().front(),
0434                       axis2dPhiExplicit);
0435     BOOST_CHECK_EQUAL(*grid2dPhiExplicit->grid().axes().back(),
0436                       axisPhi1Expected);
0437
0438     checkAllBins(dynamic_cast<GridPortalLinkT<decltype(axisPhi1Expected),
0439                                               decltype(axis2dPhiExplicit)>&>(
0440         *grid2dPhiExplicit));
0441   }
0442
0443   BOOST_TEST_CONTEXT("2D") {
0444     auto discPhi = Surface::makeShared<DiscSurface>(Transform3::Identity(),
0445                                                     30_mm, 100_mm, 45_degree);
0446
0447     Axis rBad{AxisBound, 1, 2, 5};
0448     Axis rGood{AxisBound, 30_mm, 100_mm, 5};
0449     Axis phiBad{AxisBound, 1, 2, 5};
0450     Axis phiGood{AxisBound, -45_degree, 45_degree, 5};
0451
0452     // r bad, phi bad
0453     BOOST_CHECK_THROW(GridPortalLink::make(discPhi, rBad, phiBad),
0454                       std::invalid_argument);
0455     // r bad, phi good
0456     BOOST_CHECK_THROW(GridPortalLink::make(discPhi, rBad, phiGood),
0457                       std::invalid_argument);
0458     // r good, phi bad
0459     BOOST_CHECK_THROW(GridPortalLink::make(discPhi, rGood, phiBad),
0460                       std::invalid_argument);
0461     // r good phi good
0462     auto grid = GridPortalLink::make(discPhi, rGood, phiGood);
0463     BOOST_REQUIRE_NE(grid, nullptr);
0464   }
0465 }
0466
0467 BOOST_AUTO_TEST_CASE(Plane) {
0468   // @TODO: Add plane tests
0469 }
0470
0471 BOOST_AUTO_TEST_CASE(FromTrivial) {
0472   BOOST_TEST_CONTEXT("Cylinder") {
0473     auto cyl = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
0474                                                     30_mm, 100_mm);
0475
0476     auto vol = std::make_shared<TrackingVolume>(
0477         Transform3::Identity(),
0478         std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
0479
0480     auto trivial = std::make_unique<TrivialPortalLink>(cyl, *vol);
0481     BOOST_REQUIRE(trivial);
0482
0483     BOOST_CHECK_EQUAL(trivial->resolveVolume(gctx, Vector2{1, 2}).value(),
0484                       vol.get());
0485
0486     auto gridZ = trivial->makeGrid(AxisDirection::AxisZ);
0487     BOOST_REQUIRE(gridZ);
0488
0489     BOOST_CHECK_EQUAL(gridZ->grid().axes().size(), 1);
0490     BOOST_CHECK_EQUAL(gridZ->surface().bounds(), cyl->bounds());
0491     Axis axisZExpected{AxisBound, -100_mm, 100_mm, 1};
0492     BOOST_CHECK_EQUAL(*gridZ->grid().axes().front(), axisZExpected);
0493
0494     BOOST_CHECK_EQUAL(
0495         gridZ->resolveVolume(gctx, Vector2{20_degree * 30_mm, 90_mm}).value(),
0496         vol.get());
0497
0498     auto gridRPhi = trivial->makeGrid(AxisDirection::AxisRPhi);
0499     BOOST_REQUIRE(gridRPhi);
0500
0501     BOOST_CHECK_EQUAL(gridRPhi->grid().axes().size(), 1);
0502     BOOST_CHECK_EQUAL(gridRPhi->surface().bounds(), cyl->bounds());
0503     Axis axisRPhiExpected{AxisClosed, -std::numbers::pi * 30_mm,
0504                           std::numbers::pi * 30_mm, 1};
0505     BOOST_CHECK_EQUAL(*gridRPhi->grid().axes().front(), axisRPhiExpected);
0506
0507     auto cylPhi = Surface::makeShared<CylinderSurface>(
0508         Transform3::Identity(), 30_mm, 100_mm, 30_degree);
0509
0510     auto trivialPhi = std::make_unique<TrivialPortalLink>(cylPhi, *vol);
0511     BOOST_REQUIRE(trivialPhi);
0512
0513     BOOST_CHECK_EQUAL(trivialPhi->resolveVolume(gctx, Vector2{1, 2}).value(),
0514                       vol.get());
0515
0516     auto gridRPhiSector = trivialPhi->makeGrid(AxisDirection::AxisRPhi);
0517     BOOST_REQUIRE(gridRPhiSector);
0518
0519     BOOST_CHECK_EQUAL(
0520         gridRPhiSector->resolveVolume(gctx, Vector2{20_degree * 30_mm, 90_mm})
0521             .value(),
0522         vol.get());
0523
0524     BOOST_CHECK_EQUAL(gridRPhiSector->grid().axes().size(), 1);
0525     Axis axisRPhiSectorExpected{AxisBound, -30_degree * 30_mm,
0526                                 30_degree * 30_mm, 1};
0527     BOOST_CHECK_EQUAL(*gridRPhiSector->grid().axes().front(),
0528                       axisRPhiSectorExpected);
0529   }
0530
0531   BOOST_TEST_CONTEXT("Disc") {
0532     auto disc =
0533         Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 100_mm);
0534
0535     auto vol = std::make_shared<TrackingVolume>(
0536         Transform3::Identity(),
0537         std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
0538
0539     auto trivial = std::make_unique<TrivialPortalLink>(disc, *vol);
0540     BOOST_REQUIRE(trivial);
0541
0542     // Doesn't matter which position
0543     BOOST_CHECK_EQUAL(trivial->resolveVolume(gctx, Vector2{1, 2}).value(),
0544                       vol.get());
0545
0546     auto gridR = trivial->makeGrid(AxisDirection::AxisR);
0547     BOOST_REQUIRE(gridR);
0548
0549     BOOST_CHECK_EQUAL(gridR->grid().axes().size(), 1);
0550     BOOST_CHECK_EQUAL(gridR->surface().bounds(), disc->bounds());
0551     Axis axisRExpected{AxisBound, 30_mm, 100_mm, 1};
0552     BOOST_CHECK_EQUAL(*gridR->grid().axes().front(), axisRExpected);
0553
0554     BOOST_CHECK_EQUAL(
0555         gridR->resolveVolume(gctx, Vector2{90_mm, 10_degree}).value(),
0556         vol.get());
0557
0558     auto gridPhi = trivial->makeGrid(AxisDirection::AxisPhi);
0559     BOOST_REQUIRE(gridPhi);
0560
0561     BOOST_CHECK_EQUAL(gridPhi->grid().axes().size(), 1);
0562     BOOST_CHECK_EQUAL(gridPhi->surface().bounds(), disc->bounds());
0563     Axis axisPhiExpected{AxisClosed, -std::numbers::pi, std::numbers::pi, 1};
0564     BOOST_CHECK_EQUAL(*gridPhi->grid().axes().front(), axisPhiExpected);
0565
0566     BOOST_CHECK_EQUAL(
0567         gridPhi->resolveVolume(gctx, Vector2{90_mm, 10_degree}).value(),
0568         vol.get());
0569   }
0570 }
0571
0572 BOOST_AUTO_TEST_SUITE_END()  // GridConstruction
0573
0574 BOOST_AUTO_TEST_SUITE(GridMerging)
0575
0576 BOOST_AUTO_TEST_SUITE(Merging1dCylinder)
0577
0578 BOOST_AUTO_TEST_SUITE(ZDirection)
0579
0580 BOOST_AUTO_TEST_CASE(ColinearMerge) {
0581   auto vol1 = makeDummyVolume();
0582   auto vol2 = makeDummyVolume();
0583
0584   auto cyl = Surface::makeShared<CylinderSurface>(Transform3::Identity(), 30_mm,
0585                                                   100_mm);
0586
0587   auto grid1dCyl = GridPortalLink::make(cyl, AxisDirection::AxisZ,
0588                                         Axis{AxisBound, -100_mm, 100_mm, 10});
0589   grid1dCyl->setVolume(vol1.get());
0590
0591   // Another cylinder, shifted in z
0592   auto cyl2 = Surface::makeShared<CylinderSurface>(
0593       Transform3{Translation3{Vector3::UnitZ() * 150_mm}}, 30_mm, 50_mm);
0594
0595   auto grid1dCyl2 = GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0596                                          Axis{AxisBound, -50_mm, 50_mm, 5});
0597
0598   grid1dCyl2->setVolume(vol2.get());
0599
0600   // Completely invalid
0601   BOOST_CHECK_THROW(GridPortalLink::merge(*grid1dCyl, *grid1dCyl2,
0602                                           AxisDirection::AxisPhi, *logger),
0603                     AssertionFailureException);
0604   // Invalid direction, as the cylinders are shifted in z, and can't be merged
0605   // in r x phi
0606   BOOST_CHECK_THROW(GridPortalLink::merge(*grid1dCyl, *grid1dCyl2,
0607                                           AxisDirection::AxisRPhi, *logger),
0608                     SurfaceMergingException);
0609
0610   BOOST_TEST_CONTEXT("Consistent equidistant") {
0611     auto mergedPtr = GridPortalLink::merge(*grid1dCyl, *grid1dCyl2,
0612                                            AxisDirection::AxisZ, *logger);
0613     BOOST_REQUIRE_NE(mergedPtr, nullptr);
0614     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0615     BOOST_REQUIRE_NE(merged, nullptr);
0616     BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0617     const auto& axis = *merged->grid().axes().front();
0618     BOOST_CHECK_EQUAL(axis.getMin(), -150_mm);
0619     BOOST_CHECK_EQUAL(axis.getMax(), 150_mm);
0620     BOOST_CHECK_EQUAL(axis.getNBins(), 15);
0621     BOOST_CHECK_EQUAL(axis.getType(), AxisType::Equidistant);
0622     BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0623   }
0624
0625   BOOST_TEST_CONTEXT("Inconsistent equidistant") {
0626     auto grid1dCyl2BinWidthChanged = GridPortalLink::make(
0627         cyl2, AxisDirection::AxisZ, Axis{AxisBound, -50_mm, 50_mm, 6});
0628
0629     auto mergedPtr = GridPortalLink::merge(
0630         *grid1dCyl, *grid1dCyl2BinWidthChanged, AxisDirection::AxisZ, *logger);
0631     BOOST_REQUIRE_NE(mergedPtr, nullptr);
0632     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0633     BOOST_REQUIRE_NE(merged, nullptr);
0634     BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0635     const auto& axis = *merged->grid().axes().front();
0636     BOOST_CHECK_EQUAL(axis.getMin(), -150_mm);
0637     BOOST_CHECK_EQUAL(axis.getMax(), 150_mm);
0638     BOOST_CHECK_EQUAL(axis.getNBins(), 16);
0639     BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
0640     BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0641   }
0642
0643   BOOST_TEST_CONTEXT("Right Variable") {
0644     auto gridLeft = GridPortalLink::make(cyl, AxisDirection::AxisZ,
0645                                          Axis{AxisBound, -100_mm, 100_mm, 10});
0646
0647     auto gridRight =
0648         GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0649                              Axis{AxisBound, {-50_mm, -10_mm, 10_mm, 50_mm}});
0650
0651     auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight,
0652                                            AxisDirection::AxisZ, *logger);
0653     BOOST_REQUIRE_NE(mergedPtr, nullptr);
0654     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0655     BOOST_REQUIRE_NE(merged, nullptr);
0656     BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0657     const auto& axis = *merged->grid().axes().front();
0658     BOOST_CHECK_EQUAL(axis.getMin(), -150_mm);
0659     BOOST_CHECK_EQUAL(axis.getMax(), 150_mm);
0660     BOOST_CHECK_EQUAL(axis.getNBins(), 13);
0661     BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
0662     BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0663   }
0664
0665   BOOST_TEST_CONTEXT("Left Variable") {
0666     auto gridLeft =
0667         GridPortalLink::make(cyl, AxisDirection::AxisZ,
0668                              Axis{AxisBound, {-100_mm, -80_mm, 10_mm, 100_mm}});
0669
0670     auto gridRight = GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0671                                           Axis{AxisBound, -50_mm, 50_mm, 8});
0672
0673     auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight,
0674                                            AxisDirection::AxisZ, *logger);
0675     BOOST_REQUIRE_NE(mergedPtr, nullptr);
0676     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0677     BOOST_REQUIRE_NE(merged, nullptr);
0678     BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0679     const auto& axis = *merged->grid().axes().front();
0680     BOOST_CHECK_EQUAL(axis.getMin(), -150_mm);
0681     BOOST_CHECK_EQUAL(axis.getMax(), 150_mm);
0682     BOOST_CHECK_EQUAL(axis.getNBins(), 11);
0683     BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
0684     BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0685   }
0686
0687   BOOST_TEST_CONTEXT("Both Variable") {
0688     auto gridLeft =
0689         GridPortalLink::make(cyl, AxisDirection::AxisZ,
0690                              Axis{AxisBound, {-100_mm, -80_mm, 10_mm, 100_mm}});
0691
0692     auto gridRight =
0693         GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0694                              Axis{AxisBound, {-50_mm, -10_mm, 10_mm, 50_mm}});
0695
0696     auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight,
0697                                            AxisDirection::AxisZ, *logger);
0698     BOOST_REQUIRE_NE(mergedPtr, nullptr);
0699     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0700     BOOST_REQUIRE_NE(merged, nullptr);
0701     BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0702     const auto& axis = *merged->grid().axes().front();
0703     BOOST_CHECK_EQUAL(axis.getMin(), -150_mm);
0704     BOOST_CHECK_EQUAL(axis.getMax(), 150_mm);
0705     BOOST_CHECK_EQUAL(axis.getNBins(), 6);
0706     BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
0707     BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0708   }
0709
0710   BOOST_TEST_CONTEXT("Non bound axis") {
0711     std::unique_ptr gridLeft =
0712         GridPortalLink::make(cyl, AxisDirection::AxisZ,
0713                              Axis{AxisBound, {-100_mm, -80_mm, 10_mm, 100_mm}});
0714     std::unique_ptr gridRightClosed =
0715         GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0716                              Axis{AxisClosed, {-50_mm, -10_mm, 10_mm, 50_mm}});
0717     std::unique_ptr gridRightOpen =
0718         GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0719                              Axis{AxisOpen, {-50_mm, -10_mm, 10_mm, 50_mm}});
0720
0721     auto compositeLR = PortalLinkBase::merge(
0722         copy(gridLeft), copy(gridRightClosed), AxisDirection::AxisZ, *logger);
0723     BOOST_CHECK_NE(dynamic_cast<const CompositePortalLink*>(compositeLR.get()),
0724                    nullptr);
0725     auto compositeRL = PortalLinkBase::merge(
0726         copy(gridLeft), copy(gridRightOpen), AxisDirection::AxisZ, *logger);
0727     BOOST_CHECK_NE(dynamic_cast<const CompositePortalLink*>(compositeRL.get()),
0728                    nullptr);
0729   }
0730 }
0731
0732 BOOST_AUTO_TEST_CASE(ParallelMerge) {
0733   // Merge in z direction with phi sectors
0734   auto cyl1 = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
0735                                                    30_mm, 100_mm, 35_degree);
0736
0737   auto grid1 = GridPortalLink::make(
0738       cyl1, AxisDirection::AxisRPhi,
0739
0740       Axis{AxisBound, -35_degree * 30_mm, 35_degree * 30_mm, 3});
0741   auto cyl2 = Surface::makeShared<CylinderSurface>(
0742       Transform3{Translation3{Vector3::UnitZ() * 150_mm}}, 30_mm, 50_mm,
0743       35_degree);
0744
0745   auto grid2 = GridPortalLink::make(
0746       cyl2, AxisDirection::AxisRPhi,
0747       Axis{AxisBound, -35_degree * 30_mm, 35_degree * 30_mm, 3});
0748
0749   auto merged12Ptr =
0750       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisZ, *logger);
0751   BOOST_REQUIRE_NE(merged12Ptr, nullptr);
0752   auto merged12 = dynamic_cast<const GridPortalLink*>(merged12Ptr.get());
0753   BOOST_REQUIRE_NE(merged12, nullptr);
0754
0755   BOOST_REQUIRE_EQUAL(merged12->grid().axes().size(), 2);
0756
0757   const auto& axis1 = *merged12->grid().axes().front();
0758   const auto& axis2 = *merged12->grid().axes().back();
0759   Axis axis1Expected{AxisBound, -35_degree * 30_mm, 35_degree * 30_mm, 3};
0760   BOOST_CHECK_EQUAL(axis1, axis1Expected);
0761   Axis axis2Expected{AxisBound, {-150_mm, 50_mm, 150_mm}};
0762   BOOST_CHECK_EQUAL(axis2, axis2Expected);
0763 }
0764
0765 BOOST_AUTO_TEST_SUITE_END()  // ZDirection
0766
0767 BOOST_AUTO_TEST_SUITE(RPhiDirection)
0768
0769 BOOST_AUTO_TEST_CASE(ColinearMerge) {
0770   auto cyl = Surface::makeShared<CylinderSurface>(Transform3::Identity(), 30_mm,
0771                                                   100_mm);
0772   BOOST_CHECK_THROW(GridPortalLink::make(cyl, AxisDirection::AxisRPhi,
0773                                          Axis{AxisBound, 0, 5, 5}),
0774                     std::invalid_argument);
0775
0776   auto cylNonZeroAverage = Surface::makeShared<CylinderSurface>(
0777       Transform3::Identity(), 30_mm, 100_mm, 20_degree, 45_degree);
0778
0779   BOOST_CHECK_THROW(
0780       GridPortalLink::make(
0781           cylNonZeroAverage, AxisDirection::AxisRPhi,
0782           Axis{AxisBound, -20_degree * 30_mm, 20_degree * 30_mm, 5}),
0783       std::invalid_argument);
0784
0785   BOOST_TEST_CONTEXT("Colinear merge in rPhi") {
0786     auto cylPhi1 = Surface::makeShared<CylinderSurface>(
0787         Transform3::Identity() * AngleAxis3(45_degree, Vector3::UnitZ()), 30_mm,
0788         100_mm, 20_degree, 0_degree);
0789
0790     auto cylPhi2 = Surface::makeShared<CylinderSurface>(
0791         Transform3::Identity() * AngleAxis3(105_degree, Vector3::UnitZ()),
0792         30_mm, 100_mm, 40_degree, 0_degree);
0793
0794     auto portalPhi1 = GridPortalLink::make(
0795         cylPhi1, AxisDirection::AxisRPhi,
0796         Axis{AxisBound, -20_degree * 30_mm, 20_degree * 30_mm, 5});
0797
0798     auto portalPhi2 = GridPortalLink::make(
0799         cylPhi2, AxisDirection::AxisRPhi,
0800         Axis{AxisBound, -40_degree * 30_mm, 40_degree * 30_mm, 10});
0801
0802     auto cylPhi3 = Surface::makeShared<CylinderSurface>(
0803         Transform3::Identity() * AngleAxis3(45_degree, Vector3::UnitZ()), 30_mm,
0804         100_mm, 90_degree, 0_degree);
0805
0806     auto cylPhi4 = Surface::makeShared<CylinderSurface>(
0807         Transform3::Identity() * AngleAxis3(-135_degree, Vector3::UnitZ()),
0808         30_mm, 100_mm, 90_degree, 0_degree);
0809
0810     auto portalPhi3 = GridPortalLink::make(
0811         cylPhi3, AxisDirection::AxisRPhi,
0812         Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 2});
0813
0814     auto portalPhi4 = GridPortalLink::make(
0815         cylPhi4, AxisDirection::AxisRPhi,
0816         Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 2});
0817
0818     BOOST_TEST_CONTEXT("Consistent equidistant") {
0819       auto portalMerged = GridPortalLink::merge(
0820           *portalPhi1, *portalPhi2, AxisDirection::AxisRPhi, *logger);
0821       BOOST_REQUIRE_NE(portalMerged, nullptr);
0822
0823       const auto* merged =
0824           dynamic_cast<const GridPortalLink*>(portalMerged.get());
0825       BOOST_REQUIRE_NE(merged, nullptr);
0826       BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0827       const auto& axis = *merged->grid().axes().front();
0828       BOOST_CHECK_CLOSE(axis.getMin(), -60_degree * 30_mm, 1e-9);
0829       BOOST_CHECK_CLOSE(axis.getMax(), 60_degree * 30_mm, 1e-9);
0830       BOOST_CHECK_EQUAL(axis.getNBins(), 15);
0831       BOOST_CHECK_EQUAL(axis.getType(), AxisType::Equidistant);
0832       BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0833
0834       // Test that if you merge half-circles, we get a closed axis
0835       auto portalMerged34 = GridPortalLink::merge(
0836           *portalPhi3, *portalPhi4, AxisDirection::AxisRPhi, *logger);
0837       BOOST_REQUIRE_NE(portalMerged34, nullptr);
0838
0839       const auto* merged34 =
0840           dynamic_cast<const GridPortalLink*>(portalMerged34.get());
0841       BOOST_REQUIRE_NE(merged34, nullptr);
0842       BOOST_CHECK_EQUAL(merged34->grid().axes().size(), 1);
0843       const auto& axis34 = *merged34->grid().axes().front();
0844       BOOST_CHECK_CLOSE(axis34.getMin(), -180_degree * 30_mm, 1e-9);
0845       BOOST_CHECK_CLOSE(axis34.getMax(), 180_degree * 30_mm, 1e-9);
0846       BOOST_CHECK_EQUAL(axis34.getNBins(), 4);
0847       BOOST_CHECK_EQUAL(axis34.getType(), AxisType::Equidistant);
0848       BOOST_CHECK_EQUAL(axis34.getBoundaryType(), AxisBoundaryType::Closed);
0849     }
0850
0851     BOOST_TEST_CONTEXT("Inconsistent equidistant") {
0852       auto portalPhi2Mod = GridPortalLink::make(
0853           cylPhi2, AxisDirection::AxisRPhi,
0854           Axis{AxisBound, -40_degree * 30_mm, 40_degree * 30_mm, 3});
0855
0856       auto portalMergedMod = GridPortalLink::merge(
0857           *portalPhi1, *portalPhi2Mod, AxisDirection::AxisRPhi, *logger);
0858       BOOST_REQUIRE_NE(portalMergedMod, nullptr);
0859
0860       const auto* merged12 =
0861           dynamic_cast<const GridPortalLink*>(portalMergedMod.get());
0862       BOOST_REQUIRE_NE(merged12, nullptr);
0863       BOOST_CHECK_EQUAL(merged12->grid().axes().size(), 1);
0864       const auto& axis12 = *merged12->grid().axes().front();
0865       BOOST_CHECK_CLOSE(axis12.getMin(), -60_degree * 30_mm, 1e-9);
0866       BOOST_CHECK_CLOSE(axis12.getMax(), 60_degree * 30_mm, 1e-9);
0867       BOOST_CHECK_EQUAL(axis12.getNBins(), 8);
0868       BOOST_CHECK_EQUAL(axis12.getType(), AxisType::Variable);
0869       BOOST_CHECK_EQUAL(axis12.getBoundaryType(), AxisBoundaryType::Bound);
0870
0871       std::vector<double> expected12 = {-31.4159, -17.4533, -3.49066,
0872                                         10.472,   14.6608,  18.8496,
0873                                         23.0383,  27.2271,  31.4159};
0874       CHECK_CLOSE_OR_SMALL(axis12.getBinEdges(), expected12, 1e-4, 10e-10);
0875
0876       auto portalPhi4Mod = GridPortalLink::make(
0877           cylPhi4, AxisDirection::AxisRPhi,
0878           Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 1});
0879
0880       auto portalMerged34 = GridPortalLink::merge(
0881           *portalPhi3, *portalPhi4Mod, AxisDirection::AxisRPhi, *logger);
0882       BOOST_REQUIRE_NE(portalMerged34, nullptr);
0883
0884       const auto* merged34 =
0885           dynamic_cast<const GridPortalLink*>(portalMerged34.get());
0886       BOOST_REQUIRE_NE(merged34, nullptr);
0887       BOOST_CHECK_EQUAL(merged34->grid().axes().size(), 1);
0888       const auto& axis34 = *merged34->grid().axes().front();
0889       BOOST_CHECK_CLOSE(axis34.getMin(), -180_degree * 30_mm, 1e-9);
0890       BOOST_CHECK_CLOSE(axis34.getMax(), 180_degree * 30_mm, 1e-9);
0891       BOOST_CHECK_EQUAL(axis34.getNBins(), 3);
0892       BOOST_CHECK_EQUAL(axis34.getType(), AxisType::Variable);
0893       BOOST_CHECK_EQUAL(axis34.getBoundaryType(), AxisBoundaryType::Closed);
0894
0895       // Caution: for full-azimuth cases, the ordering is preserved, you get
0896       // in what you get out. -> this can flip
0897       std::vector<double> expected34 = {-94.2478, -47.1239, 0, 94.2478};
0898       CHECK_CLOSE_OR_SMALL(axis34.getBinEdges(), expected34, 1e-4, 10e-10);
0899     }
0900
0901     BOOST_TEST_CONTEXT("Left variable") {
0902       BOOST_TEST_CONTEXT("Non-closed") {
0903         auto gridLeft =
0904             GridPortalLink::make(cylPhi1, AxisDirection::AxisRPhi,
0905                                  Axis{AxisBound,
0906                                       {-20_degree * 30_mm, -10_degree * 30_mm,
0907                                        10_degree * 30_mm, 20_degree * 30_mm}});
0908         auto gridRight = GridPortalLink::make(
0909             cylPhi2, AxisDirection::AxisRPhi,
0910             Axis{AxisBound, -40_degree * 30_mm, 40_degree * 30_mm, 3});
0911
0912         auto mergedPtr = GridPortalLink::merge(
0913             *gridLeft, *gridRight, AxisDirection::AxisRPhi, *logger);
0914         BOOST_REQUIRE_NE(mergedPtr, nullptr);
0915
0916         const auto* merged =
0917             dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0918         BOOST_REQUIRE_NE(merged, nullptr);
0919         const auto& axis = *merged->grid().axes().front();
0920         BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0921         BOOST_CHECK_CLOSE(axis.getMin(), -60_degree * 30_mm, 1e-9);
0922         BOOST_CHECK_CLOSE(axis.getMax(), 60_degree * 30_mm, 1e-9);
0923         BOOST_CHECK_EQUAL(axis.getNBins(), 6);
0924         BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
0925         BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0926
0927         std::vector<double> expected = {-31.4159, -17.4533, -3.49066, 10.472,
0928                                         15.708,   26.1799,  31.4159};
0929         CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10);
0930       }
0931
0932       BOOST_TEST_CONTEXT("Closed") {
0933         auto gridLeft = GridPortalLink::make(
0934             cylPhi4, AxisDirection::AxisRPhi,
0935             Axis{AxisBound,
0936                  {-90_degree * 30_mm, 25_degree * 30_mm, 90_degree * 30_mm}});
0937
0938         auto gridRight = GridPortalLink::make(
0939             cylPhi3, AxisDirection::AxisRPhi,
0940             Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 3});
0941
0942         auto mergedPtr = GridPortalLink::merge(
0943             *gridLeft, *gridRight, AxisDirection::AxisRPhi, *logger);
0944         BOOST_REQUIRE_NE(mergedPtr, nullptr);
0945
0946         const auto* merged =
0947             dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0948         BOOST_REQUIRE_NE(merged, nullptr);
0949         BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0950         const auto& axis = *merged->grid().axes().front();
0951         BOOST_CHECK_CLOSE(axis.getMin(), -180_degree * 30_mm, 1e-9);
0952         BOOST_CHECK_CLOSE(axis.getMax(), 180_degree * 30_mm, 1e-9);
0953         BOOST_CHECK_EQUAL(axis.getNBins(), 5);
0954         BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
0955         BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Closed);
0956
0957         // Caution: for full-azimuth cases, the ordering is preserved, you get
0958         // in what you get out. -> this can flip
0959         std::vector<double> expected = {-94.2478, -34.0339, 0,
0960                                         31.4159,  62.8319,  94.2478};
0961         CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10);
0962       }
0963     }
0964
0965     BOOST_TEST_CONTEXT("Right variable") {
0966       BOOST_TEST_CONTEXT("Non-closed") {
0967         auto gridLeft = GridPortalLink::make(
0968             cylPhi1, AxisDirection::AxisRPhi,
0969             Axis{AxisBound, -20_degree * 30_mm, 20_degree * 30_mm, 3});
0970         auto gridRight =
0971             GridPortalLink::make(cylPhi2, AxisDirection::AxisRPhi,
0972                                  Axis{AxisBound,
0973                                       {-40_degree * 30_mm, -10_degree * 30_mm,
0974                                        10_degree * 30_mm, 40_degree * 30_mm}});
0975
0976         auto mergedPtr = GridPortalLink::merge(
0977             *gridLeft, *gridRight, AxisDirection::AxisRPhi, *logger);
0978         BOOST_REQUIRE_NE(mergedPtr, nullptr);
0979
0980         const auto* merged =
0981             dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0982         BOOST_REQUIRE_NE(merged, nullptr);
0983         const auto& axis = *merged->grid().axes().front();
0984         BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0985         BOOST_CHECK_CLOSE(axis.getMin(), -60_degree * 30_mm, 1e-9);
0986         BOOST_CHECK_CLOSE(axis.getMax(), 60_degree * 30_mm, 1e-9);
0987         BOOST_CHECK_EQUAL(axis.getNBins(), 6);
0988         BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
0989         BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0990
0991         std::vector<double> expected = {-31.4159, -15.708, -5.23599, 10.472,
0992                                         17.4533,  24.4346, 31.4159};
0993         CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10);
0994       }
0995
0996       BOOST_TEST_CONTEXT("Closed") {
0997         auto gridLeft = GridPortalLink::make(
0998             cylPhi4, AxisDirection::AxisRPhi,
0999             Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 3});
1000
1001         auto gridRight = GridPortalLink::make(
1002             cylPhi3, AxisDirection::AxisRPhi,
1003             Axis{AxisBound,
1004                  {-90_degree * 30_mm, 25_degree * 30_mm, 90_degree * 30_mm}});
1005
1006         auto mergedPtr = GridPortalLink::merge(
1007             *gridLeft, *gridRight, AxisDirection::AxisRPhi, *logger);
1008         BOOST_REQUIRE_NE(mergedPtr, nullptr);
1009
1010         const auto* merged =
1011             dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1012         BOOST_REQUIRE_NE(merged, nullptr);
1013         BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1014         const auto& axis = *merged->grid().axes().front();
1015         BOOST_CHECK_CLOSE(axis.getMin(), -180_degree * 30_mm, 1e-9);
1016         BOOST_CHECK_CLOSE(axis.getMax(), 180_degree * 30_mm, 1e-9);
1017         BOOST_CHECK_EQUAL(axis.getNBins(), 5);
1018         BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
1019         BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Closed);
1020
1021         // Caution: for full-azimuth cases, the ordering is preserved, you get
1022         // in what you get out. -> this can flip
1023         std::vector<double> expected = {-94.2478, -62.8319, -31.4159,
1024                                         0,        60.2139,  94.2478};
1025         CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10);
1026       }
1027     }
1028
1029     BOOST_TEST_CONTEXT("Both variable") {
1030       BOOST_TEST_CONTEXT("Non-closed") {
1031         auto gridLeft =
1032             GridPortalLink::make(cylPhi1, AxisDirection::AxisRPhi,
1033                                  Axis{AxisBound,
1034                                       {-20_degree * 30_mm, -10_degree * 30_mm,
1035                                        10_degree * 30_mm, 20_degree * 30_mm}});
1036         auto gridRight = GridPortalLink::make(
1037             cylPhi2, AxisDirection::AxisRPhi,
1038             Axis{AxisBound,
1039                  {-40_degree * 30_mm, -5_degree * 30_mm, 40_degree * 30_mm}});
1040
1041         auto mergedPtr = GridPortalLink::merge(
1042             *gridLeft, *gridRight, AxisDirection::AxisRPhi, *logger);
1043         BOOST_REQUIRE_NE(mergedPtr, nullptr);
1044
1045         const auto* merged =
1046             dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1047         BOOST_REQUIRE_NE(merged, nullptr);
1048         const auto& axis = *merged->grid().axes().front();
1049         BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1050         BOOST_CHECK_CLOSE(axis.getMin(), -60_degree * 30_mm, 1e-9);
1051         BOOST_CHECK_CLOSE(axis.getMax(), 60_degree * 30_mm, 1e-9);
1052         BOOST_CHECK_EQUAL(axis.getNBins(), 5);
1053         BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
1054         BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
1055
1056         std::vector<double> expected = {-31.4159, -13.09,  10.472,
1057                                         15.708,   26.1799, 31.4159};
1058         CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10);
1059       }
1060
1061       BOOST_TEST_CONTEXT("Closed") {
1062         auto gridLeft = GridPortalLink::make(
1063             cylPhi4, AxisDirection::AxisRPhi,
1064             Axis{AxisBound,
1065                  {-90_degree * 30_mm, 25_degree * 30_mm, 90_degree * 30_mm}});
1066
1067         auto gridRight =
1068             GridPortalLink::make(cylPhi3, AxisDirection::AxisRPhi,
1069                                  Axis{AxisBound,
1070                                       {-90_degree * 30_mm, -10_degree * 30_mm,
1071                                        10_degree * 30_mm, 90_degree * 30_mm}});
1072
1073         auto mergedPtr = GridPortalLink::merge(
1074             *gridLeft, *gridRight, AxisDirection::AxisRPhi, *logger);
1075         BOOST_REQUIRE_NE(mergedPtr, nullptr);
1076
1077         const auto* merged =
1078             dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1079         BOOST_REQUIRE_NE(merged, nullptr);
1080         BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1081         const auto& axis = *merged->grid().axes().front();
1082         BOOST_CHECK_CLOSE(axis.getMin(), -180_degree * 30_mm, 1e-9);
1083         BOOST_CHECK_CLOSE(axis.getMax(), 180_degree * 30_mm, 1e-9);
1084         BOOST_CHECK_EQUAL(axis.getNBins(), 5);
1085         BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
1086         BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Closed);
1087
1088         // Caution: for full-azimuth cases, the ordering is preserved, you get
1089         // in what you get out. -> this can flip
1090         std::vector<double> expected = {-94.2478, -34.0339, 0,
1091                                         41.8879,  52.3599,  94.2478};
1092         CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10);
1093       }
1094     }
1095   }
1096 }
1097
1098 BOOST_AUTO_TEST_CASE(ParallelMerge) {
1099   // Merge in phi direction with z binning
1100   auto cylPhi1 = Surface::makeShared<CylinderSurface>(
1101       Transform3::Identity() * AngleAxis3(45_degree, Vector3::UnitZ()), 30_mm,
1102       100_mm, 20_degree, 0_degree);
1103
1104   auto cylPhi2 = Surface::makeShared<CylinderSurface>(
1105       Transform3::Identity() * AngleAxis3(85_degree, Vector3::UnitZ()), 30_mm,
1106       100_mm, 20_degree, 0_degree);
1107
1108   auto portalPhi1 = GridPortalLink::make(cylPhi1, AxisDirection::AxisZ,
1109                                          Axis{AxisBound, -100_mm, 100_mm, 5});
1110
1111   auto portalPhi2 = GridPortalLink::make(cylPhi2, AxisDirection::AxisZ,
1112                                          Axis{AxisBound, -100_mm, 100_mm, 5});
1113
1114   auto merged12Ptr = GridPortalLink::merge(*portalPhi1, *portalPhi2,
1115                                            AxisDirection::AxisRPhi, *logger);
1116   BOOST_REQUIRE_NE(merged12Ptr, nullptr);
1117   auto merged12 = dynamic_cast<const GridPortalLink*>(merged12Ptr.get());
1118   BOOST_REQUIRE_NE(merged12, nullptr);
1119
1120   const auto& axis1 = *merged12->grid().axes().front();
1121   const auto& axis2 = *merged12->grid().axes().back();
1122   // Phi sectors were same size, should give equidistant binning
1123   Axis axis1Expected{AxisBound, -40_degree * 30_mm, 40_degree * 30_mm, 2};
1124   BOOST_CHECK_EQUAL(axis1, axis1Expected);
1125   Axis axis2Expected{AxisBound, -100_mm, 100_mm, 5};
1126   BOOST_CHECK_EQUAL(axis2, axis2Expected);
1127 }
1128
1129 BOOST_AUTO_TEST_SUITE_END()  // RPhiDirection
1130
1131 BOOST_AUTO_TEST_SUITE_END()  // Merging1dCylinder
1132
1133 BOOST_AUTO_TEST_SUITE(Merging2dCylinder)
1134
1135 BOOST_AUTO_TEST_CASE(ZDirection) {
1136   BOOST_TEST_CONTEXT("Phi sector") {
1137     auto cyl1 = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
1138                                                      30_mm, 100_mm, 45_degree);
1139
1140     // z good, rphi good
1141     auto grid1 = GridPortalLink::make(
1142         cyl1, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5},
1143         Axis{AxisBound, -100_mm, 100_mm, 5});
1144
1145     auto trf2 = Transform3{Translation3{Vector3::UnitZ() * 150_mm}};
1146     auto cyl2 =
1147         Surface::makeShared<CylinderSurface>(trf2, 30_mm, 50_mm, 45_degree);
1148
1149     // Second grid portal with compatible phi binning
1150     auto grid2 = GridPortalLink::make(
1151         cyl2, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5},
1152         Axis{AxisBound, -50_mm, 50_mm, 5});
1153
1154     // We're merging in z direction, so the phi binnings need to be the same
1155
1156     auto mergedPtr =
1157         GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisZ, *logger);
1158     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1159     BOOST_REQUIRE_NE(mergedPtr, nullptr);
1160
1161     const auto& axis1 = *merged->grid().axes().front();
1162     const auto& axis2 = *merged->grid().axes().back();
1163
1164     BOOST_CHECK_EQUAL(axis1.getMin(), -45_degree * 30_mm);
1165     BOOST_CHECK_EQUAL(axis1.getMax(), 45_degree * 30_mm);
1166     BOOST_CHECK_EQUAL(axis1.getNBins(), 5);
1167     BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant);
1168     BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound);
1169
1170     BOOST_CHECK_EQUAL(axis2.getMin(), -150_mm);
1171     BOOST_CHECK_EQUAL(axis2.getMax(), 150_mm);
1172     BOOST_CHECK_EQUAL(axis2.getNBins(), 10);
1173     BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Variable);
1174     BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound);
1175
1176     auto grid3 = GridPortalLink::make(
1177         cyl2, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 3},
1178         Axis{AxisBound, -50_mm, 50_mm, 5});
1179
1180     auto composite = PortalLinkBase::merge(copy(grid1), copy(grid3),
1181                                            AxisDirection::AxisZ, *logger);
1182     BOOST_CHECK_NE(dynamic_cast<const CompositePortalLink*>(composite.get()),
1183                    nullptr);
1184   }
1185
1186   BOOST_TEST_CONTEXT("Check wraparound for full circle in phi") {
1187     auto cyl1 = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
1188                                                      30_mm, 100_mm, 180_degree);
1189
1190     // z good, rphi good
1191     auto grid1 = GridPortalLink::make(
1192         cyl1, Axis{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 5},
1193         Axis{AxisBound, -100_mm, 100_mm, 5});
1194
1195     auto trf2 = Transform3{Translation3{Vector3::UnitZ() * 150_mm}};
1196     auto cyl2 =
1197         Surface::makeShared<CylinderSurface>(trf2, 30_mm, 50_mm, 180_degree);
1198
1199     // Second grid portal with compatible phi binning
1200     auto grid2 = GridPortalLink::make(
1201         cyl2, Axis{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 5},
1202         Axis{AxisBound, -50_mm, 50_mm, 5});
1203
1204     auto mergedPtr =
1205         GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisZ, *logger);
1206     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1207     BOOST_REQUIRE_NE(mergedPtr, nullptr);
1208
1209     const auto& axis1 = *merged->grid().axes().front();
1210     const auto& axis2 = *merged->grid().axes().back();
1211
1212     Axis axis1Expected{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 5};
1213     BOOST_CHECK_EQUAL(axis1, axis1Expected);
1214     Axis axis2Expected{AxisBound,
1215                        {-150, -110, -70, -30, 10, 50, 70, 90, 110, 130, 150}};
1216     BOOST_CHECK_EQUAL(axis2, axis2Expected);
1217   }
1218 }
1219
1220 BOOST_AUTO_TEST_CASE(RPhiDirection) {
1221   auto cyl1 = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
1222                                                    30_mm, 100_mm, 45_degree);
1223
1224   // z good, rphi good
1225   auto grid1 = GridPortalLink::make(
1226       cyl1, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5},
1227       Axis{AxisBound, -100_mm, 100_mm, 5});
1228   BOOST_REQUIRE_NE(grid1, nullptr);
1229
1230   auto trf2 = Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}};
1231   auto cyl2 =
1232       Surface::makeShared<CylinderSurface>(trf2, 30_mm, 100_mm, 45_degree);
1233
1234   // Second grid portal with compatible phi binning
1235   auto grid2 = GridPortalLink::make(
1236       cyl2, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5},
1237       Axis{AxisBound, -100_mm, 100_mm, 5});
1238   BOOST_REQUIRE_NE(grid2, nullptr);
1239
1240   // We're merging in z direction, so the phi binnings need to be the same
1241
1242   auto mergedPtr =
1243       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisRPhi, *logger);
1244   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1245   BOOST_REQUIRE_NE(mergedPtr, nullptr);
1246
1247   const auto& axis1 = *merged->grid().axes().front();
1248   const auto& axis2 = *merged->grid().axes().back();
1249   BOOST_CHECK_CLOSE(axis1.getMin(), -90_degree * 30_mm, 1e-8);
1250   BOOST_CHECK_CLOSE(axis1.getMax(), 90_degree * 30_mm, 1e-8);
1251   BOOST_CHECK_EQUAL(axis1.getNBins(), 10);
1252   BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant);
1253   BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound);
1254   Axis axis2Expected{AxisBound, -100_mm, 100_mm, 5};
1255   BOOST_CHECK_EQUAL(axis2, axis2Expected);
1256 }
1257
1258 BOOST_AUTO_TEST_SUITE_END()  // Merging2dCylinder
1259
1260 BOOST_AUTO_TEST_SUITE(Merging1dDisc)
1261
1262 BOOST_AUTO_TEST_SUITE(RDirection)
1263
1264 BOOST_AUTO_TEST_CASE(ColinearMerge) {
1265   // Without phi sector
1266   auto disc1 =
1267       Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 100_mm);
1268
1269   auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisR,
1270                                     Axis{AxisBound, 30_mm, 100_mm, 7});
1271
1272   auto disc2 =
1273       Surface::makeShared<DiscSurface>(Transform3::Identity(), 100_mm, 150_mm);
1274
1275   auto grid2 = GridPortalLink::make(disc2, AxisDirection::AxisR,
1276                                     Axis{AxisBound, 100_mm, 150_mm, 5});
1277
1278   auto mergedPtr =
1279       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisR, *logger);
1280   BOOST_REQUIRE(mergedPtr);
1281   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1282   BOOST_REQUIRE_NE(merged, nullptr);
1283
1284   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1285   Axis axisExpected{AxisBound, 30_mm, 150_mm, 12};
1286   BOOST_CHECK_EQUAL(*merged->grid().axes().front(), axisExpected);
1287
1288   // With phi sector
1289   auto discPhi1 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1290                                                    30_mm, 100_mm, 30_degree);
1291
1292   auto discPhiGrid1 = GridPortalLink::make(discPhi1, AxisDirection::AxisR,
1293                                            Axis{AxisBound, 30_mm, 100_mm, 7});
1294
1295   auto discPhi2 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1296                                                    100_mm, 150_mm, 30_degree);
1297
1298   auto discPhiGrid2 = GridPortalLink::make(discPhi2, AxisDirection::AxisR,
1299                                            Axis{AxisBound, 100_mm, 150_mm, 5});
1300
1301   auto mergedPhiPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid2,
1302                                             AxisDirection::AxisR, *logger);
1303   BOOST_REQUIRE(mergedPhiPtr);
1304   const auto* mergedPhi =
1305       dynamic_cast<const GridPortalLink*>(mergedPhiPtr.get());
1306   BOOST_REQUIRE_NE(mergedPhi, nullptr);
1307
1308   BOOST_CHECK_EQUAL(mergedPhi->grid().axes().size(), 1);
1309   BOOST_CHECK_EQUAL(*mergedPhi->grid().axes().front(), axisExpected);
1310 }
1311
1312 BOOST_AUTO_TEST_CASE(ParallelMerge) {
1313   auto disc1 =
1314       Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 100_mm);
1315
1316   auto grid1 =
1317       GridPortalLink::make(disc1, AxisDirection::AxisPhi,
1318                            Axis{AxisClosed, -180_degree, 180_degree, 5});
1319
1320   auto disc2 =
1321       Surface::makeShared<DiscSurface>(Transform3::Identity(), 100_mm, 150_mm);
1322
1323   auto grid2 =
1324       GridPortalLink::make(disc2, AxisDirection::AxisPhi,
1325                            Axis{AxisClosed, -180_degree, 180_degree, 5});
1326
1327   auto mergedPtr =
1328       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisR, *logger);
1329   BOOST_REQUIRE(mergedPtr);
1330   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1331   BOOST_REQUIRE_NE(merged, nullptr);
1332
1333   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2);
1334   const auto& axis1 = *merged->grid().axes().front();
1335   const auto& axis2 = *merged->grid().axes().back();
1336   Axis axis1Expected{AxisBound, {30_mm, 100_mm, 150_mm}};
1337   BOOST_CHECK_EQUAL(axis1, axis1Expected);
1338   Axis axis2Expected{AxisClosed, -180_degree, 180_degree, 5};
1339   BOOST_CHECK_EQUAL(axis2, axis2Expected);
1340 }
1341
1342 BOOST_AUTO_TEST_SUITE_END()  // RDirection
1343
1344 BOOST_AUTO_TEST_SUITE(PhiDirection)
1345
1346 BOOST_AUTO_TEST_CASE(ColinearMerge) {
1347   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1348                                                 100_mm, 30_degree);
1349
1350   auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisPhi,
1351                                     Axis{AxisBound, -30_degree, 30_degree, 3});
1352
1353   auto disc2 = Surface::makeShared<DiscSurface>(
1354       Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1355       60_degree);
1356
1357   auto grid2 = GridPortalLink::make(disc2, AxisDirection::AxisPhi,
1358                                     Axis{AxisBound, -60_degree, 60_degree, 6});
1359
1360   auto mergedPtr =
1361       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
1362   BOOST_REQUIRE(mergedPtr);
1363   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1364   BOOST_REQUIRE_NE(merged, nullptr);
1365
1366   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1367   const auto& axis = *merged->grid().axes().front();
1368   BOOST_CHECK_CLOSE(axis.getMin(), -90_degree, 1e-6);
1369   BOOST_CHECK_CLOSE(axis.getMax(), 90_degree, 1e-6);
1370   BOOST_CHECK_EQUAL(axis.getNBins(), 9);
1371
1372   // Check wrapping
1373
1374   auto disc1Half = Surface::makeShared<DiscSurface>(
1375       Transform3{AngleAxis3{15_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1376       90_degree);
1377
1378   auto grid1Half =
1379       GridPortalLink::make(disc1Half, AxisDirection::AxisPhi,
1380                            Axis{AxisBound, -90_degree, 90_degree, 3});
1381
1382   auto disc2Half = Surface::makeShared<DiscSurface>(
1383       Transform3{AngleAxis3{-165_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1384       90_degree);
1385
1386   auto grid2Half =
1387       GridPortalLink::make(disc2Half, AxisDirection::AxisPhi,
1388                            Axis{AxisBound, -90_degree, 90_degree, 3});
1389
1390   auto mergedHalfPtr = GridPortalLink::merge(*grid1Half, *grid2Half,
1391                                              AxisDirection::AxisPhi, *logger);
1392   BOOST_REQUIRE(mergedHalfPtr);
1393   const auto* mergedHalf =
1394       dynamic_cast<const GridPortalLink*>(mergedHalfPtr.get());
1395   BOOST_REQUIRE_NE(mergedHalf, nullptr);
1396
1397   BOOST_CHECK_EQUAL(mergedHalf->grid().axes().size(), 1);
1398   Axis axisHalfExpected{AxisClosed, -180_degree, 180_degree, 6};
1399   BOOST_CHECK_EQUAL(axisHalfExpected, *mergedHalf->grid().axes().front());
1400 }
1401
1402 BOOST_AUTO_TEST_CASE(ParallelMerge) {
1403   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1404                                                 100_mm, 30_degree);
1405
1406   auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisR,
1407                                     Axis{AxisBound, 30_mm, 100_mm, 3});
1408
1409   auto disc2 = Surface::makeShared<DiscSurface>(
1410       Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1411       60_degree);
1412
1413   auto grid2 = GridPortalLink::make(disc2, AxisDirection::AxisR,
1414                                     Axis{AxisBound, 30_mm, 100_mm, 3});
1415
1416   auto mergedPtr =
1417       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
1418   BOOST_REQUIRE(mergedPtr);
1419   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1420   BOOST_REQUIRE_NE(merged, nullptr);
1421
1422   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2);
1423   const auto& axis1 = *merged->grid().axes().front();
1424   const auto& axis2 = *merged->grid().axes().back();
1425
1426   Axis axis1Expected{AxisBound, 30_mm, 100_mm, 3};
1427   BOOST_CHECK_EQUAL(axis1, axis1Expected);
1428
1429   BOOST_CHECK_CLOSE(axis2.getMin(), -90_degree, 1e-6);
1430   BOOST_CHECK_CLOSE(axis2.getMax(), 90_degree, 1e-6);
1431   BOOST_CHECK_EQUAL(axis2.getNBins(), 2);
1432   BOOST_CHECK_CLOSE(axis2.getBinEdges().at(1), 30_degree, 1e-6);
1433   BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Variable);
1434   BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound);
1435 }
1436
1437 BOOST_AUTO_TEST_SUITE_END()  // PhiDirection
1438
1439 BOOST_AUTO_TEST_CASE(BinFilling) {
1440   // Volumes for bin content checking
1441
1442   auto vol1 = makeDummyVolume();
1443   auto vol2 = makeDummyVolume();
1444
1445   BOOST_TEST_CONTEXT("RDirection") {
1446     auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1447                                                   60_mm, 30_degree);
1448
1449     auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisR,
1450                                       Axis{AxisBound, 30_mm, 60_mm, 2});
1451
1452     grid1->setVolume(vol1.get());
1453
1454     auto disc2 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 60_mm,
1455                                                   90_mm, 30_degree);
1456
1457     auto grid2 = GridPortalLink::make(disc2, AxisDirection::AxisR,
1458                                       Axis{AxisBound, 60_mm, 90_mm, 2});
1459
1460     grid2->setVolume(vol2.get());
1461
1462     auto mergedPtr =
1463         GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisR, *logger);
1464
1465     using merged_type =
1466         GridPortalLinkT<Axis<AxisType::Equidistant, AxisBoundaryType::Bound>>;
1467
1468     const auto* merged = dynamic_cast<const merged_type*>(mergedPtr.get());
1469     BOOST_REQUIRE(merged);
1470
1471     grid1->printContents(std::cout);
1472     grid2->printContents(std::cout);
1473     merged->printContents(std::cout);
1474
1475     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({1}), vol1.get());
1476     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({2}), vol1.get());
1477     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({3}), vol2.get());
1478     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({4}), vol2.get());
1479   }
1480
1481   BOOST_TEST_CONTEXT("PhiDirection") {
1482     auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1483                                                   100_mm, 30_degree);
1484
1485     auto grid1 =
1486         GridPortalLink::make(disc1, AxisDirection::AxisPhi,
1487                              Axis{AxisBound, -30_degree, 30_degree, 2});
1488
1489     grid1->setVolume(vol1.get());
1490
1491     auto disc2 = Surface::makeShared<DiscSurface>(
1492         Transform3{AngleAxis3{60_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1493         30_degree);
1494
1495     auto grid2 =
1496         GridPortalLink::make(disc2, AxisDirection::AxisPhi,
1497                              Axis{AxisBound, -30_degree, 30_degree, 2});
1498
1499     grid2->setVolume(vol2.get());
1500
1501     auto mergedPtr =
1502         GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
1503     BOOST_REQUIRE(mergedPtr);
1504
1505     using merged_type =
1506         GridPortalLinkT<Axis<AxisType::Equidistant, AxisBoundaryType::Bound>>;
1507
1508     const auto* merged = dynamic_cast<const merged_type*>(mergedPtr.get());
1509     BOOST_REQUIRE(merged);
1510
1511     grid1->printContents(std::cout);
1512     grid2->printContents(std::cout);
1513     merged->printContents(std::cout);
1514
1515     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({1}), vol2.get());
1516     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({2}), vol2.get());
1517     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({3}), vol1.get());
1518     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({4}), vol1.get());
1519   }
1520 }
1521
1522 BOOST_AUTO_TEST_SUITE_END()  // Merging1dDisc
1523
1524 BOOST_AUTO_TEST_SUITE(Merging2dDisc)
1525
1526 BOOST_AUTO_TEST_CASE(RDirection) {
1527   // Basic, because the parallel 1D case already tests this to some degree
1528   auto discPhi1 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1529                                                    30_mm, 100_mm, 30_degree);
1530
1531   auto discPhiGrid1 =
1532       GridPortalLink::make(discPhi1, Axis{AxisBound, 30_mm, 100_mm, 7},
1533                            Axis{AxisBound, -30_degree, 30_degree, 3});
1534
1535   auto discPhi2 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1536                                                    100_mm, 150_mm, 30_degree);
1537
1538   auto discPhiGrid2 =
1539       GridPortalLink::make(discPhi2, Axis{AxisBound, 100_mm, 150_mm, 5},
1540                            Axis{AxisBound, -30_degree, 30_degree, 3});
1541
1542   auto mergedPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid2,
1543                                          AxisDirection::AxisR, *logger);
1544   BOOST_REQUIRE(mergedPtr);
1545   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1546   BOOST_REQUIRE_NE(merged, nullptr);
1547   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2);
1548   const auto& axis1 = *merged->grid().axes().front();
1549   const auto& axis2 = *merged->grid().axes().back();
1550
1551   BOOST_CHECK_EQUAL(axis1.getMin(), 30_mm);
1552   BOOST_CHECK_EQUAL(axis1.getMax(), 150_mm);
1553   BOOST_CHECK_EQUAL(axis1.getNBins(), 12);
1554   BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant);
1555   BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound);
1556   BOOST_CHECK_EQUAL(axis2.getMin(), -30_degree);
1557   BOOST_CHECK_EQUAL(axis2.getMax(), 30_degree);
1558   BOOST_CHECK_EQUAL(axis2.getNBins(), 3);
1559   BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant);
1560 }
1561
1562 BOOST_AUTO_TEST_CASE(PhiDirection) {
1563   // Basic, because the parallel 1D case already tests this to some degree
1564   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1565                                                 100_mm, 30_degree);
1566
1567   auto grid1 = GridPortalLink::make(disc1, Axis{AxisBound, 30_mm, 100_mm, 3},
1568                                     Axis{AxisBound, -30_degree, 30_degree, 3});
1569
1570   auto disc2 = Surface::makeShared<DiscSurface>(
1571       Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1572       60_degree);
1573
1574   auto grid2 = GridPortalLink::make(disc2, Axis{AxisBound, 30_mm, 100_mm, 3},
1575                                     Axis{AxisBound, -60_degree, 60_degree, 6});
1576
1577   auto mergedPtr =
1578       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
1579   BOOST_REQUIRE(mergedPtr);
1580   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1581   BOOST_REQUIRE_NE(merged, nullptr);
1582
1583   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2);
1584   const auto& axis1 = *merged->grid().axes().front();
1585   const auto& axis2 = *merged->grid().axes().back();
1586
1587   BOOST_CHECK_EQUAL(axis1.getMin(), 30_mm);
1588   BOOST_CHECK_EQUAL(axis1.getMax(), 100_mm);
1589   BOOST_CHECK_EQUAL(axis1.getNBins(), 3);
1590   BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant);
1591   BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound);
1592   BOOST_CHECK_CLOSE(axis2.getMin(), -90_degree, 1e-6);
1593   BOOST_CHECK_CLOSE(axis2.getMax(), 90_degree, 1e-6);
1594   BOOST_CHECK_EQUAL(axis2.getNBins(), 9);
1595   BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant);
1596   BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound);
1597 }
1598
1599 BOOST_AUTO_TEST_CASE(BinFilling) {
1600   // Volumes for bin content checking
1601   // Volume shape/transform is irrelevant, only used for pointer identity
1602   auto vol1 = std::make_shared<TrackingVolume>(
1603       Transform3::Identity(),
1604       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
1605
1606   auto vol2 = std::make_shared<TrackingVolume>(
1607       Transform3::Identity(),
1608       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
1609
1610   auto fillCheckerBoard = [&](auto& grid) {
1611     auto loc = grid.numLocalBins();
1612     for (std::size_t i = 1; i <= loc[0]; ++i) {
1613       for (std::size_t j = 1; j <= loc[1]; ++j) {
1614         grid.atLocalBins({i, j}) = (i + j) % 2 == 0 ? vol1.get() : vol2.get();
1615       }
1616     }
1617   };
1618
1619   auto checkCheckerBoard = [&](const auto& grid) {
1620     auto loc = grid.numLocalBins();
1621     for (std::size_t i = 1; i <= loc[0]; ++i) {
1622       for (std::size_t j = 1; j <= loc[1]; ++j) {
1623         const auto* vol = grid.atLocalBins({i, j});
1624         if (vol != ((i + j) % 2 == 0 ? vol1.get() : vol2.get())) {
1625           BOOST_ERROR("Is not a checkerboard pattern");
1626           return;
1627         }
1628       }
1629     }
1630   };
1631
1632   BOOST_TEST_CONTEXT("RDirection") {
1633     auto discPhi1 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1634                                                      30_mm, 60_mm, 30_degree);
1635
1636     auto discPhiGrid1 =
1637         GridPortalLink::make(discPhi1, Axis{AxisBound, 30_mm, 60_mm, 2},
1638                              Axis{AxisBound, -30_degree, 30_degree, 2});
1639
1640     fillCheckerBoard(discPhiGrid1->grid());
1641     checkCheckerBoard(discPhiGrid1->grid());
1642
1643     auto discPhi2 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1644                                                      60_mm, 90_mm, 30_degree);
1645
1646     auto discPhiGrid2 =
1647         GridPortalLink::make(discPhi2, Axis{AxisBound, 60_mm, 90_mm, 2},
1648                              Axis{AxisBound, -30_degree, 30_degree, 2});
1649
1650     fillCheckerBoard(discPhiGrid2->grid());
1651     checkCheckerBoard(discPhiGrid2->grid());
1652
1653     auto mergedPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid2,
1654                                            AxisDirection::AxisR, *logger);
1655
1656     using merged_type =
1657         GridPortalLinkT<Axis<AxisType::Equidistant, AxisBoundaryType::Bound>,
1658                         Axis<AxisType::Equidistant, AxisBoundaryType::Bound>>;
1659
1660     const auto* merged = dynamic_cast<const merged_type*>(mergedPtr.get());
1661     BOOST_REQUIRE(merged);
1662     checkCheckerBoard(merged->grid());
1663
1664     // Fill a / b
1665     discPhiGrid1->setVolume(vol1.get());
1666     discPhiGrid2->setVolume(vol2.get());
1667
1668     mergedPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid2,
1669                                       AxisDirection::AxisR, *logger);
1670
1671     merged = dynamic_cast<const merged_type*>(mergedPtr.get());
1672     BOOST_REQUIRE(merged);
1673
1674     const auto* v1 = vol1.get();
1675     const auto* v2 = vol2.get();
1676
1677     std::vector<std::pair<Vector2, const TrackingVolume*>> locations = {
1678         {{40_mm, -20_degree}, v1}, {{40_mm, 20_degree}, v1},
1679         {{50_mm, -20_degree}, v1}, {{50_mm, 20_degree}, v1},
1680         {{70_mm, -20_degree}, v2}, {{70_mm, 20_degree}, v2},
1681         {{80_mm, -20_degree}, v2}, {{80_mm, 20_degree}, v2},
1682     };
1683
1684     for (const auto& [loc, vol] : locations) {
1685       BOOST_TEST_CONTEXT(loc.transpose())
1686       BOOST_CHECK_EQUAL(merged->resolveVolume(gctx, loc).value(), vol);
1687     }
1688
1689     std::vector<std::vector<const TrackingVolume*>> contents = {
1690         {v1, v1},
1691         {v1, v1},
1692         {v2, v2},
1693         {v2, v2},
1694     };
1695
1696     for (std::size_t i = 0; i < 4; ++i) {
1697       for (std::size_t j = 0; j < 2; ++j) {
1698         BOOST_CHECK_EQUAL(merged->grid().atLocalBins({i + 1, j + 1}),
1699                           contents.at(i).at(j));
1700       }
1701     }
1702   }
1703
1704   BOOST_TEST_CONTEXT("PhiDirection") {
1705     auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1706                                                   100_mm, 30_degree);
1707
1708     auto grid1 =
1709         GridPortalLink::make(disc1, Axis{AxisBound, 30_mm, 100_mm, 2},
1710                              Axis{AxisBound, -30_degree, 30_degree, 2});
1711     fillCheckerBoard(grid1->grid());
1712     checkCheckerBoard(grid1->grid());
1713
1714     auto disc2 = Surface::makeShared<DiscSurface>(
1715         Transform3{AngleAxis3{60_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1716         30_degree);
1717
1718     auto grid2 =
1719         GridPortalLink::make(disc2, Axis{AxisBound, 30_mm, 100_mm, 2},
1720                              Axis{AxisBound, -30_degree, 30_degree, 2});
1721     fillCheckerBoard(grid2->grid());
1722     checkCheckerBoard(grid2->grid());
1723
1724     auto mergedPtr =
1725         GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
1726
1727     BOOST_REQUIRE(mergedPtr);
1728
1729     using merged_type =
1730         GridPortalLinkT<Axis<AxisType::Equidistant, AxisBoundaryType::Bound>,
1731                         Axis<AxisType::Equidistant, AxisBoundaryType::Bound>>;
1732
1733     const auto* merged = dynamic_cast<const merged_type*>(mergedPtr.get());
1734     BOOST_REQUIRE(merged);
1735
1736     checkCheckerBoard(merged->grid());
1737
1738     // Fill a / b
1739     grid1->setVolume(vol1.get());
1740     grid2->setVolume(vol2.get());
1741
1742     mergedPtr =
1743         GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
1744     merged = dynamic_cast<const merged_type*>(mergedPtr.get());
1745     BOOST_REQUIRE(merged);
1746
1747     const auto* v1 = vol1.get();
1748     const auto* v2 = vol2.get();
1749
1750     std::vector<std::pair<Vector2, const TrackingVolume*>> locations = {
1751         {{40_mm, -50_degree}, v2}, {{40_mm, -10_degree}, v2},
1752         {{50_mm, -50_degree}, v2}, {{50_mm, -10_degree}, v2},
1753         {{40_mm, 10_degree}, v1},  {{50_mm, 50_degree}, v1},
1754         {{50_mm, 10_degree}, v1},  {{50_mm, 50_degree}, v1},
1755     };
1756
1757     for (const auto& [loc, vol] : locations) {
1758       BOOST_TEST_CONTEXT(loc.transpose())
1759       BOOST_CHECK_EQUAL(merged->resolveVolume(gctx, loc).value(), vol);
1760     }
1761
1762     std::vector<std::vector<const TrackingVolume*>> contents = {
1763         {v2, v2, v1, v1},
1764         {v2, v2, v1, v1},
1765     };
1766
1767     for (std::size_t i = 0; i < 2; ++i) {
1768       for (std::size_t j = 0; j < 4; ++j) {
1769         BOOST_CHECK_EQUAL(merged->grid().atLocalBins({i + 1, j + 1}),
1770                           contents.at(i).at(j));
1771       }
1772     }
1773   }
1774 }
1775
1776 BOOST_AUTO_TEST_SUITE_END()  // Merging2dDisc
1777
1778 BOOST_AUTO_TEST_SUITE(MergingMixedDisc)
1779
1780 BOOST_AUTO_TEST_CASE(RDirection) {
1781   auto discPhi1 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1782                                                    30_mm, 100_mm, 30_degree);
1783
1784   auto discPhiGrid1 =
1785       GridPortalLink::make(discPhi1, Axis{AxisBound, 30_mm, 100_mm, 7},
1786                            Axis{AxisBound, -30_degree, 30_degree, 3});
1787
1788   auto discPhi2 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1789                                                    100_mm, 150_mm, 30_degree);
1790
1791   auto discPhiGrid21dPhi =
1792       GridPortalLink::make(discPhi2, AxisDirection::AxisPhi,
1793                            Axis{AxisBound, -30_degree, 30_degree, 3});
1794
1795   auto merged12PhiPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid21dPhi,
1796                                               AxisDirection::AxisR, *logger);
1797   BOOST_REQUIRE(merged12PhiPtr);
1798   const auto* merged12Phi =
1799       dynamic_cast<const GridPortalLink*>(merged12PhiPtr.get());
1800   BOOST_REQUIRE_NE(merged12Phi, nullptr);
1801
1802   auto merged21PhiPtr = GridPortalLink::merge(*discPhiGrid21dPhi, *discPhiGrid1,
1803                                               AxisDirection::AxisR, *logger);
1804   BOOST_REQUIRE(merged21PhiPtr);
1805   const auto* merged21Phi =
1806       dynamic_cast<const GridPortalLink*>(merged21PhiPtr.get());
1807   BOOST_REQUIRE_NE(merged21Phi, nullptr);
1808
1809   BOOST_CHECK_EQUAL(merged12Phi->grid(), merged21Phi->grid());
1810
1811   BOOST_CHECK_EQUAL(merged12Phi->grid().axes().size(), 2);
1812   const auto& axis1 = *merged12Phi->grid().axes().front();
1813   const auto& axis2 = *merged12Phi->grid().axes().back();
1814
1815   BOOST_CHECK_EQUAL(axis1.getMin(), 30_mm);
1816   BOOST_CHECK_EQUAL(axis1.getMax(), 150_mm);
1817   BOOST_CHECK_EQUAL(axis1.getNBins(), 8);
1818   BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Variable);
1819   BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound);
1820   BOOST_CHECK_EQUAL(axis2.getMin(), -30_degree);
1821   BOOST_CHECK_EQUAL(axis2.getMax(), 30_degree);
1822   BOOST_CHECK_EQUAL(axis2.getNBins(), 3);
1823   BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant);
1824   BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound);
1825
1826   auto discPhiGrid21dR = GridPortalLink::make(
1827       discPhi2, AxisDirection::AxisR, Axis{AxisBound, 100_mm, 150_mm, 5});
1828
1829   auto merged12RPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid21dR,
1830                                             AxisDirection::AxisR, *logger);
1831   BOOST_REQUIRE(merged12RPtr);
1832   const auto* merged12R =
1833       dynamic_cast<const GridPortalLink*>(merged12RPtr.get());
1834   BOOST_REQUIRE_NE(merged12R, nullptr);
1835
1836   auto merged21RPtr = GridPortalLink::merge(*discPhiGrid21dR, *discPhiGrid1,
1837                                             AxisDirection::AxisR, *logger);
1838   BOOST_REQUIRE(merged21RPtr);
1839   const auto* merged21R =
1840       dynamic_cast<const GridPortalLink*>(merged21RPtr.get());
1841   BOOST_REQUIRE_NE(merged21R, nullptr);
1842   BOOST_CHECK_EQUAL(merged12R->grid(), merged21R->grid());
1843
1844   BOOST_CHECK_EQUAL(merged12R->grid().axes().size(), 2);
1845   const auto& axis1R = *merged12R->grid().axes().front();
1846   const auto& axis2R = *merged12R->grid().axes().back();
1847
1848   BOOST_CHECK_EQUAL(axis1R.getMin(), 30_mm);
1849   BOOST_CHECK_EQUAL(axis1R.getMax(), 150_mm);
1850   BOOST_CHECK_EQUAL(axis1R.getNBins(), 12);
1851   BOOST_CHECK_EQUAL(axis1R.getType(), AxisType::Equidistant);
1852   BOOST_CHECK_EQUAL(axis1R.getBoundaryType(), AxisBoundaryType::Bound);
1853   BOOST_CHECK_EQUAL(axis2R.getMin(), -30_degree);
1854   BOOST_CHECK_EQUAL(axis2R.getMax(), 30_degree);
1855   BOOST_CHECK_EQUAL(axis2R.getNBins(), 3);
1856   BOOST_CHECK_EQUAL(axis2R.getType(), AxisType::Equidistant);
1857   BOOST_CHECK_EQUAL(axis2R.getBoundaryType(), AxisBoundaryType::Bound);
1858 }
1859
1860 BOOST_AUTO_TEST_CASE(PhiDirection) {
1861   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1862                                                 100_mm, 30_degree);
1863
1864   auto grid1 = GridPortalLink::make(disc1, Axis{AxisBound, 30_mm, 100_mm, 3},
1865                                     Axis{AxisBound, -30_degree, 30_degree, 3});
1866
1867   auto disc2 = Surface::makeShared<DiscSurface>(
1868       Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1869       60_degree);
1870
1871   auto grid21dPhi =
1872       GridPortalLink::make(disc2, AxisDirection::AxisPhi,
1873
1874                            Axis{AxisBound, -60_degree, 60_degree, 6});
1875
1876   auto merged12PhiPtr = GridPortalLink::merge(*grid1, *grid21dPhi,
1877                                               AxisDirection::AxisPhi, *logger);
1878   BOOST_REQUIRE(merged12PhiPtr);
1879   const auto* merged12Phi =
1880       dynamic_cast<const GridPortalLink*>(merged12PhiPtr.get());
1881   BOOST_REQUIRE_NE(merged12Phi, nullptr);
1882
1883   BOOST_CHECK_EQUAL(merged12Phi->grid().axes().size(), 2);
1884   const auto& axis1 = *merged12Phi->grid().axes().front();
1885   const auto& axis2 = *merged12Phi->grid().axes().back();
1886
1887   BOOST_CHECK_EQUAL(axis1.getMin(), 30_mm);
1888   BOOST_CHECK_EQUAL(axis1.getMax(), 100_mm);
1889   BOOST_CHECK_EQUAL(axis1.getNBins(), 3);
1890   BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant);
1891   BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound);
1892   BOOST_CHECK_CLOSE(axis2.getMin(), -90_degree, 1e-6);
1893   BOOST_CHECK_CLOSE(axis2.getMax(), 90_degree, 1e-6);
1894   BOOST_CHECK_EQUAL(axis2.getNBins(), 9);
1895   BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant);
1896   BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound);
1897
1898   auto grid21dR = GridPortalLink::make(disc2, AxisDirection::AxisR,
1899                                        Axis{AxisBound, 30_mm, 100_mm, 3});
1900
1901   auto merged12RPtr =
1902       GridPortalLink::merge(*grid1, *grid21dR, AxisDirection::AxisPhi, *logger);
1903   BOOST_REQUIRE(merged12RPtr);
1904   const auto* merged12R =
1905       dynamic_cast<const GridPortalLink*>(merged12RPtr.get());
1906   BOOST_REQUIRE_NE(merged12R, nullptr);
1907
1908   BOOST_CHECK_EQUAL(merged12R->grid().axes().size(), 2);
1909   const auto& axis1R = *merged12R->grid().axes().front();
1910   const auto& axis2R = *merged12R->grid().axes().back();
1911
1912   BOOST_CHECK_EQUAL(axis1R.getMin(), 30_mm);
1913   BOOST_CHECK_EQUAL(axis1R.getMax(), 100_mm);
1914   BOOST_CHECK_EQUAL(axis1R.getNBins(), 3);
1915   BOOST_CHECK_EQUAL(axis1R.getType(), AxisType::Equidistant);
1916   BOOST_CHECK_EQUAL(axis1R.getBoundaryType(), AxisBoundaryType::Bound);
1917   BOOST_CHECK_CLOSE(axis2R.getMin(), -90_degree, 1e-6);
1918   BOOST_CHECK_CLOSE(axis2R.getMax(), 90_degree, 1e-6);
1919   BOOST_CHECK_EQUAL(axis2R.getNBins(), 4);
1920   BOOST_CHECK_EQUAL(axis2R.getType(), AxisType::Variable);
1921   BOOST_CHECK_EQUAL(axis2R.getBoundaryType(), AxisBoundaryType::Bound);
1922 }
1923
1924 BOOST_AUTO_TEST_SUITE_END()  // MergingMixedDisc
1925
1926 BOOST_AUTO_TEST_SUITE(MergingCrossDisc)
1927
1928 BOOST_AUTO_TEST_CASE(RDirection) {
1929   // Volumes for bin content checking
1930   auto vol1 = makeDummyVolume();
1931   auto vol2 = makeDummyVolume();
1932   auto vol3 = makeDummyVolume();
1933   auto vol4 = makeDummyVolume();
1934
1935   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1936                                                 100_mm, 30_degree);
1937
1938   auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisR,
1939                                     Axis{AxisBound, 30_mm, 100_mm, 2});
1940   grid1->grid().atLocalBins({1}) = vol1.get();
1941   grid1->grid().atLocalBins({2}) = vol2.get();
1942
1943   auto disc2 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 100_mm,
1944                                                 150_mm, 30_degree);
1945
1946   auto grid2 = GridPortalLink::make(disc2, AxisDirection::AxisPhi,
1947                                     Axis{AxisBound, -30_degree, 30_degree, 2});
1948
1949   grid2->grid().atLocalBins({1}) = vol3.get();
1950   grid2->grid().atLocalBins({2}) = vol4.get();
1951
1952   auto mergedPtr =
1953       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisR, *logger);
1954
1955   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1956   BOOST_REQUIRE_NE(merged, nullptr);
1957
1958   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2);
1959   const auto& axis1 = *merged->grid().axes().front();
1960   const auto& axis2 = *merged->grid().axes().back();
1961   Axis axisExpectedR{AxisBound, {30_mm, 65_mm, 100_mm, 150_mm}};
1962   Axis axisExpectedPhi{AxisBound, -30_degree, 30_degree, 2};
1963   BOOST_CHECK_EQUAL(axis1, axisExpectedR);
1964   BOOST_CHECK_EQUAL(axis2, axisExpectedPhi);
1965
1966   std::vector<std::pair<Vector2, TrackingVolume*>> locations = {
1967       {{40_mm, -15_degree}, vol1.get()},  {{40_mm, 15_degree}, vol1.get()},
1968       {{90_mm, -15_degree}, vol2.get()},  {{90_mm, 15_degree}, vol2.get()},
1969
1970       {{110_mm, -15_degree}, vol3.get()}, {{110_mm, 15_degree}, vol4.get()},
1971       {{140_mm, -15_degree}, vol3.get()}, {{140_mm, 15_degree}, vol4.get()},
1972   };
1973
1974   for (const auto& [loc, vol] : locations) {
1975     BOOST_TEST_CONTEXT(loc.transpose())
1976     BOOST_CHECK_EQUAL(merged->resolveVolume(gctx, loc).value(), vol);
1977   }
1978
1979   grid1->printContents(std::cout);
1980   grid2->printContents(std::cout);
1981   merged->printContents(std::cout);
1982 }
1983
1984 BOOST_AUTO_TEST_CASE(PhiDirection) {
1985   // Volumes for bin content checking
1986
1987   auto vol1 = makeDummyVolume();
1988   auto vol2 = makeDummyVolume();
1989   auto vol3 = makeDummyVolume();
1990   auto vol4 = makeDummyVolume();
1991
1992   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1993                                                 100_mm, 30_degree);
1994
1995   auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisR,
1996                                     Axis{AxisBound, 30_mm, 100_mm, 2});
1997
1998   grid1->grid().atLocalBins({1}) = vol1.get();
1999   grid1->grid().atLocalBins({2}) = vol2.get();
2000
2001   auto disc2 = Surface::makeShared<DiscSurface>(
2002       Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
2003       60_degree);
2004
2005   auto grid2 = GridPortalLink::make(disc2, AxisDirection::AxisPhi,
2006                                     Axis{AxisBound, -60_degree, 60_degree, 2});
2007
2008   grid2->grid().atLocalBins({1}) = vol3.get();
2009   grid2->grid().atLocalBins({2}) = vol4.get();
2010
2011   auto mergedPtr =
2012       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
2013
2014   using merged_type =
2015       GridPortalLinkT<Axis<AxisType::Equidistant, AxisBoundaryType::Bound>,
2016                       Axis<AxisType::Equidistant, AxisBoundaryType::Bound>>;
2017
2018   const auto* merged = dynamic_cast<const merged_type*>(mergedPtr.get());
2019   BOOST_REQUIRE_NE(merged, nullptr);
2020
2021   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2);
2022   const auto& axis1 = *merged->grid().axes().front();
2023   const auto& axis2 = *merged->grid().axes().back();
2024   Axis axisExpectedR{AxisBound, 30_mm, 100_mm, 2};
2025   BOOST_CHECK_EQUAL(axis1, axisExpectedR);
2026
2027   BOOST_CHECK_CLOSE(axis2.getMin(), -90_degree, 1e-6);
2028   BOOST_CHECK_CLOSE(axis2.getMax(), 90_degree, 1e-6);
2029   BOOST_CHECK_EQUAL(axis2.getNBins(), 3);
2030   BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant);
2031   BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound);
2032
2033   std::vector<std::pair<Vector2, TrackingVolume*>> locations = {
2034       {{40_mm, 45_degree}, vol1.get()},  {{40_mm, 0_degree}, vol4.get()},
2035       {{40_mm, -80_degree}, vol3.get()}, {{90_mm, 45_degree}, vol2.get()},
2036       {{90_mm, 0_degree}, vol4.get()},   {{90_mm, -80_degree}, vol3.get()},
2037   };
2038
2039   grid1->printContents(std::cout);
2040   grid2->printContents(std::cout);
2041   merged->printContents(std::cout);
2042
2043   for (const auto& [loc, vol] : locations) {
2044     BOOST_TEST_CONTEXT((Vector2{loc[0], loc[1] / 1_degree}.transpose()))
2045     BOOST_CHECK_EQUAL(merged->resolveVolume(gctx, loc).value(), vol);
2046   }
2047 }
2048
2049 BOOST_AUTO_TEST_SUITE_END()  // MergeCrossDisc
2050
2051 BOOST_AUTO_TEST_SUITE_END()  // GridMerging
2052
2053 BOOST_AUTO_TEST_CASE(CompositeConstruction) {
2054   // Arbitrary volumes for testing only
2055   auto vol1 = std::make_shared<TrackingVolume>(
2056       Transform3::Identity(),
2057       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2058
2059   auto vol2 = std::make_shared<TrackingVolume>(
2060       Transform3::Identity(),
2061       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2062
2063   auto vol3 = std::make_shared<TrackingVolume>(
2064       Transform3::Identity(),
2065       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2066
2067   auto disc1 =
2068       Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 60_mm);
2069
2070   auto trivial1 = std::make_unique<TrivialPortalLink>(disc1, *vol1);
2071
2072   auto disc2 =
2073       Surface::makeShared<DiscSurface>(Transform3::Identity(), 60_mm, 90_mm);
2074   auto trivial2 = std::make_unique<TrivialPortalLink>(disc2, *vol2);
2075
2076   auto composite = std::make_unique<CompositePortalLink>(
2077       copy(trivial1), copy(trivial2), AxisDirection::AxisR);
2078
2079   auto compositeCopy = std::make_unique<CompositePortalLink>(
2080       copy(trivial1), copy(trivial2), AxisDirection::AxisR);
2081
2082   BOOST_CHECK_EQUAL(
2083       composite->resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(),
2084       vol1.get());
2085   BOOST_CHECK_EQUAL(
2086       composite->resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(),
2087       vol2.get());
2088
2089   BOOST_CHECK_EQUAL(composite->depth(), 1);
2090
2091   // Test exception on different surface types
2092   auto cyl = Surface::makeShared<CylinderSurface>(Transform3::Identity(), 30_mm,
2093                                                   40_mm);
2094   auto trivialCyl = std::make_unique<TrivialPortalLink>(cyl, *vol3);
2095   BOOST_CHECK_THROW(CompositePortalLink(copy(trivial1), copy(trivialCyl),
2096                                         AxisDirection::AxisR),
2097                     std::invalid_argument);
2098
2099   auto disc3 =
2100       Surface::makeShared<DiscSurface>(Transform3::Identity(), 90_mm, 120_mm);
2101   auto trivial3 = std::make_unique<TrivialPortalLink>(disc3, *vol3);
2102
2103   // Test exception on un-mergable surfaces
2104   BOOST_CHECK_THROW(
2105       CompositePortalLink(copy(trivial1), copy(trivial3), AxisDirection::AxisR),
2106       SurfaceMergingException);
2107
2108   // Composite with a composite (this should work regardless of flattening)
2109
2110   CompositePortalLink composite2(std::move(composite), copy(trivial3),
2111                                  AxisDirection::AxisR, false);
2112
2113   BOOST_CHECK_EQUAL(
2114       composite2.resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(),
2115       vol1.get());
2116   BOOST_CHECK_EQUAL(
2117       composite2.resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(),
2118       vol2.get());
2119   BOOST_CHECK_EQUAL(
2120       composite2.resolveVolume(gctx, Vector2{100_mm, 0_degree}).value(),
2121       vol3.get());
2122
2123   // Two without flattening
2124   BOOST_CHECK_EQUAL(composite2.depth(), 2);
2125
2126   CompositePortalLink composite2Flat(std::move(compositeCopy), copy(trivial3),
2127                                      AxisDirection::AxisR, true);
2128
2129   // One because of flattening
2130   BOOST_CHECK_EQUAL(composite2Flat.depth(), 1);
2131
2132   BOOST_CHECK_EQUAL(
2133       composite2Flat.resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(),
2134       vol1.get());
2135   BOOST_CHECK_EQUAL(
2136       composite2Flat.resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(),
2137       vol2.get());
2138   BOOST_CHECK_EQUAL(
2139       composite2Flat.resolveVolume(gctx, Vector2{100_mm, 0_degree}).value(),
2140       vol3.get());
2141 }
2142
2143 BOOST_AUTO_TEST_SUITE(PortalMerging)
2144
2145 BOOST_DATA_TEST_CASE(TrivialTrivial,
2146                      (boost::unit_test::data::make(0, -135, 180, 45) *
2147                       boost::unit_test::data::make(Vector3{0_mm, 0_mm, 0_mm},
2148                                                    Vector3{20_mm, 0_mm, 0_mm},
2149                                                    Vector3{0_mm, 20_mm, 0_mm},
2150                                                    Vector3{20_mm, 20_mm, 0_mm},
2151                                                    Vector3{0_mm, 0_mm, 20_mm})),
2152                      angle, offset) {
2153   Transform3 base =
2154       AngleAxis3(angle * 1_degree, Vector3::UnitX()) * Translation3(offset);
2155
2156   BOOST_TEST_CONTEXT("RDirection") {
2157     auto vol1 = std::make_shared<TrackingVolume>(
2158         base, std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2159
2160     auto vol2 = std::make_shared<TrackingVolume>(
2161         base, std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2162
2163     auto vol3 = std::make_shared<TrackingVolume>(
2164         base, std::make_shared<CylinderVolumeBounds>(40_mm, 50_mm, 100_mm));
2165
2166     auto disc1 =
2167         Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 60_mm);
2168
2169     auto disc2 =
2170         Surface::makeShared<DiscSurface>(Transform3::Identity(), 60_mm, 90_mm);
2171
2172     auto disc3 =
2173         Surface::makeShared<DiscSurface>(Transform3::Identity(), 90_mm, 120_mm);
2174
2175     auto trivial1 = std::make_unique<TrivialPortalLink>(disc1, *vol1);
2176     BOOST_REQUIRE(trivial1);
2177     auto trivial2 = std::make_unique<TrivialPortalLink>(disc2, *vol2);
2178     BOOST_REQUIRE(trivial2);
2179     auto trivial3 = std::make_unique<TrivialPortalLink>(disc3, *vol3);
2180     BOOST_REQUIRE(trivial3);
2181
2182     auto grid1 = trivial1->makeGrid(AxisDirection::AxisR);
2183     auto compGridTrivial = PortalLinkBase::merge(
2184         std::move(grid1), std::make_unique<TrivialPortalLink>(*trivial2),
2185         AxisDirection::AxisR, *logger);
2186     BOOST_REQUIRE(compGridTrivial);
2187     BOOST_CHECK_EQUAL(dynamic_cast<CompositePortalLink&>(*compGridTrivial)
2188                           .makeGrid(gctx, *logger),
2189                       nullptr);
2190
2191     auto composite =
2192         PortalLinkBase::merge(std::move(trivial1), std::move(trivial2),
2193                               AxisDirection::AxisR, *logger);
2194     BOOST_REQUIRE(composite);
2195
2196     auto grid12 =
2197         dynamic_cast<CompositePortalLink&>(*composite).makeGrid(gctx, *logger);
2198     BOOST_REQUIRE(grid12);
2199
2200     BOOST_CHECK_EQUAL(
2201         grid12->resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(),
2202         vol1.get());
2203
2204     BOOST_CHECK_EQUAL(
2205         grid12->resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(),
2206         vol2.get());
2207
2208     composite = PortalLinkBase::merge(std::move(composite), std::move(trivial3),
2209                                       AxisDirection::AxisR, *logger);
2210     BOOST_REQUIRE(composite);
2211
2212     auto grid123 =
2213         dynamic_cast<CompositePortalLink&>(*composite).makeGrid(gctx, *logger);
2214     BOOST_REQUIRE(grid123);
2215
2216     BOOST_CHECK_EQUAL(
2217         grid123->resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(),
2218         vol1.get());
2219
2220     BOOST_CHECK_EQUAL(
2221         grid123->resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(),
2222         vol2.get());
2223
2224     BOOST_CHECK_EQUAL(
2225         grid123->resolveVolume(gctx, Vector2{100_mm, 0_degree}).value(),
2226         vol3.get());
2227   }
2228
2229   BOOST_TEST_CONTEXT("ZDirection") {
2230     auto vol1 = std::make_shared<TrackingVolume>(
2231         base, std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2232
2233     auto vol2 = std::make_shared<TrackingVolume>(
2234         base * Translation3{Vector3::UnitZ() * 200},
2235         std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2236
2237     auto vol3 = std::make_shared<TrackingVolume>(
2238         base * Translation3{Vector3::UnitZ() * 400},
2239         std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2240
2241     auto cyl1 = Surface::makeShared<CylinderSurface>(base, 40_mm, 100_mm);
2242
2243     auto cyl2 = Surface::makeShared<CylinderSurface>(
2244         base * Translation3{Vector3::UnitZ() * 200}, 40_mm, 100_mm);
2245
2246     auto cyl3 = Surface::makeShared<CylinderSurface>(
2247         base * Translation3{Vector3::UnitZ() * 400}, 40_mm, 100_mm);
2248
2249     auto trivial1 = std::make_unique<TrivialPortalLink>(cyl1, *vol1);
2250     BOOST_REQUIRE(trivial1);
2251     auto trivial2 = std::make_unique<TrivialPortalLink>(cyl2, *vol2);
2252     BOOST_REQUIRE(trivial2);
2253     auto trivial3 = std::make_unique<TrivialPortalLink>(cyl3, *vol3);
2254     BOOST_REQUIRE(trivial3);
2255
2256     auto grid1 = trivial1->makeGrid(AxisDirection::AxisZ);
2257     auto compGridTrivial = PortalLinkBase::merge(
2258         std::move(grid1), std::make_unique<TrivialPortalLink>(*trivial2),
2259         AxisDirection::AxisZ, *logger);
2260     BOOST_REQUIRE(compGridTrivial);
2261     BOOST_CHECK_EQUAL(dynamic_cast<CompositePortalLink&>(*compGridTrivial)
2262                           .makeGrid(gctx, *logger),
2263                       nullptr);
2264
2265     auto composite =
2266         PortalLinkBase::merge(std::move(trivial1), std::move(trivial2),
2267                               AxisDirection::AxisZ, *logger);
2268     BOOST_REQUIRE(composite);
2269
2270     auto grid12 =
2271         dynamic_cast<CompositePortalLink&>(*composite).makeGrid(gctx, *logger);
2272     BOOST_REQUIRE(grid12);
2273
2274     BOOST_CHECK_EQUAL(
2275         grid12->resolveVolume(gctx, Vector2{40_mm, -40_mm}).value(),
2276         vol1.get());
2277
2278     BOOST_CHECK_EQUAL(
2279         grid12->resolveVolume(gctx, Vector2{40_mm, 40_mm}).value(), vol2.get());
2280
2281     composite = PortalLinkBase::merge(std::move(composite), std::move(trivial3),
2282                                       AxisDirection::AxisZ, *logger);
2283     BOOST_REQUIRE(composite);
2284
2285     auto grid123 =
2286         dynamic_cast<CompositePortalLink&>(*composite).makeGrid(gctx, *logger);
2287     BOOST_REQUIRE(grid123);
2288
2289     BOOST_CHECK_EQUAL(
2290         grid123->resolveVolume(gctx, Vector2{40_mm, -110_mm}).value(),
2291         vol1.get());
2292
2293     BOOST_CHECK_EQUAL(
2294         grid123->resolveVolume(gctx, Vector2{40_mm, -10_mm}).value(),
2295         vol2.get());
2296
2297     BOOST_CHECK_EQUAL(
2298         grid123->resolveVolume(gctx, Vector2{40_mm, 190_mm}).value(),
2299         vol3.get());
2300   }
2301 }
2302
2303 BOOST_AUTO_TEST_CASE(TrivialGridR) {
2304   auto vol1 = std::make_shared<TrackingVolume>(
2305       Transform3::Identity(),
2306       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2307
2308   auto vol2 = std::make_shared<TrackingVolume>(
2309       Transform3::Identity(),
2310       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2311
2312   auto disc1 =
2313       Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 60_mm);
2314
2315   auto disc2 =
2316       Surface::makeShared<DiscSurface>(Transform3::Identity(), 60_mm, 90_mm);
2317
2318   auto trivial = std::make_unique<TrivialPortalLink>(disc2, *vol2);
2319   BOOST_REQUIRE(trivial);
2320
2321   auto gridPhi = GridPortalLink::make(
2322       disc1, AxisDirection::AxisPhi,
2323       Axis{AxisClosed, -std::numbers::pi, std::numbers::pi, 2});
2324   gridPhi->setVolume(vol1.get());
2325
2326   auto gridR = GridPortalLink::make(disc1, AxisDirection::AxisR,
2327                                     Axis{AxisBound, 30_mm, 60_mm, 2});
2328   gridR->setVolume(vol1.get());
2329
2330   BOOST_TEST_CONTEXT("Colinear") {
2331     auto merged = PortalLinkBase::merge(copy(trivial), copy(gridR),
2332                                         AxisDirection::AxisR, *logger);
2333     BOOST_REQUIRE(merged);
2334     BOOST_CHECK_NE(dynamic_cast<CompositePortalLink*>(merged.get()), nullptr);
2335   }
2336
2337   BOOST_TEST_CONTEXT("Orthogonal") {
2338     auto merged = PortalLinkBase::merge(copy(gridPhi), copy(trivial),
2339                                         AxisDirection::AxisR, *logger);
2340     BOOST_REQUIRE(merged);
2341     BOOST_CHECK_NE(dynamic_cast<CompositePortalLink*>(merged.get()), nullptr);
2342   }
2343 }
2344
2345 BOOST_AUTO_TEST_CASE(TrivialGridPhi) {
2346   auto vol1 = std::make_shared<TrackingVolume>(
2347       Transform3::Identity(),
2348       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2349
2350   auto vol2 = std::make_shared<TrackingVolume>(
2351       Transform3::Identity(),
2352       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2353
2354   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
2355                                                 100_mm, 30_degree);
2356
2357   auto disc2 = Surface::makeShared<DiscSurface>(
2358       Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
2359       60_degree);
2360
2361   auto trivial = std::make_unique<TrivialPortalLink>(disc2, *vol2);
2362   BOOST_REQUIRE(trivial);
2363
2364   auto gridPhi = GridPortalLink::make(
2365       disc1, AxisDirection::AxisPhi, Axis{AxisBound, -30_degree, 30_degree, 2});
2366   gridPhi->setVolume(vol1.get());
2367
2368   auto gridR = GridPortalLink::make(disc1, AxisDirection::AxisR,
2369                                     Axis{AxisBound, 30_mm, 100_mm, 2});
2370   gridR->setVolume(vol1.get());
2371
2372   BOOST_TEST_CONTEXT("Colinear") {
2373     auto merged = PortalLinkBase::merge(copy(trivial), copy(gridPhi),
2374                                         AxisDirection::AxisPhi, *logger);
2375     BOOST_REQUIRE(merged);
2376     BOOST_CHECK_NE(dynamic_cast<CompositePortalLink*>(merged.get()), nullptr);
2377   }
2378
2379   BOOST_TEST_CONTEXT("Orthogonal") {
2380     auto merged = PortalLinkBase::merge(copy(gridR), copy(trivial),
2381                                         AxisDirection::AxisPhi, *logger);
2382     BOOST_REQUIRE(merged);
2383     BOOST_CHECK_NE(dynamic_cast<CompositePortalLink*>(merged.get()), nullptr);
2384   }
2385 }
2386
2387 BOOST_AUTO_TEST_CASE(CompositeOther) {
2388   auto vol1 = std::make_shared<TrackingVolume>(
2389       Transform3::Identity(),
2390       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2391
2392   auto vol2 = std::make_shared<TrackingVolume>(
2393       Transform3::Identity(),
2394       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2395
2396   auto vol3 = std::make_shared<TrackingVolume>(
2397       Transform3::Identity(),
2398       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2399
2400   auto disc1 =
2401       Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 60_mm);
2402   auto disc2 =
2403       Surface::makeShared<DiscSurface>(Transform3::Identity(), 60_mm, 90_mm);
2404   auto disc3 =
2405       Surface::makeShared<DiscSurface>(Transform3::Identity(), 90_mm, 120_mm);
2406
2407   auto grid1 = GridPortalLink::make(disc1, *vol1, AxisDirection::AxisR);
2408   auto trivial2 = std::make_unique<TrivialPortalLink>(disc2, *vol2);
2409
2410   auto composite12 = std::make_unique<CompositePortalLink>(
2411       std::move(grid1), std::move(trivial2), AxisDirection::AxisR);
2412
2413   BOOST_CHECK_EQUAL(composite12->depth(), 1);
2414   BOOST_CHECK_EQUAL(composite12->size(), 2);
2415
2416   auto trivial3 = std::make_unique<TrivialPortalLink>(disc3, *vol3);
2417
2418   auto composite123Ptr =
2419       PortalLinkBase::merge(std::move(composite12), std::move(trivial3),
2420                             AxisDirection::AxisR, *logger);
2421
2422   const auto* composite123 =
2423       dynamic_cast<const CompositePortalLink*>(composite123Ptr.get());
2424   BOOST_REQUIRE(composite123);
2425
2426   BOOST_CHECK_EQUAL(composite123->depth(), 1);
2427
2428   BOOST_CHECK_EQUAL(
2429       composite123->resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(),
2430       vol1.get());
2431   BOOST_CHECK_EQUAL(
2432       composite123->resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(),
2433       vol2.get());
2434   BOOST_CHECK_EQUAL(
2435       composite123->resolveVolume(gctx, Vector2{100_mm, 0_degree}).value(),
2436       vol3.get());
2437
2438   BOOST_CHECK_EQUAL(composite123->size(), 3);
2439 }
2440
2441 BOOST_AUTO_TEST_SUITE_END()  // PortalMerging
2442
2443 BOOST_AUTO_TEST_SUITE_END()  // Geometry
2444
2445 }  // namespace Acts::Test