cpu/material/material_maps.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 // Detray include(s)
0010 #include "detray/definitions/indexing.hpp"
0011 #include "detray/geometry/mask.hpp"
0012 #include "detray/geometry/shapes.hpp"
0013 #include "detray/material/concepts.hpp"
0014 #include "detray/material/material_map.hpp"
0015
0016 // Detray test include(s)
0017 #include "detray/test/framework/types.hpp"
0018
0019 // GTest include(s)
0020 #include <gtest/gtest.h>
0021
0022 using namespace detray;
0023 using namespace detray::axis;
0024
0025 using test_algebra = test::algebra;
0026 using scalar = test::scalar;
0027
0028 using material_t =
0029     typename material_grid_factory<test_algebra>::bin_type::entry_type;
0030
0031 namespace {
0032
0033 material_grid_factory<test_algebra> mat_map_factory{};
0034
0035 }  // anonymous namespace
0036
0037 /// Unittest: Test the construction of an annulus shaped material map
0038 GTEST_TEST(detray_material, annulus_map) {
0039   constexpr scalar minR{7.2f * unit<scalar>::mm};
0040   constexpr scalar maxR{12.0f * unit<scalar>::mm};
0041   constexpr scalar minPhi{0.74195f};
0042   constexpr scalar maxPhi{1.33970f};
0043
0044   mask<annulus2D, test_algebra> ann2{0u,     minR, maxR, minPhi,
0045                                      maxPhi, 0.f,  -2.f, 2.f};
0046
0047   auto annulus_map = mat_map_factory.new_grid(ann2, {10u, 20u});
0048
0049   static_assert(concepts::material_map<decltype(annulus_map)>);
0050   static_assert(!concepts::homogeneous_material<decltype(annulus_map)>);
0051   static_assert(concepts::surface_material<decltype(annulus_map)>);
0052   static_assert(!concepts::volume_material<decltype(annulus_map)>);
0053
0054   EXPECT_EQ(annulus_map.dim, 2u);
0055   EXPECT_EQ(annulus_map.nbins(), 200u);
0056
0057   auto r_axis = annulus_map.get_axis<label::e_r>();
0058   EXPECT_EQ(r_axis.nbins(), 10u);
0059   EXPECT_EQ(r_axis.min(), minR);
0060   EXPECT_EQ(r_axis.max(), maxR);
0061
0062   auto phi_axis = annulus_map.get_axis<label::e_phi>();
0063   EXPECT_EQ(phi_axis.nbins(), 20u);
0064   EXPECT_EQ(phi_axis.min(), -minPhi);
0065   EXPECT_EQ(phi_axis.max(), maxPhi);
0066
0067   // Add some material entries
0068   scalar thickness = 2.f * unit<scalar>::mm;
0069   for (dindex gbin = 0; gbin < annulus_map.nbins(); ++gbin) {
0070     annulus_map.template populate<replace<>>(
0071         gbin, material_t(silicon_tml<scalar>{}, thickness));
0072     thickness += 1.f * unit<scalar>::mm;
0073   }
0074
0075   EXPECT_EQ(annulus_map.at(0, 0),
0076             material_t(silicon_tml<scalar>{}, 2.f * unit<scalar>::mm));
0077   EXPECT_EQ(annulus_map.at(1, 0),
0078             material_t(silicon_tml<scalar>{}, 3.f * unit<scalar>::mm));
0079   EXPECT_EQ(annulus_map.at(22, 0),
0080             material_t(silicon_tml<scalar>{}, 24.f * unit<scalar>::mm));
0081   EXPECT_EQ(annulus_map.at(199, 0),
0082             material_t(silicon_tml<scalar>{}, 201.f * unit<scalar>::mm));
0083   EXPECT_FALSE(annulus_map.at(199, 0) ==
0084                material_t(gold<scalar>{}, 201.f * unit<scalar>::mm));
0085 }
0086
0087 /// Unittest: Test the construction of a cylindrical material map
0088 GTEST_TEST(detray_material, cylinder_map) {
0089   constexpr scalar r{3.f * unit<scalar>::mm};
0090   constexpr scalar hz{4.f * unit<scalar>::mm};
0091
0092   mask<cylinder2D, test_algebra> cyl{0u, r, -hz, hz};
0093
0094   auto cylinder_map = mat_map_factory.new_grid(cyl, {10u, 20u});
0095
0096   static_assert(concepts::material_map<decltype(cylinder_map)>);
0097   static_assert(!concepts::homogeneous_material<decltype(cylinder_map)>);
0098   static_assert(concepts::surface_material<decltype(cylinder_map)>);
0099   static_assert(!concepts::volume_material<decltype(cylinder_map)>);
0100
0101   EXPECT_EQ(cylinder_map.dim, 2u);
0102   EXPECT_EQ(cylinder_map.nbins(), 200u);
0103
0104   auto rphi_axis = cylinder_map.get_axis<label::e_rphi>();
0105   EXPECT_EQ(rphi_axis.nbins(), 10u);
0106   EXPECT_EQ(rphi_axis.min(), -constant<scalar>::pi);
0107   EXPECT_EQ(rphi_axis.max(), constant<scalar>::pi);
0108
0109   auto z_axis = cylinder_map.get_axis<label::e_cyl_z>();
0110   EXPECT_EQ(z_axis.nbins(), 20u);
0111   EXPECT_EQ(z_axis.min(), -hz);
0112   EXPECT_EQ(z_axis.max(), hz);
0113
0114   scalar thickness = 2.f * unit<scalar>::mm;
0115   for (dindex gbin = 0; gbin < cylinder_map.nbins(); ++gbin) {
0116     cylinder_map.template populate<replace<>>(
0117         gbin, material_t(vacuum<scalar>{}, thickness));
0118     thickness += 1.f * unit<scalar>::mm;
0119   }
0120
0121   EXPECT_EQ(cylinder_map.at(0, 0),
0122             material_t(vacuum<scalar>{}, 2.f * unit<scalar>::mm));
0123   EXPECT_EQ(cylinder_map.at(1, 0),
0124             material_t(vacuum<scalar>{}, 3.f * unit<scalar>::mm));
0125   EXPECT_EQ(cylinder_map.at(22, 0),
0126             material_t(vacuum<scalar>{}, 24.f * unit<scalar>::mm));
0127   EXPECT_EQ(cylinder_map.at(199, 0),
0128             material_t(vacuum<scalar>{}, 201.f * unit<scalar>::mm));
0129   EXPECT_FALSE(cylinder_map.at(199, 0) ==
0130                material_t(gold<scalar>{}, 201.f * unit<scalar>::mm));
0131 }
0132
0133 /// Unittest: Test the construction of a rectangular material map
0134 GTEST_TEST(detray_material, rectangle_map) {
0135   constexpr scalar hx{1.f * unit<scalar>::mm};
0136   constexpr scalar hy{9.3f * unit<scalar>::mm};
0137
0138   mask<rectangle2D, test_algebra> r2{0u, hx, hy};
0139
0140   auto rectangle_map = mat_map_factory.new_grid(r2, {10u, 20u});
0141
0142   static_assert(concepts::material_map<decltype(rectangle_map)>);
0143   static_assert(!concepts::homogeneous_material<decltype(rectangle_map)>);
0144   static_assert(concepts::surface_material<decltype(rectangle_map)>);
0145   static_assert(!concepts::volume_material<decltype(rectangle_map)>);
0146
0147   EXPECT_EQ(rectangle_map.dim, 2u);
0148   EXPECT_EQ(rectangle_map.nbins(), 200u);
0149
0150   auto x_axis = rectangle_map.get_axis<label::e_x>();
0151   EXPECT_EQ(x_axis.nbins(), 10u);
0152   EXPECT_EQ(x_axis.min(), -hx);
0153   EXPECT_EQ(x_axis.max(), hx);
0154
0155   auto y_axis = rectangle_map.get_axis<label::e_y>();
0156   EXPECT_EQ(y_axis.nbins(), 20u);
0157   EXPECT_EQ(y_axis.min(), -hy);
0158   EXPECT_EQ(y_axis.max(), hy);
0159
0160   scalar thickness = 2.f * unit<scalar>::mm;
0161   for (dindex gbin = 0; gbin < rectangle_map.nbins(); ++gbin) {
0162     rectangle_map.template populate<replace<>>(
0163         gbin, material_t(oxygen_gas<scalar>{}, thickness));
0164     thickness += 1.f * unit<scalar>::mm;
0165   }
0166
0167   EXPECT_EQ(rectangle_map.at(0, 0),
0168             material_t(oxygen_gas<scalar>{}, 2.f * unit<scalar>::mm));
0169   EXPECT_EQ(rectangle_map.at(1, 0),
0170             material_t(oxygen_gas<scalar>{}, 3.f * unit<scalar>::mm));
0171   EXPECT_EQ(rectangle_map.at(22, 0),
0172             material_t(oxygen_gas<scalar>{}, 24.f * unit<scalar>::mm));
0173   EXPECT_EQ(rectangle_map.at(199, 0),
0174             material_t(oxygen_gas<scalar>{}, 201.f * unit<scalar>::mm));
0175   EXPECT_FALSE(rectangle_map.at(199, 0) ==
0176                material_t(gold<scalar>{}, 201.f * unit<scalar>::mm));
0177 }
0178
0179 /// Unittest: Test the construction of a ring shaped material map
0180 GTEST_TEST(detray_material, disc_map) {
0181   constexpr scalar inner_r{0.f * unit<scalar>::mm};
0182   constexpr scalar outer_r{3.5f * unit<scalar>::mm};
0183
0184   mask<ring2D, test_algebra> r2{0u, inner_r, outer_r};
0185
0186   auto disc_map = mat_map_factory.new_grid(r2, {10u, 20u});
0187
0188   static_assert(concepts::material_map<decltype(disc_map)>);
0189   static_assert(!concepts::homogeneous_material<decltype(disc_map)>);
0190   static_assert(concepts::surface_material<decltype(disc_map)>);
0191   static_assert(!concepts::volume_material<decltype(disc_map)>);
0192
0193   EXPECT_EQ(disc_map.dim, 2u);
0194   EXPECT_EQ(disc_map.nbins(), 200u);
0195
0196   auto r_axis = disc_map.get_axis<label::e_r>();
0197   EXPECT_EQ(r_axis.nbins(), 10u);
0198   EXPECT_EQ(r_axis.min(), inner_r);
0199   EXPECT_EQ(r_axis.max(), outer_r);
0200
0201   auto phi_axis = disc_map.get_axis<label::e_phi>();
0202   EXPECT_EQ(phi_axis.nbins(), 20u);
0203   EXPECT_EQ(phi_axis.min(), -constant<scalar>::pi);
0204   EXPECT_EQ(phi_axis.max(), constant<scalar>::pi);
0205
0206   scalar thickness = 2.f * unit<scalar>::mm;
0207   for (dindex gbin = 0; gbin < disc_map.nbins(); ++gbin) {
0208     disc_map.template populate<replace<>>(
0209         gbin, material_t(aluminium<scalar>{}, thickness));
0210     thickness += 1.f * unit<scalar>::mm;
0211   }
0212
0213   EXPECT_EQ(disc_map.at(0, 0),
0214             material_t(aluminium<scalar>{}, 2.f * unit<scalar>::mm));
0215   EXPECT_EQ(disc_map.at(1, 0),
0216             material_t(aluminium<scalar>{}, 3.f * unit<scalar>::mm));
0217   EXPECT_EQ(disc_map.at(22, 0),
0218             material_t(aluminium<scalar>{}, 24.f * unit<scalar>::mm));
0219   EXPECT_EQ(disc_map.at(199, 0),
0220             material_t(aluminium<scalar>{}, 201.f * unit<scalar>::mm));
0221   EXPECT_FALSE(disc_map.at(199, 0) ==
0222                material_t(gold<scalar>{}, 201.f * unit<scalar>::mm));
0223 }
0224
0225 /// Unittest: Test the construction of a trapezoid material map
0226 GTEST_TEST(detray_material, trapezoid_map) {
0227   constexpr scalar hx_miny{1.f * unit<scalar>::mm};
0228   constexpr scalar hx_maxy{3.f * unit<scalar>::mm};
0229   constexpr scalar hy{2.f * unit<scalar>::mm};
0230   constexpr scalar divisor{1.f / (2.f * hy)};
0231
0232   mask<trapezoid2D, test_algebra> t2{0u, hx_miny, hx_maxy, hy, divisor};
0233
0234   auto trapezoid_map = mat_map_factory.new_grid(t2, {10u, 20u});
0235
0236   static_assert(concepts::material_map<decltype(trapezoid_map)>);
0237   static_assert(!concepts::homogeneous_material<decltype(trapezoid_map)>);
0238   static_assert(concepts::surface_material<decltype(trapezoid_map)>);
0239   static_assert(!concepts::volume_material<decltype(trapezoid_map)>);
0240
0241   EXPECT_EQ(trapezoid_map.dim, 2u);
0242   EXPECT_EQ(trapezoid_map.nbins(), 200u);
0243
0244   auto x_axis = trapezoid_map.get_axis<label::e_x>();
0245   EXPECT_EQ(x_axis.nbins(), 10u);
0246   EXPECT_EQ(x_axis.min(), -hx_maxy);
0247   EXPECT_EQ(x_axis.max(), hx_maxy);
0248
0249   auto y_axis = trapezoid_map.get_axis<label::e_y>();
0250   EXPECT_EQ(y_axis.nbins(), 20u);
0251   EXPECT_EQ(y_axis.min(), -hy);
0252   EXPECT_EQ(y_axis.max(), hy);
0253
0254   scalar thickness = 2.f * unit<scalar>::mm;
0255   for (dindex gbin = 0; gbin < trapezoid_map.nbins(); ++gbin) {
0256     trapezoid_map.template populate<replace<>>(
0257         gbin, material_t(gold<scalar>{}, thickness));
0258     thickness += 1.f * unit<scalar>::mm;
0259   }
0260
0261   EXPECT_EQ(trapezoid_map.at(0, 0),
0262             material_t(gold<scalar>{}, 2.f * unit<scalar>::mm));
0263   EXPECT_EQ(trapezoid_map.at(1, 0),
0264             material_t(gold<scalar>{}, 3.f * unit<scalar>::mm));
0265   EXPECT_EQ(trapezoid_map.at(22, 0),
0266             material_t(gold<scalar>{}, 24.f * unit<scalar>::mm));
0267   EXPECT_EQ(trapezoid_map.at(199, 0),
0268             material_t(gold<scalar>{}, 201.f * unit<scalar>::mm));
0269   EXPECT_FALSE(trapezoid_map.at(199, 0) ==
0270                material_t(aluminium<scalar>{}, 201.f * unit<scalar>::mm));
0271 }
0272
0273 /// Unittest: Test the material grid comparisons
0274 GTEST_TEST(detray_material, material_grid_comparison) {
0275   /** Allows to create a regular grid to check the equality operator
0276    * grids can differ in:
0277    * - type (will never be compared)
0278    * - bins/axes
0279    * - entries (i.e. material data)
0280    */
0281   auto createGrid = [](const scalar hx, const scalar hy, unsigned int bx,
0282                        unsigned int by, bool distort_entries) {
0283     mask<rectangle2D, test_algebra> r2{0u, hx, hy};
0284     auto material_grid = mat_map_factory.new_grid(r2, {bx, by});
0285
0286     // Fill the material grid with some data
0287     scalar thickness = 2.f * unit<scalar>::mm;
0288     for (dindex gbin = 0; gbin < material_grid.nbins(); ++gbin) {
0289       material_grid.template populate<replace<>>(
0290           gbin, material_t(oxygen_gas<scalar>{}, thickness));
0291       thickness += 1.f * unit<scalar>::mm;
0292       if (distort_entries) {
0293         thickness += 1.f * unit<scalar>::mm;
0294       }
0295     }
0296     // Return it
0297     return material_grid;
0298   };
0299
0300   // Two equal grids
0301   auto grid_ref = createGrid(10.f, 20.f, 10u, 20u, false);
0302   auto grid_eq = createGrid(10.f, 20.f, 10u, 20u, false);
0303   EXPECT_EQ(grid_ref, grid_eq);
0304
0305   // One grid with different size
0306   auto grid_neq_size = createGrid(11.f, 21.f, 10u, 20u, false);
0307   EXPECT_NE(grid_ref, grid_neq_size);
0308
0309   // One grid with different binning
0310   auto grid_neq_bins = createGrid(10.f, 20.f, 11u, 21u, false);
0311   EXPECT_NE(grid_ref, grid_neq_bins);
0312
0313   // One grid with different entries
0314   auto grid_neq_entries = createGrid(10.f, 20.f, 10u, 20u, true);
0315   EXPECT_NE(grid_ref, grid_neq_entries);
0316 }