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/CuboidVolumeBounds.hpp"
0019 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
0020 #include "Acts/Geometry/GridPortalLink.hpp"
0021 #include "Acts/Geometry/TrackingVolume.hpp"
0022 #include "Acts/Geometry/TrivialPortalLink.hpp"
0023 #include "Acts/Surfaces/CylinderSurface.hpp"
0024 #include "Acts/Surfaces/RadialBounds.hpp"
0025 #include "Acts/Surfaces/RectangleBounds.hpp"
0026 #include "Acts/Surfaces/SurfaceMergingException.hpp"
0027 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0028 #include "Acts/Utilities/AxisDefinitions.hpp"
0029 #include "Acts/Utilities/ThrowAssert.hpp"
0030
0031 #include <cstdio>
0032 #include <iostream>
0033 #include <memory>
0034 #include <numbers>
0035 #include <stdexcept>
0036
0037 using namespace Acts::UnitLiterals;
0038
0039 namespace Acts::Test {
0040
0041 auto logger = Acts::getDefaultLogger("UnitTests", Acts::Logging::VERBOSE);
0042
0043 struct Fixture {
0044   Logging::Level m_level;
0045   Fixture() {
0046     m_level = Acts::Logging::getFailureThreshold();
0047     Acts::Logging::setFailureThreshold(Acts::Logging::FATAL);
0048   }
0049
0050   ~Fixture() { Acts::Logging::setFailureThreshold(m_level); }
0051 };
0052
0053 std::shared_ptr<TrackingVolume> makeDummyVolume() {
0054   return std::make_shared<TrackingVolume>(
0055       Transform3::Identity(),
0056       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
0057 }
0058
0059 GeometryContext gctx;
0060
0061 template <typename T>
0062 std::unique_ptr<T> copy(const std::unique_ptr<T>& p) {
0063   return std::make_unique<T>(*p);
0064 }
0065
0066 template <typename link_t>
0067 void visitBins(const link_t& link,
0068                const std::function<void(const TrackingVolume*)>& func) {
0069   auto& grid = link.grid();
0070   auto loc = grid.numLocalBins();
0071   if constexpr (std::decay_t<decltype(grid)>::DIM == 1) {
0072     for (std::size_t i = 1; i <= loc[0]; i++) {
0073       func(grid.atLocalBins({i}));
0074     }
0075   } else {
0076     for (std::size_t i = 1; i <= loc[0]; i++) {
0077       for (std::size_t j = 1; j <= loc[1]; j++) {
0078         func(grid.atLocalBins({i, j}));
0079       }
0080     }
0081   }
0082 }
0083
0084 BOOST_FIXTURE_TEST_SUITE(Geometry, Fixture)
0085
0086 BOOST_AUTO_TEST_SUITE(GridConstruction)
0087
0088 BOOST_AUTO_TEST_CASE(Cylinder) {
0089   BOOST_TEST_CONTEXT("1D") {
0090     auto cyl = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
0091                                                     30_mm, 100_mm);
0092
0093     // Volume for bin testing
0094     auto vol = std::make_shared<TrackingVolume>(
0095         Transform3::Identity(),
0096         std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
0097
0098     // Incompatible binning
0099     BOOST_CHECK_THROW(GridPortalLink::make(cyl, AxisDirection::AxisZ,
0100                                            Axis{AxisBound, 0, 5, 5}),
0101                       std::invalid_argument);
0102
0103     auto grid1dCyl = GridPortalLink::make(cyl, AxisDirection::AxisZ,
0104                                           Axis{AxisBound, -100_mm, 100_mm, 10});
0105     BOOST_REQUIRE(grid1dCyl);
0106     grid1dCyl->setVolume(vol.get());
0107
0108     // Throws because non-closed axis
0109     BOOST_CHECK_THROW(GridPortalLink::make(cyl, AxisDirection::AxisRPhi,
0110                                            Axis{AxisBound, -180_degree * 30_mm,
0111                                                 180_degree * 30_mm, 10}),
0112                       std::invalid_argument);
0113
0114     auto grid1dCylRPhi = GridPortalLink::make(
0115         cyl, AxisDirection::AxisRPhi,
0116         Axis{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 10});
0117     BOOST_REQUIRE_NE(grid1dCylRPhi, nullptr);
0118     grid1dCylRPhi->setVolume(vol.get());
0119
0120     Axis axisExpected{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 10};
0121     BOOST_CHECK_EQUAL(grid1dCylRPhi->grid().axes().size(), 1);
0122     const auto& axis = *grid1dCylRPhi->grid().axes().front();
0123     BOOST_CHECK_EQUAL(axis, axisExpected);
0124
0125     // Another cylinder, shifted in z
0126     auto cyl2 = Surface::makeShared<CylinderSurface>(
0127         Transform3{Translation3{Vector3::UnitZ() * 150_mm}}, 30_mm, 50_mm);
0128
0129     auto grid1dCyl2 = GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0130                                            Axis{AxisBound, -50_mm, 50_mm, 5});
0131
0132     // Test exception on cylinder with non-zero average phi
0133     auto cylNonZeroAverage = Surface::makeShared<CylinderSurface>(
0134         Transform3::Identity(), 30_mm, 100_mm, 20_degree, 45_degree);
0135     BOOST_CHECK_THROW(
0136         GridPortalLink::make(cylNonZeroAverage, AxisDirection::AxisZ,
0137                              Axis{AxisBound, -100_mm, 100_mm, 10}),
0138         std::invalid_argument);
0139
0140     auto checkAllBins = [&](const auto& link) {
0141       visitBins(link, [&](const TrackingVolume* content) {
0142         BOOST_CHECK_EQUAL(content, vol.get());
0143       });
0144     };
0145
0146     checkAllBins(*grid1dCyl);
0147     checkAllBins(*grid1dCylRPhi);
0148
0149     // Extend to a 2D grid with auto phi binning
0150
0151     auto grid2dCyl1 = grid1dCyl->extendTo2d(nullptr);
0152     BOOST_REQUIRE(grid2dCyl1);
0153     BOOST_CHECK_EQUAL(grid2dCyl1->grid().axes().size(), 2);
0154     BOOST_CHECK_EQUAL(grid2dCyl1->surface().bounds(), cyl->bounds());
0155     const auto* axis1 = grid2dCyl1->grid().axes().front();
0156     const auto* axis2 = grid2dCyl1->grid().axes().back();
0157
0158     Axis axis1Expected{AxisClosed, -std::numbers::pi * 30_mm,
0159                        std::numbers::pi * 30_mm, 1};
0160     BOOST_CHECK_EQUAL(*axis1, axis1Expected);
0161     Axis axis2Expected{AxisBound, -100_mm, 100_mm, 10};
0162     BOOST_CHECK_EQUAL(*axis2, axis2Expected);
0163
0164     auto& concrete = dynamic_cast<
0165         GridPortalLinkT<decltype(axis1Expected), decltype(axis2Expected)>&>(
0166         *grid2dCyl1);
0167
0168     checkAllBins(concrete);
0169
0170     Axis axis1Explicit{AxisClosed, -std::numbers::pi * 30_mm,
0171                        std::numbers::pi * 30_mm, 13};
0172     auto grid2dCyl1Explicit = grid1dCyl->extendTo2d(&axis1Explicit);
0173     BOOST_REQUIRE(grid2dCyl1Explicit);
0174     BOOST_CHECK_EQUAL(grid2dCyl1Explicit->grid().axes().size(), 2);
0175     axis1 = grid2dCyl1Explicit->grid().axes().front();
0176     axis2 = grid2dCyl1Explicit->grid().axes().back();
0177
0178     BOOST_CHECK_EQUAL(*axis1, axis1Explicit);
0179     BOOST_CHECK_EQUAL(*axis2, axis2Expected);
0180
0181     auto& concrete2 = dynamic_cast<
0182         GridPortalLinkT<decltype(axis1Explicit), decltype(axis2Expected)>&>(
0183         *grid2dCyl1Explicit);
0184
0185     checkAllBins(concrete2);
0186
0187     auto cylPhi = Surface::makeShared<CylinderSurface>(
0188         Transform3::Identity(), 30_mm, 100_mm, 45_degree);
0189     std::unique_ptr<GridPortalLink> grid1dCylPhi = GridPortalLink::make(
0190         cylPhi, AxisDirection::AxisZ, Axis{AxisBound, -100_mm, 100_mm, 10});
0191
0192     grid1dCylPhi->setVolume(vol.get());
0193
0194     // Check that phi sector portal does not accept closed axis
0195     BOOST_CHECK_THROW(GridPortalLink::make(cylPhi, AxisDirection::AxisRPhi,
0196                                            Axis{AxisClosed, -45_degree * 30_mm,
0197                                                 45_degree * 30_mm, 10}),
0198                       std::invalid_argument);
0199
0200     auto grid2dCylPhi = grid1dCylPhi->extendTo2d(nullptr);
0201     BOOST_CHECK_EQUAL(grid2dCylPhi->grid().axes().size(), 2);
0202     BOOST_CHECK_EQUAL(grid2dCylPhi->surface().bounds(), cylPhi->bounds());
0203     const auto* axis1Phi = grid2dCylPhi->grid().axes().front();
0204     const auto* axis2Phi = grid2dCylPhi->grid().axes().back();
0205
0206     Axis axis1PhiExpected{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 1};
0207     BOOST_CHECK_EQUAL(*axis1Phi, axis1PhiExpected);
0208     Axis axis2PhiExpected{AxisBound, -100_mm, 100_mm, 10};
0209     BOOST_CHECK_EQUAL(*axis2Phi, axis2PhiExpected);
0210
0211     auto& concrete3 =
0212         dynamic_cast<GridPortalLinkT<decltype(axis1PhiExpected),
0213                                      decltype(axis2PhiExpected)>&>(
0214             *grid2dCylPhi);
0215
0216     checkAllBins(concrete3);
0217
0218     Axis axis1PhiExplicit{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 13};
0219     auto grid2dCylPhiExplicit = grid1dCylPhi->extendTo2d(&axis1PhiExplicit);
0220     BOOST_REQUIRE(grid2dCylPhiExplicit);
0221     BOOST_CHECK_EQUAL(grid2dCylPhiExplicit->grid().axes().size(), 2);
0222     axis1Phi = grid2dCylPhiExplicit->grid().axes().front();
0223     axis2Phi = grid2dCylPhiExplicit->grid().axes().back();
0224     BOOST_CHECK_EQUAL(*axis1Phi, axis1PhiExplicit);
0225     BOOST_CHECK_EQUAL(*axis2Phi, axis2PhiExpected);
0226
0227     auto& concrete4 =
0228         dynamic_cast<GridPortalLinkT<decltype(axis1PhiExplicit),
0229                                      decltype(axis2PhiExpected)>&>(
0230             *grid2dCylPhiExplicit);
0231
0232     checkAllBins(concrete4);
0233   }
0234
0235   BOOST_TEST_CONTEXT("2D") {
0236     auto cyl = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
0237                                                     30_mm, 100_mm, 45_degree);
0238
0239     // z bad, rphi bad
0240     BOOST_CHECK_THROW(GridPortalLink::make(cyl, Axis{AxisBound, 1, 2, 5},
0241                                            Axis{AxisBound, 3_mm, 4_mm, 5}),
0242                       std::invalid_argument);
0243
0244     // z good, rphi bad
0245     BOOST_CHECK_THROW(GridPortalLink::make(cyl, Axis{AxisBound, 3_mm, 4_mm, 5},
0246                                            Axis{AxisBound, -100_mm, 100_m, 5}),
0247                       std::invalid_argument);
0248
0249     // z bad, rphi good
0250     BOOST_CHECK_THROW(
0251         GridPortalLink::make(
0252             cyl, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5},
0253             Axis{AxisBound, -80_mm, 100_mm, 5}),
0254         std::invalid_argument);
0255
0256     auto grid1 = GridPortalLink::make(
0257         cyl, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5},
0258         Axis{AxisBound, -100_mm, 100_mm, 5});
0259
0260     auto cylFull = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
0261                                                         30_mm, 100_mm);
0262
0263     // Throws because non-closed axis
0264     BOOST_CHECK_THROW(GridPortalLink::make(cylFull,
0265                                            Axis{AxisBound, -180_degree * 30_mm,
0266                                                 180_degree * 30_mm, 5},
0267                                            Axis{AxisBound, -100_mm, 100_mm, 5}),
0268                       std::invalid_argument);
0269
0270     auto gridFull = GridPortalLink::make(
0271         cylFull, Axis{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 5},
0272         Axis{AxisBound, -100_mm, 100_mm, 5});
0273
0274     BOOST_CHECK_EQUAL(gridFull->grid().axes().size(), 2);
0275     BOOST_CHECK_EQUAL(gridFull->grid().axes().size(), 2);
0276     Axis axisFullExpected{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm,
0277                           5};
0278     const auto& axisFull = *gridFull->grid().axes().front();
0279     BOOST_CHECK_EQUAL(axisFull, axisFullExpected);
0280   }
0281 }
0282
0283 BOOST_AUTO_TEST_CASE(Disc) {
0284   using enum AxisType;
0285   using enum AxisBoundaryType;
0286   BOOST_TEST_CONTEXT("1D") {
0287     auto disc1 =
0288         Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 100_mm);
0289
0290     // Volume for bin testing
0291     auto vol = std::make_shared<TrackingVolume>(
0292         Transform3::Identity(),
0293         std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
0294
0295     BOOST_CHECK_THROW(GridPortalLink::make(disc1, AxisDirection::AxisZ,
0296                                            Axis{AxisBound, 30_mm, 100_mm, 3}),
0297                       std::invalid_argument);
0298
0299     // Check exception for full disc and non-closed phi axis
0300     BOOST_CHECK_THROW(
0301         GridPortalLink::make(disc1, AxisDirection::AxisPhi,
0302                              Axis{AxisBound, -180_degree, 180_degree, 3}),
0303         std::invalid_argument);
0304
0305     auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisR,
0306                                       Axis{AxisBound, 30_mm, 100_mm, 3});
0307     BOOST_REQUIRE_NE(grid1, nullptr);
0308     BOOST_CHECK_EQUAL(grid1->grid().axes().size(), 1);
0309     const auto& axis = *grid1->grid().axes().front();
0310     Axis axis1Expected{AxisBound, 30_mm, 100_mm, 3};
0311     BOOST_CHECK_EQUAL(axis, axis1Expected);
0312
0313     grid1->setVolume(vol.get());
0314
0315     auto discPhi = Surface::makeShared<DiscSurface>(Transform3::Identity(),
0316                                                     30_mm, 100_mm, 45_degree);
0317
0318     // Check thet disc with phi sector does not accept closed axis
0319     BOOST_CHECK_THROW(
0320         GridPortalLink::make(discPhi, AxisDirection::AxisPhi,
0321                              Axis{AxisClosed, -45_degree, 45_degree, 3}),
0322         std::invalid_argument);
0323
0324     auto gridPhi =
0325         GridPortalLink::make(discPhi, AxisDirection::AxisPhi,
0326                              Axis{AxisBound, -45_degree, 45_degree, 3});
0327     BOOST_REQUIRE_NE(gridPhi, nullptr);
0328     gridPhi->setVolume(vol.get());
0329
0330     // Test exception on disc with non-zero average phi
0331     auto discNonZeroAverage = Surface::makeShared<DiscSurface>(
0332         Transform3::Identity(),
0333         std::make_shared<RadialBounds>(30_mm, 100_mm, 45_degree, 75_degree));
0334     BOOST_CHECK_THROW(
0335         GridPortalLink::make(discNonZeroAverage, AxisDirection::AxisR,
0336                              Axis{AxisBound, 30_mm, 100_mm, 3}),
0337         std::invalid_argument);
0338
0339     BOOST_CHECK_EQUAL(gridPhi->grid().axes().size(), 1);
0340     const auto& axisPhi = *gridPhi->grid().axes().front();
0341     Axis axisPhi1Expected{AxisBound, -45_degree, 45_degree, 3};
0342     BOOST_CHECK_EQUAL(axisPhi, axisPhi1Expected);
0343
0344     auto checkAllBins = [&](const auto& grid) {
0345       visitBins(grid, [&](const TrackingVolume* content) {
0346         BOOST_CHECK_EQUAL(content, vol.get());
0347       });
0348     };
0349
0350     checkAllBins(*grid1);
0351     checkAllBins(*gridPhi);
0352
0353     // Test making 2D grids from the 1D ones
0354     auto grid2d = grid1->extendTo2d(nullptr);
0355     BOOST_REQUIRE(grid2d);
0356     BOOST_CHECK_EQUAL(grid2d->grid().axes().size(), 2);
0357     const auto* axis1 = grid2d->grid().axes().front();
0358     BOOST_CHECK_EQUAL(*axis1, axis1Expected);
0359     const auto* axis2 = grid2d->grid().axes().back();
0360     BOOST_CHECK_CLOSE(axis2->getMin(), -180_degree, 1e-6);
0361     BOOST_CHECK_CLOSE(axis2->getMax(), 180_degree, 1e-6);
0362     BOOST_CHECK_EQUAL(axis2->getNBins(), 1);
0363     BOOST_CHECK_EQUAL(axis2->getType(), AxisType::Equidistant);
0364     BOOST_CHECK_EQUAL(axis2->getBoundaryType(), AxisBoundaryType::Closed);
0365
0366     checkAllBins(
0367         dynamic_cast<GridPortalLinkT<decltype(axisPhi1Expected),
0368                                      Axis<Equidistant, Closed>>&>(*grid2d));
0369
0370     Axis axis2Explicit{AxisClosed, -180_degree, 180_degree, 3};
0371     auto grid2dExplicit = grid1->extendTo2d(&axis2Explicit);
0372     BOOST_REQUIRE(grid2dExplicit);
0373     BOOST_CHECK_EQUAL(grid2dExplicit->grid().axes().size(), 2);
0374     axis1 = grid2dExplicit->grid().axes().front();
0375     axis2 = grid2dExplicit->grid().axes().back();
0376     BOOST_CHECK_EQUAL(*axis1, axis1Expected);
0377     BOOST_CHECK_EQUAL(*axis2, axis2Explicit);
0378
0379     checkAllBins(dynamic_cast<GridPortalLinkT<decltype(axisPhi1Expected),
0380                                               decltype(axis2Explicit)>&>(
0381         *grid2dExplicit));
0382
0383     auto gridPhiBinnedInR = GridPortalLink::make(
0384         discPhi, AxisDirection::AxisR, Axis{AxisBound, 30_mm, 100_mm, 3});
0385     gridPhiBinnedInR->setVolume(vol.get());
0386     auto grid2dPhiNonClosed = gridPhiBinnedInR->extendTo2d(nullptr);
0387     BOOST_REQUIRE(grid2dPhiNonClosed);
0388     BOOST_CHECK_EQUAL(grid2dPhiNonClosed->grid().axes().size(), 2);
0389     Axis gridPhiBinnedInRExpected{AxisBound, 30_mm, 100_mm, 3};
0390     BOOST_CHECK_EQUAL(*grid2dPhiNonClosed->grid().axes().front(),
0391                       gridPhiBinnedInRExpected);
0392     const auto* axisPhiNonClosed = grid2dPhiNonClosed->grid().axes().back();
0393     BOOST_CHECK_CLOSE(axisPhiNonClosed->getMin(), -45_degree, 1e-6);
0394     BOOST_CHECK_CLOSE(axisPhiNonClosed->getMax(), 45_degree, 1e-6);
0395     BOOST_CHECK_EQUAL(axisPhiNonClosed->getNBins(), 1);
0396     BOOST_CHECK_EQUAL(axisPhiNonClosed->getType(), AxisType::Equidistant);
0397     BOOST_CHECK_EQUAL(axisPhiNonClosed->getBoundaryType(),
0398                       AxisBoundaryType::Bound);
0399
0400     checkAllBins(
0401         dynamic_cast<GridPortalLinkT<decltype(gridPhiBinnedInRExpected),
0402                                      Axis<Equidistant, Bound>>&>(
0403             *grid2dPhiNonClosed));
0404
0405     Axis axisPhiNonClosedExplicit{AxisBound, -45_degree, 45_degree, 3};
0406     auto grid2dPhiNonClosedExplicit =
0407         gridPhiBinnedInR->extendTo2d(&axisPhiNonClosedExplicit);
0408     BOOST_REQUIRE(grid2dPhiNonClosedExplicit);
0409     BOOST_CHECK_EQUAL(grid2dPhiNonClosedExplicit->grid().axes().size(), 2);
0410     axisPhiNonClosed = grid2dPhiNonClosedExplicit->grid().axes().back();
0411     BOOST_CHECK_EQUAL(*axisPhiNonClosed, axisPhiNonClosedExplicit);
0412     BOOST_CHECK_EQUAL(*grid2dPhiNonClosedExplicit->grid().axes().front(),
0413                       gridPhiBinnedInRExpected);
0414
0415     checkAllBins(
0416         dynamic_cast<GridPortalLinkT<decltype(gridPhiBinnedInRExpected),
0417                                      decltype(axisPhiNonClosedExplicit)>&>(
0418             *grid2dPhiNonClosedExplicit));
0419
0420     auto grid2dPhi = gridPhi->extendTo2d(nullptr);
0421     BOOST_REQUIRE(grid2dPhi);
0422     BOOST_CHECK_EQUAL(grid2dPhi->grid().axes().size(), 2);
0423     Axis axis2dPhiExpected{AxisBound, 30_mm, 100_mm, 1};
0424     BOOST_CHECK_EQUAL(*grid2dPhi->grid().axes().front(), axis2dPhiExpected);
0425     BOOST_CHECK_EQUAL(*grid2dPhi->grid().axes().back(), axisPhi1Expected);
0426
0427     checkAllBins(
0428         dynamic_cast<GridPortalLinkT<decltype(axis2dPhiExpected),
0429                                      decltype(axisPhi1Expected)>&>(*grid2dPhi));
0430
0431     Axis axis2dPhiExplicit{AxisBound, 30_mm, 100_mm, 3};
0432     auto grid2dPhiExplicit = gridPhi->extendTo2d(&axis2dPhiExplicit);
0433     BOOST_REQUIRE(grid2dPhiExplicit);
0434     BOOST_CHECK_EQUAL(grid2dPhiExplicit->grid().axes().size(), 2);
0435     BOOST_CHECK_EQUAL(*grid2dPhiExplicit->grid().axes().front(),
0436                       axis2dPhiExplicit);
0437     BOOST_CHECK_EQUAL(*grid2dPhiExplicit->grid().axes().back(),
0438                       axisPhi1Expected);
0439
0440     checkAllBins(dynamic_cast<GridPortalLinkT<decltype(axisPhi1Expected),
0441                                               decltype(axis2dPhiExplicit)>&>(
0442         *grid2dPhiExplicit));
0443   }
0444
0445   BOOST_TEST_CONTEXT("2D") {
0446     auto discPhi = Surface::makeShared<DiscSurface>(Transform3::Identity(),
0447                                                     30_mm, 100_mm, 45_degree);
0448
0449     Axis rBad{AxisBound, 1, 2, 5};
0450     Axis rGood{AxisBound, 30_mm, 100_mm, 5};
0451     Axis phiBad{AxisBound, 1, 2, 5};
0452     Axis phiGood{AxisBound, -45_degree, 45_degree, 5};
0453
0454     // r bad, phi bad
0455     BOOST_CHECK_THROW(GridPortalLink::make(discPhi, rBad, phiBad),
0456                       std::invalid_argument);
0457     // r bad, phi good
0458     BOOST_CHECK_THROW(GridPortalLink::make(discPhi, rBad, phiGood),
0459                       std::invalid_argument);
0460     // r good, phi bad
0461     BOOST_CHECK_THROW(GridPortalLink::make(discPhi, rGood, phiBad),
0462                       std::invalid_argument);
0463     // r good phi good
0464     auto grid = GridPortalLink::make(discPhi, rGood, phiGood);
0465     BOOST_REQUIRE_NE(grid, nullptr);
0466   }
0467 }
0468
0469 BOOST_AUTO_TEST_CASE(Plane) {
0470   using enum AxisType;
0471   using enum AxisBoundaryType;
0472   BOOST_TEST_CONTEXT("1D") {
0473     // Initialize surfaces
0474     auto rBounds = std::make_shared<const RectangleBounds>(30_mm, 100_mm);
0475     auto planeX =
0476         Surface::makeShared<PlaneSurface>(Transform3::Identity(), rBounds);
0477     auto planeY =
0478         Surface::makeShared<PlaneSurface>(Transform3::Identity(), rBounds);
0479
0480     // Volume for bin testing
0481     auto vol = std::make_shared<TrackingVolume>(
0482         Transform3::Identity(),
0483         std::make_shared<CuboidVolumeBounds>(30_mm, 40_mm, 100_mm));
0484
0485     // Initialize grids
0486     auto gridX = GridPortalLink::make(planeX, AxisDirection::AxisX,
0487                                       Axis{AxisBound, -30_mm, 30_mm, 3});
0488     auto gridY = GridPortalLink::make(planeY, AxisDirection::AxisY,
0489                                       Axis{AxisBound, -100_mm, 100_mm, 3});
0490
0491     // Check grid X
0492     BOOST_REQUIRE_NE(gridX, nullptr);
0493     BOOST_CHECK_EQUAL(gridX->grid().axes().size(), 1);
0494     const auto& axisX = *gridX->grid().axes().front();
0495     Axis axisXExpected{AxisBound, -30_mm, 30_mm, 3};
0496     BOOST_CHECK_EQUAL(axisX, axisXExpected);
0497
0498     // Check grid Y
0499     BOOST_REQUIRE_NE(gridY, nullptr);
0500     BOOST_CHECK_EQUAL(gridY->grid().axes().size(), 1);
0501     const auto& axisY = *gridY->grid().axes().front();
0502     Axis axisYExpected{AxisBound, -100_mm, 100_mm, 3};
0503     BOOST_CHECK_EQUAL(axisY, axisYExpected);
0504
0505     // Check gridX/gridX bin content
0506     auto checkAllBins = [&](const auto& grid) {
0507       visitBins(grid, [&](const TrackingVolume* content) {
0508         BOOST_CHECK_EQUAL(content, vol.get());
0509       });
0510     };
0511
0512     gridX->setVolume(vol.get());
0513     checkAllBins(*gridX);
0514
0515     gridY->setVolume(vol.get());
0516     checkAllBins(*gridY);
0517
0518     // Test making 2D grids from the 1D ones
0519
0520     // Test grid X
0521     auto gridX2d = gridX->extendTo2d(nullptr);
0522     BOOST_REQUIRE(gridX2d);
0523     BOOST_CHECK_EQUAL(gridX2d->grid().axes().size(), 2);
0524     const auto* axisX1 = gridX2d->grid().axes().front();
0525     BOOST_CHECK_EQUAL(*axisX1, axisXExpected);
0526     const auto* axisX2 = gridX2d->grid().axes().back();
0527     BOOST_CHECK_EQUAL(axisX2->getMin(), -100_mm);
0528     BOOST_CHECK_EQUAL(axisX2->getMax(), 100_mm);
0529     BOOST_CHECK_EQUAL(axisX2->getNBins(), 1);
0530     BOOST_CHECK_EQUAL(axisX2->getType(), AxisType::Equidistant);
0531     BOOST_CHECK_EQUAL(axisX2->getBoundaryType(), AxisBoundaryType::Bound);
0532
0533     // Test grid X explicit
0534     Axis axisYExplicit{AxisClosed, -100_mm, 100_mm, 3};
0535     auto gridX2dExplicit = gridX->extendTo2d(&axisYExplicit);
0536     BOOST_REQUIRE(gridX2dExplicit);
0537     BOOST_CHECK_EQUAL(gridX2dExplicit->grid().axes().size(), 2);
0538     axisX1 = gridX2dExplicit->grid().axes().front();
0539     axisX2 = gridX2dExplicit->grid().axes().back();
0540     BOOST_CHECK_EQUAL(*axisX1, axisXExpected);
0541     BOOST_CHECK_EQUAL(*axisX2, axisYExplicit);
0542
0543     checkAllBins(
0544         dynamic_cast<
0545             GridPortalLinkT<decltype(axisXExpected), decltype(axisYExplicit)>&>(
0546             *gridX2dExplicit));
0547
0548     // Test grid Y
0549     auto gridY2d = gridY->extendTo2d(nullptr);
0550     BOOST_REQUIRE(gridY2d);
0551     BOOST_CHECK_EQUAL(gridY2d->grid().axes().size(), 2);
0552     const auto* axisY1 = gridY2d->grid().axes().front();
0553     BOOST_CHECK_EQUAL(axisY1->getMin(), -30_mm);
0554     BOOST_CHECK_EQUAL(axisY1->getMax(), 30_mm);
0555     BOOST_CHECK_EQUAL(axisY1->getNBins(), 1);
0556     BOOST_CHECK_EQUAL(axisY1->getType(), AxisType::Equidistant);
0557     BOOST_CHECK_EQUAL(axisY1->getBoundaryType(), AxisBoundaryType::Bound);
0558     const auto* axisY2 = gridY2d->grid().axes().back();
0559     BOOST_CHECK_EQUAL(*axisY2, axisYExpected);
0560
0561     // Test grid Y explicit
0562     Axis axisXExplicit{AxisClosed, -30_mm, 30_mm, 3};
0563     auto gridY2dExplicit = gridY->extendTo2d(&axisXExplicit);
0564     BOOST_REQUIRE(gridY2dExplicit);
0565     BOOST_CHECK_EQUAL(gridY2dExplicit->grid().axes().size(), 2);
0566     axisY1 = gridY2dExplicit->grid().axes().front();
0567     axisY2 = gridY2dExplicit->grid().axes().back();
0568     BOOST_CHECK_EQUAL(*axisY1, axisXExplicit);
0569     BOOST_CHECK_EQUAL(*axisY2, axisYExpected);
0570
0571     checkAllBins(
0572         dynamic_cast<
0573             GridPortalLinkT<decltype(axisXExplicit), decltype(axisYExpected)>&>(
0574             *gridY2dExplicit));
0575   }
0576   BOOST_TEST_CONTEXT("2D") {
0577     auto rBounds = std::make_shared<const RectangleBounds>(30_mm, 100_mm);
0578
0579     auto plane =
0580         Surface::makeShared<PlaneSurface>(Transform3::Identity(), rBounds);
0581
0582     Axis xAxis{AxisBound, -30_mm, 30_mm, 5};
0583     Axis yAxis{AxisBound, -100_mm, 100_mm, 5};
0584
0585     auto grid = GridPortalLink::make(plane, xAxis, yAxis);
0586     BOOST_REQUIRE_NE(grid, nullptr);
0587   }
0588 }
0589
0590 BOOST_AUTO_TEST_CASE(FromTrivial) {
0591   BOOST_TEST_CONTEXT("Cylinder") {
0592     auto cyl = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
0593                                                     30_mm, 100_mm);
0594
0595     auto vol = std::make_shared<TrackingVolume>(
0596         Transform3::Identity(),
0597         std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
0598
0599     auto trivial = std::make_unique<TrivialPortalLink>(cyl, *vol);
0600     BOOST_REQUIRE(trivial);
0601
0602     BOOST_CHECK_EQUAL(trivial->resolveVolume(gctx, Vector2{1, 2}).value(),
0603                       vol.get());
0604
0605     auto gridZ = trivial->makeGrid(AxisDirection::AxisZ);
0606     BOOST_REQUIRE(gridZ);
0607
0608     BOOST_CHECK_EQUAL(gridZ->grid().axes().size(), 1);
0609     BOOST_CHECK_EQUAL(gridZ->surface().bounds(), cyl->bounds());
0610     Axis axisZExpected{AxisBound, -100_mm, 100_mm, 1};
0611     BOOST_CHECK_EQUAL(*gridZ->grid().axes().front(), axisZExpected);
0612
0613     BOOST_CHECK_EQUAL(
0614         gridZ->resolveVolume(gctx, Vector2{20_degree * 30_mm, 90_mm}).value(),
0615         vol.get());
0616
0617     auto gridRPhi = trivial->makeGrid(AxisDirection::AxisRPhi);
0618     BOOST_REQUIRE(gridRPhi);
0619
0620     BOOST_CHECK_EQUAL(gridRPhi->grid().axes().size(), 1);
0621     BOOST_CHECK_EQUAL(gridRPhi->surface().bounds(), cyl->bounds());
0622     Axis axisRPhiExpected{AxisClosed, -std::numbers::pi * 30_mm,
0623                           std::numbers::pi * 30_mm, 1};
0624     BOOST_CHECK_EQUAL(*gridRPhi->grid().axes().front(), axisRPhiExpected);
0625
0626     auto cylPhi = Surface::makeShared<CylinderSurface>(
0627         Transform3::Identity(), 30_mm, 100_mm, 30_degree);
0628
0629     auto trivialPhi = std::make_unique<TrivialPortalLink>(cylPhi, *vol);
0630     BOOST_REQUIRE(trivialPhi);
0631
0632     BOOST_CHECK_EQUAL(trivialPhi->resolveVolume(gctx, Vector2{1, 2}).value(),
0633                       vol.get());
0634
0635     auto gridRPhiSector = trivialPhi->makeGrid(AxisDirection::AxisRPhi);
0636     BOOST_REQUIRE(gridRPhiSector);
0637
0638     BOOST_CHECK_EQUAL(
0639         gridRPhiSector->resolveVolume(gctx, Vector2{20_degree * 30_mm, 90_mm})
0640             .value(),
0641         vol.get());
0642
0643     BOOST_CHECK_EQUAL(gridRPhiSector->grid().axes().size(), 1);
0644     Axis axisRPhiSectorExpected{AxisBound, -30_degree * 30_mm,
0645                                 30_degree * 30_mm, 1};
0646     BOOST_CHECK_EQUAL(*gridRPhiSector->grid().axes().front(),
0647                       axisRPhiSectorExpected);
0648   }
0649
0650   BOOST_TEST_CONTEXT("Disc") {
0651     auto disc =
0652         Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 100_mm);
0653
0654     auto vol = std::make_shared<TrackingVolume>(
0655         Transform3::Identity(),
0656         std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
0657
0658     auto trivial = std::make_unique<TrivialPortalLink>(disc, *vol);
0659     BOOST_REQUIRE(trivial);
0660
0661     // Doesn't matter which position
0662     BOOST_CHECK_EQUAL(trivial->resolveVolume(gctx, Vector2{1, 2}).value(),
0663                       vol.get());
0664
0665     auto gridR = trivial->makeGrid(AxisDirection::AxisR);
0666     BOOST_REQUIRE(gridR);
0667
0668     BOOST_CHECK_EQUAL(gridR->grid().axes().size(), 1);
0669     BOOST_CHECK_EQUAL(gridR->surface().bounds(), disc->bounds());
0670     Axis axisRExpected{AxisBound, 30_mm, 100_mm, 1};
0671     BOOST_CHECK_EQUAL(*gridR->grid().axes().front(), axisRExpected);
0672
0673     BOOST_CHECK_EQUAL(
0674         gridR->resolveVolume(gctx, Vector2{90_mm, 10_degree}).value(),
0675         vol.get());
0676
0677     auto gridPhi = trivial->makeGrid(AxisDirection::AxisPhi);
0678     BOOST_REQUIRE(gridPhi);
0679
0680     BOOST_CHECK_EQUAL(gridPhi->grid().axes().size(), 1);
0681     BOOST_CHECK_EQUAL(gridPhi->surface().bounds(), disc->bounds());
0682     Axis axisPhiExpected{AxisClosed, -std::numbers::pi, std::numbers::pi, 1};
0683     BOOST_CHECK_EQUAL(*gridPhi->grid().axes().front(), axisPhiExpected);
0684
0685     BOOST_CHECK_EQUAL(
0686         gridPhi->resolveVolume(gctx, Vector2{90_mm, 10_degree}).value(),
0687         vol.get());
0688   }
0689   BOOST_TEST_CONTEXT("Plane") {
0690     auto rBounds = std::make_shared<const RectangleBounds>(30_mm, 100_mm);
0691
0692     auto plane =
0693         Surface::makeShared<PlaneSurface>(Transform3::Identity(), rBounds);
0694
0695     auto vol = std::make_shared<TrackingVolume>(
0696         Transform3::Identity(),
0697         std::make_shared<CuboidVolumeBounds>(30_mm, 40_mm, 100_mm));
0698
0699     auto trivial = std::make_unique<TrivialPortalLink>(plane, *vol);
0700     BOOST_REQUIRE(trivial);
0701
0702     // Doesn't matter which position
0703     BOOST_CHECK_EQUAL(
0704         trivial->resolveVolume(gctx, Vector2{10_mm, 20_mm}).value(), vol.get());
0705
0706     auto gridX = trivial->makeGrid(AxisDirection::AxisX);
0707     BOOST_REQUIRE(gridX);
0708
0709     BOOST_CHECK_EQUAL(gridX->grid().axes().size(), 1);
0710     BOOST_CHECK_EQUAL(gridX->surface().bounds(), plane->bounds());
0711     Axis axisXExpected{AxisBound, -30_mm, 30_mm, 1};
0712     BOOST_CHECK_EQUAL(*gridX->grid().axes().front(), axisXExpected);
0713
0714     BOOST_CHECK_EQUAL(gridX->resolveVolume(gctx, Vector2{20_mm, 10_mm}).value(),
0715                       vol.get());
0716
0717     auto gridY = trivial->makeGrid(AxisDirection::AxisY);
0718     BOOST_REQUIRE(gridY);
0719
0720     BOOST_CHECK_EQUAL(gridY->grid().axes().size(), 1);
0721     BOOST_CHECK_EQUAL(gridY->surface().bounds(), plane->bounds());
0722     Axis axisYExpected{AxisBound, -100_mm, 100_mm, 1};
0723     BOOST_CHECK_EQUAL(*gridY->grid().axes().front(), axisYExpected);
0724
0725     BOOST_CHECK_EQUAL(gridY->resolveVolume(gctx, Vector2{15_mm, 20_mm}).value(),
0726                       vol.get());
0727   }
0728 }
0729
0730 BOOST_AUTO_TEST_SUITE_END()  // GridConstruction
0731
0732 BOOST_AUTO_TEST_SUITE(GridMerging)
0733
0734 BOOST_AUTO_TEST_SUITE(Merging1dCylinder)
0735
0736 BOOST_AUTO_TEST_SUITE(ZDirection)
0737
0738 BOOST_AUTO_TEST_CASE(ColinearMerge) {
0739   auto vol1 = makeDummyVolume();
0740   auto vol2 = makeDummyVolume();
0741
0742   auto cyl = Surface::makeShared<CylinderSurface>(Transform3::Identity(), 30_mm,
0743                                                   100_mm);
0744
0745   auto grid1dCyl = GridPortalLink::make(cyl, AxisDirection::AxisZ,
0746                                         Axis{AxisBound, -100_mm, 100_mm, 10});
0747   grid1dCyl->setVolume(vol1.get());
0748
0749   // Another cylinder, shifted in z
0750   auto cyl2 = Surface::makeShared<CylinderSurface>(
0751       Transform3{Translation3{Vector3::UnitZ() * 150_mm}}, 30_mm, 50_mm);
0752
0753   auto grid1dCyl2 = GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0754                                          Axis{AxisBound, -50_mm, 50_mm, 5});
0755
0756   grid1dCyl2->setVolume(vol2.get());
0757
0758   // Completely invalid
0759   BOOST_CHECK_THROW(GridPortalLink::merge(*grid1dCyl, *grid1dCyl2,
0760                                           AxisDirection::AxisPhi, *logger),
0761                     AssertionFailureException);
0762   // Invalid direction, as the cylinders are shifted in z, and can't be merged
0763   // in r x phi
0764   BOOST_CHECK_THROW(GridPortalLink::merge(*grid1dCyl, *grid1dCyl2,
0765                                           AxisDirection::AxisRPhi, *logger),
0766                     SurfaceMergingException);
0767
0768   BOOST_TEST_CONTEXT("Consistent equidistant") {
0769     auto mergedPtr = GridPortalLink::merge(*grid1dCyl, *grid1dCyl2,
0770                                            AxisDirection::AxisZ, *logger);
0771     BOOST_REQUIRE_NE(mergedPtr, nullptr);
0772     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0773     BOOST_REQUIRE_NE(merged, nullptr);
0774     BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0775     const auto& axis = *merged->grid().axes().front();
0776     BOOST_CHECK_EQUAL(axis.getMin(), -150_mm);
0777     BOOST_CHECK_EQUAL(axis.getMax(), 150_mm);
0778     BOOST_CHECK_EQUAL(axis.getNBins(), 15);
0779     BOOST_CHECK_EQUAL(axis.getType(), AxisType::Equidistant);
0780     BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0781   }
0782
0783   BOOST_TEST_CONTEXT("Inconsistent equidistant") {
0784     auto grid1dCyl2BinWidthChanged = GridPortalLink::make(
0785         cyl2, AxisDirection::AxisZ, Axis{AxisBound, -50_mm, 50_mm, 6});
0786
0787     auto mergedPtr = GridPortalLink::merge(
0788         *grid1dCyl, *grid1dCyl2BinWidthChanged, AxisDirection::AxisZ, *logger);
0789     BOOST_REQUIRE_NE(mergedPtr, nullptr);
0790     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0791     BOOST_REQUIRE_NE(merged, nullptr);
0792     BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0793     const auto& axis = *merged->grid().axes().front();
0794     BOOST_CHECK_EQUAL(axis.getMin(), -150_mm);
0795     BOOST_CHECK_EQUAL(axis.getMax(), 150_mm);
0796     BOOST_CHECK_EQUAL(axis.getNBins(), 16);
0797     BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
0798     BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0799   }
0800
0801   BOOST_TEST_CONTEXT("Right Variable") {
0802     auto gridLeft = GridPortalLink::make(cyl, AxisDirection::AxisZ,
0803                                          Axis{AxisBound, -100_mm, 100_mm, 10});
0804
0805     auto gridRight =
0806         GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0807                              Axis{AxisBound, {-50_mm, -10_mm, 10_mm, 50_mm}});
0808
0809     auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight,
0810                                            AxisDirection::AxisZ, *logger);
0811     BOOST_REQUIRE_NE(mergedPtr, nullptr);
0812     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0813     BOOST_REQUIRE_NE(merged, nullptr);
0814     BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0815     const auto& axis = *merged->grid().axes().front();
0816     BOOST_CHECK_EQUAL(axis.getMin(), -150_mm);
0817     BOOST_CHECK_EQUAL(axis.getMax(), 150_mm);
0818     BOOST_CHECK_EQUAL(axis.getNBins(), 13);
0819     BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
0820     BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0821   }
0822
0823   BOOST_TEST_CONTEXT("Left Variable") {
0824     auto gridLeft =
0825         GridPortalLink::make(cyl, AxisDirection::AxisZ,
0826                              Axis{AxisBound, {-100_mm, -80_mm, 10_mm, 100_mm}});
0827
0828     auto gridRight = GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0829                                           Axis{AxisBound, -50_mm, 50_mm, 8});
0830
0831     auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight,
0832                                            AxisDirection::AxisZ, *logger);
0833     BOOST_REQUIRE_NE(mergedPtr, nullptr);
0834     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0835     BOOST_REQUIRE_NE(merged, nullptr);
0836     BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0837     const auto& axis = *merged->grid().axes().front();
0838     BOOST_CHECK_EQUAL(axis.getMin(), -150_mm);
0839     BOOST_CHECK_EQUAL(axis.getMax(), 150_mm);
0840     BOOST_CHECK_EQUAL(axis.getNBins(), 11);
0841     BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
0842     BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0843   }
0844
0845   BOOST_TEST_CONTEXT("Both Variable") {
0846     auto gridLeft =
0847         GridPortalLink::make(cyl, AxisDirection::AxisZ,
0848                              Axis{AxisBound, {-100_mm, -80_mm, 10_mm, 100_mm}});
0849
0850     auto gridRight =
0851         GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0852                              Axis{AxisBound, {-50_mm, -10_mm, 10_mm, 50_mm}});
0853
0854     auto mergedPtr = GridPortalLink::merge(*gridLeft, *gridRight,
0855                                            AxisDirection::AxisZ, *logger);
0856     BOOST_REQUIRE_NE(mergedPtr, nullptr);
0857     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
0858     BOOST_REQUIRE_NE(merged, nullptr);
0859     BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0860     const auto& axis = *merged->grid().axes().front();
0861     BOOST_CHECK_EQUAL(axis.getMin(), -150_mm);
0862     BOOST_CHECK_EQUAL(axis.getMax(), 150_mm);
0863     BOOST_CHECK_EQUAL(axis.getNBins(), 6);
0864     BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
0865     BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0866   }
0867
0868   BOOST_TEST_CONTEXT("Non bound axis") {
0869     std::unique_ptr gridLeft =
0870         GridPortalLink::make(cyl, AxisDirection::AxisZ,
0871                              Axis{AxisBound, {-100_mm, -80_mm, 10_mm, 100_mm}});
0872     std::unique_ptr gridRightClosed =
0873         GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0874                              Axis{AxisClosed, {-50_mm, -10_mm, 10_mm, 50_mm}});
0875     std::unique_ptr gridRightOpen =
0876         GridPortalLink::make(cyl2, AxisDirection::AxisZ,
0877                              Axis{AxisOpen, {-50_mm, -10_mm, 10_mm, 50_mm}});
0878
0879     auto compositeLR = PortalLinkBase::merge(
0880         copy(gridLeft), copy(gridRightClosed), AxisDirection::AxisZ, *logger);
0881     BOOST_CHECK_NE(dynamic_cast<const CompositePortalLink*>(compositeLR.get()),
0882                    nullptr);
0883     auto compositeRL = PortalLinkBase::merge(
0884         copy(gridLeft), copy(gridRightOpen), AxisDirection::AxisZ, *logger);
0885     BOOST_CHECK_NE(dynamic_cast<const CompositePortalLink*>(compositeRL.get()),
0886                    nullptr);
0887   }
0888 }
0889
0890 BOOST_AUTO_TEST_CASE(ParallelMerge) {
0891   // Merge in z direction with phi sectors
0892   auto cyl1 = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
0893                                                    30_mm, 100_mm, 35_degree);
0894
0895   auto grid1 = GridPortalLink::make(
0896       cyl1, AxisDirection::AxisRPhi,
0897
0898       Axis{AxisBound, -35_degree * 30_mm, 35_degree * 30_mm, 3});
0899   auto cyl2 = Surface::makeShared<CylinderSurface>(
0900       Transform3{Translation3{Vector3::UnitZ() * 150_mm}}, 30_mm, 50_mm,
0901       35_degree);
0902
0903   auto grid2 = GridPortalLink::make(
0904       cyl2, AxisDirection::AxisRPhi,
0905       Axis{AxisBound, -35_degree * 30_mm, 35_degree * 30_mm, 3});
0906
0907   auto merged12Ptr =
0908       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisZ, *logger);
0909   BOOST_REQUIRE_NE(merged12Ptr, nullptr);
0910   auto merged12 = dynamic_cast<const GridPortalLink*>(merged12Ptr.get());
0911   BOOST_REQUIRE_NE(merged12, nullptr);
0912
0913   BOOST_REQUIRE_EQUAL(merged12->grid().axes().size(), 2);
0914
0915   const auto& axis1 = *merged12->grid().axes().front();
0916   const auto& axis2 = *merged12->grid().axes().back();
0917   Axis axis1Expected{AxisBound, -35_degree * 30_mm, 35_degree * 30_mm, 3};
0918   BOOST_CHECK_EQUAL(axis1, axis1Expected);
0919   Axis axis2Expected{AxisBound, {-150_mm, 50_mm, 150_mm}};
0920   BOOST_CHECK_EQUAL(axis2, axis2Expected);
0921 }
0922
0923 BOOST_AUTO_TEST_SUITE_END()  // ZDirection
0924
0925 BOOST_AUTO_TEST_SUITE(RPhiDirection)
0926
0927 BOOST_AUTO_TEST_CASE(ColinearMerge) {
0928   auto cyl = Surface::makeShared<CylinderSurface>(Transform3::Identity(), 30_mm,
0929                                                   100_mm);
0930   BOOST_CHECK_THROW(GridPortalLink::make(cyl, AxisDirection::AxisRPhi,
0931                                          Axis{AxisBound, 0, 5, 5}),
0932                     std::invalid_argument);
0933
0934   auto cylNonZeroAverage = Surface::makeShared<CylinderSurface>(
0935       Transform3::Identity(), 30_mm, 100_mm, 20_degree, 45_degree);
0936
0937   BOOST_CHECK_THROW(
0938       GridPortalLink::make(
0939           cylNonZeroAverage, AxisDirection::AxisRPhi,
0940           Axis{AxisBound, -20_degree * 30_mm, 20_degree * 30_mm, 5}),
0941       std::invalid_argument);
0942
0943   BOOST_TEST_CONTEXT("Colinear merge in rPhi") {
0944     auto cylPhi1 = Surface::makeShared<CylinderSurface>(
0945         Transform3::Identity() * AngleAxis3(45_degree, Vector3::UnitZ()), 30_mm,
0946         100_mm, 20_degree, 0_degree);
0947
0948     auto cylPhi2 = Surface::makeShared<CylinderSurface>(
0949         Transform3::Identity() * AngleAxis3(105_degree, Vector3::UnitZ()),
0950         30_mm, 100_mm, 40_degree, 0_degree);
0951
0952     auto portalPhi1 = GridPortalLink::make(
0953         cylPhi1, AxisDirection::AxisRPhi,
0954         Axis{AxisBound, -20_degree * 30_mm, 20_degree * 30_mm, 5});
0955
0956     auto portalPhi2 = GridPortalLink::make(
0957         cylPhi2, AxisDirection::AxisRPhi,
0958         Axis{AxisBound, -40_degree * 30_mm, 40_degree * 30_mm, 10});
0959
0960     auto cylPhi3 = Surface::makeShared<CylinderSurface>(
0961         Transform3::Identity() * AngleAxis3(45_degree, Vector3::UnitZ()), 30_mm,
0962         100_mm, 90_degree, 0_degree);
0963
0964     auto cylPhi4 = Surface::makeShared<CylinderSurface>(
0965         Transform3::Identity() * AngleAxis3(-135_degree, Vector3::UnitZ()),
0966         30_mm, 100_mm, 90_degree, 0_degree);
0967
0968     auto portalPhi3 = GridPortalLink::make(
0969         cylPhi3, AxisDirection::AxisRPhi,
0970         Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 2});
0971
0972     auto portalPhi4 = GridPortalLink::make(
0973         cylPhi4, AxisDirection::AxisRPhi,
0974         Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 2});
0975
0976     BOOST_TEST_CONTEXT("Consistent equidistant") {
0977       auto portalMerged = GridPortalLink::merge(
0978           *portalPhi1, *portalPhi2, AxisDirection::AxisRPhi, *logger);
0979       BOOST_REQUIRE_NE(portalMerged, nullptr);
0980
0981       const auto* merged =
0982           dynamic_cast<const GridPortalLink*>(portalMerged.get());
0983       BOOST_REQUIRE_NE(merged, nullptr);
0984       BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
0985       const auto& axis = *merged->grid().axes().front();
0986       BOOST_CHECK_CLOSE(axis.getMin(), -60_degree * 30_mm, 1e-9);
0987       BOOST_CHECK_CLOSE(axis.getMax(), 60_degree * 30_mm, 1e-9);
0988       BOOST_CHECK_EQUAL(axis.getNBins(), 15);
0989       BOOST_CHECK_EQUAL(axis.getType(), AxisType::Equidistant);
0990       BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
0991
0992       // Test that if you merge half-circles, we get a closed axis
0993       auto portalMerged34 = GridPortalLink::merge(
0994           *portalPhi3, *portalPhi4, AxisDirection::AxisRPhi, *logger);
0995       BOOST_REQUIRE_NE(portalMerged34, nullptr);
0996
0997       const auto* merged34 =
0998           dynamic_cast<const GridPortalLink*>(portalMerged34.get());
0999       BOOST_REQUIRE_NE(merged34, nullptr);
1000       BOOST_CHECK_EQUAL(merged34->grid().axes().size(), 1);
1001       const auto& axis34 = *merged34->grid().axes().front();
1002       BOOST_CHECK_CLOSE(axis34.getMin(), -180_degree * 30_mm, 1e-9);
1003       BOOST_CHECK_CLOSE(axis34.getMax(), 180_degree * 30_mm, 1e-9);
1004       BOOST_CHECK_EQUAL(axis34.getNBins(), 4);
1005       BOOST_CHECK_EQUAL(axis34.getType(), AxisType::Equidistant);
1006       BOOST_CHECK_EQUAL(axis34.getBoundaryType(), AxisBoundaryType::Closed);
1007     }
1008
1009     BOOST_TEST_CONTEXT("Inconsistent equidistant") {
1010       auto portalPhi2Mod = GridPortalLink::make(
1011           cylPhi2, AxisDirection::AxisRPhi,
1012           Axis{AxisBound, -40_degree * 30_mm, 40_degree * 30_mm, 3});
1013
1014       auto portalMergedMod = GridPortalLink::merge(
1015           *portalPhi1, *portalPhi2Mod, AxisDirection::AxisRPhi, *logger);
1016       BOOST_REQUIRE_NE(portalMergedMod, nullptr);
1017
1018       const auto* merged12 =
1019           dynamic_cast<const GridPortalLink*>(portalMergedMod.get());
1020       BOOST_REQUIRE_NE(merged12, nullptr);
1021       BOOST_CHECK_EQUAL(merged12->grid().axes().size(), 1);
1022       const auto& axis12 = *merged12->grid().axes().front();
1023       BOOST_CHECK_CLOSE(axis12.getMin(), -60_degree * 30_mm, 1e-9);
1024       BOOST_CHECK_CLOSE(axis12.getMax(), 60_degree * 30_mm, 1e-9);
1025       BOOST_CHECK_EQUAL(axis12.getNBins(), 8);
1026       BOOST_CHECK_EQUAL(axis12.getType(), AxisType::Variable);
1027       BOOST_CHECK_EQUAL(axis12.getBoundaryType(), AxisBoundaryType::Bound);
1028
1029       std::vector<double> expected12 = {-31.4159, -17.4533, -3.49066,
1030                                         10.472,   14.6608,  18.8496,
1031                                         23.0383,  27.2271,  31.4159};
1032       CHECK_CLOSE_OR_SMALL(axis12.getBinEdges(), expected12, 1e-4, 10e-10);
1033
1034       auto portalPhi4Mod = GridPortalLink::make(
1035           cylPhi4, AxisDirection::AxisRPhi,
1036           Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 1});
1037
1038       auto portalMerged34 = GridPortalLink::merge(
1039           *portalPhi3, *portalPhi4Mod, AxisDirection::AxisRPhi, *logger);
1040       BOOST_REQUIRE_NE(portalMerged34, nullptr);
1041
1042       const auto* merged34 =
1043           dynamic_cast<const GridPortalLink*>(portalMerged34.get());
1044       BOOST_REQUIRE_NE(merged34, nullptr);
1045       BOOST_CHECK_EQUAL(merged34->grid().axes().size(), 1);
1046       const auto& axis34 = *merged34->grid().axes().front();
1047       BOOST_CHECK_CLOSE(axis34.getMin(), -180_degree * 30_mm, 1e-9);
1048       BOOST_CHECK_CLOSE(axis34.getMax(), 180_degree * 30_mm, 1e-9);
1049       BOOST_CHECK_EQUAL(axis34.getNBins(), 3);
1050       BOOST_CHECK_EQUAL(axis34.getType(), AxisType::Variable);
1051       BOOST_CHECK_EQUAL(axis34.getBoundaryType(), AxisBoundaryType::Closed);
1052
1053       // Caution: for full-azimuth cases, the ordering is preserved, you get
1054       // in what you get out. -> this can flip
1055       std::vector<double> expected34 = {-94.2478, -47.1239, 0, 94.2478};
1056       CHECK_CLOSE_OR_SMALL(axis34.getBinEdges(), expected34, 1e-4, 10e-10);
1057     }
1058
1059     BOOST_TEST_CONTEXT("Left variable") {
1060       BOOST_TEST_CONTEXT("Non-closed") {
1061         auto gridLeft =
1062             GridPortalLink::make(cylPhi1, AxisDirection::AxisRPhi,
1063                                  Axis{AxisBound,
1064                                       {-20_degree * 30_mm, -10_degree * 30_mm,
1065                                        10_degree * 30_mm, 20_degree * 30_mm}});
1066         auto gridRight = GridPortalLink::make(
1067             cylPhi2, AxisDirection::AxisRPhi,
1068             Axis{AxisBound, -40_degree * 30_mm, 40_degree * 30_mm, 3});
1069
1070         auto mergedPtr = GridPortalLink::merge(
1071             *gridLeft, *gridRight, AxisDirection::AxisRPhi, *logger);
1072         BOOST_REQUIRE_NE(mergedPtr, nullptr);
1073
1074         const auto* merged =
1075             dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1076         BOOST_REQUIRE_NE(merged, nullptr);
1077         const auto& axis = *merged->grid().axes().front();
1078         BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1079         BOOST_CHECK_CLOSE(axis.getMin(), -60_degree * 30_mm, 1e-9);
1080         BOOST_CHECK_CLOSE(axis.getMax(), 60_degree * 30_mm, 1e-9);
1081         BOOST_CHECK_EQUAL(axis.getNBins(), 6);
1082         BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
1083         BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
1084
1085         std::vector<double> expected = {-31.4159, -17.4533, -3.49066, 10.472,
1086                                         15.708,   26.1799,  31.4159};
1087         CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10);
1088       }
1089
1090       BOOST_TEST_CONTEXT("Closed") {
1091         auto gridLeft = GridPortalLink::make(
1092             cylPhi4, AxisDirection::AxisRPhi,
1093             Axis{AxisBound,
1094                  {-90_degree * 30_mm, 25_degree * 30_mm, 90_degree * 30_mm}});
1095
1096         auto gridRight = GridPortalLink::make(
1097             cylPhi3, AxisDirection::AxisRPhi,
1098             Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 3});
1099
1100         auto mergedPtr = GridPortalLink::merge(
1101             *gridLeft, *gridRight, AxisDirection::AxisRPhi, *logger);
1102         BOOST_REQUIRE_NE(mergedPtr, nullptr);
1103
1104         const auto* merged =
1105             dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1106         BOOST_REQUIRE_NE(merged, nullptr);
1107         BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1108         const auto& axis = *merged->grid().axes().front();
1109         BOOST_CHECK_CLOSE(axis.getMin(), -180_degree * 30_mm, 1e-9);
1110         BOOST_CHECK_CLOSE(axis.getMax(), 180_degree * 30_mm, 1e-9);
1111         BOOST_CHECK_EQUAL(axis.getNBins(), 5);
1112         BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
1113         BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Closed);
1114
1115         // Caution: for full-azimuth cases, the ordering is preserved, you get
1116         // in what you get out. -> this can flip
1117         std::vector<double> expected = {-94.2478, -34.0339, 0,
1118                                         31.4159,  62.8319,  94.2478};
1119         CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10);
1120       }
1121     }
1122
1123     BOOST_TEST_CONTEXT("Right variable") {
1124       BOOST_TEST_CONTEXT("Non-closed") {
1125         auto gridLeft = GridPortalLink::make(
1126             cylPhi1, AxisDirection::AxisRPhi,
1127             Axis{AxisBound, -20_degree * 30_mm, 20_degree * 30_mm, 3});
1128         auto gridRight =
1129             GridPortalLink::make(cylPhi2, AxisDirection::AxisRPhi,
1130                                  Axis{AxisBound,
1131                                       {-40_degree * 30_mm, -10_degree * 30_mm,
1132                                        10_degree * 30_mm, 40_degree * 30_mm}});
1133
1134         auto mergedPtr = GridPortalLink::merge(
1135             *gridLeft, *gridRight, AxisDirection::AxisRPhi, *logger);
1136         BOOST_REQUIRE_NE(mergedPtr, nullptr);
1137
1138         const auto* merged =
1139             dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1140         BOOST_REQUIRE_NE(merged, nullptr);
1141         const auto& axis = *merged->grid().axes().front();
1142         BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1143         BOOST_CHECK_CLOSE(axis.getMin(), -60_degree * 30_mm, 1e-9);
1144         BOOST_CHECK_CLOSE(axis.getMax(), 60_degree * 30_mm, 1e-9);
1145         BOOST_CHECK_EQUAL(axis.getNBins(), 6);
1146         BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
1147         BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
1148
1149         std::vector<double> expected = {-31.4159, -15.708, -5.23599, 10.472,
1150                                         17.4533,  24.4346, 31.4159};
1151         CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10);
1152       }
1153
1154       BOOST_TEST_CONTEXT("Closed") {
1155         auto gridLeft = GridPortalLink::make(
1156             cylPhi4, AxisDirection::AxisRPhi,
1157             Axis{AxisBound, -90_degree * 30_mm, 90_degree * 30_mm, 3});
1158
1159         auto gridRight = GridPortalLink::make(
1160             cylPhi3, AxisDirection::AxisRPhi,
1161             Axis{AxisBound,
1162                  {-90_degree * 30_mm, 25_degree * 30_mm, 90_degree * 30_mm}});
1163
1164         auto mergedPtr = GridPortalLink::merge(
1165             *gridLeft, *gridRight, AxisDirection::AxisRPhi, *logger);
1166         BOOST_REQUIRE_NE(mergedPtr, nullptr);
1167
1168         const auto* merged =
1169             dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1170         BOOST_REQUIRE_NE(merged, nullptr);
1171         BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1172         const auto& axis = *merged->grid().axes().front();
1173         BOOST_CHECK_CLOSE(axis.getMin(), -180_degree * 30_mm, 1e-9);
1174         BOOST_CHECK_CLOSE(axis.getMax(), 180_degree * 30_mm, 1e-9);
1175         BOOST_CHECK_EQUAL(axis.getNBins(), 5);
1176         BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
1177         BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Closed);
1178
1179         // Caution: for full-azimuth cases, the ordering is preserved, you get
1180         // in what you get out. -> this can flip
1181         std::vector<double> expected = {-94.2478, -62.8319, -31.4159,
1182                                         0,        60.2139,  94.2478};
1183         CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10);
1184       }
1185     }
1186
1187     BOOST_TEST_CONTEXT("Both variable") {
1188       BOOST_TEST_CONTEXT("Non-closed") {
1189         auto gridLeft =
1190             GridPortalLink::make(cylPhi1, AxisDirection::AxisRPhi,
1191                                  Axis{AxisBound,
1192                                       {-20_degree * 30_mm, -10_degree * 30_mm,
1193                                        10_degree * 30_mm, 20_degree * 30_mm}});
1194         auto gridRight = GridPortalLink::make(
1195             cylPhi2, AxisDirection::AxisRPhi,
1196             Axis{AxisBound,
1197                  {-40_degree * 30_mm, -5_degree * 30_mm, 40_degree * 30_mm}});
1198
1199         auto mergedPtr = GridPortalLink::merge(
1200             *gridLeft, *gridRight, AxisDirection::AxisRPhi, *logger);
1201         BOOST_REQUIRE_NE(mergedPtr, nullptr);
1202
1203         const auto* merged =
1204             dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1205         BOOST_REQUIRE_NE(merged, nullptr);
1206         const auto& axis = *merged->grid().axes().front();
1207         BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1208         BOOST_CHECK_CLOSE(axis.getMin(), -60_degree * 30_mm, 1e-9);
1209         BOOST_CHECK_CLOSE(axis.getMax(), 60_degree * 30_mm, 1e-9);
1210         BOOST_CHECK_EQUAL(axis.getNBins(), 5);
1211         BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
1212         BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Bound);
1213
1214         std::vector<double> expected = {-31.4159, -13.09,  10.472,
1215                                         15.708,   26.1799, 31.4159};
1216         CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10);
1217       }
1218
1219       BOOST_TEST_CONTEXT("Closed") {
1220         auto gridLeft = GridPortalLink::make(
1221             cylPhi4, AxisDirection::AxisRPhi,
1222             Axis{AxisBound,
1223                  {-90_degree * 30_mm, 25_degree * 30_mm, 90_degree * 30_mm}});
1224
1225         auto gridRight =
1226             GridPortalLink::make(cylPhi3, AxisDirection::AxisRPhi,
1227                                  Axis{AxisBound,
1228                                       {-90_degree * 30_mm, -10_degree * 30_mm,
1229                                        10_degree * 30_mm, 90_degree * 30_mm}});
1230
1231         auto mergedPtr = GridPortalLink::merge(
1232             *gridLeft, *gridRight, AxisDirection::AxisRPhi, *logger);
1233         BOOST_REQUIRE_NE(mergedPtr, nullptr);
1234
1235         const auto* merged =
1236             dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1237         BOOST_REQUIRE_NE(merged, nullptr);
1238         BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1239         const auto& axis = *merged->grid().axes().front();
1240         BOOST_CHECK_CLOSE(axis.getMin(), -180_degree * 30_mm, 1e-9);
1241         BOOST_CHECK_CLOSE(axis.getMax(), 180_degree * 30_mm, 1e-9);
1242         BOOST_CHECK_EQUAL(axis.getNBins(), 5);
1243         BOOST_CHECK_EQUAL(axis.getType(), AxisType::Variable);
1244         BOOST_CHECK_EQUAL(axis.getBoundaryType(), AxisBoundaryType::Closed);
1245
1246         // Caution: for full-azimuth cases, the ordering is preserved, you get
1247         // in what you get out. -> this can flip
1248         std::vector<double> expected = {-94.2478, -34.0339, 0,
1249                                         41.8879,  52.3599,  94.2478};
1250         CHECK_CLOSE_OR_SMALL(axis.getBinEdges(), expected, 1e-4, 10e-10);
1251       }
1252     }
1253   }
1254 }
1255
1256 BOOST_AUTO_TEST_CASE(ParallelMerge) {
1257   // Merge in phi direction with z binning
1258   auto cylPhi1 = Surface::makeShared<CylinderSurface>(
1259       Transform3::Identity() * AngleAxis3(45_degree, Vector3::UnitZ()), 30_mm,
1260       100_mm, 20_degree, 0_degree);
1261
1262   auto cylPhi2 = Surface::makeShared<CylinderSurface>(
1263       Transform3::Identity() * AngleAxis3(85_degree, Vector3::UnitZ()), 30_mm,
1264       100_mm, 20_degree, 0_degree);
1265
1266   auto portalPhi1 = GridPortalLink::make(cylPhi1, AxisDirection::AxisZ,
1267                                          Axis{AxisBound, -100_mm, 100_mm, 5});
1268
1269   auto portalPhi2 = GridPortalLink::make(cylPhi2, AxisDirection::AxisZ,
1270                                          Axis{AxisBound, -100_mm, 100_mm, 5});
1271
1272   auto merged12Ptr = GridPortalLink::merge(*portalPhi1, *portalPhi2,
1273                                            AxisDirection::AxisRPhi, *logger);
1274   BOOST_REQUIRE_NE(merged12Ptr, nullptr);
1275   auto merged12 = dynamic_cast<const GridPortalLink*>(merged12Ptr.get());
1276   BOOST_REQUIRE_NE(merged12, nullptr);
1277
1278   const auto& axis1 = *merged12->grid().axes().front();
1279   const auto& axis2 = *merged12->grid().axes().back();
1280   // Phi sectors were same size, should give equidistant binning
1281   Axis axis1Expected{AxisBound, -40_degree * 30_mm, 40_degree * 30_mm, 2};
1282   BOOST_CHECK_EQUAL(axis1, axis1Expected);
1283   Axis axis2Expected{AxisBound, -100_mm, 100_mm, 5};
1284   BOOST_CHECK_EQUAL(axis2, axis2Expected);
1285 }
1286
1287 BOOST_AUTO_TEST_SUITE_END()  // RPhiDirection
1288
1289 BOOST_AUTO_TEST_SUITE_END()  // Merging1dCylinder
1290
1291 BOOST_AUTO_TEST_SUITE(Merging2dCylinder)
1292
1293 BOOST_AUTO_TEST_CASE(ZDirection) {
1294   BOOST_TEST_CONTEXT("Phi sector") {
1295     auto cyl1 = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
1296                                                      30_mm, 100_mm, 45_degree);
1297
1298     // z good, rphi good
1299     auto grid1 = GridPortalLink::make(
1300         cyl1, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5},
1301         Axis{AxisBound, -100_mm, 100_mm, 5});
1302
1303     auto trf2 = Transform3{Translation3{Vector3::UnitZ() * 150_mm}};
1304     auto cyl2 =
1305         Surface::makeShared<CylinderSurface>(trf2, 30_mm, 50_mm, 45_degree);
1306
1307     // Second grid portal with compatible phi binning
1308     auto grid2 = GridPortalLink::make(
1309         cyl2, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5},
1310         Axis{AxisBound, -50_mm, 50_mm, 5});
1311
1312     // We're merging in z direction, so the phi binnings need to be the same
1313
1314     auto mergedPtr =
1315         GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisZ, *logger);
1316     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1317     BOOST_REQUIRE_NE(mergedPtr, nullptr);
1318
1319     const auto& axis1 = *merged->grid().axes().front();
1320     const auto& axis2 = *merged->grid().axes().back();
1321
1322     BOOST_CHECK_EQUAL(axis1.getMin(), -45_degree * 30_mm);
1323     BOOST_CHECK_EQUAL(axis1.getMax(), 45_degree * 30_mm);
1324     BOOST_CHECK_EQUAL(axis1.getNBins(), 5);
1325     BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant);
1326     BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound);
1327
1328     BOOST_CHECK_EQUAL(axis2.getMin(), -150_mm);
1329     BOOST_CHECK_EQUAL(axis2.getMax(), 150_mm);
1330     BOOST_CHECK_EQUAL(axis2.getNBins(), 10);
1331     BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Variable);
1332     BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound);
1333
1334     auto grid3 = GridPortalLink::make(
1335         cyl2, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 3},
1336         Axis{AxisBound, -50_mm, 50_mm, 5});
1337
1338     auto composite = PortalLinkBase::merge(copy(grid1), copy(grid3),
1339                                            AxisDirection::AxisZ, *logger);
1340     BOOST_CHECK_NE(dynamic_cast<const CompositePortalLink*>(composite.get()),
1341                    nullptr);
1342   }
1343
1344   BOOST_TEST_CONTEXT("Check wraparound for full circle in phi") {
1345     auto cyl1 = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
1346                                                      30_mm, 100_mm, 180_degree);
1347
1348     // z good, rphi good
1349     auto grid1 = GridPortalLink::make(
1350         cyl1, Axis{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 5},
1351         Axis{AxisBound, -100_mm, 100_mm, 5});
1352
1353     auto trf2 = Transform3{Translation3{Vector3::UnitZ() * 150_mm}};
1354     auto cyl2 =
1355         Surface::makeShared<CylinderSurface>(trf2, 30_mm, 50_mm, 180_degree);
1356
1357     // Second grid portal with compatible phi binning
1358     auto grid2 = GridPortalLink::make(
1359         cyl2, Axis{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 5},
1360         Axis{AxisBound, -50_mm, 50_mm, 5});
1361
1362     auto mergedPtr =
1363         GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisZ, *logger);
1364     const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1365     BOOST_REQUIRE_NE(mergedPtr, nullptr);
1366
1367     const auto& axis1 = *merged->grid().axes().front();
1368     const auto& axis2 = *merged->grid().axes().back();
1369
1370     Axis axis1Expected{AxisClosed, -180_degree * 30_mm, 180_degree * 30_mm, 5};
1371     BOOST_CHECK_EQUAL(axis1, axis1Expected);
1372     Axis axis2Expected{AxisBound,
1373                        {-150, -110, -70, -30, 10, 50, 70, 90, 110, 130, 150}};
1374     BOOST_CHECK_EQUAL(axis2, axis2Expected);
1375   }
1376 }
1377
1378 BOOST_AUTO_TEST_CASE(RPhiDirection) {
1379   auto cyl1 = Surface::makeShared<CylinderSurface>(Transform3::Identity(),
1380                                                    30_mm, 100_mm, 45_degree);
1381
1382   // z good, rphi good
1383   auto grid1 = GridPortalLink::make(
1384       cyl1, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5},
1385       Axis{AxisBound, -100_mm, 100_mm, 5});
1386   BOOST_REQUIRE_NE(grid1, nullptr);
1387
1388   auto trf2 = Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}};
1389   auto cyl2 =
1390       Surface::makeShared<CylinderSurface>(trf2, 30_mm, 100_mm, 45_degree);
1391
1392   // Second grid portal with compatible phi binning
1393   auto grid2 = GridPortalLink::make(
1394       cyl2, Axis{AxisBound, -45_degree * 30_mm, 45_degree * 30_mm, 5},
1395       Axis{AxisBound, -100_mm, 100_mm, 5});
1396   BOOST_REQUIRE_NE(grid2, nullptr);
1397
1398   // We're merging in z direction, so the phi binnings need to be the same
1399
1400   auto mergedPtr =
1401       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisRPhi, *logger);
1402   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1403   BOOST_REQUIRE_NE(mergedPtr, nullptr);
1404
1405   const auto& axis1 = *merged->grid().axes().front();
1406   const auto& axis2 = *merged->grid().axes().back();
1407   BOOST_CHECK_CLOSE(axis1.getMin(), -90_degree * 30_mm, 1e-8);
1408   BOOST_CHECK_CLOSE(axis1.getMax(), 90_degree * 30_mm, 1e-8);
1409   BOOST_CHECK_EQUAL(axis1.getNBins(), 10);
1410   BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant);
1411   BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound);
1412   Axis axis2Expected{AxisBound, -100_mm, 100_mm, 5};
1413   BOOST_CHECK_EQUAL(axis2, axis2Expected);
1414 }
1415
1416 BOOST_AUTO_TEST_SUITE_END()  // Merging2dCylinder
1417
1418 BOOST_AUTO_TEST_SUITE(Merging1dDisc)
1419
1420 BOOST_AUTO_TEST_SUITE(RDirection)
1421
1422 BOOST_AUTO_TEST_CASE(ColinearMerge) {
1423   // Without phi sector
1424   auto disc1 =
1425       Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 100_mm);
1426
1427   auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisR,
1428                                     Axis{AxisBound, 30_mm, 100_mm, 7});
1429
1430   auto disc2 =
1431       Surface::makeShared<DiscSurface>(Transform3::Identity(), 100_mm, 150_mm);
1432
1433   auto grid2 = GridPortalLink::make(disc2, AxisDirection::AxisR,
1434                                     Axis{AxisBound, 100_mm, 150_mm, 5});
1435
1436   auto mergedPtr =
1437       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisR, *logger);
1438   BOOST_REQUIRE(mergedPtr);
1439   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1440   BOOST_REQUIRE_NE(merged, nullptr);
1441
1442   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1443   Axis axisExpected{AxisBound, 30_mm, 150_mm, 12};
1444   BOOST_CHECK_EQUAL(*merged->grid().axes().front(), axisExpected);
1445
1446   // With phi sector
1447   auto discPhi1 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1448                                                    30_mm, 100_mm, 30_degree);
1449
1450   auto discPhiGrid1 = GridPortalLink::make(discPhi1, AxisDirection::AxisR,
1451                                            Axis{AxisBound, 30_mm, 100_mm, 7});
1452
1453   auto discPhi2 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1454                                                    100_mm, 150_mm, 30_degree);
1455
1456   auto discPhiGrid2 = GridPortalLink::make(discPhi2, AxisDirection::AxisR,
1457                                            Axis{AxisBound, 100_mm, 150_mm, 5});
1458
1459   auto mergedPhiPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid2,
1460                                             AxisDirection::AxisR, *logger);
1461   BOOST_REQUIRE(mergedPhiPtr);
1462   const auto* mergedPhi =
1463       dynamic_cast<const GridPortalLink*>(mergedPhiPtr.get());
1464   BOOST_REQUIRE_NE(mergedPhi, nullptr);
1465
1466   BOOST_CHECK_EQUAL(mergedPhi->grid().axes().size(), 1);
1467   BOOST_CHECK_EQUAL(*mergedPhi->grid().axes().front(), axisExpected);
1468 }
1469
1470 BOOST_AUTO_TEST_CASE(ParallelMerge) {
1471   auto disc1 =
1472       Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 100_mm);
1473
1474   auto grid1 =
1475       GridPortalLink::make(disc1, AxisDirection::AxisPhi,
1476                            Axis{AxisClosed, -180_degree, 180_degree, 5});
1477
1478   auto disc2 =
1479       Surface::makeShared<DiscSurface>(Transform3::Identity(), 100_mm, 150_mm);
1480
1481   auto grid2 =
1482       GridPortalLink::make(disc2, AxisDirection::AxisPhi,
1483                            Axis{AxisClosed, -180_degree, 180_degree, 5});
1484
1485   auto mergedPtr =
1486       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisR, *logger);
1487   BOOST_REQUIRE(mergedPtr);
1488   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1489   BOOST_REQUIRE_NE(merged, nullptr);
1490
1491   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2);
1492   const auto& axis1 = *merged->grid().axes().front();
1493   const auto& axis2 = *merged->grid().axes().back();
1494   Axis axis1Expected{AxisBound, {30_mm, 100_mm, 150_mm}};
1495   BOOST_CHECK_EQUAL(axis1, axis1Expected);
1496   Axis axis2Expected{AxisClosed, -180_degree, 180_degree, 5};
1497   BOOST_CHECK_EQUAL(axis2, axis2Expected);
1498 }
1499
1500 BOOST_AUTO_TEST_SUITE_END()  // RDirection
1501
1502 BOOST_AUTO_TEST_SUITE(PhiDirection)
1503
1504 BOOST_AUTO_TEST_CASE(ColinearMerge) {
1505   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1506                                                 100_mm, 30_degree);
1507
1508   auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisPhi,
1509                                     Axis{AxisBound, -30_degree, 30_degree, 3});
1510
1511   auto disc2 = Surface::makeShared<DiscSurface>(
1512       Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1513       60_degree);
1514
1515   auto grid2 = GridPortalLink::make(disc2, AxisDirection::AxisPhi,
1516                                     Axis{AxisBound, -60_degree, 60_degree, 6});
1517
1518   auto mergedPtr =
1519       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
1520   BOOST_REQUIRE(mergedPtr);
1521   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1522   BOOST_REQUIRE_NE(merged, nullptr);
1523
1524   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 1);
1525   const auto& axis = *merged->grid().axes().front();
1526   BOOST_CHECK_CLOSE(axis.getMin(), -90_degree, 1e-6);
1527   BOOST_CHECK_CLOSE(axis.getMax(), 90_degree, 1e-6);
1528   BOOST_CHECK_EQUAL(axis.getNBins(), 9);
1529
1530   // Check wrapping
1531
1532   auto disc1Half = Surface::makeShared<DiscSurface>(
1533       Transform3{AngleAxis3{15_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1534       90_degree);
1535
1536   auto grid1Half =
1537       GridPortalLink::make(disc1Half, AxisDirection::AxisPhi,
1538                            Axis{AxisBound, -90_degree, 90_degree, 3});
1539
1540   auto disc2Half = Surface::makeShared<DiscSurface>(
1541       Transform3{AngleAxis3{-165_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1542       90_degree);
1543
1544   auto grid2Half =
1545       GridPortalLink::make(disc2Half, AxisDirection::AxisPhi,
1546                            Axis{AxisBound, -90_degree, 90_degree, 3});
1547
1548   auto mergedHalfPtr = GridPortalLink::merge(*grid1Half, *grid2Half,
1549                                              AxisDirection::AxisPhi, *logger);
1550   BOOST_REQUIRE(mergedHalfPtr);
1551   const auto* mergedHalf =
1552       dynamic_cast<const GridPortalLink*>(mergedHalfPtr.get());
1553   BOOST_REQUIRE_NE(mergedHalf, nullptr);
1554
1555   BOOST_CHECK_EQUAL(mergedHalf->grid().axes().size(), 1);
1556   Axis axisHalfExpected{AxisClosed, -180_degree, 180_degree, 6};
1557   BOOST_CHECK_EQUAL(axisHalfExpected, *mergedHalf->grid().axes().front());
1558 }
1559
1560 BOOST_AUTO_TEST_CASE(ParallelMerge) {
1561   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1562                                                 100_mm, 30_degree);
1563
1564   auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisR,
1565                                     Axis{AxisBound, 30_mm, 100_mm, 3});
1566
1567   auto disc2 = Surface::makeShared<DiscSurface>(
1568       Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1569       60_degree);
1570
1571   auto grid2 = GridPortalLink::make(disc2, AxisDirection::AxisR,
1572                                     Axis{AxisBound, 30_mm, 100_mm, 3});
1573
1574   auto mergedPtr =
1575       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
1576   BOOST_REQUIRE(mergedPtr);
1577   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1578   BOOST_REQUIRE_NE(merged, nullptr);
1579
1580   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2);
1581   const auto& axis1 = *merged->grid().axes().front();
1582   const auto& axis2 = *merged->grid().axes().back();
1583
1584   Axis axis1Expected{AxisBound, 30_mm, 100_mm, 3};
1585   BOOST_CHECK_EQUAL(axis1, axis1Expected);
1586
1587   BOOST_CHECK_CLOSE(axis2.getMin(), -90_degree, 1e-6);
1588   BOOST_CHECK_CLOSE(axis2.getMax(), 90_degree, 1e-6);
1589   BOOST_CHECK_EQUAL(axis2.getNBins(), 2);
1590   BOOST_CHECK_CLOSE(axis2.getBinEdges().at(1), 30_degree, 1e-6);
1591   BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Variable);
1592   BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound);
1593 }
1594
1595 BOOST_AUTO_TEST_SUITE_END()  // PhiDirection
1596
1597 BOOST_AUTO_TEST_CASE(BinFilling) {
1598   // Volumes for bin content checking
1599
1600   auto vol1 = makeDummyVolume();
1601   auto vol2 = makeDummyVolume();
1602
1603   BOOST_TEST_CONTEXT("RDirection") {
1604     auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1605                                                   60_mm, 30_degree);
1606
1607     auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisR,
1608                                       Axis{AxisBound, 30_mm, 60_mm, 2});
1609
1610     grid1->setVolume(vol1.get());
1611
1612     auto disc2 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 60_mm,
1613                                                   90_mm, 30_degree);
1614
1615     auto grid2 = GridPortalLink::make(disc2, AxisDirection::AxisR,
1616                                       Axis{AxisBound, 60_mm, 90_mm, 2});
1617
1618     grid2->setVolume(vol2.get());
1619
1620     auto mergedPtr =
1621         GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisR, *logger);
1622
1623     using merged_type =
1624         GridPortalLinkT<Axis<AxisType::Equidistant, AxisBoundaryType::Bound>>;
1625
1626     const auto* merged = dynamic_cast<const merged_type*>(mergedPtr.get());
1627     BOOST_REQUIRE(merged);
1628
1629     grid1->printContents(std::cout);
1630     grid2->printContents(std::cout);
1631     merged->printContents(std::cout);
1632
1633     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({1}), vol1.get());
1634     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({2}), vol1.get());
1635     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({3}), vol2.get());
1636     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({4}), vol2.get());
1637   }
1638
1639   BOOST_TEST_CONTEXT("PhiDirection") {
1640     auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1641                                                   100_mm, 30_degree);
1642
1643     auto grid1 =
1644         GridPortalLink::make(disc1, AxisDirection::AxisPhi,
1645                              Axis{AxisBound, -30_degree, 30_degree, 2});
1646
1647     grid1->setVolume(vol1.get());
1648
1649     auto disc2 = Surface::makeShared<DiscSurface>(
1650         Transform3{AngleAxis3{60_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1651         30_degree);
1652
1653     auto grid2 =
1654         GridPortalLink::make(disc2, AxisDirection::AxisPhi,
1655                              Axis{AxisBound, -30_degree, 30_degree, 2});
1656
1657     grid2->setVolume(vol2.get());
1658
1659     auto mergedPtr =
1660         GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
1661     BOOST_REQUIRE(mergedPtr);
1662
1663     using merged_type =
1664         GridPortalLinkT<Axis<AxisType::Equidistant, AxisBoundaryType::Bound>>;
1665
1666     const auto* merged = dynamic_cast<const merged_type*>(mergedPtr.get());
1667     BOOST_REQUIRE(merged);
1668
1669     grid1->printContents(std::cout);
1670     grid2->printContents(std::cout);
1671     merged->printContents(std::cout);
1672
1673     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({1}), vol2.get());
1674     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({2}), vol2.get());
1675     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({3}), vol1.get());
1676     BOOST_CHECK_EQUAL(merged->grid().atLocalBins({4}), vol1.get());
1677   }
1678 }
1679
1680 BOOST_AUTO_TEST_SUITE_END()  // Merging1dDisc
1681
1682 BOOST_AUTO_TEST_SUITE(Merging2dDisc)
1683
1684 BOOST_AUTO_TEST_CASE(RDirection) {
1685   // Basic, because the parallel 1D case already tests this to some degree
1686   auto discPhi1 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1687                                                    30_mm, 100_mm, 30_degree);
1688
1689   auto discPhiGrid1 =
1690       GridPortalLink::make(discPhi1, Axis{AxisBound, 30_mm, 100_mm, 7},
1691                            Axis{AxisBound, -30_degree, 30_degree, 3});
1692
1693   auto discPhi2 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1694                                                    100_mm, 150_mm, 30_degree);
1695
1696   auto discPhiGrid2 =
1697       GridPortalLink::make(discPhi2, Axis{AxisBound, 100_mm, 150_mm, 5},
1698                            Axis{AxisBound, -30_degree, 30_degree, 3});
1699
1700   auto mergedPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid2,
1701                                          AxisDirection::AxisR, *logger);
1702   BOOST_REQUIRE(mergedPtr);
1703   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1704   BOOST_REQUIRE_NE(merged, nullptr);
1705   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2);
1706   const auto& axis1 = *merged->grid().axes().front();
1707   const auto& axis2 = *merged->grid().axes().back();
1708
1709   BOOST_CHECK_EQUAL(axis1.getMin(), 30_mm);
1710   BOOST_CHECK_EQUAL(axis1.getMax(), 150_mm);
1711   BOOST_CHECK_EQUAL(axis1.getNBins(), 12);
1712   BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant);
1713   BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound);
1714   BOOST_CHECK_EQUAL(axis2.getMin(), -30_degree);
1715   BOOST_CHECK_EQUAL(axis2.getMax(), 30_degree);
1716   BOOST_CHECK_EQUAL(axis2.getNBins(), 3);
1717   BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant);
1718 }
1719
1720 BOOST_AUTO_TEST_CASE(PhiDirection) {
1721   // Basic, because the parallel 1D case already tests this to some degree
1722   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1723                                                 100_mm, 30_degree);
1724
1725   auto grid1 = GridPortalLink::make(disc1, Axis{AxisBound, 30_mm, 100_mm, 3},
1726                                     Axis{AxisBound, -30_degree, 30_degree, 3});
1727
1728   auto disc2 = Surface::makeShared<DiscSurface>(
1729       Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1730       60_degree);
1731
1732   auto grid2 = GridPortalLink::make(disc2, Axis{AxisBound, 30_mm, 100_mm, 3},
1733                                     Axis{AxisBound, -60_degree, 60_degree, 6});
1734
1735   auto mergedPtr =
1736       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
1737   BOOST_REQUIRE(mergedPtr);
1738   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
1739   BOOST_REQUIRE_NE(merged, nullptr);
1740
1741   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2);
1742   const auto& axis1 = *merged->grid().axes().front();
1743   const auto& axis2 = *merged->grid().axes().back();
1744
1745   BOOST_CHECK_EQUAL(axis1.getMin(), 30_mm);
1746   BOOST_CHECK_EQUAL(axis1.getMax(), 100_mm);
1747   BOOST_CHECK_EQUAL(axis1.getNBins(), 3);
1748   BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant);
1749   BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound);
1750   BOOST_CHECK_CLOSE(axis2.getMin(), -90_degree, 1e-6);
1751   BOOST_CHECK_CLOSE(axis2.getMax(), 90_degree, 1e-6);
1752   BOOST_CHECK_EQUAL(axis2.getNBins(), 9);
1753   BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant);
1754   BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound);
1755 }
1756
1757 BOOST_AUTO_TEST_CASE(BinFilling) {
1758   // Volumes for bin content checking
1759   // Volume shape/transform is irrelevant, only used for pointer identity
1760   auto vol1 = std::make_shared<TrackingVolume>(
1761       Transform3::Identity(),
1762       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
1763
1764   auto vol2 = std::make_shared<TrackingVolume>(
1765       Transform3::Identity(),
1766       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
1767
1768   auto fillCheckerBoard = [&](auto& grid) {
1769     auto loc = grid.numLocalBins();
1770     for (std::size_t i = 1; i <= loc[0]; ++i) {
1771       for (std::size_t j = 1; j <= loc[1]; ++j) {
1772         grid.atLocalBins({i, j}) = (i + j) % 2 == 0 ? vol1.get() : vol2.get();
1773       }
1774     }
1775   };
1776
1777   auto checkCheckerBoard = [&](const auto& grid) {
1778     auto loc = grid.numLocalBins();
1779     for (std::size_t i = 1; i <= loc[0]; ++i) {
1780       for (std::size_t j = 1; j <= loc[1]; ++j) {
1781         const auto* vol = grid.atLocalBins({i, j});
1782         if (vol != ((i + j) % 2 == 0 ? vol1.get() : vol2.get())) {
1783           BOOST_ERROR("Is not a checkerboard pattern");
1784           return;
1785         }
1786       }
1787     }
1788   };
1789
1790   BOOST_TEST_CONTEXT("RDirection") {
1791     auto discPhi1 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1792                                                      30_mm, 60_mm, 30_degree);
1793
1794     auto discPhiGrid1 =
1795         GridPortalLink::make(discPhi1, Axis{AxisBound, 30_mm, 60_mm, 2},
1796                              Axis{AxisBound, -30_degree, 30_degree, 2});
1797
1798     fillCheckerBoard(discPhiGrid1->grid());
1799     checkCheckerBoard(discPhiGrid1->grid());
1800
1801     auto discPhi2 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1802                                                      60_mm, 90_mm, 30_degree);
1803
1804     auto discPhiGrid2 =
1805         GridPortalLink::make(discPhi2, Axis{AxisBound, 60_mm, 90_mm, 2},
1806                              Axis{AxisBound, -30_degree, 30_degree, 2});
1807
1808     fillCheckerBoard(discPhiGrid2->grid());
1809     checkCheckerBoard(discPhiGrid2->grid());
1810
1811     auto mergedPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid2,
1812                                            AxisDirection::AxisR, *logger);
1813
1814     using merged_type =
1815         GridPortalLinkT<Axis<AxisType::Equidistant, AxisBoundaryType::Bound>,
1816                         Axis<AxisType::Equidistant, AxisBoundaryType::Bound>>;
1817
1818     const auto* merged = dynamic_cast<const merged_type*>(mergedPtr.get());
1819     BOOST_REQUIRE(merged);
1820     checkCheckerBoard(merged->grid());
1821
1822     // Fill a / b
1823     discPhiGrid1->setVolume(vol1.get());
1824     discPhiGrid2->setVolume(vol2.get());
1825
1826     mergedPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid2,
1827                                       AxisDirection::AxisR, *logger);
1828
1829     merged = dynamic_cast<const merged_type*>(mergedPtr.get());
1830     BOOST_REQUIRE(merged);
1831
1832     const auto* v1 = vol1.get();
1833     const auto* v2 = vol2.get();
1834
1835     std::vector<std::pair<Vector2, const TrackingVolume*>> locations = {
1836         {{40_mm, -20_degree}, v1}, {{40_mm, 20_degree}, v1},
1837         {{50_mm, -20_degree}, v1}, {{50_mm, 20_degree}, v1},
1838         {{70_mm, -20_degree}, v2}, {{70_mm, 20_degree}, v2},
1839         {{80_mm, -20_degree}, v2}, {{80_mm, 20_degree}, v2},
1840     };
1841
1842     for (const auto& [loc, vol] : locations) {
1843       BOOST_TEST_CONTEXT(loc.transpose())
1844       BOOST_CHECK_EQUAL(merged->resolveVolume(gctx, loc).value(), vol);
1845     }
1846
1847     std::vector<std::vector<const TrackingVolume*>> contents = {
1848         {v1, v1},
1849         {v1, v1},
1850         {v2, v2},
1851         {v2, v2},
1852     };
1853
1854     for (std::size_t i = 0; i < 4; ++i) {
1855       for (std::size_t j = 0; j < 2; ++j) {
1856         BOOST_CHECK_EQUAL(merged->grid().atLocalBins({i + 1, j + 1}),
1857                           contents.at(i).at(j));
1858       }
1859     }
1860   }
1861
1862   BOOST_TEST_CONTEXT("PhiDirection") {
1863     auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
1864                                                   100_mm, 30_degree);
1865
1866     auto grid1 =
1867         GridPortalLink::make(disc1, Axis{AxisBound, 30_mm, 100_mm, 2},
1868                              Axis{AxisBound, -30_degree, 30_degree, 2});
1869     fillCheckerBoard(grid1->grid());
1870     checkCheckerBoard(grid1->grid());
1871
1872     auto disc2 = Surface::makeShared<DiscSurface>(
1873         Transform3{AngleAxis3{60_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
1874         30_degree);
1875
1876     auto grid2 =
1877         GridPortalLink::make(disc2, Axis{AxisBound, 30_mm, 100_mm, 2},
1878                              Axis{AxisBound, -30_degree, 30_degree, 2});
1879     fillCheckerBoard(grid2->grid());
1880     checkCheckerBoard(grid2->grid());
1881
1882     auto mergedPtr =
1883         GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
1884
1885     BOOST_REQUIRE(mergedPtr);
1886
1887     using merged_type =
1888         GridPortalLinkT<Axis<AxisType::Equidistant, AxisBoundaryType::Bound>,
1889                         Axis<AxisType::Equidistant, AxisBoundaryType::Bound>>;
1890
1891     const auto* merged = dynamic_cast<const merged_type*>(mergedPtr.get());
1892     BOOST_REQUIRE(merged);
1893
1894     checkCheckerBoard(merged->grid());
1895
1896     // Fill a / b
1897     grid1->setVolume(vol1.get());
1898     grid2->setVolume(vol2.get());
1899
1900     mergedPtr =
1901         GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
1902     merged = dynamic_cast<const merged_type*>(mergedPtr.get());
1903     BOOST_REQUIRE(merged);
1904
1905     const auto* v1 = vol1.get();
1906     const auto* v2 = vol2.get();
1907
1908     std::vector<std::pair<Vector2, const TrackingVolume*>> locations = {
1909         {{40_mm, -50_degree}, v2}, {{40_mm, -10_degree}, v2},
1910         {{50_mm, -50_degree}, v2}, {{50_mm, -10_degree}, v2},
1911         {{40_mm, 10_degree}, v1},  {{50_mm, 50_degree}, v1},
1912         {{50_mm, 10_degree}, v1},  {{50_mm, 50_degree}, v1},
1913     };
1914
1915     for (const auto& [loc, vol] : locations) {
1916       BOOST_TEST_CONTEXT(loc.transpose())
1917       BOOST_CHECK_EQUAL(merged->resolveVolume(gctx, loc).value(), vol);
1918     }
1919
1920     std::vector<std::vector<const TrackingVolume*>> contents = {
1921         {v2, v2, v1, v1},
1922         {v2, v2, v1, v1},
1923     };
1924
1925     for (std::size_t i = 0; i < 2; ++i) {
1926       for (std::size_t j = 0; j < 4; ++j) {
1927         BOOST_CHECK_EQUAL(merged->grid().atLocalBins({i + 1, j + 1}),
1928                           contents.at(i).at(j));
1929       }
1930     }
1931   }
1932 }
1933
1934 BOOST_AUTO_TEST_SUITE_END()  // Merging2dDisc
1935
1936 BOOST_AUTO_TEST_SUITE(MergingMixedDisc)
1937
1938 BOOST_AUTO_TEST_CASE(RDirection) {
1939   auto discPhi1 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1940                                                    30_mm, 100_mm, 30_degree);
1941
1942   auto discPhiGrid1 =
1943       GridPortalLink::make(discPhi1, Axis{AxisBound, 30_mm, 100_mm, 7},
1944                            Axis{AxisBound, -30_degree, 30_degree, 3});
1945
1946   auto discPhi2 = Surface::makeShared<DiscSurface>(Transform3::Identity(),
1947                                                    100_mm, 150_mm, 30_degree);
1948
1949   auto discPhiGrid21dPhi =
1950       GridPortalLink::make(discPhi2, AxisDirection::AxisPhi,
1951                            Axis{AxisBound, -30_degree, 30_degree, 3});
1952
1953   auto merged12PhiPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid21dPhi,
1954                                               AxisDirection::AxisR, *logger);
1955   BOOST_REQUIRE(merged12PhiPtr);
1956   const auto* merged12Phi =
1957       dynamic_cast<const GridPortalLink*>(merged12PhiPtr.get());
1958   BOOST_REQUIRE_NE(merged12Phi, nullptr);
1959
1960   auto merged21PhiPtr = GridPortalLink::merge(*discPhiGrid21dPhi, *discPhiGrid1,
1961                                               AxisDirection::AxisR, *logger);
1962   BOOST_REQUIRE(merged21PhiPtr);
1963   const auto* merged21Phi =
1964       dynamic_cast<const GridPortalLink*>(merged21PhiPtr.get());
1965   BOOST_REQUIRE_NE(merged21Phi, nullptr);
1966
1967   BOOST_CHECK_EQUAL(merged12Phi->grid(), merged21Phi->grid());
1968
1969   BOOST_CHECK_EQUAL(merged12Phi->grid().axes().size(), 2);
1970   const auto& axis1 = *merged12Phi->grid().axes().front();
1971   const auto& axis2 = *merged12Phi->grid().axes().back();
1972
1973   BOOST_CHECK_EQUAL(axis1.getMin(), 30_mm);
1974   BOOST_CHECK_EQUAL(axis1.getMax(), 150_mm);
1975   BOOST_CHECK_EQUAL(axis1.getNBins(), 8);
1976   BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Variable);
1977   BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound);
1978   BOOST_CHECK_EQUAL(axis2.getMin(), -30_degree);
1979   BOOST_CHECK_EQUAL(axis2.getMax(), 30_degree);
1980   BOOST_CHECK_EQUAL(axis2.getNBins(), 3);
1981   BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant);
1982   BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound);
1983
1984   auto discPhiGrid21dR = GridPortalLink::make(
1985       discPhi2, AxisDirection::AxisR, Axis{AxisBound, 100_mm, 150_mm, 5});
1986
1987   auto merged12RPtr = GridPortalLink::merge(*discPhiGrid1, *discPhiGrid21dR,
1988                                             AxisDirection::AxisR, *logger);
1989   BOOST_REQUIRE(merged12RPtr);
1990   const auto* merged12R =
1991       dynamic_cast<const GridPortalLink*>(merged12RPtr.get());
1992   BOOST_REQUIRE_NE(merged12R, nullptr);
1993
1994   auto merged21RPtr = GridPortalLink::merge(*discPhiGrid21dR, *discPhiGrid1,
1995                                             AxisDirection::AxisR, *logger);
1996   BOOST_REQUIRE(merged21RPtr);
1997   const auto* merged21R =
1998       dynamic_cast<const GridPortalLink*>(merged21RPtr.get());
1999   BOOST_REQUIRE_NE(merged21R, nullptr);
2000   BOOST_CHECK_EQUAL(merged12R->grid(), merged21R->grid());
2001
2002   BOOST_CHECK_EQUAL(merged12R->grid().axes().size(), 2);
2003   const auto& axis1R = *merged12R->grid().axes().front();
2004   const auto& axis2R = *merged12R->grid().axes().back();
2005
2006   BOOST_CHECK_EQUAL(axis1R.getMin(), 30_mm);
2007   BOOST_CHECK_EQUAL(axis1R.getMax(), 150_mm);
2008   BOOST_CHECK_EQUAL(axis1R.getNBins(), 12);
2009   BOOST_CHECK_EQUAL(axis1R.getType(), AxisType::Equidistant);
2010   BOOST_CHECK_EQUAL(axis1R.getBoundaryType(), AxisBoundaryType::Bound);
2011   BOOST_CHECK_EQUAL(axis2R.getMin(), -30_degree);
2012   BOOST_CHECK_EQUAL(axis2R.getMax(), 30_degree);
2013   BOOST_CHECK_EQUAL(axis2R.getNBins(), 3);
2014   BOOST_CHECK_EQUAL(axis2R.getType(), AxisType::Equidistant);
2015   BOOST_CHECK_EQUAL(axis2R.getBoundaryType(), AxisBoundaryType::Bound);
2016 }
2017
2018 BOOST_AUTO_TEST_CASE(PhiDirection) {
2019   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
2020                                                 100_mm, 30_degree);
2021
2022   auto grid1 = GridPortalLink::make(disc1, Axis{AxisBound, 30_mm, 100_mm, 3},
2023                                     Axis{AxisBound, -30_degree, 30_degree, 3});
2024
2025   auto disc2 = Surface::makeShared<DiscSurface>(
2026       Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
2027       60_degree);
2028
2029   auto grid21dPhi =
2030       GridPortalLink::make(disc2, AxisDirection::AxisPhi,
2031
2032                            Axis{AxisBound, -60_degree, 60_degree, 6});
2033
2034   auto merged12PhiPtr = GridPortalLink::merge(*grid1, *grid21dPhi,
2035                                               AxisDirection::AxisPhi, *logger);
2036   BOOST_REQUIRE(merged12PhiPtr);
2037   const auto* merged12Phi =
2038       dynamic_cast<const GridPortalLink*>(merged12PhiPtr.get());
2039   BOOST_REQUIRE_NE(merged12Phi, nullptr);
2040
2041   BOOST_CHECK_EQUAL(merged12Phi->grid().axes().size(), 2);
2042   const auto& axis1 = *merged12Phi->grid().axes().front();
2043   const auto& axis2 = *merged12Phi->grid().axes().back();
2044
2045   BOOST_CHECK_EQUAL(axis1.getMin(), 30_mm);
2046   BOOST_CHECK_EQUAL(axis1.getMax(), 100_mm);
2047   BOOST_CHECK_EQUAL(axis1.getNBins(), 3);
2048   BOOST_CHECK_EQUAL(axis1.getType(), AxisType::Equidistant);
2049   BOOST_CHECK_EQUAL(axis1.getBoundaryType(), AxisBoundaryType::Bound);
2050   BOOST_CHECK_CLOSE(axis2.getMin(), -90_degree, 1e-6);
2051   BOOST_CHECK_CLOSE(axis2.getMax(), 90_degree, 1e-6);
2052   BOOST_CHECK_EQUAL(axis2.getNBins(), 9);
2053   BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant);
2054   BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound);
2055
2056   auto grid21dR = GridPortalLink::make(disc2, AxisDirection::AxisR,
2057                                        Axis{AxisBound, 30_mm, 100_mm, 3});
2058
2059   auto merged12RPtr =
2060       GridPortalLink::merge(*grid1, *grid21dR, AxisDirection::AxisPhi, *logger);
2061   BOOST_REQUIRE(merged12RPtr);
2062   const auto* merged12R =
2063       dynamic_cast<const GridPortalLink*>(merged12RPtr.get());
2064   BOOST_REQUIRE_NE(merged12R, nullptr);
2065
2066   BOOST_CHECK_EQUAL(merged12R->grid().axes().size(), 2);
2067   const auto& axis1R = *merged12R->grid().axes().front();
2068   const auto& axis2R = *merged12R->grid().axes().back();
2069
2070   BOOST_CHECK_EQUAL(axis1R.getMin(), 30_mm);
2071   BOOST_CHECK_EQUAL(axis1R.getMax(), 100_mm);
2072   BOOST_CHECK_EQUAL(axis1R.getNBins(), 3);
2073   BOOST_CHECK_EQUAL(axis1R.getType(), AxisType::Equidistant);
2074   BOOST_CHECK_EQUAL(axis1R.getBoundaryType(), AxisBoundaryType::Bound);
2075   BOOST_CHECK_CLOSE(axis2R.getMin(), -90_degree, 1e-6);
2076   BOOST_CHECK_CLOSE(axis2R.getMax(), 90_degree, 1e-6);
2077   BOOST_CHECK_EQUAL(axis2R.getNBins(), 4);
2078   BOOST_CHECK_EQUAL(axis2R.getType(), AxisType::Variable);
2079   BOOST_CHECK_EQUAL(axis2R.getBoundaryType(), AxisBoundaryType::Bound);
2080 }
2081
2082 BOOST_AUTO_TEST_SUITE_END()  // MergingMixedDisc
2083
2084 BOOST_AUTO_TEST_SUITE(MergingCrossDisc)
2085
2086 BOOST_AUTO_TEST_CASE(RDirection) {
2087   // Volumes for bin content checking
2088   auto vol1 = makeDummyVolume();
2089   auto vol2 = makeDummyVolume();
2090   auto vol3 = makeDummyVolume();
2091   auto vol4 = makeDummyVolume();
2092
2093   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
2094                                                 100_mm, 30_degree);
2095
2096   auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisR,
2097                                     Axis{AxisBound, 30_mm, 100_mm, 2});
2098   grid1->grid().atLocalBins({1}) = vol1.get();
2099   grid1->grid().atLocalBins({2}) = vol2.get();
2100
2101   auto disc2 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 100_mm,
2102                                                 150_mm, 30_degree);
2103
2104   auto grid2 = GridPortalLink::make(disc2, AxisDirection::AxisPhi,
2105                                     Axis{AxisBound, -30_degree, 30_degree, 2});
2106
2107   grid2->grid().atLocalBins({1}) = vol3.get();
2108   grid2->grid().atLocalBins({2}) = vol4.get();
2109
2110   auto mergedPtr =
2111       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisR, *logger);
2112
2113   const auto* merged = dynamic_cast<const GridPortalLink*>(mergedPtr.get());
2114   BOOST_REQUIRE_NE(merged, nullptr);
2115
2116   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2);
2117   const auto& axis1 = *merged->grid().axes().front();
2118   const auto& axis2 = *merged->grid().axes().back();
2119   Axis axisExpectedR{AxisBound, {30_mm, 65_mm, 100_mm, 150_mm}};
2120   Axis axisExpectedPhi{AxisBound, -30_degree, 30_degree, 2};
2121   BOOST_CHECK_EQUAL(axis1, axisExpectedR);
2122   BOOST_CHECK_EQUAL(axis2, axisExpectedPhi);
2123
2124   std::vector<std::pair<Vector2, TrackingVolume*>> locations = {
2125       {{40_mm, -15_degree}, vol1.get()},  {{40_mm, 15_degree}, vol1.get()},
2126       {{90_mm, -15_degree}, vol2.get()},  {{90_mm, 15_degree}, vol2.get()},
2127
2128       {{110_mm, -15_degree}, vol3.get()}, {{110_mm, 15_degree}, vol4.get()},
2129       {{140_mm, -15_degree}, vol3.get()}, {{140_mm, 15_degree}, vol4.get()},
2130   };
2131
2132   for (const auto& [loc, vol] : locations) {
2133     BOOST_TEST_CONTEXT(loc.transpose())
2134     BOOST_CHECK_EQUAL(merged->resolveVolume(gctx, loc).value(), vol);
2135   }
2136
2137   grid1->printContents(std::cout);
2138   grid2->printContents(std::cout);
2139   merged->printContents(std::cout);
2140 }
2141
2142 BOOST_AUTO_TEST_CASE(PhiDirection) {
2143   // Volumes for bin content checking
2144
2145   auto vol1 = makeDummyVolume();
2146   auto vol2 = makeDummyVolume();
2147   auto vol3 = makeDummyVolume();
2148   auto vol4 = makeDummyVolume();
2149
2150   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
2151                                                 100_mm, 30_degree);
2152
2153   auto grid1 = GridPortalLink::make(disc1, AxisDirection::AxisR,
2154                                     Axis{AxisBound, 30_mm, 100_mm, 2});
2155
2156   grid1->grid().atLocalBins({1}) = vol1.get();
2157   grid1->grid().atLocalBins({2}) = vol2.get();
2158
2159   auto disc2 = Surface::makeShared<DiscSurface>(
2160       Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
2161       60_degree);
2162
2163   auto grid2 = GridPortalLink::make(disc2, AxisDirection::AxisPhi,
2164                                     Axis{AxisBound, -60_degree, 60_degree, 2});
2165
2166   grid2->grid().atLocalBins({1}) = vol3.get();
2167   grid2->grid().atLocalBins({2}) = vol4.get();
2168
2169   auto mergedPtr =
2170       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisPhi, *logger);
2171
2172   using merged_type =
2173       GridPortalLinkT<Axis<AxisType::Equidistant, AxisBoundaryType::Bound>,
2174                       Axis<AxisType::Equidistant, AxisBoundaryType::Bound>>;
2175
2176   const auto* merged = dynamic_cast<const merged_type*>(mergedPtr.get());
2177   BOOST_REQUIRE_NE(merged, nullptr);
2178
2179   BOOST_CHECK_EQUAL(merged->grid().axes().size(), 2);
2180   const auto& axis1 = *merged->grid().axes().front();
2181   const auto& axis2 = *merged->grid().axes().back();
2182   Axis axisExpectedR{AxisBound, 30_mm, 100_mm, 2};
2183   BOOST_CHECK_EQUAL(axis1, axisExpectedR);
2184
2185   BOOST_CHECK_CLOSE(axis2.getMin(), -90_degree, 1e-6);
2186   BOOST_CHECK_CLOSE(axis2.getMax(), 90_degree, 1e-6);
2187   BOOST_CHECK_EQUAL(axis2.getNBins(), 3);
2188   BOOST_CHECK_EQUAL(axis2.getType(), AxisType::Equidistant);
2189   BOOST_CHECK_EQUAL(axis2.getBoundaryType(), AxisBoundaryType::Bound);
2190
2191   std::vector<std::pair<Vector2, TrackingVolume*>> locations = {
2192       {{40_mm, 45_degree}, vol1.get()},  {{40_mm, 0_degree}, vol4.get()},
2193       {{40_mm, -80_degree}, vol3.get()}, {{90_mm, 45_degree}, vol2.get()},
2194       {{90_mm, 0_degree}, vol4.get()},   {{90_mm, -80_degree}, vol3.get()},
2195   };
2196
2197   grid1->printContents(std::cout);
2198   grid2->printContents(std::cout);
2199   merged->printContents(std::cout);
2200
2201   for (const auto& [loc, vol] : locations) {
2202     BOOST_TEST_CONTEXT((Vector2{loc[0], loc[1] / 1_degree}.transpose()))
2203     BOOST_CHECK_EQUAL(merged->resolveVolume(gctx, loc).value(), vol);
2204   }
2205 }
2206
2207 BOOST_AUTO_TEST_SUITE_END()  // MergeCrossDisc
2208
2209 BOOST_AUTO_TEST_SUITE(Merging1dPlane)
2210
2211 BOOST_AUTO_TEST_CASE(ColinearMerge) {
2212   auto rBounds1 = std::make_shared<const RectangleBounds>(30_mm, 100_mm);
2213
2214   auto plane1 =
2215       Surface::makeShared<PlaneSurface>(Transform3::Identity(), rBounds1);
2216
2217   auto gridX1 = GridPortalLink::make(plane1, AxisDirection::AxisX,
2218                                      Axis{AxisBound, -30_mm, 30_mm, 6});
2219   auto gridY1 = GridPortalLink::make(plane1, AxisDirection::AxisY,
2220                                      Axis{AxisBound, -100_mm, 100_mm, 10});
2221
2222   auto rBounds2 = std::make_shared<const RectangleBounds>(80_mm, 100_mm);
2223   Translation3 offsetX{110_mm, 0., 0.};
2224   auto planeX2 = Surface::makeShared<PlaneSurface>(
2225       Transform3::Identity() * offsetX, rBounds2);
2226
2227   auto rBounds3 = std::make_shared<const RectangleBounds>(30_mm, 20_mm);
2228   Translation3 offsetY{0, 120_mm, 0.};
2229   auto planeY2 = Surface::makeShared<PlaneSurface>(
2230       Transform3::Identity() * offsetY, rBounds3);
2231
2232   auto gridX2 = GridPortalLink::make(planeX2, AxisDirection::AxisX,
2233                                      Axis{AxisBound, -80_mm, 80_mm, 16});
2234   auto gridY2 = GridPortalLink::make(planeY2, AxisDirection::AxisY,
2235                                      Axis{AxisBound, -20_mm, 20_mm, 2});
2236
2237   auto mergedPtrX =
2238       GridPortalLink::merge(*gridX1, *gridX2, AxisDirection::AxisX, *logger);
2239   BOOST_REQUIRE(mergedPtrX);
2240   const auto* mergedX = dynamic_cast<const GridPortalLink*>(mergedPtrX.get());
2241   BOOST_REQUIRE_NE(mergedX, nullptr);
2242
2243   BOOST_CHECK_EQUAL(mergedX->grid().axes().size(), 1);
2244   Axis axisXExpected{AxisBound, -110_mm, 110_mm, 22};
2245   BOOST_CHECK_EQUAL(*mergedX->grid().axes().front(), axisXExpected);
2246
2247   auto mergedPtrY =
2248       GridPortalLink::merge(*gridY1, *gridY2, AxisDirection::AxisY, *logger);
2249   BOOST_REQUIRE(mergedPtrY);
2250   const auto* mergedY = dynamic_cast<const GridPortalLink*>(mergedPtrY.get());
2251   BOOST_REQUIRE_NE(mergedY, nullptr);
2252
2253   BOOST_CHECK_EQUAL(mergedY->grid().axes().size(), 1);
2254   Axis axisYExpected{AxisBound, -120_mm, 120_mm, 12};
2255   BOOST_CHECK_EQUAL(*mergedY->grid().axes().front(), axisYExpected);
2256 }
2257
2258 BOOST_AUTO_TEST_CASE(ParallelMerge) {
2259   auto rBounds = std::make_shared<const RectangleBounds>(30_mm, 100_mm);
2260   auto plane1 =
2261       Surface::makeShared<PlaneSurface>(Transform3::Identity(), rBounds);
2262
2263   auto grid1X = GridPortalLink::make(plane1, AxisDirection::AxisX,
2264                                      Axis{AxisBound, -30_mm, 30_mm, 6});
2265
2266   auto grid1Y = GridPortalLink::make(plane1, AxisDirection::AxisY,
2267                                      Axis{AxisBound, -100_mm, 100_mm, 5});
2268
2269   Translation3 offsetX{60_mm, 0, 0.};
2270   auto plane2 = Surface::makeShared<PlaneSurface>(
2271       Transform3::Identity() * offsetX, rBounds);
2272   auto grid2 = GridPortalLink::make(plane2, AxisDirection::AxisY,
2273                                     Axis{AxisBound, -100_mm, 100_mm, 5});
2274
2275   Translation3 offsetY{0, 200_mm, 0.};
2276   auto plane3 = Surface::makeShared<PlaneSurface>(
2277       Transform3::Identity() * offsetY, rBounds);
2278   auto grid3 = GridPortalLink::make(plane3, AxisDirection::AxisX,
2279                                     Axis{AxisBound, -30_mm, 30_mm, 6});
2280
2281   auto mergedPtrX =
2282       GridPortalLink::merge(*grid1Y, *grid2, AxisDirection::AxisX, *logger);
2283   BOOST_REQUIRE(mergedPtrX);
2284   const auto* mergedX = dynamic_cast<const GridPortalLink*>(mergedPtrX.get());
2285   BOOST_REQUIRE_NE(mergedX, nullptr);
2286
2287   BOOST_CHECK_EQUAL(mergedX->grid().axes().size(), 2);
2288   const auto& axisX1 = *mergedX->grid().axes().front();
2289   const auto& axisX2 = *mergedX->grid().axes().back();
2290   Axis axisX1Expected{AxisBound, -60_mm, 60_mm, 2};
2291   BOOST_CHECK_EQUAL(axisX1, axisX1Expected);
2292   Axis axisX2Expected{AxisBound, -100_mm, 100_mm, 5};
2293   BOOST_CHECK_EQUAL(axisX2, axisX2Expected);
2294
2295   auto mergedPtrY =
2296       GridPortalLink::merge(*grid1X, *grid3, AxisDirection::AxisY, *logger);
2297   BOOST_REQUIRE(mergedPtrY);
2298   const auto* mergedY = dynamic_cast<const GridPortalLink*>(mergedPtrY.get());
2299   BOOST_REQUIRE_NE(mergedY, nullptr);
2300
2301   BOOST_CHECK_EQUAL(mergedY->grid().axes().size(), 2);
2302   const auto& axisY1 = *mergedY->grid().axes().front();
2303   const auto& axisY2 = *mergedY->grid().axes().back();
2304   Axis axisY1Expected{AxisBound, -30_mm, 30_mm, 6};
2305   BOOST_CHECK_EQUAL(axisY1, axisY1Expected);
2306   Axis axisY2Expected{AxisBound, -200_mm, 200_mm, 2};
2307   BOOST_CHECK_EQUAL(axisY2, axisY2Expected);
2308 }
2309
2310 BOOST_AUTO_TEST_CASE(BinFilling) {
2311   // Volumes for bin content checking
2312   auto vol1 = std::make_shared<TrackingVolume>(
2313       Transform3::Identity(),
2314       std::make_shared<CuboidVolumeBounds>(30_mm, 40_mm, 100_mm));
2315   auto vol2 = std::make_shared<TrackingVolume>(
2316       Transform3::Identity(),
2317       std::make_shared<CuboidVolumeBounds>(30_mm, 40_mm, 100_mm));
2318
2319   auto rBounds = std::make_shared<const RectangleBounds>(30_mm, 100_mm);
2320   auto plane1 =
2321       Surface::makeShared<PlaneSurface>(Transform3::Identity(), rBounds);
2322
2323   auto grid1X = GridPortalLink::make(plane1, AxisDirection::AxisX,
2324                                      Axis{AxisBound, -30_mm, 30_mm, 2});
2325   grid1X->setVolume(vol1.get());
2326
2327   auto grid1Y = GridPortalLink::make(plane1, AxisDirection::AxisY,
2328                                      Axis{AxisBound, -100_mm, 100_mm, 2});
2329   grid1Y->setVolume(vol1.get());
2330
2331   Translation3 offsetX{60_mm, 0., 0.};
2332   auto plane2 = Surface::makeShared<PlaneSurface>(
2333       Transform3::Identity() * offsetX, rBounds);
2334   auto grid2 = GridPortalLink::make(plane2, AxisDirection::AxisX,
2335                                     Axis{AxisBound, -30_mm, 30_mm, 2});
2336   grid2->setVolume(vol2.get());
2337
2338   Translation3 offsetY{0., 200_mm, 0.};
2339   auto plane3 = Surface::makeShared<PlaneSurface>(
2340       Transform3::Identity() * offsetY, rBounds);
2341   auto grid3 = GridPortalLink::make(plane3, AxisDirection::AxisY,
2342                                     Axis{AxisBound, -100_mm, 100_mm, 2});
2343   grid3->setVolume(vol2.get());
2344
2345   auto mergedPtrX =
2346       GridPortalLink::merge(*grid1X, *grid2, AxisDirection::AxisX, *logger);
2347
2348   auto mergedPtrY =
2349       GridPortalLink::merge(*grid1Y, *grid3, AxisDirection::AxisY, *logger);
2350
2351   using merged_type =
2352       GridPortalLinkT<Axis<AxisType::Equidistant, AxisBoundaryType::Bound>>;
2353
2354   const auto* mergedX = dynamic_cast<const merged_type*>(mergedPtrX.get());
2355   BOOST_REQUIRE(mergedX);
2356
2357   grid1X->printContents(std::cout);
2358   grid2->printContents(std::cout);
2359   mergedX->printContents(std::cout);
2360
2361   BOOST_CHECK_EQUAL(mergedX->grid().atLocalBins({1}), vol1.get());
2362   BOOST_CHECK_EQUAL(mergedX->grid().atLocalBins({2}), vol1.get());
2363   BOOST_CHECK_EQUAL(mergedX->grid().atLocalBins({3}), vol2.get());
2364   BOOST_CHECK_EQUAL(mergedX->grid().atLocalBins({4}), vol2.get());
2365
2366   const auto* mergedY = dynamic_cast<const merged_type*>(mergedPtrX.get());
2367   BOOST_REQUIRE(mergedY);
2368
2369   BOOST_CHECK_EQUAL(mergedY->grid().atLocalBins({1}), vol1.get());
2370   BOOST_CHECK_EQUAL(mergedY->grid().atLocalBins({2}), vol1.get());
2371   BOOST_CHECK_EQUAL(mergedY->grid().atLocalBins({3}), vol2.get());
2372   BOOST_CHECK_EQUAL(mergedY->grid().atLocalBins({4}), vol2.get());
2373
2374   grid1X->printContents(std::cout);
2375   grid2->printContents(std::cout);
2376   grid3->printContents(std::cout);
2377   mergedX->printContents(std::cout);
2378   mergedY->printContents(std::cout);
2379 }
2380
2381 BOOST_AUTO_TEST_SUITE_END()  // Merging1dPlane
2382
2383 BOOST_AUTO_TEST_SUITE(Merging2dPlane)
2384
2385 BOOST_AUTO_TEST_CASE(XYDirection) {
2386   // Basic, because the parallel 1D case already tests this to some degree
2387   auto rBounds = std::make_shared<const RectangleBounds>(30_mm, 100_mm);
2388   auto plane1 =
2389       Surface::makeShared<PlaneSurface>(Transform3::Identity(), rBounds);
2390
2391   auto grid1 = GridPortalLink::make(plane1, Axis{AxisBound, -30_mm, 30_mm, 10},
2392                                     Axis{AxisBound, -100_mm, 100_mm, 3});
2393
2394   Translation3 offsetX{60_mm, 0., 0.};
2395   auto plane2 = Surface::makeShared<PlaneSurface>(
2396       Transform3::Identity() * offsetX, rBounds);
2397   auto grid2 = GridPortalLink::make(plane2, Axis{AxisBound, -30_mm, 30_mm, 10},
2398                                     Axis{AxisBound, -100_mm, 100_mm, 3});
2399
2400   Translation3 offsetY{0., 200_mm, 0.};
2401   auto plane3 = Surface::makeShared<PlaneSurface>(
2402       Transform3::Identity() * offsetY, rBounds);
2403   auto grid3 = GridPortalLink::make(plane3, Axis{AxisBound, -30_mm, 30_mm, 10},
2404                                     Axis{AxisBound, -100_mm, 100_mm, 3});
2405
2406   auto mergedPtrX =
2407       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisX, *logger);
2408
2409   auto mergedPtrY =
2410       GridPortalLink::merge(*grid1, *grid3, AxisDirection::AxisY, *logger);
2411
2412   BOOST_REQUIRE(mergedPtrX);
2413   const auto* mergedX = dynamic_cast<const GridPortalLink*>(mergedPtrX.get());
2414   BOOST_REQUIRE_NE(mergedX, nullptr);
2415   BOOST_CHECK_EQUAL(mergedX->grid().axes().size(), 2);
2416   const auto& axisX1 = *mergedX->grid().axes().front();
2417   const auto& axisX2 = *mergedX->grid().axes().back();
2418
2419   BOOST_CHECK_EQUAL(axisX1.getMin(), -60_mm);
2420   BOOST_CHECK_EQUAL(axisX1.getMax(), 60_mm);
2421   BOOST_CHECK_EQUAL(axisX1.getNBins(), 20);
2422   BOOST_CHECK_EQUAL(axisX1.getType(), AxisType::Equidistant);
2423   BOOST_CHECK_EQUAL(axisX1.getBoundaryType(), AxisBoundaryType::Bound);
2424   BOOST_CHECK_EQUAL(axisX2.getMin(), -100_mm);
2425   BOOST_CHECK_EQUAL(axisX2.getMax(), 100_mm);
2426   BOOST_CHECK_EQUAL(axisX2.getNBins(), 3);
2427   BOOST_CHECK_EQUAL(axisX2.getType(), AxisType::Equidistant);
2428
2429   BOOST_REQUIRE(mergedPtrY);
2430   const auto* mergedY = dynamic_cast<const GridPortalLink*>(mergedPtrY.get());
2431   BOOST_REQUIRE_NE(mergedY, nullptr);
2432   BOOST_CHECK_EQUAL(mergedY->grid().axes().size(), 2);
2433   const auto& axisY1 = *mergedY->grid().axes().front();
2434   const auto& axisY2 = *mergedY->grid().axes().back();
2435
2436   BOOST_CHECK_EQUAL(axisY1.getMin(), -30_mm);
2437   BOOST_CHECK_EQUAL(axisY1.getMax(), 30_mm);
2438   BOOST_CHECK_EQUAL(axisY1.getNBins(), 10);
2439   BOOST_CHECK_EQUAL(axisY1.getType(), AxisType::Equidistant);
2440   BOOST_CHECK_EQUAL(axisY1.getBoundaryType(), AxisBoundaryType::Bound);
2441   BOOST_CHECK_EQUAL(axisY2.getMin(), -200_mm);
2442   BOOST_CHECK_EQUAL(axisY2.getMax(), 200_mm);
2443   BOOST_CHECK_EQUAL(axisY2.getNBins(), 6);
2444   BOOST_CHECK_EQUAL(axisY2.getType(), AxisType::Equidistant);
2445 }
2446
2447 BOOST_AUTO_TEST_CASE(BinFilling) {
2448   // Volumes for bin content checking
2449   // Volume shape/transform is irrelevant, only used for pointer identity
2450   auto vol1 = std::make_shared<TrackingVolume>(
2451       Transform3::Identity(),
2452       std::make_shared<CuboidVolumeBounds>(30_mm, 40_mm, 100_mm));
2453
2454   auto vol2 = std::make_shared<TrackingVolume>(
2455       Transform3::Identity(),
2456       std::make_shared<CuboidVolumeBounds>(30_mm, 40_mm, 100_mm));
2457
2458   auto fillCheckerBoard = [&](auto& grid) {
2459     auto loc = grid.numLocalBins();
2460     for (std::size_t i = 1; i <= loc[0]; ++i) {
2461       for (std::size_t j = 1; j <= loc[1]; ++j) {
2462         grid.atLocalBins({i, j}) = (i + j) % 2 == 0 ? vol1.get() : vol2.get();
2463       }
2464     }
2465   };
2466
2467   auto checkCheckerBoard = [&](const auto& grid) {
2468     auto loc = grid.numLocalBins();
2469     for (std::size_t i = 1; i <= loc[0]; ++i) {
2470       for (std::size_t j = 1; j <= loc[1]; ++j) {
2471         const auto* vol = grid.atLocalBins({i, j});
2472         if (vol != ((i + j) % 2 == 0 ? vol1.get() : vol2.get())) {
2473           BOOST_ERROR("Is not a checkerboard pattern");
2474           return;
2475         }
2476       }
2477     }
2478   };
2479
2480   // Basic, because the parallel 1D case already tests this to some degree
2481   auto rBounds = std::make_shared<const RectangleBounds>(30_mm, 100_mm);
2482   auto plane1 =
2483       Surface::makeShared<PlaneSurface>(Transform3::Identity(), rBounds);
2484
2485   auto grid1 = GridPortalLink::make(plane1, Axis{AxisBound, -30_mm, 30_mm, 2},
2486                                     Axis{AxisBound, -100_mm, 100_mm, 2});
2487
2488   Translation3 offsetX{60_mm, 0., 0.};
2489   auto plane2 = Surface::makeShared<PlaneSurface>(
2490       Transform3::Identity() * offsetX, rBounds);
2491   auto grid2 = GridPortalLink::make(plane2, Axis{AxisBound, -30_mm, 30_mm, 2},
2492                                     Axis{AxisBound, -100_mm, 100_mm, 2});
2493
2494   Translation3 offsetY{0., 200_mm, 0.};
2495   auto plane3 = Surface::makeShared<PlaneSurface>(
2496       Transform3::Identity() * offsetY, rBounds);
2497   auto grid3 = GridPortalLink::make(plane3, Axis{AxisBound, -30_mm, 30_mm, 2},
2498                                     Axis{AxisBound, -100_mm, 100_mm, 2});
2499
2500   fillCheckerBoard(grid1->grid());
2501   checkCheckerBoard(grid1->grid());
2502
2503   fillCheckerBoard(grid2->grid());
2504   checkCheckerBoard(grid2->grid());
2505
2506   fillCheckerBoard(grid3->grid());
2507   checkCheckerBoard(grid3->grid());
2508
2509   auto mergedPtrX =
2510       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisX, *logger);
2511   auto mergedPtrY =
2512       GridPortalLink::merge(*grid1, *grid3, AxisDirection::AxisY, *logger);
2513
2514   using merged_type =
2515       GridPortalLinkT<Axis<AxisType::Equidistant, AxisBoundaryType::Bound>,
2516                       Axis<AxisType::Equidistant, AxisBoundaryType::Bound>>;
2517
2518   const auto* mergedX = dynamic_cast<const merged_type*>(mergedPtrX.get());
2519   BOOST_REQUIRE(mergedX);
2520   checkCheckerBoard(mergedX->grid());
2521
2522   const auto* mergedY = dynamic_cast<const merged_type*>(mergedPtrY.get());
2523   BOOST_REQUIRE(mergedY);
2524   checkCheckerBoard(mergedY->grid());
2525
2526   // Fill a / b
2527   grid1->setVolume(vol1.get());
2528   grid2->setVolume(vol2.get());
2529   grid3->setVolume(vol2.get());
2530
2531   mergedPtrX =
2532       GridPortalLink::merge(*grid1, *grid2, AxisDirection::AxisX, *logger);
2533   mergedPtrY =
2534       GridPortalLink::merge(*grid1, *grid3, AxisDirection::AxisY, *logger);
2535
2536   mergedX = dynamic_cast<const merged_type*>(mergedPtrX.get());
2537   BOOST_REQUIRE(mergedX);
2538
2539   mergedY = dynamic_cast<const merged_type*>(mergedPtrY.get());
2540   BOOST_REQUIRE(mergedY);
2541
2542   const auto* v1 = vol1.get();
2543   const auto* v2 = vol2.get();
2544
2545   std::vector<std::pair<Vector2, const TrackingVolume*>> locationsX = {
2546       {{-45_mm, 0_mm}, v1},
2547       {{-15_mm, 0_mm}, v1},
2548       {{15_mm, 0_mm}, v2},
2549       {{45_mm, 0_mm}, v2},
2550   };
2551   std::vector<std::pair<Vector2, const TrackingVolume*>> locationsY = {
2552       {{0_mm, -150_mm}, v1},
2553       {{0_mm, -50_mm}, v1},
2554       {{0_mm, 50_mm}, v2},
2555       {{0_mm, 150_mm}, v2},
2556   };
2557
2558   for (const auto& [loc, vol] : locationsX) {
2559     BOOST_TEST_CONTEXT(loc.transpose())
2560     BOOST_CHECK_EQUAL(mergedX->resolveVolume(gctx, loc).value(), vol);
2561   }
2562   for (const auto& [loc, vol] : locationsY) {
2563     BOOST_TEST_CONTEXT(loc.transpose())
2564     BOOST_CHECK_EQUAL(mergedY->resolveVolume(gctx, loc).value(), vol);
2565   }
2566
2567   std::vector<std::vector<const TrackingVolume*>> contentsX = {
2568       {v1, v1},
2569       {v1, v1},
2570       {v2, v2},
2571       {v2, v2},
2572   };
2573   std::vector<std::vector<const TrackingVolume*>> contentsY = {
2574       {v1, v1, v2, v2},
2575       {v1, v1, v2, v2},
2576   };
2577
2578   for (std::size_t i = 0; i < 4; ++i) {
2579     for (std::size_t j = 0; j < 2; ++j) {
2580       BOOST_CHECK_EQUAL(mergedX->grid().atLocalBins({i + 1, j + 1}),
2581                         contentsX.at(i).at(j));
2582     }
2583   }
2584   for (std::size_t i = 0; i < 2; ++i) {
2585     for (std::size_t j = 0; j < 4; ++j) {
2586       BOOST_CHECK_EQUAL(mergedY->grid().atLocalBins({i + 1, j + 1}),
2587                         contentsY.at(i).at(j));
2588     }
2589   }
2590 }
2591
2592 BOOST_AUTO_TEST_SUITE_END()  // Merging2dPlane
2593
2594 BOOST_AUTO_TEST_SUITE_END()  // GridMerging
2595
2596 BOOST_AUTO_TEST_CASE(CompositeConstruction) {
2597   // Arbitrary volumes for testing only
2598   auto vol1 = std::make_shared<TrackingVolume>(
2599       Transform3::Identity(),
2600       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2601
2602   auto vol2 = std::make_shared<TrackingVolume>(
2603       Transform3::Identity(),
2604       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2605
2606   auto vol3 = std::make_shared<TrackingVolume>(
2607       Transform3::Identity(),
2608       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2609
2610   auto disc1 =
2611       Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 60_mm);
2612
2613   auto trivial1 = std::make_unique<TrivialPortalLink>(disc1, *vol1);
2614
2615   auto disc2 =
2616       Surface::makeShared<DiscSurface>(Transform3::Identity(), 60_mm, 90_mm);
2617   auto trivial2 = std::make_unique<TrivialPortalLink>(disc2, *vol2);
2618
2619   auto composite = std::make_unique<CompositePortalLink>(
2620       copy(trivial1), copy(trivial2), AxisDirection::AxisR);
2621
2622   auto compositeCopy = std::make_unique<CompositePortalLink>(
2623       copy(trivial1), copy(trivial2), AxisDirection::AxisR);
2624
2625   BOOST_CHECK_EQUAL(
2626       composite->resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(),
2627       vol1.get());
2628   BOOST_CHECK_EQUAL(
2629       composite->resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(),
2630       vol2.get());
2631
2632   BOOST_CHECK_EQUAL(composite->depth(), 1);
2633
2634   // Test exception on different surface types
2635   auto cyl = Surface::makeShared<CylinderSurface>(Transform3::Identity(), 30_mm,
2636                                                   40_mm);
2637   auto trivialCyl = std::make_unique<TrivialPortalLink>(cyl, *vol3);
2638   BOOST_CHECK_THROW(CompositePortalLink(copy(trivial1), copy(trivialCyl),
2639                                         AxisDirection::AxisR),
2640                     std::invalid_argument);
2641
2642   auto disc3 =
2643       Surface::makeShared<DiscSurface>(Transform3::Identity(), 90_mm, 120_mm);
2644   auto trivial3 = std::make_unique<TrivialPortalLink>(disc3, *vol3);
2645
2646   // Test exception on un-mergable surfaces
2647   BOOST_CHECK_THROW(
2648       CompositePortalLink(copy(trivial1), copy(trivial3), AxisDirection::AxisR),
2649       SurfaceMergingException);
2650
2651   // Composite with a composite (this should work regardless of flattening)
2652
2653   CompositePortalLink composite2(std::move(composite), copy(trivial3),
2654                                  AxisDirection::AxisR, false);
2655
2656   BOOST_CHECK_EQUAL(
2657       composite2.resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(),
2658       vol1.get());
2659   BOOST_CHECK_EQUAL(
2660       composite2.resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(),
2661       vol2.get());
2662   BOOST_CHECK_EQUAL(
2663       composite2.resolveVolume(gctx, Vector2{100_mm, 0_degree}).value(),
2664       vol3.get());
2665
2666   // Two without flattening
2667   BOOST_CHECK_EQUAL(composite2.depth(), 2);
2668
2669   CompositePortalLink composite2Flat(std::move(compositeCopy), copy(trivial3),
2670                                      AxisDirection::AxisR, true);
2671
2672   // One because of flattening
2673   BOOST_CHECK_EQUAL(composite2Flat.depth(), 1);
2674
2675   BOOST_CHECK_EQUAL(
2676       composite2Flat.resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(),
2677       vol1.get());
2678   BOOST_CHECK_EQUAL(
2679       composite2Flat.resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(),
2680       vol2.get());
2681   BOOST_CHECK_EQUAL(
2682       composite2Flat.resolveVolume(gctx, Vector2{100_mm, 0_degree}).value(),
2683       vol3.get());
2684 }
2685
2686 BOOST_AUTO_TEST_SUITE(PortalMerging)
2687
2688 BOOST_DATA_TEST_CASE(TrivialTrivial,
2689                      (boost::unit_test::data::make(0, -135, 180, 45) *
2690                       boost::unit_test::data::make(Vector3{0_mm, 0_mm, 0_mm},
2691                                                    Vector3{20_mm, 0_mm, 0_mm},
2692                                                    Vector3{0_mm, 20_mm, 0_mm},
2693                                                    Vector3{20_mm, 20_mm, 0_mm},
2694                                                    Vector3{0_mm, 0_mm, 20_mm})),
2695                      angle, offset) {
2696   Transform3 base =
2697       AngleAxis3(angle * 1_degree, Vector3::UnitX()) * Translation3(offset);
2698
2699   BOOST_TEST_CONTEXT("RDirection") {
2700     auto vol1 = std::make_shared<TrackingVolume>(
2701         base, std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2702
2703     auto vol2 = std::make_shared<TrackingVolume>(
2704         base, std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2705
2706     auto vol3 = std::make_shared<TrackingVolume>(
2707         base, std::make_shared<CylinderVolumeBounds>(40_mm, 50_mm, 100_mm));
2708
2709     auto disc1 =
2710         Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 60_mm);
2711
2712     auto disc2 =
2713         Surface::makeShared<DiscSurface>(Transform3::Identity(), 60_mm, 90_mm);
2714
2715     auto disc3 =
2716         Surface::makeShared<DiscSurface>(Transform3::Identity(), 90_mm, 120_mm);
2717
2718     auto trivial1 = std::make_unique<TrivialPortalLink>(disc1, *vol1);
2719     BOOST_REQUIRE(trivial1);
2720     auto trivial2 = std::make_unique<TrivialPortalLink>(disc2, *vol2);
2721     BOOST_REQUIRE(trivial2);
2722     auto trivial3 = std::make_unique<TrivialPortalLink>(disc3, *vol3);
2723     BOOST_REQUIRE(trivial3);
2724
2725     auto grid1 = trivial1->makeGrid(AxisDirection::AxisR);
2726     auto compGridTrivial = PortalLinkBase::merge(
2727         std::move(grid1), std::make_unique<TrivialPortalLink>(*trivial2),
2728         AxisDirection::AxisR, *logger);
2729     BOOST_REQUIRE(compGridTrivial);
2730     BOOST_CHECK_EQUAL(dynamic_cast<CompositePortalLink&>(*compGridTrivial)
2731                           .makeGrid(gctx, *logger),
2732                       nullptr);
2733
2734     auto composite =
2735         PortalLinkBase::merge(std::move(trivial1), std::move(trivial2),
2736                               AxisDirection::AxisR, *logger);
2737     BOOST_REQUIRE(composite);
2738
2739     auto grid12 =
2740         dynamic_cast<CompositePortalLink&>(*composite).makeGrid(gctx, *logger);
2741     BOOST_REQUIRE(grid12);
2742
2743     BOOST_CHECK_EQUAL(
2744         grid12->resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(),
2745         vol1.get());
2746
2747     BOOST_CHECK_EQUAL(
2748         grid12->resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(),
2749         vol2.get());
2750
2751     composite = PortalLinkBase::merge(std::move(composite), std::move(trivial3),
2752                                       AxisDirection::AxisR, *logger);
2753     BOOST_REQUIRE(composite);
2754
2755     auto grid123 =
2756         dynamic_cast<CompositePortalLink&>(*composite).makeGrid(gctx, *logger);
2757     BOOST_REQUIRE(grid123);
2758
2759     BOOST_CHECK_EQUAL(
2760         grid123->resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(),
2761         vol1.get());
2762
2763     BOOST_CHECK_EQUAL(
2764         grid123->resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(),
2765         vol2.get());
2766
2767     BOOST_CHECK_EQUAL(
2768         grid123->resolveVolume(gctx, Vector2{100_mm, 0_degree}).value(),
2769         vol3.get());
2770   }
2771
2772   BOOST_TEST_CONTEXT("ZDirection") {
2773     auto vol1 = std::make_shared<TrackingVolume>(
2774         base, std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2775
2776     auto vol2 = std::make_shared<TrackingVolume>(
2777         base * Translation3{Vector3::UnitZ() * 200},
2778         std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2779
2780     auto vol3 = std::make_shared<TrackingVolume>(
2781         base * Translation3{Vector3::UnitZ() * 400},
2782         std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2783
2784     auto cyl1 = Surface::makeShared<CylinderSurface>(base, 40_mm, 100_mm);
2785
2786     auto cyl2 = Surface::makeShared<CylinderSurface>(
2787         base * Translation3{Vector3::UnitZ() * 200}, 40_mm, 100_mm);
2788
2789     auto cyl3 = Surface::makeShared<CylinderSurface>(
2790         base * Translation3{Vector3::UnitZ() * 400}, 40_mm, 100_mm);
2791
2792     auto trivial1 = std::make_unique<TrivialPortalLink>(cyl1, *vol1);
2793     BOOST_REQUIRE(trivial1);
2794     auto trivial2 = std::make_unique<TrivialPortalLink>(cyl2, *vol2);
2795     BOOST_REQUIRE(trivial2);
2796     auto trivial3 = std::make_unique<TrivialPortalLink>(cyl3, *vol3);
2797     BOOST_REQUIRE(trivial3);
2798
2799     auto grid1 = trivial1->makeGrid(AxisDirection::AxisZ);
2800     auto compGridTrivial = PortalLinkBase::merge(
2801         std::move(grid1), std::make_unique<TrivialPortalLink>(*trivial2),
2802         AxisDirection::AxisZ, *logger);
2803     BOOST_REQUIRE(compGridTrivial);
2804     BOOST_CHECK_EQUAL(dynamic_cast<CompositePortalLink&>(*compGridTrivial)
2805                           .makeGrid(gctx, *logger),
2806                       nullptr);
2807
2808     auto composite =
2809         PortalLinkBase::merge(std::move(trivial1), std::move(trivial2),
2810                               AxisDirection::AxisZ, *logger);
2811     BOOST_REQUIRE(composite);
2812
2813     auto grid12 =
2814         dynamic_cast<CompositePortalLink&>(*composite).makeGrid(gctx, *logger);
2815     BOOST_REQUIRE(grid12);
2816
2817     BOOST_CHECK_EQUAL(
2818         grid12->resolveVolume(gctx, Vector2{40_mm, -40_mm}).value(),
2819         vol1.get());
2820
2821     BOOST_CHECK_EQUAL(
2822         grid12->resolveVolume(gctx, Vector2{40_mm, 40_mm}).value(), vol2.get());
2823
2824     composite = PortalLinkBase::merge(std::move(composite), std::move(trivial3),
2825                                       AxisDirection::AxisZ, *logger);
2826     BOOST_REQUIRE(composite);
2827
2828     auto grid123 =
2829         dynamic_cast<CompositePortalLink&>(*composite).makeGrid(gctx, *logger);
2830     BOOST_REQUIRE(grid123);
2831
2832     BOOST_CHECK_EQUAL(
2833         grid123->resolveVolume(gctx, Vector2{40_mm, -110_mm}).value(),
2834         vol1.get());
2835
2836     BOOST_CHECK_EQUAL(
2837         grid123->resolveVolume(gctx, Vector2{40_mm, -10_mm}).value(),
2838         vol2.get());
2839
2840     BOOST_CHECK_EQUAL(
2841         grid123->resolveVolume(gctx, Vector2{40_mm, 190_mm}).value(),
2842         vol3.get());
2843   }
2844   BOOST_TEST_CONTEXT("PlaneXDirection") {
2845     auto vol1 = std::make_shared<TrackingVolume>(
2846         Transform3::Identity(),
2847         std::make_shared<CuboidVolumeBounds>(30_mm, 40_mm, 100_mm));
2848
2849     auto vol2 = std::make_shared<TrackingVolume>(
2850         Transform3::Identity() * Translation3{Vector3::UnitX() * 60},
2851         std::make_shared<CuboidVolumeBounds>(30_mm, 40_mm, 100_mm));
2852
2853     auto vol3 = std::make_shared<TrackingVolume>(
2854         Transform3::Identity() * Translation3{Vector3::UnitX() * 120},
2855         std::make_shared<CuboidVolumeBounds>(30_mm, 40_mm, 100_mm));
2856
2857     auto rBounds = std::make_shared<const RectangleBounds>(30_mm, 40_mm);
2858     auto plane1 =
2859         Surface::makeShared<PlaneSurface>(Transform3::Identity(), rBounds);
2860
2861     auto plane2 = Surface::makeShared<PlaneSurface>(
2862         Transform3::Identity() * Translation3{Vector3::UnitX() * 60}, rBounds);
2863
2864     auto plane3 = Surface::makeShared<PlaneSurface>(
2865         Transform3::Identity() * Translation3{Vector3::UnitX() * 120}, rBounds);
2866
2867     auto trivial1 = std::make_unique<TrivialPortalLink>(plane1, *vol1);
2868     BOOST_REQUIRE(trivial1);
2869     auto trivial2 = std::make_unique<TrivialPortalLink>(plane2, *vol2);
2870     BOOST_REQUIRE(trivial2);
2871     auto trivial3 = std::make_unique<TrivialPortalLink>(plane3, *vol3);
2872     BOOST_REQUIRE(trivial3);
2873
2874     auto grid1 = trivial1->makeGrid(AxisDirection::AxisX);
2875     auto compGridTrivial = PortalLinkBase::merge(
2876         std::move(grid1), std::make_unique<TrivialPortalLink>(*trivial2),
2877         AxisDirection::AxisX, *logger);
2878     BOOST_REQUIRE(compGridTrivial);
2879     BOOST_CHECK_EQUAL(dynamic_cast<CompositePortalLink&>(*compGridTrivial)
2880                           .makeGrid(gctx, *logger),
2881                       nullptr);
2882
2883     auto composite =
2884         PortalLinkBase::merge(std::move(trivial1), std::move(trivial2),
2885                               AxisDirection::AxisX, *logger);
2886     BOOST_REQUIRE(composite);
2887
2888     auto grid12 =
2889         dynamic_cast<CompositePortalLink&>(*composite).makeGrid(gctx, *logger);
2890     BOOST_REQUIRE(grid12);
2891
2892     BOOST_CHECK_EQUAL(grid12->resolveVolume(gctx, Vector2{-30_mm, 0}).value(),
2893                       vol1.get());
2894
2895     BOOST_CHECK_EQUAL(grid12->resolveVolume(gctx, Vector2{30_mm, 0}).value(),
2896                       vol2.get());
2897
2898     composite = PortalLinkBase::merge(std::move(composite), std::move(trivial3),
2899                                       AxisDirection::AxisX, *logger);
2900     BOOST_REQUIRE(composite);
2901
2902     auto grid123 =
2903         dynamic_cast<CompositePortalLink&>(*composite).makeGrid(gctx, *logger);
2904     BOOST_REQUIRE(grid123);
2905
2906     BOOST_CHECK_EQUAL(
2907         grid123->resolveVolume(gctx, Vector2{-80_mm, 0_mm}).value(),
2908         vol1.get());
2909
2910     BOOST_CHECK_EQUAL(grid123->resolveVolume(gctx, Vector2{0_mm, 0_mm}).value(),
2911                       vol2.get());
2912
2913     BOOST_CHECK_EQUAL(
2914         grid123->resolveVolume(gctx, Vector2{80_mm, 0_mm}).value(), vol3.get());
2915   }
2916   BOOST_TEST_CONTEXT("PlaneYDirection") {
2917     auto vol1 = std::make_shared<TrackingVolume>(
2918         Transform3::Identity(),
2919         std::make_shared<CuboidVolumeBounds>(30_mm, 40_mm, 100_mm));
2920
2921     auto vol2 = std::make_shared<TrackingVolume>(
2922         Transform3::Identity() * Translation3{Vector3::UnitY() * 80},
2923         std::make_shared<CuboidVolumeBounds>(30_mm, 40_mm, 100_mm));
2924
2925     auto vol3 = std::make_shared<TrackingVolume>(
2926         Transform3::Identity() * Translation3{Vector3::UnitY() * 160},
2927         std::make_shared<CuboidVolumeBounds>(30_mm, 40_mm, 100_mm));
2928
2929     auto rBounds = std::make_shared<const RectangleBounds>(30_mm, 40_mm);
2930     auto plane1 =
2931         Surface::makeShared<PlaneSurface>(Transform3::Identity(), rBounds);
2932
2933     auto plane2 = Surface::makeShared<PlaneSurface>(
2934         Transform3::Identity() * Translation3{Vector3::UnitY() * 80}, rBounds);
2935
2936     auto plane3 = Surface::makeShared<PlaneSurface>(
2937         Transform3::Identity() * Translation3{Vector3::UnitY() * 160}, rBounds);
2938
2939     auto trivial1 = std::make_unique<TrivialPortalLink>(plane1, *vol1);
2940     BOOST_REQUIRE(trivial1);
2941     auto trivial2 = std::make_unique<TrivialPortalLink>(plane2, *vol2);
2942     BOOST_REQUIRE(trivial2);
2943     auto trivial3 = std::make_unique<TrivialPortalLink>(plane3, *vol3);
2944     BOOST_REQUIRE(trivial3);
2945
2946     auto grid1 = trivial1->makeGrid(AxisDirection::AxisY);
2947     auto compGridTrivial = PortalLinkBase::merge(
2948         std::move(grid1), std::make_unique<TrivialPortalLink>(*trivial2),
2949         AxisDirection::AxisY, *logger);
2950     BOOST_REQUIRE(compGridTrivial);
2951     BOOST_CHECK_EQUAL(dynamic_cast<CompositePortalLink&>(*compGridTrivial)
2952                           .makeGrid(gctx, *logger),
2953                       nullptr);
2954
2955     auto composite =
2956         PortalLinkBase::merge(std::move(trivial1), std::move(trivial2),
2957                               AxisDirection::AxisY, *logger);
2958     BOOST_REQUIRE(composite);
2959
2960     auto grid12 =
2961         dynamic_cast<CompositePortalLink&>(*composite).makeGrid(gctx, *logger);
2962     BOOST_REQUIRE(grid12);
2963
2964     BOOST_CHECK_EQUAL(
2965         grid12->resolveVolume(gctx, Vector2{0_mm, -40_mm}).value(), vol1.get());
2966
2967     BOOST_CHECK_EQUAL(grid12->resolveVolume(gctx, Vector2{0_mm, 40_mm}).value(),
2968                       vol2.get());
2969
2970     composite = PortalLinkBase::merge(std::move(composite), std::move(trivial3),
2971                                       AxisDirection::AxisY, *logger);
2972     BOOST_REQUIRE(composite);
2973
2974     auto grid123 =
2975         dynamic_cast<CompositePortalLink&>(*composite).makeGrid(gctx, *logger);
2976     BOOST_REQUIRE(grid123);
2977
2978     BOOST_CHECK_EQUAL(
2979         grid123->resolveVolume(gctx, Vector2{0_mm, -110_mm}).value(),
2980         vol1.get());
2981
2982     BOOST_CHECK_EQUAL(
2983         grid123->resolveVolume(gctx, Vector2{0_mm, -10_mm}).value(),
2984         vol2.get());
2985
2986     BOOST_CHECK_EQUAL(
2987         grid123->resolveVolume(gctx, Vector2{0_mm, 110_mm}).value(),
2988         vol3.get());
2989   }
2990 }
2991
2992 BOOST_AUTO_TEST_CASE(TrivialGridR) {
2993   auto vol1 = std::make_shared<TrackingVolume>(
2994       Transform3::Identity(),
2995       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
2996
2997   auto vol2 = std::make_shared<TrackingVolume>(
2998       Transform3::Identity(),
2999       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
3000
3001   auto disc1 =
3002       Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 60_mm);
3003
3004   auto disc2 =
3005       Surface::makeShared<DiscSurface>(Transform3::Identity(), 60_mm, 90_mm);
3006
3007   auto trivial = std::make_unique<TrivialPortalLink>(disc2, *vol2);
3008   BOOST_REQUIRE(trivial);
3009
3010   auto gridPhi = GridPortalLink::make(
3011       disc1, AxisDirection::AxisPhi,
3012       Axis{AxisClosed, -std::numbers::pi, std::numbers::pi, 2});
3013   gridPhi->setVolume(vol1.get());
3014
3015   auto gridR = GridPortalLink::make(disc1, AxisDirection::AxisR,
3016                                     Axis{AxisBound, 30_mm, 60_mm, 2});
3017   gridR->setVolume(vol1.get());
3018
3019   BOOST_TEST_CONTEXT("Colinear") {
3020     auto merged = PortalLinkBase::merge(copy(trivial), copy(gridR),
3021                                         AxisDirection::AxisR, *logger);
3022     BOOST_REQUIRE(merged);
3023     BOOST_CHECK_NE(dynamic_cast<CompositePortalLink*>(merged.get()), nullptr);
3024   }
3025
3026   BOOST_TEST_CONTEXT("Orthogonal") {
3027     auto merged = PortalLinkBase::merge(copy(gridPhi), copy(trivial),
3028                                         AxisDirection::AxisR, *logger);
3029     BOOST_REQUIRE(merged);
3030     BOOST_CHECK_NE(dynamic_cast<CompositePortalLink*>(merged.get()), nullptr);
3031   }
3032 }
3033
3034 BOOST_AUTO_TEST_CASE(TrivialGridPhi) {
3035   auto vol1 = std::make_shared<TrackingVolume>(
3036       Transform3::Identity(),
3037       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
3038
3039   auto vol2 = std::make_shared<TrackingVolume>(
3040       Transform3::Identity(),
3041       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
3042
3043   auto disc1 = Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm,
3044                                                 100_mm, 30_degree);
3045
3046   auto disc2 = Surface::makeShared<DiscSurface>(
3047       Transform3{AngleAxis3{90_degree, Vector3::UnitZ()}}, 30_mm, 100_mm,
3048       60_degree);
3049
3050   auto trivial = std::make_unique<TrivialPortalLink>(disc2, *vol2);
3051   BOOST_REQUIRE(trivial);
3052
3053   auto gridPhi = GridPortalLink::make(
3054       disc1, AxisDirection::AxisPhi, Axis{AxisBound, -30_degree, 30_degree, 2});
3055   gridPhi->setVolume(vol1.get());
3056
3057   auto gridR = GridPortalLink::make(disc1, AxisDirection::AxisR,
3058                                     Axis{AxisBound, 30_mm, 100_mm, 2});
3059   gridR->setVolume(vol1.get());
3060
3061   BOOST_TEST_CONTEXT("Colinear") {
3062     auto merged = PortalLinkBase::merge(copy(trivial), copy(gridPhi),
3063                                         AxisDirection::AxisPhi, *logger);
3064     BOOST_REQUIRE(merged);
3065     BOOST_CHECK_NE(dynamic_cast<CompositePortalLink*>(merged.get()), nullptr);
3066   }
3067
3068   BOOST_TEST_CONTEXT("Orthogonal") {
3069     auto merged = PortalLinkBase::merge(copy(gridR), copy(trivial),
3070                                         AxisDirection::AxisPhi, *logger);
3071     BOOST_REQUIRE(merged);
3072     BOOST_CHECK_NE(dynamic_cast<CompositePortalLink*>(merged.get()), nullptr);
3073   }
3074 }
3075
3076 BOOST_AUTO_TEST_CASE(TrivialGridXY) {
3077   auto vol1 = std::make_shared<TrackingVolume>(
3078       Transform3::Identity(),
3079       std::make_shared<CuboidVolumeBounds>(30_mm, 30_mm, 30_mm));
3080
3081   auto vol2 = std::make_shared<TrackingVolume>(
3082       Transform3::Identity() * Translation3{Vector3::UnitX() * 60},
3083       std::make_shared<CuboidVolumeBounds>(30_mm, 30_mm, 30_mm));
3084
3085   auto rBounds = std::make_shared<const RectangleBounds>(30_mm, 30_mm);
3086   auto plane1 =
3087       Surface::makeShared<PlaneSurface>(Transform3::Identity(), rBounds);
3088
3089   auto plane2 = Surface::makeShared<PlaneSurface>(
3090       Transform3::Identity() * Translation3{Vector3::UnitX() * 60}, rBounds);
3091
3092   auto gridX = GridPortalLink::make(plane1, AxisDirection::AxisX,
3093                                     Axis{AxisBound, -30_mm, 30_mm, 2});
3094   gridX->setVolume(vol1.get());
3095
3096   auto gridY = GridPortalLink::make(plane1, AxisDirection::AxisY,
3097                                     Axis{AxisBound, -30_mm, 30_mm, 2});
3098   gridY->setVolume(vol1.get());
3099
3100   auto trivial = std::make_unique<TrivialPortalLink>(plane2, *vol2);
3101   BOOST_REQUIRE(trivial);
3102
3103   BOOST_TEST_CONTEXT("Colinear") {
3104     auto merged = PortalLinkBase::merge(copy(trivial), copy(gridX),
3105                                         AxisDirection::AxisX, *logger);
3106     BOOST_REQUIRE(merged);
3107     BOOST_CHECK_NE(dynamic_cast<CompositePortalLink*>(merged.get()), nullptr);
3108   }
3109
3110   BOOST_TEST_CONTEXT("Orthogonal") {
3111     auto merged = PortalLinkBase::merge(copy(gridY), copy(trivial),
3112                                         AxisDirection::AxisX, *logger);
3113     BOOST_REQUIRE(merged);
3114     BOOST_CHECK_NE(dynamic_cast<CompositePortalLink*>(merged.get()), nullptr);
3115   }
3116 }
3117
3118 BOOST_AUTO_TEST_CASE(CompositeOther) {
3119   auto vol1 = std::make_shared<TrackingVolume>(
3120       Transform3::Identity(),
3121       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
3122
3123   auto vol2 = std::make_shared<TrackingVolume>(
3124       Transform3::Identity(),
3125       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
3126
3127   auto vol3 = std::make_shared<TrackingVolume>(
3128       Transform3::Identity(),
3129       std::make_shared<CylinderVolumeBounds>(30_mm, 40_mm, 100_mm));
3130
3131   auto disc1 =
3132       Surface::makeShared<DiscSurface>(Transform3::Identity(), 30_mm, 60_mm);
3133   auto disc2 =
3134       Surface::makeShared<DiscSurface>(Transform3::Identity(), 60_mm, 90_mm);
3135   auto disc3 =
3136       Surface::makeShared<DiscSurface>(Transform3::Identity(), 90_mm, 120_mm);
3137
3138   auto grid1 = GridPortalLink::make(disc1, *vol1, AxisDirection::AxisR);
3139   auto trivial2 = std::make_unique<TrivialPortalLink>(disc2, *vol2);
3140
3141   auto composite12 = std::make_unique<CompositePortalLink>(
3142       std::move(grid1), std::move(trivial2), AxisDirection::AxisR);
3143
3144   BOOST_CHECK_EQUAL(composite12->depth(), 1);
3145   BOOST_CHECK_EQUAL(composite12->size(), 2);
3146
3147   auto trivial3 = std::make_unique<TrivialPortalLink>(disc3, *vol3);
3148
3149   auto composite123Ptr =
3150       PortalLinkBase::merge(std::move(composite12), std::move(trivial3),
3151                             AxisDirection::AxisR, *logger);
3152
3153   const auto* composite123 =
3154       dynamic_cast<const CompositePortalLink*>(composite123Ptr.get());
3155   BOOST_REQUIRE(composite123);
3156
3157   BOOST_CHECK_EQUAL(composite123->depth(), 1);
3158
3159   BOOST_CHECK_EQUAL(
3160       composite123->resolveVolume(gctx, Vector2{40_mm, 0_degree}).value(),
3161       vol1.get());
3162   BOOST_CHECK_EQUAL(
3163       composite123->resolveVolume(gctx, Vector2{70_mm, 0_degree}).value(),
3164       vol2.get());
3165   BOOST_CHECK_EQUAL(
3166       composite123->resolveVolume(gctx, Vector2{100_mm, 0_degree}).value(),
3167       vol3.get());
3168
3169   BOOST_CHECK_EQUAL(composite123->size(), 3);
3170 }
3171
3172 BOOST_AUTO_TEST_SUITE_END()  // PortalMerging
3173
3174 BOOST_AUTO_TEST_SUITE_END()  // Geometry
3175
3176 }  // namespace Acts::Test