cpu/builders/material_map_builder.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/builders/material_map_builder.hpp"
0011
0012 #include "detray/builders/cuboid_portal_generator.hpp"
0013 #include "detray/builders/detector_builder.hpp"
0014 #include "detray/builders/material_map_factory.hpp"
0015 #include "detray/builders/surface_factory.hpp"
0016 #include "detray/builders/volume_builder.hpp"
0017 #include "detray/core/detector.hpp"
0018 #include "detray/definitions/indexing.hpp"
0019 #include "detray/utils/ranges.hpp"
0020
0021 // Detray test include(s)
0022 #include "detray/test/framework/types.hpp"
0023
0024 // Vecmem include(s)
0025 #include <vecmem/memory/host_memory_resource.hpp>
0026
0027 // GTest include(s)
0028 #include <gtest/gtest.h>
0029
0030 // System include(s)
0031 #include <limits>
0032 #include <memory>
0033 #include <vector>
0034
0035 using namespace detray;
0036
0037 namespace {
0038
0039 using point3 = test::point3;
0040
0041 using metadata_t = test::default_metadata;
0042 using detector_t = detector<metadata_t>;
0043 using mat_id = typename detector_t::material::id;
0044 using bin_index_t = axis::multi_bin<2u>;
0045
0046 /// Add generate input material for material maps
0047 template <typename material_factory_t, concepts::scalar scalar_t>
0048 auto add_material_data(const material_factory_t& mat_factory, mat_id id,
0049                        std::size_t sf_index, scalar_t t,
0050                        material<scalar_t> mat = silicon<scalar_t>()) {
0051   typename material_factory_t::element_type::data_type mat_data{sf_index};
0052   std::vector<bin_index_t> m_bins{};
0053   std::vector<std::size_t> n_bins{5u, 10u};
0054   std::vector<std::vector<scalar_t>> axis_spans = {};
0055
0056   // Add material for every bin
0057   for (auto [i, j] : detray::views::cartesian_product{
0058            detray::views::iota{0u, 5u}, detray::views::iota{0u, 10u}}) {
0059     m_bins.push_back({i, j});
0060     mat_data.append(t, mat);
0061     t += 0.25f * unit<scalar_t>::mm;
0062   }
0063
0064   mat_factory->add_material(id, std::move(mat_data), std::move(n_bins),
0065                             std::move(axis_spans), std::move(m_bins));
0066 }
0067
0068 }  // anonymous namespace
0069
0070 /// Integration test: material builder as volume builder decorator
0071 GTEST_TEST(detray_builders, decorator_material_map_builder) {
0072   using test_algebra = typename detector_t::algebra_type;
0073   using scalar = dscalar<test_algebra>;
0074   using transform3 = dtransform3D<test_algebra>;
0075   using mask_id = typename detector_t::masks::id;
0076
0077   using pt_cylinder_factory_t =
0078       surface_factory<detector_t, concentric_cylinder2D>;
0079   using rectangle_factory = surface_factory<detector_t, rectangle2D>;
0080   using trapezoid_factory = surface_factory<detector_t, trapezoid2D>;
0081   using cylinder_factory = surface_factory<detector_t, cylinder2D>;
0082
0083   using mat_factory_t = material_map_factory<detector_t, bin_index_t>;
0084
0085   vecmem::host_memory_resource host_mr;
0086   detector_t d(host_mr);
0087   auto geo_ctx = typename detector_t::geometry_context{};
0088
0089   auto vbuilder =
0090       std::make_unique<volume_builder<detector_t>>(volume_id::e_cylinder);
0091   auto mat_builder = material_map_builder<detector_t>{std::move(vbuilder)};
0092
0093   EXPECT_TRUE(d.volumes().empty());
0094
0095   // Add some portals first
0096   auto pt_cyl_factory = std::make_unique<pt_cylinder_factory_t>();
0097
0098   pt_cyl_factory->push_back({surface_id::e_portal,
0099                              transform3(point3{0.f, 0.f, 0.f}), 1u,
0100                              std::vector<scalar>{10.f, -1500.f, 1500.f}});
0101   pt_cyl_factory->push_back({surface_id::e_portal,
0102                              transform3(point3{0.f, 0.f, 0.f}), 2u,
0103                              std::vector<scalar>{20.f, -1500.f, 1500.f}});
0104
0105   // Then some passive and sensitive surfaces
0106   auto rect_factory = std::make_unique<rectangle_factory>();
0107
0108   typename rectangle_factory::sf_data_collection rect_sf_data;
0109   rect_sf_data.emplace_back(surface_id::e_sensitive,
0110                             transform3(point3{0.f, 0.f, -10.f}), 0u,
0111                             std::vector<scalar>{10.f, 8.f});
0112   rect_sf_data.emplace_back(surface_id::e_sensitive,
0113                             transform3(point3{0.f, 0.f, -20.f}), 0u,
0114                             std::vector<scalar>{10.f, 8.f});
0115   rect_sf_data.emplace_back(surface_id::e_sensitive,
0116                             transform3(point3{0.f, 0.f, -30.f}), 0u,
0117                             std::vector<scalar>{10.f, 8.f});
0118   rect_factory->push_back(std::move(rect_sf_data));
0119
0120   auto trpz_factory = std::make_unique<trapezoid_factory>();
0121
0122   trpz_factory->push_back({surface_id::e_sensitive,
0123                            transform3(point3{0.f, 0.f, 1000.f}), 0u,
0124                            std::vector<scalar>{1.f, 3.f, 2.f, 0.25f}});
0125
0126   auto cyl_factory = std::make_unique<cylinder_factory>();
0127
0128   cyl_factory->push_back({surface_id::e_passive,
0129                           transform3(point3{0.f, 0.f, 0.f}), 0u,
0130                           std::vector<scalar>{5.f, -1300.f, 1300.f}});
0131
0132   // Now add the material for each surface
0133   auto mat_pt_cyl_factory =
0134       std::make_shared<mat_factory_t>(std::move(pt_cyl_factory));
0135   scalar t{1.f * unit<scalar>::mm};
0136   add_material_data(mat_pt_cyl_factory, mat_id::e_concentric_cylinder2D_map, 0u,
0137                     t, silicon<scalar>());
0138   add_material_data(mat_pt_cyl_factory, mat_id::e_concentric_cylinder2D_map, 1u,
0139                     t, silicon<scalar>());
0140
0141   auto mat_rect_factory =
0142       std::make_shared<mat_factory_t>(std::move(rect_factory));
0143   t = 1.f * unit<scalar>::mm;
0144   add_material_data(mat_rect_factory, mat_id::e_rectangle2D_map, 2u, t,
0145                     tungsten<scalar>());
0146   // No material for surface with index 3
0147   t = 3.f * unit<scalar>::mm;
0148   add_material_data(mat_rect_factory, mat_id::e_rectangle2D_map, 4u, t,
0149                     tungsten<scalar>());
0150
0151   auto mat_trpz_factory =
0152       std::make_shared<mat_factory_t>(std::move(trpz_factory));
0153   t = 1.f * unit<scalar>::mm;
0154   add_material_data(mat_trpz_factory, mat_id::e_trapezoid2D_map, 5u, t,
0155                     tungsten<scalar>());
0156
0157   auto mat_cyl_factory =
0158       std::make_shared<mat_factory_t>(std::move(cyl_factory));
0159   t = 1.5f * unit<scalar>::mm;
0160   add_material_data(mat_cyl_factory, mat_id::e_cylinder2D_map, 6u, t,
0161                     gold<scalar>());
0162
0163   // Add surfaces and material to detector
0164   mat_builder.add_surfaces(mat_pt_cyl_factory, geo_ctx);
0165   mat_builder.add_surfaces(mat_rect_factory, geo_ctx);
0166   mat_builder.add_surfaces(mat_trpz_factory, geo_ctx);
0167   mat_builder.add_surfaces(mat_cyl_factory, geo_ctx);
0168
0169   // Add the volume to the detector
0170   mat_builder.build(d);
0171
0172   //
0173   // check results
0174   //
0175   const auto& vol = d.volumes().back();
0176   EXPECT_TRUE(d.volumes().size() == 1u);
0177   EXPECT_EQ(vol.index(), 0u);
0178   EXPECT_EQ(vol.id(), volume_id::e_cylinder);
0179
0180   EXPECT_EQ(d.surfaces().size(), 7u);
0181   EXPECT_EQ(d.transform_store().size(), 8u);
0182   EXPECT_EQ(d.mask_store().template size<mask_id::e_concentric_cylinder2D>(),
0183             2u);
0184   EXPECT_EQ(d.mask_store().template size<mask_id::e_ring2D>(), 0u);
0185   EXPECT_EQ(d.mask_store().template size<mask_id::e_cylinder2D>(), 1u);
0186   EXPECT_EQ(d.mask_store().template size<mask_id::e_rectangle2D>(), 3u);
0187   EXPECT_EQ(d.mask_store().template size<mask_id::e_trapezoid2D>(), 1u);
0188
0189   EXPECT_EQ(d.material_store().template size<mat_id::e_material_slab>(), 0u);
0190   EXPECT_EQ(d.material_store().template size<mat_id::e_material_rod>(), 0u);
0191   EXPECT_EQ(d.material_store().template size<mat_id::e_ring2D_map>(), 0u);
0192   EXPECT_EQ(d.material_store().template size<mat_id::e_annulus2D_map>(), 0u);
0193   EXPECT_EQ(d.material_store().template size<mat_id::e_cylinder2D_map>(), 1u);
0194   EXPECT_EQ(
0195       d.material_store().template size<mat_id::e_concentric_cylinder2D_map>(),
0196       2u);
0197   // Rectangle and trapezoid surfaces have the same grid geometry
0198   EXPECT_EQ(d.material_store().template size<mat_id::e_rectangle2D_map>(), 3u);
0199   EXPECT_EQ(d.material_store().template size<mat_id::e_trapezoid2D_map>(), 3u);
0200
0201   // Check the material links
0202   std::size_t pt_cyl_idx{0u};
0203   std::size_t cyl_idx{0u};
0204   std::size_t cart_idx{0u};
0205   for (auto& sf_desc : d.surfaces()) {
0206     const auto& mat_link = sf_desc.material();
0207     switch (mat_link.id()) {
0208       case mat_id::e_cylinder2D_map: {
0209         EXPECT_EQ(mat_link.index(), cyl_idx++) << sf_desc;
0210         break;
0211       }
0212       case mat_id::e_concentric_cylinder2D_map: {
0213         EXPECT_EQ(mat_link.index(), pt_cyl_idx++) << sf_desc;
0214         break;
0215       }
0216       case mat_id::e_rectangle2D_map: {
0217         EXPECT_EQ(mat_link.index(), cart_idx++) << sf_desc;
0218         break;
0219       }
0220       case mat_id::e_none: {
0221         // No material on surface 3
0222         EXPECT_TRUE(sf_desc.index() == 3u) << "No material on: " << sf_desc;
0223         EXPECT_TRUE(mat_link.id() == mat_id::e_none) << sf_desc;
0224         EXPECT_TRUE(mat_link.is_invalid()) << sf_desc;
0225         break;
0226       }
0227       default: {
0228         EXPECT_TRUE(false) << sf_desc;
0229       }
0230     }
0231   }
0232
0233   // Check the material map content
0234   for (auto cyl_mat_grid :
0235        d.material_store().template get<mat_id::e_concentric_cylinder2D_map>()) {
0236     EXPECT_EQ(cyl_mat_grid.nbins(), 50u);
0237     EXPECT_EQ(cyl_mat_grid.size(), 50u);
0238
0239     auto r_axis = cyl_mat_grid.template get_axis<axis::label::e_rphi>();
0240     EXPECT_EQ(r_axis.nbins(), 5u);
0241     auto z_axis = cyl_mat_grid.template get_axis<axis::label::e_cyl_z>();
0242     EXPECT_EQ(z_axis.nbins(), 10u);
0243
0244     for (const auto& mat_slab : cyl_mat_grid.all()) {
0245       EXPECT_TRUE(mat_slab.get_material() == silicon<scalar>() ||
0246                   mat_slab.get_material() == gold<scalar>());
0247     }
0248   }
0249
0250   for (auto cart_mat_grid :
0251        d.material_store().template get<mat_id::e_rectangle2D_map>()) {
0252     EXPECT_EQ(cart_mat_grid.nbins(), 50u);
0253     EXPECT_EQ(cart_mat_grid.size(), 50u);
0254
0255     auto x_axis = cart_mat_grid.template get_axis<axis::label::e_x>();
0256     EXPECT_EQ(x_axis.nbins(), 5u);
0257     auto y_axis = cart_mat_grid.template get_axis<axis::label::e_y>();
0258     EXPECT_EQ(y_axis.nbins(), 10u);
0259
0260     for (const auto& mat_slab : cart_mat_grid.all()) {
0261       EXPECT_TRUE(mat_slab.get_material() == tungsten<scalar>());
0262     }
0263   }
0264 }