Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-10-22 07:53:19

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/unit_test.hpp>
0010 
0011 #include "Acts/Material/GridSurfaceMaterial.hpp"
0012 #include "Acts/Material/GridSurfaceMaterialFactory.hpp"
0013 #include "Acts/Material/Material.hpp"
0014 #include "Acts/Material/MaterialSlab.hpp"
0015 #include "Acts/Utilities/GridAxisGenerators.hpp"
0016 #include "Acts/Utilities/ProtoAxis.hpp"
0017 #include "Acts/Utilities/VectorHelpers.hpp"
0018 
0019 #include <numbers>
0020 #include <vector>
0021 
0022 using namespace Acts;
0023 
0024 // this is a global access to the x coordinate
0025 class GlobalAccessX final : public GridAccess::IGlobalToGridLocal {
0026  public:
0027   std::array<double, 1u> g2X(const Vector3& global) const {
0028     return {global.x()};
0029   }
0030 };
0031 
0032 class LocalAccessX final : public GridAccess::IBoundToGridLocal {
0033  public:
0034   std::array<double, 1u> l2X(const Vector2& local) const { return {local.x()}; }
0035 };
0036 
0037 class GlobalAccessPhi final : public GridAccess::IGlobalToGridLocal {
0038  public:
0039   std::array<double, 1u> g2Phi(const Vector3& global) const {
0040     return {std::atan2(global.y(), global.x())};
0041   }
0042 };
0043 
0044 class LocalAccessPhi final : public GridAccess::IBoundToGridLocal {
0045  public:
0046   std::array<double, 1u> l2Phi(const Vector2& local) const {
0047     return {std::atan2(local.y(), local.x())};
0048   }
0049 };
0050 
0051 class GlobalAccessXY final : public GridAccess::IGlobalToGridLocal {
0052  public:
0053   std::array<double, 2u> g2XY(const Vector3& global) const {
0054     return {global.x(), global.y()};
0055   }
0056 };
0057 
0058 class LocalAccessXY final : public GridAccess::IBoundToGridLocal {
0059  public:
0060   std::array<double, 2u> l2XY(const Vector2& local) const {
0061     return {local.x(), local.y()};
0062   }
0063 };
0064 
0065 class GlobalToZPhi final : public GridAccess::IGlobalToGridLocal {
0066  public:
0067   double zShift = 0.;
0068 
0069   explicit GlobalToZPhi(double shift) : zShift(shift) {}
0070 
0071   std::array<double, 2u> g2ZPhi(const Vector3& global) const {
0072     return {global.z() + zShift, VectorHelpers::phi(global)};
0073   }
0074 };
0075 
0076 // Local on cylinder surface is rPhi, z
0077 class LocalToZPhi final : public GridAccess::IBoundToGridLocal {
0078  public:
0079   double radius = 1.;
0080 
0081   explicit LocalToZPhi(double r) : radius(r) {}
0082 
0083   std::array<double, 2u> l2ZPhi(const Vector2& local) const {
0084     return {local[1u], local[0u] / radius};
0085   }
0086 };
0087 
0088 namespace ActsTests {
0089 
0090 BOOST_AUTO_TEST_SUITE(MaterialSuite)
0091 
0092 // This test covers some wrongly configured cases
0093 BOOST_AUTO_TEST_CASE(GridIndexedMaterial_invalid_bound2Grid_Unconnected) {
0094   std::vector<MaterialSlab> material;
0095 
0096   using EqBound = GridAxisGenerators::EqBound;
0097   using EqGrid = EqBound::grid_type<std::size_t>;
0098 
0099   EqBound eqBound{{0., 5.}, 5};
0100   EqGrid eqGrid{eqBound()};
0101 
0102   IndexedSurfaceMaterial<EqGrid>::BoundToGridLocalDelegate bToX;
0103 
0104   auto globalX = std::make_unique<const GlobalAccessX>();
0105   IndexedSurfaceMaterial<EqGrid>::GlobalToGridLocalDelegate gToX;
0106   gToX.connect<&GlobalAccessX::g2X>(std::move(globalX));
0107 
0108   BOOST_CHECK_THROW(
0109       auto ism = IndexedSurfaceMaterial<EqGrid>(
0110           std::move(eqGrid), IndexedMaterialAccessor{std::move(material)},
0111           std::move(bToX), std::move(gToX)),
0112       std::invalid_argument);
0113 }
0114 
0115 // This test covers some wrongly configured cases
0116 BOOST_AUTO_TEST_CASE(GridIndexedMaterial_invalid_global2Grid_Unconnected) {
0117   std::vector<MaterialSlab> material;
0118 
0119   using EqBound = GridAxisGenerators::EqBound;
0120   using EqGrid = EqBound::grid_type<std::size_t>;
0121 
0122   EqBound eqBound{{0., 5.}, 5};
0123   EqGrid eqGrid{eqBound()};
0124 
0125   auto localX = std::make_unique<const LocalAccessX>();
0126   IndexedSurfaceMaterial<EqGrid>::BoundToGridLocalDelegate bToX;
0127   bToX.connect<&LocalAccessX::l2X>(std::move(localX));
0128 
0129   IndexedSurfaceMaterial<EqGrid>::GlobalToGridLocalDelegate gToX;
0130 
0131   BOOST_CHECK_THROW(
0132       auto ism = IndexedSurfaceMaterial<EqGrid>(
0133           std::move(eqGrid), IndexedMaterialAccessor{std::move(material)},
0134           std::move(bToX), std::move(gToX)),
0135       std::invalid_argument);
0136 }
0137 
0138 // This test covers the locally indexed grid material in 1D
0139 BOOST_AUTO_TEST_CASE(GridMaterial1D) {
0140   std::vector<MaterialSlab> material;
0141   material.emplace_back(Material::Vacuum(), 0.0);  // vacuum
0142   material.emplace_back(Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0),
0143                         1.0);
0144   material.emplace_back(
0145       Material::fromMolarDensity(11.0, 12.0, 13.0, 14.0, 15.0), 2.0);
0146   material.emplace_back(
0147       Material::fromMolarDensity(21.0, 22.0, 23.0, 24.0, 25.0), 3.0);
0148   material.emplace_back(
0149       Material::fromMolarDensity(31.0, 32.0, 33.0, 34.0, 35.0), 4.0);
0150 
0151   // Bound, equidistant axis
0152   ProtoAxis pAxisX(AxisBoundaryType::Bound, 0.0, 5.0, 5);
0153 
0154   auto localX = std::make_unique<const LocalAccessX>();
0155   GridAccess::BoundToGridLocal1DimDelegate bToX;
0156   bToX.connect<&LocalAccessX::l2X>(std::move(localX));
0157 
0158   auto globalX = std::make_unique<const GlobalAccessX>();
0159   GridAccess::GlobalToGridLocal1DimDelegate gToX;
0160   gToX.connect<&GlobalAccessX::g2X>(std::move(globalX));
0161 
0162   auto ismX = GridSurfaceMaterialFactory::create(pAxisX, GridMaterialAccessor{},
0163                                                  std::move(bToX),
0164                                                  std::move(gToX), material);
0165 
0166   BOOST_CHECK(ismX != nullptr);
0167 
0168   // Global access test
0169   Vector3 g0(0.5, 0., 0.);
0170   Vector3 g1(1.5, 0., 0.);
0171   Vector3 g2(2.5, 0., 0.);
0172   Vector3 g3(3.5, 0., 0.);
0173   Vector3 g4(4.5, 0., 0.);
0174 
0175   const MaterialSlab& mg0 = ismX->materialSlab(g0);
0176   const MaterialSlab& mg1 = ismX->materialSlab(g1);
0177   const MaterialSlab& mg2 = ismX->materialSlab(g2);
0178   const MaterialSlab& mg3 = ismX->materialSlab(g3);
0179   const MaterialSlab& mg4 = ismX->materialSlab(g4);
0180 
0181   BOOST_CHECK(mg0.material().isVacuum());
0182   BOOST_CHECK_EQUAL(mg1.material().X0(), 1.);
0183   BOOST_CHECK_EQUAL(mg2.material().X0(), 11.);
0184   BOOST_CHECK_EQUAL(mg3.material().X0(), 21.);
0185   BOOST_CHECK_EQUAL(mg4.material().X0(), 31.);
0186 
0187   // Try the same with Closed access
0188   // Bound, equidistant axis
0189   ProtoAxis pAxisPhi(AxisBoundaryType::Closed, -std::numbers::pi,
0190                      std::numbers::pi, 8);
0191 
0192   auto localPhi = std::make_unique<const LocalAccessPhi>();
0193   GridAccess::BoundToGridLocal1DimDelegate bToPhi;
0194   bToPhi.connect<&LocalAccessPhi::l2Phi>(std::move(localPhi));
0195 
0196   auto globalPhi = std::make_unique<const GlobalAccessPhi>();
0197   GridAccess::GlobalToGridLocal1DimDelegate gToPhi;
0198   gToPhi.connect<&GlobalAccessPhi::g2Phi>(std::move(globalPhi));
0199 
0200   std::vector<MaterialSlab> materialPhi;
0201   materialPhi.emplace_back(Material::Vacuum(), 0.0);  // vacuum
0202   materialPhi.emplace_back(Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0),
0203                            1.0);
0204   materialPhi.emplace_back(Material::Vacuum(), 0.0);  // vacuum
0205   materialPhi.emplace_back(
0206       Material::fromMolarDensity(11.0, 12.0, 13.0, 14.0, 15.0), 2.0);
0207   materialPhi.emplace_back(Material::Vacuum(), 0.0);  // vacuum
0208   materialPhi.emplace_back(
0209       Material::fromMolarDensity(21.0, 22.0, 23.0, 24.0, 25.0), 3.0);
0210   materialPhi.emplace_back(Material::Vacuum(), 0.0);  // vacuum
0211   materialPhi.emplace_back(
0212       Material::fromMolarDensity(31.0, 32.0, 33.0, 34.0, 35.0), 4.0);
0213 
0214   auto ismPhi = GridSurfaceMaterialFactory::create(
0215       pAxisPhi, GridMaterialAccessor{}, std::move(bToPhi), std::move(gToPhi),
0216       materialPhi);
0217 
0218   BOOST_CHECK(ismPhi != nullptr);
0219 
0220   for (std::size_t i = 0; i < 8; ++i) {
0221     double alpha = -std::numbers::pi + (i + 0.5) * std::numbers::pi / 4.;
0222     Vector2 query{std::cos(alpha), std::sin(alpha)};
0223     const MaterialSlab& m = ismPhi->materialSlab(query);
0224     if (i % 2 == 0) {
0225       BOOST_CHECK(m.material().isVacuum());
0226     } else {
0227       BOOST_CHECK_EQUAL(m.material().X0(), materialPhi[i].material().X0());
0228     }
0229   }
0230 }
0231 
0232 // This test covers the locally indexed grid material in 1D
0233 BOOST_AUTO_TEST_CASE(GridMaterial2D) {
0234   std::vector<std::vector<MaterialSlab>> material2x3;
0235   // This is a material matrix 2 bins in x and 3 bins in y
0236   std::vector<MaterialSlab> materialRow0;
0237   materialRow0.emplace_back(Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0),
0238                             1.0);
0239   materialRow0.emplace_back(
0240       Material::fromMolarDensity(11.0, 12.0, 13.0, 14.0, 15.0), 2.0);
0241   materialRow0.emplace_back(
0242       Material::fromMolarDensity(21.0, 22.0, 23.0, 24.0, 25.0), 3.0);
0243   std::vector<MaterialSlab> materialRow1;
0244   materialRow1.emplace_back(Material::fromMolarDensity(2.0, 2.0, 3.0, 4.0, 5.0),
0245                             1.0);
0246   materialRow1.emplace_back(
0247       Material::fromMolarDensity(12.0, 12.0, 13.0, 14.0, 15.0), 2.0);
0248   materialRow1.emplace_back(
0249       Material::fromMolarDensity(22.0, 22.0, 23.0, 24.0, 25.0), 3.0);
0250   // This gives a row major matrix
0251   material2x3.push_back(std::move(materialRow0));
0252   material2x3.push_back(std::move(materialRow1));
0253 
0254   BOOST_CHECK(material2x3[0][0].material().X0() == 1.);
0255   BOOST_CHECK(material2x3[0][1].material().X0() == 11.);
0256   BOOST_CHECK(material2x3[0][2].material().X0() == 21.);
0257   BOOST_CHECK(material2x3[1][0].material().X0() == 2.);
0258   BOOST_CHECK(material2x3[1][1].material().X0() == 12.);
0259   BOOST_CHECK(material2x3[1][2].material().X0() == 22.);
0260 
0261   ProtoAxis pAxisX(AxisBoundaryType::Bound, -1.0, 1.0, 2);
0262   ProtoAxis pAxisY(AxisBoundaryType::Bound, -1.5, 1.5, 3);
0263 
0264   std::vector<std::vector<MaterialSlab>> materialXY = material2x3;
0265 
0266   auto localXY = std::make_unique<const LocalAccessXY>();
0267   GridAccess::BoundToGridLocal2DimDelegate bToXY;
0268   bToXY.connect<&LocalAccessXY::l2XY>(std::move(localXY));
0269 
0270   auto globalXY = std::make_unique<const GlobalAccessXY>();
0271   GridAccess::GlobalToGridLocal2DimDelegate gToXY;
0272   gToXY.connect<&GlobalAccessXY::g2XY>(std::move(globalXY));
0273 
0274   auto ismXY = GridSurfaceMaterialFactory::create(
0275       pAxisX, pAxisY, GridMaterialAccessor{}, std::move(bToXY),
0276       std::move(gToXY), materialXY);
0277 
0278   BOOST_CHECK(ismXY != nullptr);
0279 
0280   // Global access test
0281   Vector3 g00(-0.5, -1.5, 0.);
0282   Vector3 g01(-0.5, 0., 0.);
0283   Vector3 g02(-0.5, 1.5, 0.);
0284   Vector3 g10(0.5, -1.5, 0.);
0285   Vector3 g11(0.5, 0., 0.);
0286   Vector3 g12(0.5, 1.5, 0.);
0287 
0288   const MaterialSlab& mg00 = ismXY->materialSlab(g00);
0289   const MaterialSlab& mg01 = ismXY->materialSlab(g01);
0290   const MaterialSlab& mg02 = ismXY->materialSlab(g02);
0291   const MaterialSlab& mg10 = ismXY->materialSlab(g10);
0292   const MaterialSlab& mg11 = ismXY->materialSlab(g11);
0293   const MaterialSlab& mg12 = ismXY->materialSlab(g12);
0294 
0295   BOOST_CHECK_EQUAL(mg00.material().X0(), 1.);
0296   BOOST_CHECK_EQUAL(mg01.material().X0(), 11.);
0297   BOOST_CHECK_EQUAL(mg02.material().X0(), 21.);
0298   BOOST_CHECK_EQUAL(mg10.material().X0(), 2.);
0299   BOOST_CHECK_EQUAL(mg11.material().X0(), 12.);
0300   BOOST_CHECK_EQUAL(mg12.material().X0(), 22.);
0301 
0302   // Let's try a ZPhi model as well
0303   auto materialZPhi = material2x3;
0304 
0305   ProtoAxis pAxisZ(AxisBoundaryType::Bound, -1.0, 1.0, 2);
0306   ProtoAxis pAxisPhi(AxisBoundaryType::Closed, -std::numbers::pi,
0307                      std::numbers::pi, 3);
0308 
0309   auto localZPhi = std::make_unique<const LocalToZPhi>(1.);
0310   GridAccess::BoundToGridLocal2DimDelegate bToZPhi;
0311   bToZPhi.connect<&LocalToZPhi::l2ZPhi>(std::move(localZPhi));
0312 
0313   auto globalZPhi = std::make_unique<const GlobalToZPhi>(0.);
0314   GridAccess::GlobalToGridLocal2DimDelegate gToZPhi;
0315   gToZPhi.connect<&GlobalToZPhi::g2ZPhi>(std::move(globalZPhi));
0316 
0317   auto ismZPhi = GridSurfaceMaterialFactory::create(
0318       pAxisZ, pAxisPhi, GridMaterialAccessor{}, std::move(bToZPhi),
0319       std::move(gToZPhi), materialZPhi);
0320 
0321   BOOST_CHECK(ismZPhi != nullptr);
0322 
0323   // Local access test - trick here is also to
0324   // see BoundtoGridLocal switches from r*phi, z -> phi,z
0325   Vector2 l00(-0.5 * std::numbers::pi, -0.5);
0326   Vector2 l01(0., -0.5);
0327   Vector2 l02(0.5 * std::numbers::pi, -0.5);
0328   Vector2 l10(-0.5 * std::numbers::pi, 0.5);
0329   Vector2 l11(0., 0.5);
0330   Vector2 l12(0.5 * std::numbers::pi, 0.5);
0331 
0332   const MaterialSlab& ml00 = ismZPhi->materialSlab(l00);
0333   const MaterialSlab& ml01 = ismZPhi->materialSlab(l01);
0334   const MaterialSlab& ml02 = ismZPhi->materialSlab(l02);
0335   const MaterialSlab& ml10 = ismZPhi->materialSlab(l10);
0336   const MaterialSlab& ml11 = ismZPhi->materialSlab(l11);
0337   const MaterialSlab& ml12 = ismZPhi->materialSlab(l12);
0338 
0339   BOOST_CHECK_EQUAL(ml00.material().X0(), 1.);
0340   BOOST_CHECK_EQUAL(ml01.material().X0(), 11.);
0341   BOOST_CHECK_EQUAL(ml02.material().X0(), 21.);
0342   BOOST_CHECK_EQUAL(ml10.material().X0(), 2.);
0343   BOOST_CHECK_EQUAL(ml11.material().X0(), 12.);
0344   BOOST_CHECK_EQUAL(ml12.material().X0(), 22.);
0345 
0346   // Test the closed character of the phi axis
0347   Vector2 l03(1.05 * std::numbers::pi, -0.5);
0348   const MaterialSlab& ml03 = ismZPhi->materialSlab(l03);
0349   BOOST_CHECK(ml03.material().X0() == ml00.material().X0());
0350 }
0351 
0352 // This test covers the locally indexed grid material in 1D
0353 BOOST_AUTO_TEST_CASE(GridIndexedMaterial1D) {
0354   std::vector<MaterialSlab> material;
0355   material.emplace_back(Material::Vacuum(), 0.0);  // vacuum
0356   material.emplace_back(Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0),
0357                         1.0);
0358   material.emplace_back(
0359       Material::fromMolarDensity(11.0, 12.0, 13.0, 14.0, 15.0), 2.0);
0360   material.emplace_back(
0361       Material::fromMolarDensity(21.0, 22.0, 23.0, 24.0, 25.0), 3.0);
0362 
0363   using EqBound = GridAxisGenerators::EqBound;
0364   using EqGrid = EqBound::grid_type<std::size_t>;
0365   using Point = EqGrid::point_t;
0366 
0367   EqBound eqBound{{0., 5.}, 5};
0368   EqGrid eqGrid{eqBound()};
0369 
0370   eqGrid.atPosition(Point{0.5}) = 1u;  // material 1
0371   eqGrid.atPosition(Point{1.5}) = 0u;  // vacuum
0372   eqGrid.atPosition(Point{2.5}) = 2u;  // material 2
0373   eqGrid.atPosition(Point{3.5}) = 2u;  // material 2
0374   eqGrid.atPosition(Point{4.5}) = 3u;  // material 3
0375 
0376   auto localX = std::make_unique<const LocalAccessX>();
0377   IndexedSurfaceMaterial<EqGrid>::BoundToGridLocalDelegate bToX;
0378   bToX.connect<&LocalAccessX::l2X>(std::move(localX));
0379 
0380   auto globalX = std::make_unique<const GlobalAccessX>();
0381   IndexedSurfaceMaterial<EqGrid>::GlobalToGridLocalDelegate gToX;
0382   gToX.connect<&GlobalAccessX::g2X>(std::move(globalX));
0383 
0384   IndexedSurfaceMaterial<EqGrid> ism(
0385       std::move(eqGrid), IndexedMaterialAccessor{std::move(material)},
0386       std::move(bToX), std::move(gToX));
0387 
0388   // Global access test
0389   Vector3 g0(0.5, 0., 0.);
0390   Vector3 g1(1.5, 0., 0.);
0391   Vector3 g2(2.5, 0., 0.);
0392   Vector3 g3(3.5, 0., 0.);
0393   Vector3 g4(4.5, 0., 0.);
0394 
0395   const MaterialSlab& mg0 = ism.materialSlab(g0);
0396   const MaterialSlab& mg1 = ism.materialSlab(g1);
0397   const MaterialSlab& mg2 = ism.materialSlab(g2);
0398   const MaterialSlab& mg3 = ism.materialSlab(g3);
0399   const MaterialSlab& mg4 = ism.materialSlab(g4);
0400 
0401   BOOST_CHECK_EQUAL(mg0.material().X0(), 1.);
0402   BOOST_CHECK(mg1.material().isVacuum());
0403   BOOST_CHECK_EQUAL(mg2.material().X0(), 11.);
0404   BOOST_CHECK_EQUAL(mg3.material().X0(), 11.);
0405   BOOST_CHECK_EQUAL(mg4.material().X0(), 21.);
0406 
0407   // Local access test
0408   Vector2 l0(0.5, 0.);
0409   Vector2 l1(1.5, 0.);
0410   Vector2 l2(2.5, 0.);
0411   Vector2 l3(3.5, 0.);
0412   Vector2 l4(4.5, 0.);
0413 
0414   const MaterialSlab& ml0 = ism.materialSlab(l0);
0415   const MaterialSlab& ml1 = ism.materialSlab(l1);
0416   const MaterialSlab& ml2 = ism.materialSlab(l2);
0417   const MaterialSlab& ml3 = ism.materialSlab(l3);
0418   const MaterialSlab& ml4 = ism.materialSlab(l4);
0419 
0420   BOOST_CHECK_EQUAL(ml0.material().X0(), 1.);
0421   BOOST_CHECK(ml1.material().isVacuum());
0422   BOOST_CHECK_EQUAL(ml2.material().X0(), 11.);
0423   BOOST_CHECK_EQUAL(ml3.material().X0(), 11.);
0424   BOOST_CHECK_EQUAL(ml4.material().X0(), 21.);
0425 
0426   // Now scale it - and access again
0427   ism.scale(2.);
0428   const MaterialSlab& sml0 = ism.materialSlab(l0);
0429   const MaterialSlab& sml1 = ism.materialSlab(l1);
0430   const MaterialSlab& sml2 = ism.materialSlab(l2);
0431   const MaterialSlab& sml3 = ism.materialSlab(l3);
0432   const MaterialSlab& sml4 = ism.materialSlab(l4);
0433 
0434   BOOST_CHECK_EQUAL(sml0.thickness(), 2.);
0435   BOOST_CHECK(sml1.material().isVacuum());
0436   BOOST_CHECK_EQUAL(sml2.thickness(), 4.);
0437   BOOST_CHECK_EQUAL(sml3.thickness(), 4.);
0438   BOOST_CHECK_EQUAL(sml4.thickness(), 6.);
0439 
0440   // Now test with the protoAxis creation method
0441 
0442   std::vector<MaterialSlab> materialStorage;
0443   materialStorage.emplace_back(Material::Vacuum(), 0.0);  // vacuum
0444   materialStorage.emplace_back(
0445       Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0), 1.0);
0446   materialStorage.emplace_back(
0447       Material::fromMolarDensity(11.0, 12.0, 13.0, 14.0, 15.0), 2.0);
0448   materialStorage.emplace_back(
0449       Material::fromMolarDensity(21.0, 22.0, 23.0, 24.0, 25.0), 3.0);
0450 
0451   std::vector<std::size_t> indexPayload = {0u, 1u, 2u, 3u, 0u, 3u, 2u, 1u, 0u};
0452 
0453   auto indexedAccessor = IndexedMaterialAccessor{std::move(materialStorage)};
0454 
0455   // An X proto axis
0456   ProtoAxis pAxisX(AxisBoundaryType::Bound, 0.0, 9.0, 9);
0457 
0458   auto localXidx = std::make_unique<const LocalAccessX>();
0459   IndexedSurfaceMaterial<EqGrid>::BoundToGridLocalDelegate bToXidx;
0460   bToXidx.connect<&LocalAccessX::l2X>(std::move(localXidx));
0461 
0462   auto globalXidx = std::make_unique<const GlobalAccessX>();
0463   IndexedSurfaceMaterial<EqGrid>::GlobalToGridLocalDelegate gToXidx;
0464   gToXidx.connect<&GlobalAccessX::g2X>(std::move(globalXidx));
0465 
0466   auto ismXidx = GridSurfaceMaterialFactory::create(
0467       pAxisX, std::move(indexedAccessor), std::move(bToXidx),
0468       std::move(gToXidx), indexPayload);
0469 
0470   // Check construction
0471   BOOST_CHECK(ismXidx != nullptr);
0472   // The vacuum (==0) indexed entries
0473   BOOST_CHECK(ismXidx->materialSlab(Vector3{0.5, 0., 0.}).isVacuum());
0474   BOOST_CHECK(ismXidx->materialSlab(Vector3{4.5, 0., 0.}).isVacuum());
0475   BOOST_CHECK(ismXidx->materialSlab(Vector3{8.5, 0., 0.}).isVacuum());
0476   // The material 1 (==1) indexed entries
0477   BOOST_CHECK_EQUAL(ismXidx->materialSlab(Vector3{1.5, 0., 0.}).material().X0(),
0478                     1.);
0479   BOOST_CHECK_EQUAL(ismXidx->materialSlab(Vector3{7.5, 0., 0.}).material().X0(),
0480                     1.);
0481   // The material 2 (==2) indexed entries
0482   BOOST_CHECK_EQUAL(ismXidx->materialSlab(Vector3{2.5, 0., 0.}).material().X0(),
0483                     11.);
0484   BOOST_CHECK_EQUAL(ismXidx->materialSlab(Vector3{6.5, 0., 0.}).material().X0(),
0485                     11.);
0486   // The material 3 (==3) indexed entries
0487   BOOST_CHECK_EQUAL(ismXidx->materialSlab(Vector3{3.5, 0., 0.}).material().X0(),
0488                     21.);
0489   BOOST_CHECK_EQUAL(ismXidx->materialSlab(Vector3{5.5, 0., 0.}).material().X0(),
0490                     21.);
0491 }
0492 
0493 // This test covers the locally indexed grid material in 2D
0494 BOOST_AUTO_TEST_CASE(GridIndexedMaterial2D) {
0495   std::vector<MaterialSlab> material;
0496   material.emplace_back(Material::Vacuum(), 1.0);  // vacuum
0497   material.emplace_back(Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0),
0498                         1.0);
0499   material.emplace_back(
0500       Material::fromMolarDensity(11.0, 12.0, 13.0, 14.0, 15.0), 1.0);
0501   material.emplace_back(
0502       Material::fromMolarDensity(21.0, 22.0, 23.0, 24.0, 25.0), 1.0);
0503 
0504   //  Test (1) with explicit grid creation
0505   std::vector<MaterialSlab> materialT1 = material;
0506   using EqBoundEqClosed = GridAxisGenerators::EqBoundEqClosed;
0507   using EqEqGrid = EqBoundEqClosed::grid_type<std::size_t>;
0508   using Point = EqEqGrid::point_t;
0509 
0510   // 2 bins in z, 4 bins in phi
0511   EqBoundEqClosed eqeqBound{
0512       {-1., 1.}, 2, {-std::numbers::pi, std::numbers::pi}, 4};
0513   EqEqGrid eqeqGrid{eqeqBound()};
0514 
0515   eqeqGrid.atPosition(Point{-0.5, -std::numbers::pi * 0.75}) =
0516       1u;                                                          // material 1
0517   eqeqGrid.atPosition(Point{-0.5, -std::numbers::pi / 4.}) = 1u;   // material 1
0518   eqeqGrid.atPosition(Point{-0.5, std::numbers::pi / 4.}) = 0u;    // vacuum
0519   eqeqGrid.atPosition(Point{-0.5, std::numbers::pi * 0.75}) = 2u;  // material 2
0520 
0521   eqeqGrid.atPosition(Point{0.5, -std::numbers::pi * 0.75}) = 0u;  // vacuum
0522   eqeqGrid.atPosition(Point{0.5, -std::numbers::pi / 4.}) = 3u;    // material 3
0523   eqeqGrid.atPosition(Point{0.5, std::numbers::pi / 4.}) = 3u;     // material 3
0524   eqeqGrid.atPosition(Point{0.5, std::numbers::pi * 0.75}) = 0u;   // vacuum
0525 
0526   // With radius 20
0527   auto boundToGridT1 = std::make_unique<const LocalToZPhi>(20.);
0528   IndexedSurfaceMaterial<EqEqGrid>::BoundToGridLocalDelegate bToZPhiT1;
0529   bToZPhiT1.connect<&LocalToZPhi::l2ZPhi>(std::move(boundToGridT1));
0530 
0531   // With z shift 10
0532   auto globalToGridT1 = std::make_unique<const GlobalToZPhi>(10.);
0533   IndexedSurfaceMaterial<EqEqGrid>::GlobalToGridLocalDelegate gToZphiT1;
0534   gToZphiT1.connect<&GlobalToZPhi::g2ZPhi>(std::move(globalToGridT1));
0535 
0536   // Create the indexed material grid
0537   IndexedSurfaceMaterial<EqEqGrid> ism(
0538       std::move(eqeqGrid), IndexedMaterialAccessor{std::move(materialT1)},
0539       std::move(bToZPhiT1), std::move(gToZphiT1));
0540 
0541   // Test with proto grid greation
0542   auto materialT2 = material;
0543 
0544   auto boundToGridT2 = std::make_unique<const LocalToZPhi>(20.);
0545   GridAccess::BoundToGridLocal2DimDelegate bToZPhiT2;
0546   bToZPhiT2.connect<&LocalToZPhi::l2ZPhi>(std::move(boundToGridT2));
0547 
0548   auto globalToGridT2 = std::make_unique<const GlobalToZPhi>(10.);
0549   GridAccess::GlobalToGridLocal2DimDelegate gToZphiT2;
0550   gToZphiT2.connect<&GlobalToZPhi::g2ZPhi>(std::move(globalToGridT2));
0551 
0552   ProtoAxis pAxisZ(AxisBoundaryType::Bound, -1.0, 1.0, 2);
0553   ProtoAxis pAxisPhi(AxisBoundaryType::Closed, -std::numbers::pi,
0554                      std::numbers::pi, 4);
0555 
0556   std::vector<std::vector<std::size_t>> indexPayload = {
0557       std::vector<std::size_t>{1u, 1u, 0u, 2u},
0558       std::vector<std::size_t>{0u, 3u, 3u, 0u}};
0559 
0560   auto ismZPhi = GridSurfaceMaterialFactory::create(
0561       pAxisZ, pAxisPhi, IndexedMaterialAccessor{std::move(materialT2)},
0562       std::move(bToZPhiT2), std::move(gToZphiT2), indexPayload);
0563 
0564   // Global access test, both should give material 1
0565   Vector3 g0(-0.5, -0.5, -10.5);
0566   const MaterialSlab& mg0T1 = ism.materialSlab(g0);
0567   const MaterialSlab& mg0T2 = ismZPhi->materialSlab(g0);
0568   BOOST_CHECK_EQUAL(mg0T1.material().X0(), 1.);
0569   BOOST_CHECK_EQUAL(mg0T2.material().X0(), 1.);
0570 
0571   Vector3 g1(0.5, -0.5, -11.5);  // checking out of bound access
0572   const MaterialSlab& mg1T1 = ism.materialSlab(g1);
0573   const MaterialSlab& mg1T2 = ismZPhi->materialSlab(g1);
0574   BOOST_CHECK_EQUAL(mg1T1.material().X0(), 1.);
0575   BOOST_CHECK_EQUAL(mg1T2.material().X0(), 1.);
0576 
0577   Vector3 g2(0.5, 0.5, -10.5);
0578   const MaterialSlab& mg2T1 = ism.materialSlab(g2);
0579   const MaterialSlab& mg2T2 = ismZPhi->materialSlab(g2);
0580   BOOST_CHECK(mg2T1.material().isVacuum());  // vacuum
0581   BOOST_CHECK(mg2T2.material().isVacuum());  // vacuum
0582 
0583   Vector3 g3(0.5, 0.5,
0584              -9.5);  // should be material 3, same phi but different z
0585   const MaterialSlab& mg3T1 = ism.materialSlab(g3);
0586   const MaterialSlab& mg3T2 = ismZPhi->materialSlab(g3);
0587   BOOST_CHECK_EQUAL(mg3T1.material().X0(), 21.);
0588   BOOST_CHECK_EQUAL(mg3T2.material().X0(), 21.);
0589 }
0590 
0591 // This test covers the globally indexed grid material with non-shared material
0592 BOOST_AUTO_TEST_CASE(GridGloballyIndexedMaterialNonShared) {
0593   auto material = std::make_shared<std::vector<MaterialSlab>>();
0594 
0595   material->emplace_back(Material::Vacuum(), 0.0);  // vacuum
0596   material->emplace_back(Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0),
0597                          1.0);
0598   material->emplace_back(
0599       Material::fromMolarDensity(11.0, 12.0, 13.0, 14.0, 15.0), 2.0);
0600   material->emplace_back(
0601       Material::fromMolarDensity(21.0, 22.0, 23.0, 24.0, 25.0), 3.0);
0602   material->emplace_back(
0603       Material::fromMolarDensity(31.0, 22.0, 23.0, 24.0, 25.0), 4.0);
0604 
0605   using EqBound = GridAxisGenerators::EqBound;
0606   using EqGrid = EqBound::grid_type<std::size_t>;
0607   using Point = EqGrid::point_t;
0608 
0609   EqBound eqBound{{0., 5.}, 5};
0610   EqGrid eqGrid{eqBound()};
0611 
0612   eqGrid.atPosition(Point{0.5}) = 1u;  // material 1
0613   eqGrid.atPosition(Point{1.5}) = 0u;  // vacuum
0614   eqGrid.atPosition(Point{2.5}) = 2u;  // material 2
0615   eqGrid.atPosition(Point{3.5}) = 2u;  // material 2
0616   eqGrid.atPosition(Point{4.5}) = 3u;  // material 3
0617 
0618   auto localX = std::make_unique<const LocalAccessX>();
0619   IndexedSurfaceMaterial<EqGrid>::BoundToGridLocalDelegate bToX;
0620   bToX.connect<&LocalAccessX::l2X>(std::move(localX));
0621 
0622   auto globalX = std::make_unique<const GlobalAccessX>();
0623   IndexedSurfaceMaterial<EqGrid>::GlobalToGridLocalDelegate gToX;
0624   gToX.connect<&GlobalAccessX::g2X>(std::move(globalX));
0625 
0626   GloballyIndexedSurfaceMaterial<EqGrid> ism(
0627       std::move(eqGrid), GloballyIndexedMaterialAccessor{material, false},
0628       std::move(bToX), std::move(gToX));
0629 
0630   // Local access test
0631   Vector2 l0(0.5, 0.);
0632   Vector2 l1(1.5, 0.);
0633   Vector2 l2(2.5, 0.);
0634   Vector2 l3(3.5, 0.);
0635   Vector2 l4(4.5, 0.);
0636 
0637   const MaterialSlab& ml0 = ism.materialSlab(l0);
0638   const MaterialSlab& ml1 = ism.materialSlab(l1);
0639   const MaterialSlab& ml2 = ism.materialSlab(l2);
0640   const MaterialSlab& ml3 = ism.materialSlab(l3);
0641   const MaterialSlab& ml4 = ism.materialSlab(l4);
0642 
0643   BOOST_CHECK_EQUAL(ml0.material().X0(), 1.);
0644   BOOST_CHECK(ml1.material().isVacuum());
0645   BOOST_CHECK_EQUAL(ml2.material().X0(), 11.);
0646   BOOST_CHECK_EQUAL(ml3.material().X0(), 11.);
0647   BOOST_CHECK_EQUAL(ml4.material().X0(), 21.);
0648 
0649   EqBound eqBound1{{0., 5.}, 1};
0650   EqGrid eqGrid1{eqBound1()};
0651 
0652   auto localX1 = std::make_unique<const LocalAccessX>();
0653   IndexedSurfaceMaterial<EqGrid>::BoundToGridLocalDelegate bToX1;
0654   bToX1.connect<&LocalAccessX::l2X>(std::move(localX1));
0655 
0656   auto globalX1 = std::make_unique<const GlobalAccessX>();
0657   IndexedSurfaceMaterial<EqGrid>::GlobalToGridLocalDelegate gToX1;
0658   gToX1.connect<&GlobalAccessX::g2X>(std::move(globalX1));
0659 
0660   eqGrid1.atPosition(Point{2.5}) = 4u;  // material 4
0661 
0662   GloballyIndexedSurfaceMaterial<EqGrid> ism1(
0663       std::move(eqGrid1), GloballyIndexedMaterialAccessor{material, false},
0664       std::move(bToX1), std::move(gToX1));
0665 
0666   Vector2 l0g1(2.5, 0.);
0667   const MaterialSlab& ml0g1 = ism1.materialSlab(l0g1);
0668   BOOST_CHECK_EQUAL(ml0g1.material().X0(), 31.);
0669 
0670   // Scale
0671   ism1.scale(2.);
0672   const MaterialSlab& sml0g1 = ism1.materialSlab(l0g1);
0673   BOOST_CHECK_EQUAL(sml0g1.thickness(), 8.);
0674 
0675   // First one stays unscaled
0676   const MaterialSlab& sml0 = ism.materialSlab(l0);
0677   BOOST_CHECK_EQUAL(sml0.thickness(), 1.);
0678 }
0679 
0680 // This test covers the globally indexed grid material with shared
0681 BOOST_AUTO_TEST_CASE(GridGloballyIndexedMaterialShared) {
0682   auto material = std::make_shared<std::vector<MaterialSlab>>();
0683 
0684   material->emplace_back(Material::Vacuum(), 0.0);  // vacuum
0685   material->emplace_back(Material::fromMolarDensity(1.0, 2.0, 3.0, 4.0, 5.0),
0686                          1.0);
0687 
0688   using EqBound = GridAxisGenerators::EqBound;
0689   using EqGrid = EqBound::grid_type<std::size_t>;
0690   using Point = EqGrid::point_t;
0691 
0692   EqBound eqBound0{{0., 5.}, 1};
0693   EqGrid eqGrid0{eqBound0()};
0694 
0695   eqGrid0.atPosition(Point{2.5}) = 1u;  // material 1
0696   auto localX0 = std::make_unique<const LocalAccessX>();
0697   IndexedSurfaceMaterial<EqGrid>::BoundToGridLocalDelegate bToX0;
0698   bToX0.connect<&LocalAccessX::l2X>(std::move(localX0));
0699 
0700   auto globalX0 = std::make_unique<const GlobalAccessX>();
0701   IndexedSurfaceMaterial<EqGrid>::GlobalToGridLocalDelegate gToX0;
0702   gToX0.connect<&GlobalAccessX::g2X>(std::move(globalX0));
0703 
0704   GloballyIndexedSurfaceMaterial<EqGrid> ism0(
0705       std::move(eqGrid0), GloballyIndexedMaterialAccessor{material, true},
0706       std::move(bToX0), std::move(gToX0));
0707 
0708   EqBound eqBound1{{0., 5.}, 1};
0709   EqGrid eqGrid1{eqBound1()};
0710 
0711   eqGrid1.atPosition(Point{2.5}) = 1u;  // material 1
0712   auto localX1 = std::make_unique<const LocalAccessX>();
0713   IndexedSurfaceMaterial<EqGrid>::BoundToGridLocalDelegate bToX1;
0714   bToX1.connect<&LocalAccessX::l2X>(std::move(localX1));
0715 
0716   auto globalX1 = std::make_unique<const GlobalAccessX>();
0717   IndexedSurfaceMaterial<EqGrid>::GlobalToGridLocalDelegate gToX1;
0718   gToX1.connect<&GlobalAccessX::g2X>(std::move(globalX1));
0719 
0720   GloballyIndexedSurfaceMaterial<EqGrid> ism1(
0721       std::move(eqGrid1), GloballyIndexedMaterialAccessor{material, true},
0722       std::move(bToX1), std::move(gToX1));
0723 
0724   Vector2 l0(2.5, 0.);
0725 
0726   // check grid material 0
0727   const MaterialSlab& ml0 = ism0.materialSlab(l0);
0728   BOOST_CHECK_EQUAL(ml0.material().X0(), 1.);
0729 
0730   const MaterialSlab& ml0g1 = ism1.materialSlab(l0);
0731   BOOST_CHECK_EQUAL(ml0g1.material().X0(), 1.);
0732 
0733   // scaling shared material should throw a std::invalid_argument
0734   BOOST_CHECK_THROW(ism1.scale(2.), std::invalid_argument);
0735 }
0736 
0737 // This test covers the grid material (non-indexed accessor)
0738 //
0739 // In this setup, the material is not indexed, but filled directly
0740 // into the grid structure.
0741 BOOST_AUTO_TEST_CASE(GridSurfaceMaterialTests) {
0742   using EqBound = GridAxisGenerators::EqBound;
0743   using EqGrid = EqBound::grid_type<MaterialSlab>;
0744   using Point = EqGrid::point_t;
0745 
0746   EqBound eqBound{{0., 5.}, 5};
0747   EqGrid eqGrid{eqBound()};
0748 
0749   eqGrid.atPosition(Point{0.5}) = MaterialSlab::Vacuum(0.0);
0750   eqGrid.atPosition(Point{1.5}) = MaterialSlab::Vacuum(1.0);
0751   eqGrid.atPosition(Point{2.5}) = MaterialSlab::Vacuum(2.0);
0752   eqGrid.atPosition(Point{3.5}) = MaterialSlab::Vacuum(3.0);
0753   eqGrid.atPosition(Point{4.5}) = MaterialSlab::Vacuum(4.0);
0754 
0755   auto localX = std::make_unique<const LocalAccessX>();
0756   IndexedSurfaceMaterial<EqGrid>::BoundToGridLocalDelegate bToX;
0757   bToX.connect<&LocalAccessX::l2X>(std::move(localX));
0758 
0759   auto globalX = std::make_unique<const GlobalAccessX>();
0760   IndexedSurfaceMaterial<EqGrid>::GlobalToGridLocalDelegate gToX;
0761   gToX.connect<&GlobalAccessX::g2X>(std::move(globalX));
0762 
0763   GridSurfaceMaterial<EqGrid> gsm(std::move(eqGrid), GridMaterialAccessor{},
0764                                   std::move(bToX), std::move(gToX));
0765 
0766   // Global access test
0767   Vector3 g0(0.5, 0., 0.);
0768   Vector3 g1(1.5, 0., 0.);
0769   Vector3 g2(2.5, 0., 0.);
0770   Vector3 g3(3.5, 0., 0.);
0771   Vector3 g4(4.5, 0., 0.);
0772 
0773   const MaterialSlab& mg0 = gsm.materialSlab(g0);
0774   const MaterialSlab& mg1 = gsm.materialSlab(g1);
0775   const MaterialSlab& mg2 = gsm.materialSlab(g2);
0776   const MaterialSlab& mg3 = gsm.materialSlab(g3);
0777   const MaterialSlab& mg4 = gsm.materialSlab(g4);
0778 
0779   BOOST_CHECK_EQUAL(mg0.thickness(), 0.);
0780   BOOST_CHECK_EQUAL(mg1.thickness(), 1.);
0781   BOOST_CHECK_EQUAL(mg2.thickness(), 2.);
0782   BOOST_CHECK_EQUAL(mg3.thickness(), 3.);
0783   BOOST_CHECK_EQUAL(mg4.thickness(), 4.);
0784 
0785   // Local access test
0786   Vector2 l0(0.5, 0.);
0787   Vector2 l1(1.5, 0.);
0788   Vector2 l2(2.5, 0.);
0789   Vector2 l3(3.5, 0.);
0790   Vector2 l4(4.5, 0.);
0791 
0792   const MaterialSlab& ml0 = gsm.materialSlab(l0);
0793   const MaterialSlab& ml1 = gsm.materialSlab(l1);
0794   const MaterialSlab& ml2 = gsm.materialSlab(l2);
0795   const MaterialSlab& ml3 = gsm.materialSlab(l3);
0796   const MaterialSlab& ml4 = gsm.materialSlab(l4);
0797 
0798   BOOST_CHECK_EQUAL(ml0.thickness(), 0.);
0799   BOOST_CHECK_EQUAL(ml1.thickness(), 1.);
0800   BOOST_CHECK_EQUAL(ml2.thickness(), 2.);
0801   BOOST_CHECK_EQUAL(ml3.thickness(), 3.);
0802   BOOST_CHECK_EQUAL(ml4.thickness(), 4.);
0803 
0804   // Now scale it - and access again
0805   gsm.scale(2.);
0806 
0807   const MaterialSlab& sml0 = gsm.materialSlab(l0);
0808   const MaterialSlab& sml1 = gsm.materialSlab(l1);
0809   const MaterialSlab& sml2 = gsm.materialSlab(l2);
0810   const MaterialSlab& sml3 = gsm.materialSlab(l3);
0811   const MaterialSlab& sml4 = gsm.materialSlab(l4);
0812 
0813   BOOST_CHECK_EQUAL(sml0.thickness(), 0.);
0814   BOOST_CHECK_EQUAL(sml1.thickness(), 2.);
0815   BOOST_CHECK_EQUAL(sml2.thickness(), 4.);
0816   BOOST_CHECK_EQUAL(sml3.thickness(), 6.);
0817   BOOST_CHECK_EQUAL(sml4.thickness(), 8.);
0818 }
0819 
0820 BOOST_AUTO_TEST_SUITE_END()
0821 
0822 }  // namespace ActsTests