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