Core/Material/MaterialGridHelperTests.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/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Material/AccumulatedVolumeMaterial.hpp"
0013 #include "Acts/Material/Material.hpp"
0014 #include "Acts/Material/MaterialGridHelper.hpp"
0015 #include "Acts/Material/MaterialSlab.hpp"
0016 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0017 #include "Acts/Utilities/Axis.hpp"
0018 #include "Acts/Utilities/AxisDefinitions.hpp"
0019 #include "Acts/Utilities/BinUtility.hpp"
0020 #include "Acts/Utilities/BinningType.hpp"
0021 #include "Acts/Utilities/Grid.hpp"
0022
0023 #include <cmath>
0024 #include <functional>
0025 #include <memory>
0026 #include <numbers>
0027 #include <utility>
0028 #include <vector>
0029
0030 namespace Acts::Test {
0031
0032 using EAxis = Acts::Axis<AxisType::Equidistant>;
0033 using Grid2D = Acts::Grid<Acts::AccumulatedVolumeMaterial, EAxis, EAxis>;
0034 using Grid3D = Acts::Grid<Acts::AccumulatedVolumeMaterial, EAxis, EAxis, EAxis>;
0035 using MaterialGrid2D =
0036     Acts::Grid<Acts::Material::ParametersVector, EAxis, EAxis>;
0037 using MaterialGrid3D =
0038     Acts::Grid<Acts::Material::ParametersVector, EAxis, EAxis, EAxis>;
0039
0040 /// @brief Various test for the Material in the case of a Cuboid volume and 2D
0041 /// Grid
0042 BOOST_AUTO_TEST_CASE(Square_Grid_test) {
0043   BinUtility bu(7, -3., 3., open, AxisDirection::AxisX);
0044   bu += BinUtility(3, -2., 2., open, AxisDirection::AxisY);
0045   auto bd = bu.binningData();
0046   std::function<Acts::Vector2(Acts::Vector3)> transfoGlobalToLocal;
0047
0048   Grid2D Grid = createGrid2D(bu, transfoGlobalToLocal);
0049
0050   // Test Global To Local transform
0051   Acts::Vector3 pos(1., 2., 3.);
0052   Acts::Vector2 pos_2d(1., 2.);
0053   BOOST_CHECK_EQUAL(pos_2d, transfoGlobalToLocal(pos));
0054
0055   // Test Grid
0056   BOOST_CHECK_EQUAL(Grid.dimensions(), 2);
0057
0058   BOOST_CHECK_EQUAL(Grid.numLocalBins()[0], bd[0].bins());
0059   BOOST_CHECK_EQUAL(Grid.numLocalBins()[1], bd[1].bins());
0060
0061   BOOST_CHECK_EQUAL(Grid.minPosition()[0], bd[0].min);
0062   BOOST_CHECK_EQUAL(Grid.minPosition()[1], bd[1].min);
0063
0064   float max1 = bd[0].max + std::abs(bd[0].max - bd[0].min) / (bd[0].bins() - 1);
0065   float max2 = bd[1].max + std::abs(bd[1].max - bd[1].min) / (bd[1].bins() - 1);
0066
0067   BOOST_CHECK_EQUAL(Grid.maxPosition()[0], max1);
0068   BOOST_CHECK_EQUAL(Grid.maxPosition()[1], max2);
0069
0070   // Test pos to index
0071   Grid2D::index_t index1 = {1, 1};
0072   Grid2D::index_t index2 = {7, 2};
0073   Grid2D::index_t index3 = {1, 3};
0074
0075   Acts::Vector3 pos1 = {-2.6, -1.5, -0.7};
0076   Acts::Vector3 pos2 = {2.8, 0, 0.2};
0077   Acts::Vector3 pos3 = {-2.7, 1.8, 0.8};
0078
0079   for (int i = 0; i < 2; i++) {
0080     BOOST_CHECK_EQUAL(
0081         Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(pos1))[i],
0082         index1[i]);
0083     BOOST_CHECK_EQUAL(
0084         Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(pos2))[i],
0085         index2[i]);
0086     BOOST_CHECK_EQUAL(
0087         Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(pos3))[i],
0088         index3[i]);
0089   }
0090   // Test material mapping
0091
0092   std::vector<Acts::Vector3> vectPos1;
0093   vectPos1.push_back(pos1);
0094   std::vector<Acts::Vector3> vectPos2;
0095   vectPos2.push_back(pos2);
0096   std::vector<Acts::Vector3> vectPos3;
0097   vectPos3.push_back(pos3);
0098
0099   std::vector<std::pair<MaterialSlab, std::vector<Vector3>>> matRecord;
0100   Material mat1 = Material::fromMolarDensity(1., 2., 3., 4., 5.);
0101   Material mat2 = Material::fromMolarDensity(6., 7., 8., 9., 10.);
0102   Material vacuum;
0103
0104   MaterialSlab matprop1(mat1, 1);
0105   MaterialSlab matprop2(mat2, 1);
0106
0107   matRecord.clear();
0108   matRecord.push_back(std::make_pair(matprop1, vectPos1));
0109   matRecord.push_back(std::make_pair(matprop2, vectPos2));
0110
0111   // Walk over each property
0112   for (const auto& rm : matRecord) {
0113     // Walk over each point associated with the properties
0114     for (const auto& point : rm.second) {
0115       // Search for fitting grid point and accumulate
0116       Acts::Grid2D::index_t index =
0117           Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(point));
0118       Grid.atLocalBins(index).accumulate(rm.first);
0119     }
0120   }
0121
0122   MaterialGrid2D matMap = mapMaterialPoints(Grid);
0123
0124   CHECK_CLOSE_REL(matMap.atLocalBins(index1), mat1.parameters(), 1e-4);
0125   CHECK_CLOSE_REL(matMap.atLocalBins(index2), mat2.parameters(), 1e-4);
0126   BOOST_CHECK_EQUAL(matMap.atLocalBins(index3), vacuum.parameters());
0127 }
0128
0129 /// @brief Various test for the Material in the case of a Cylindrical volume
0130 /// with a 2D grid
0131 BOOST_AUTO_TEST_CASE(PhiZ_Grid_test) {
0132   BinUtility bu(2, -2., 2., open, AxisDirection::AxisZ);
0133   bu += BinUtility(3, -std::numbers::pi, std::numbers::pi, closed,
0134                    AxisDirection::AxisPhi);
0135   auto bd = bu.binningData();
0136   std::function<Acts::Vector2(Acts::Vector3)> transfoGlobalToLocal;
0137
0138   Grid2D Grid = createGrid2D(bu, transfoGlobalToLocal);
0139
0140   // Test Global To Local transform
0141   Acts::Vector3 pos(1., 2., 3.);
0142
0143   CHECK_CLOSE_REL(transfoGlobalToLocal(pos)[1], atan2(2, 1), 1e-4);
0144   CHECK_CLOSE_REL(transfoGlobalToLocal(pos)[0], 3, 1e-4);
0145
0146   // Test Grid
0147   BOOST_CHECK_EQUAL(Grid.dimensions(), 2);
0148
0149   BOOST_CHECK_EQUAL(Grid.numLocalBins()[0], bd[0].bins());
0150   BOOST_CHECK_EQUAL(Grid.numLocalBins()[1], bd[1].bins());
0151
0152   BOOST_CHECK_EQUAL(Grid.minPosition()[0], bd[0].min);
0153   BOOST_CHECK_EQUAL(Grid.minPosition()[1], bd[1].min);
0154
0155   float max1 = bd[0].max + std::abs(bd[0].max - bd[0].min) / (bd[0].bins() - 1);
0156   float max2 = bd[1].max + std::abs(bd[1].max - bd[1].min) / (bd[1].bins() - 1);
0157
0158   BOOST_CHECK_EQUAL(Grid.maxPosition()[0], max1);
0159   BOOST_CHECK_EQUAL(Grid.maxPosition()[1], max2);
0160
0161   // Test pos to index
0162   Grid2D::index_t index1 = {1, 1};
0163   Grid2D::index_t index2 = {1, 2};
0164   Grid2D::index_t index3 = {2, 3};
0165
0166   Acts::Vector3 pos1 = {-0.2, -1, -1};
0167   Acts::Vector3 pos2 = {3.6, 0., -1.5};
0168   Acts::Vector3 pos3 = {-1, 0.3, 0.8};
0169
0170   for (int i = 0; i < 2; i++) {
0171     BOOST_CHECK_EQUAL(
0172         Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(pos1))[i],
0173         index1[i]);
0174     BOOST_CHECK_EQUAL(
0175         Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(pos2))[i],
0176         index2[i]);
0177     BOOST_CHECK_EQUAL(
0178         Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(pos3))[i],
0179         index3[i]);
0180   }
0181
0182   // Test material mapping
0183   std::vector<Acts::Vector3> vectPos1;
0184   vectPos1.push_back(pos1);
0185   std::vector<Acts::Vector3> vectPos2;
0186   vectPos2.push_back(pos2);
0187   std::vector<Acts::Vector3> vectPos3;
0188   vectPos3.push_back(pos3);
0189
0190   std::vector<std::pair<MaterialSlab, std::vector<Vector3>>> matRecord;
0191   Material mat1 = Material::fromMolarDensity(1., 2., 3., 4., 5.);
0192   Material mat2 = Material::fromMolarDensity(6., 7., 8., 9., 10.);
0193   Material vacuum;
0194
0195   MaterialSlab matprop1(mat1, 1);
0196   MaterialSlab matprop2(mat2, 1);
0197
0198   matRecord.clear();
0199   matRecord.push_back(std::make_pair(matprop1, vectPos1));
0200   matRecord.push_back(std::make_pair(matprop2, vectPos2));
0201
0202   // Walk over each property
0203   for (const auto& rm : matRecord) {
0204     // Walk over each point associated with the properties
0205     for (const auto& point : rm.second) {
0206       // Search for fitting grid point and accumulate
0207       Acts::Grid2D::index_t index =
0208           Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(point));
0209       Grid.atLocalBins(index).accumulate(rm.first);
0210     }
0211   }
0212
0213   MaterialGrid2D matMap = mapMaterialPoints(Grid);
0214
0215   CHECK_CLOSE_REL(matMap.atLocalBins(index1), mat1.parameters(), 1e-4);
0216   CHECK_CLOSE_REL(matMap.atLocalBins(index2), mat2.parameters(), 1e-4);
0217   BOOST_CHECK_EQUAL(matMap.atLocalBins(index3), vacuum.parameters());
0218 }
0219
0220 /// @brief Various test for the Material in the case of a Cuboid volume
0221 BOOST_AUTO_TEST_CASE(Cubic_Grid_test) {
0222   BinUtility bu(7, -3., 3., open, AxisDirection::AxisX);
0223   bu += BinUtility(3, -2., 2., open, AxisDirection::AxisY);
0224   bu += BinUtility(2, -1., 1., open, AxisDirection::AxisZ);
0225   auto bd = bu.binningData();
0226   std::function<Acts::Vector3(Acts::Vector3)> transfoGlobalToLocal;
0227
0228   Grid3D Grid = createGrid3D(bu, transfoGlobalToLocal);
0229
0230   // Test Global To Local transform
0231   Acts::Vector3 pos(1., 2., 3.);
0232   BOOST_CHECK_EQUAL(pos, transfoGlobalToLocal(pos));
0233
0234   // Test Grid
0235   BOOST_CHECK_EQUAL(Grid.dimensions(), 3);
0236
0237   BOOST_CHECK_EQUAL(Grid.numLocalBins()[0], bd[0].bins());
0238   BOOST_CHECK_EQUAL(Grid.numLocalBins()[1], bd[1].bins());
0239   BOOST_CHECK_EQUAL(Grid.numLocalBins()[2], bd[2].bins());
0240
0241   BOOST_CHECK_EQUAL(Grid.minPosition()[0], bd[0].min);
0242   BOOST_CHECK_EQUAL(Grid.minPosition()[1], bd[1].min);
0243   BOOST_CHECK_EQUAL(Grid.minPosition()[2], bd[2].min);
0244
0245   float max1 = bd[0].max + std::abs(bd[0].max - bd[0].min) / (bd[0].bins() - 1);
0246   float max2 = bd[1].max + std::abs(bd[1].max - bd[1].min) / (bd[1].bins() - 1);
0247   float max3 = bd[2].max + std::abs(bd[2].max - bd[2].min) / (bd[2].bins() - 1);
0248
0249   BOOST_CHECK_EQUAL(Grid.maxPosition()[0], max1);
0250   BOOST_CHECK_EQUAL(Grid.maxPosition()[1], max2);
0251   BOOST_CHECK_EQUAL(Grid.maxPosition()[2], max3);
0252
0253   // Test pos to index
0254   Grid3D::index_t index1 = {1, 1, 1};
0255   Grid3D::index_t index2 = {7, 2, 2};
0256   Grid3D::index_t index3 = {1, 3, 2};
0257
0258   Acts::Vector3 pos1 = {-2.6, -1.5, -0.7};
0259   Acts::Vector3 pos2 = {2.8, 0, 0.2};
0260   Acts::Vector3 pos3 = {-2.7, 1.8, 0.8};
0261
0262   for (int i = 0; i < 3; i++) {
0263     BOOST_CHECK_EQUAL(
0264         Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(pos1))[i],
0265         index1[i]);
0266     BOOST_CHECK_EQUAL(
0267         Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(pos2))[i],
0268         index2[i]);
0269     BOOST_CHECK_EQUAL(
0270         Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(pos3))[i],
0271         index3[i]);
0272   }
0273   // Test material mapping
0274   std::vector<Acts::Vector3> vectPos1;
0275   vectPos1.push_back(pos1);
0276   std::vector<Acts::Vector3> vectPos2;
0277   vectPos2.push_back(pos2);
0278   std::vector<Acts::Vector3> vectPos3;
0279   vectPos3.push_back(pos3);
0280
0281   std::vector<std::pair<MaterialSlab, std::vector<Vector3>>> matRecord;
0282   Material mat1 = Material::fromMolarDensity(1., 2., 3., 4., 5.);
0283   Material mat2 = Material::fromMolarDensity(6., 7., 8., 9., 10.);
0284   Material vacuum;
0285
0286   MaterialSlab matprop1(mat1, 1);
0287   MaterialSlab matprop2(mat2, 1);
0288
0289   matRecord.clear();
0290   matRecord.push_back(std::make_pair(matprop1, vectPos1));
0291   matRecord.push_back(std::make_pair(matprop2, vectPos2));
0292
0293   // Walk over each property
0294   for (const auto& rm : matRecord) {
0295     // Walk over each point associated with the properties
0296     for (const auto& point : rm.second) {
0297       // Search for fitting grid point and accumulate
0298       Acts::Grid3D::index_t index =
0299           Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(point));
0300       Grid.atLocalBins(index).accumulate(rm.first);
0301     }
0302   }
0303
0304   MaterialGrid3D matMap = mapMaterialPoints(Grid);
0305
0306   CHECK_CLOSE_REL(matMap.atLocalBins(index1), mat1.parameters(), 1e-4);
0307   CHECK_CLOSE_REL(matMap.atLocalBins(index2), mat2.parameters(), 1e-4);
0308   BOOST_CHECK_EQUAL(matMap.atLocalBins(index3), vacuum.parameters());
0309 }
0310
0311 /// @brief Various test for the Material in the case of a Cylindrical volume
0312 BOOST_AUTO_TEST_CASE(Cylindrical_Grid_test) {
0313   BinUtility bu(4, 1., 4., open, AxisDirection::AxisR);
0314   bu += BinUtility(3, -std::numbers::pi, std::numbers::pi, closed,
0315                    AxisDirection::AxisPhi);
0316   bu += BinUtility(2, -2., 2., open, AxisDirection::AxisZ);
0317   auto bd = bu.binningData();
0318   std::function<Acts::Vector3(Acts::Vector3)> transfoGlobalToLocal;
0319
0320   Grid3D Grid = createGrid3D(bu, transfoGlobalToLocal);
0321
0322   // Test Global To Local transform
0323   Acts::Vector3 pos(1., 2., 3.);
0324
0325   CHECK_CLOSE_REL(transfoGlobalToLocal(pos)[0], sqrt(5), 1e-4);
0326   CHECK_CLOSE_REL(transfoGlobalToLocal(pos)[1], atan2(2, 1), 1e-4);
0327   CHECK_CLOSE_REL(transfoGlobalToLocal(pos)[2], 3, 1e-4);
0328
0329   // Test Grid
0330   BOOST_CHECK_EQUAL(Grid.dimensions(), 3);
0331
0332   BOOST_CHECK_EQUAL(Grid.numLocalBins()[0], bd[0].bins());
0333   BOOST_CHECK_EQUAL(Grid.numLocalBins()[1], bd[1].bins());
0334   BOOST_CHECK_EQUAL(Grid.numLocalBins()[2], bd[2].bins());
0335
0336   BOOST_CHECK_EQUAL(Grid.minPosition()[0], bd[0].min);
0337   BOOST_CHECK_EQUAL(Grid.minPosition()[1], bd[1].min);
0338   BOOST_CHECK_EQUAL(Grid.minPosition()[2], bd[2].min);
0339
0340   float max1 = bd[0].max + std::abs(bd[0].max - bd[0].min) / (bd[0].bins() - 1);
0341   float max2 = bd[1].max + std::abs(bd[1].max - bd[1].min) / (bd[1].bins() - 1);
0342   float max3 = bd[2].max + std::abs(bd[2].max - bd[2].min) / (bd[2].bins() - 1);
0343
0344   BOOST_CHECK_EQUAL(Grid.maxPosition()[0], max1);
0345   BOOST_CHECK_EQUAL(Grid.maxPosition()[1], max2);
0346   BOOST_CHECK_EQUAL(Grid.maxPosition()[2], max3);
0347
0348   // Test pos to index
0349   Grid3D::index_t index1 = {1, 1, 1};
0350   Grid3D::index_t index2 = {4, 2, 1};
0351   Grid3D::index_t index3 = {1, 3, 2};
0352
0353   Acts::Vector3 pos1 = {-0.2, -1, -1};
0354   Acts::Vector3 pos2 = {3.6, 0., -1.5};
0355   Acts::Vector3 pos3 = {-1, 0.3, 0.8};
0356
0357   for (int i = 0; i < 3; i++) {
0358     BOOST_CHECK_EQUAL(
0359         Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(pos1))[i],
0360         index1[i]);
0361     BOOST_CHECK_EQUAL(
0362         Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(pos2))[i],
0363         index2[i]);
0364     BOOST_CHECK_EQUAL(
0365         Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(pos3))[i],
0366         index3[i]);
0367   }
0368
0369   // Test material mapping
0370   std::vector<Acts::Vector3> vectPos1;
0371   vectPos1.push_back(pos1);
0372   std::vector<Acts::Vector3> vectPos2;
0373   vectPos2.push_back(pos2);
0374   std::vector<Acts::Vector3> vectPos3;
0375   vectPos3.push_back(pos3);
0376
0377   std::vector<std::pair<MaterialSlab, std::vector<Vector3>>> matRecord;
0378   Material mat1 = Material::fromMolarDensity(1., 2., 3., 4., 5.);
0379   Material mat2 = Material::fromMolarDensity(6., 7., 8., 9., 10.);
0380   Material vacuum;
0381
0382   MaterialSlab matprop1(mat1, 1);
0383   MaterialSlab matprop2(mat2, 1);
0384
0385   matRecord.clear();
0386   matRecord.push_back(std::make_pair(matprop1, vectPos1));
0387   matRecord.push_back(std::make_pair(matprop2, vectPos2));
0388
0389   // Walk over each property
0390   for (const auto& rm : matRecord) {
0391     // Walk over each point associated with the properties
0392     for (const auto& point : rm.second) {
0393       // Search for fitting grid point and accumulate
0394       Acts::Grid3D::index_t index =
0395           Grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(point));
0396       Grid.atLocalBins(index).accumulate(rm.first);
0397     }
0398   }
0399
0400   MaterialGrid3D matMap = mapMaterialPoints(Grid);
0401
0402   CHECK_CLOSE_REL(matMap.atLocalBins(index1), mat1.parameters(), 1e-4);
0403   CHECK_CLOSE_REL(matMap.atLocalBins(index2), mat2.parameters(), 1e-4);
0404   BOOST_CHECK_EQUAL(matMap.atLocalBins(index3), vacuum.parameters());
0405 }
0406
0407 }  // namespace Acts::Test