Detray/src/DetrayMaterialConverter.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 "Acts/Plugins/Detray/DetrayMaterialConverter.hpp"
0010
0011 #include "Acts/Detector/Detector.hpp"
0012 #include "Acts/Material/BinnedSurfaceMaterial.hpp"
0013 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
0014 #include "Acts/Material/ProtoSurfaceMaterial.hpp"
0015 #include "Acts/Plugins/Detray/DetrayConversionUtils.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Utilities/BinUtility.hpp"
0018 #include "Acts/Utilities/BinningType.hpp"
0019 #include "Acts/Utilities/Helpers.hpp"
0020
0021 #include <numbers>
0022 #include <stdexcept>
0023
0024 namespace {
0025
0026 struct MaterialSurfaceSelector {
0027   std::vector<const Acts::Surface*> surfaces = {};
0028
0029   /// @param surface is the test surface
0030   void operator()(const Acts::Surface* surface) {
0031     if (surface->surfaceMaterial() != nullptr &&
0032         !rangeContainsValue(surfaces, surface)) {
0033       surfaces.push_back(surface);
0034     }
0035   }
0036 };
0037
0038 }  // namespace
0039
0040 detray::io::material_slab_payload
0041 Acts::DetrayMaterialConverter::convertMaterialSlab(
0042     const MaterialSlab& materialSlab) {
0043   detray::io::material_slab_payload slab;
0044   // Fill the material parameters and the thickness
0045   const auto& material = materialSlab.material();
0046   slab.thickness = materialSlab.thickness();
0047   slab.mat = detray::io::material_payload{
0048       {material.X0(), material.L0(), material.Ar(), material.Z(),
0049        material.massDensity(), material.molarDensity(), 0.}};
0050   slab.type = detray::io::material_id::slab;
0051   return slab;
0052 }
0053
0054 detray::io::detector_homogeneous_material_payload
0055 Acts::DetrayMaterialConverter::convertHomogeneousSurfaceMaterial(
0056     const DetrayConversionUtils::Cache& cCache,
0057     const Experimental::Detector& detector, const Logger& logger) {
0058   detray::io::detector_homogeneous_material_payload materialPayload;
0059
0060   for (const auto volume : detector.volumes()) {
0061     auto volumeIndex = cCache.volumeLinks.find(volume->geometryId());
0062     if (volumeIndex != cCache.volumeLinks.end()) {
0063       // The volume material payload & its link
0064       detray::io::material_volume_payload volumePayload;
0065       detray::io::single_link_payload volumeLink;
0066       volumeLink.link = volumeIndex->second;
0067       volumePayload.volume_link = volumeLink;
0068       // Now run through surfaces and portals to find the material
0069       MaterialSurfaceSelector selector;
0070       volume->visitSurfaces(selector);
0071       ACTS_DEBUG("DetrayMaterialConverter: found "
0072                  << selector.surfaces.size()
0073                  << " surfaces/portals with material in volume "
0074                  << volume->name());
0075       for (const auto surface : selector.surfaces) {
0076         const auto* surfaceMaterial = surface->surfaceMaterial();
0077         auto homogeneousMaterial =
0078             dynamic_cast<const HomogeneousSurfaceMaterial*>(surfaceMaterial);
0079         if (homogeneousMaterial != nullptr) {
0080           // Convert the material slab
0081           auto materialSlab = homogeneousMaterial->materialSlab();
0082           detray::io::material_slab_payload slabPayload =
0083               convertMaterialSlab(materialSlab);
0084           // Find the surfaces to assign
0085           auto vIndex = cCache.volumeIndex(volume);
0086           auto localSurfaceLinks = cCache.localSurfaceLinks.find(vIndex);
0087           if (localSurfaceLinks != cCache.localSurfaceLinks.end()) {
0088             // Find the surface link
0089             auto surfaceIndices =
0090                 localSurfaceLinks->second.equal_range(surface->geometryId());
0091             // Loop over the equal range and fill one grid each, this is needed
0092             // as the initial portal could be split into multiple surfaces
0093             for (auto itr = surfaceIndices.first; itr != surfaceIndices.second;
0094                  ++itr) {
0095               // Make an identified link copy for every matching surface
0096               slabPayload.surface.link = itr->second;
0097               volumePayload.mat_slabs.push_back(slabPayload);
0098             }
0099           } else {
0100             ACTS_WARNING(
0101                 "DetrayMaterialConverter: no local surface links found");
0102           }
0103         }
0104       }
0105       materialPayload.volumes.push_back(volumePayload);
0106     } else {
0107       ACTS_WARNING("DetrayMaterialConverter: volume " << volume->name()
0108                                                       << " not found in cache");
0109     }
0110   }
0111
0112   return materialPayload;
0113 }
0114
0115 detray::io::grid_payload<detray::io::material_slab_payload,
0116                          detray::io::material_id>
0117 Acts::DetrayMaterialConverter::convertGridSurfaceMaterial(
0118     const ISurfaceMaterial& material, const Logger& logger) {
0119   detray::io::grid_payload<detray::io::material_slab_payload,
0120                            detray::io::material_id>
0121       materialGrid;
0122
0123   // Check the material types
0124   // (1) homogeneous -> skip
0125   auto homogeneousMaterial =
0126       dynamic_cast<const HomogeneousSurfaceMaterial*>(&material);
0127   if (homogeneousMaterial != nullptr) {
0128     ACTS_DEBUG(
0129         "DetrayMaterialConverter: found homogeneous surface material, ignored "
0130         "as this should be handled by the homogeneous material conversion.");
0131     return materialGrid;
0132   }
0133   // (2) - binned material -> convert into grid structure
0134   auto binnedMaterial = dynamic_cast<const BinnedSurfaceMaterial*>(&material);
0135   if (binnedMaterial != nullptr) {
0136     ACTS_VERBOSE("DetrayMaterialConverter: found binned surface material");
0137
0138     // BinUtility modifications
0139     bool swapped = false;
0140     // Get the bin utility (make a copy as we may modify it)
0141     // Detray expects 2-dimensional grid, currently supported are
0142     // x-y, r-phi, phi-z
0143     BinUtility bUtility = binnedMaterial->binUtility();
0144     // Turn the bin value into a 2D grid
0145     if (bUtility.dimensions() == 1u) {
0146       if (bUtility.binningData()[0u].binvalue == AxisDirection::AxisX) {
0147         // Turn to X-Y
0148         bUtility += BinUtility(1u, std::numeric_limits<float>::lowest(),
0149                                std::numeric_limits<float>::max(),
0150                                BinningOption::closed, AxisDirection::AxisY);
0151       } else if (bUtility.binningData()[0u].binvalue == AxisDirection::AxisY) {
0152         // Turn to X-Y
0153         BinUtility nbUtility(1u, std::numeric_limits<float>::lowest(),
0154                              std::numeric_limits<float>::max(),
0155                              BinningOption::closed, AxisDirection::AxisX);
0156         nbUtility += bUtility;
0157         bUtility = std::move(nbUtility);
0158         swapped = true;
0159       } else if (bUtility.binningData()[0u].binvalue == AxisDirection::AxisR) {
0160         // Turn to R-Phi
0161         bUtility += BinUtility(1u, -std::numbers::pi, std::numbers::pi, closed,
0162                                AxisDirection::AxisPhi);
0163       } else if (bUtility.binningData()[0u].binvalue == AxisDirection::AxisZ) {
0164         // Turn to Phi-Z - swap needed
0165         BinUtility nbUtility(1u, -std::numbers::pi, std::numbers::pi, closed,
0166                              AxisDirection::AxisPhi);
0167         nbUtility += bUtility;
0168         bUtility = std::move(nbUtility);
0169         swapped = true;
0170       } else {
0171         std::invalid_argument("Unsupported binning for Detray");
0172       }
0173     } else if (bUtility.dimensions() == 2u &&
0174                bUtility.binningData()[0u].binvalue == AxisDirection::AxisZ &&
0175                bUtility.binningData()[1u].binvalue == AxisDirection::AxisPhi) {
0176       BinUtility nbUtility(bUtility.binningData()[1u]);
0177       nbUtility += BinUtility{bUtility.binningData()[0u]};
0178       bUtility = std::move(nbUtility);
0179       swapped = true;
0180     }
0181
0182     AxisDirection bVal0 = bUtility.binningData()[0u].binvalue;
0183     AxisDirection bVal1 = bUtility.binningData()[1u].binvalue;
0184
0185     // Translate into grid index type
0186     detray::io::material_id gridIndexType = detray::io::material_id::unknown;
0187     if (bVal0 == AxisDirection::AxisR && bVal1 == AxisDirection::AxisPhi) {
0188       gridIndexType = detray::io::material_id::ring2_map;
0189     } else if (bVal0 == AxisDirection::AxisPhi &&
0190                bVal1 == AxisDirection::AxisZ) {
0191       gridIndexType = detray::io::material_id::concentric_cylinder2_map;
0192     } else if (bVal0 == AxisDirection::AxisX && bVal1 == AxisDirection::AxisY) {
0193       gridIndexType = detray::io::material_id::rectangle2_map;
0194     } else {
0195       std::runtime_error(
0196           "DetrayMaterialConverter: Unsupported binning for Detray");
0197     }
0198
0199     detray::io::typed_link_payload<detray::io::material_id> linkPayload{
0200         gridIndexType, 0u};
0201     materialGrid.grid_link = linkPayload;
0202
0203     // Now convert the (modified) bin utility
0204     for (const auto& bData : bUtility.binningData()) {
0205       auto axisPayload = DetrayConversionUtils::convertBinningData(bData);
0206       materialGrid.axes.push_back(axisPayload);
0207     }
0208
0209     // Convert the material slabs from the material matrix
0210     auto materialMatrix = binnedMaterial->fullMaterial();
0211     for (std::size_t ib1 = 0; ib1 < materialMatrix.size(); ++ib1) {
0212       for (std::size_t ib0 = 0; ib0 < materialMatrix[0u].size(); ++ib0) {
0213         // Translate into a local bin
0214         std::size_t lb0 = swapped ? ib1 : ib0;
0215         std::size_t lb1 = swapped ? ib0 : ib1;
0216         detray::io::material_slab_payload slab =
0217             convertMaterialSlab(materialMatrix[ib1][ib0]);
0218         detray::io::grid_bin_payload<detray::io::material_slab_payload> slabBin{
0219             {static_cast<unsigned int>(lb0), static_cast<unsigned int>(lb1)},
0220             {slab}};
0221         // Fill into the grid
0222         materialGrid.bins.push_back(slabBin);
0223       }
0224     }
0225     return materialGrid;
0226   }
0227
0228   if (dynamic_cast<const Acts::ProtoSurfaceMaterial*>(&material) != nullptr ||
0229       dynamic_cast<const Acts::ProtoGridSurfaceMaterial*>(&material) !=
0230           nullptr) {
0231     ACTS_WARNING(
0232         "DetrayMaterialConverter: ProtoSurfaceMaterial and "
0233         "ProtoGridSurfaceMaterial are not being translated, consider to switch "
0234         "material conversion off.");
0235     return materialGrid;
0236   }
0237
0238   throw std::invalid_argument(
0239       "DetrayMaterialConverter: unknown surface material type detected.");
0240 }
0241
0242 detray::io::detector_grids_payload<detray::io::material_slab_payload,
0243                                    detray::io::material_id>
0244 Acts::DetrayMaterialConverter::convertGridSurfaceMaterial(
0245     const DetrayConversionUtils::Cache& cCache,
0246     const Experimental::Detector& detector, const Logger& logger) {
0247   // The material grid payload
0248   detray::io::detector_grids_payload<detray::io::material_slab_payload,
0249                                      detray::io::material_id>
0250       materialGrids;
0251
0252   using DetrayMaterialGrid =
0253       detray::io::grid_payload<detray::io::material_slab_payload,
0254                                detray::io::material_id>;
0255
0256   // Loop over the volumes in order to assign the right volume links
0257   for (const auto& volume : detector.volumes()) {
0258     // Per volume surface selector
0259     MaterialSurfaceSelector selector;
0260     volume->visitSurfaces(selector);
0261     ACTS_VERBOSE("DetrayMaterialConverter: found "
0262                  << selector.surfaces.size()
0263                  << " surfaces/portals with material in volume "
0264                  << volume->name());
0265     // Find the voluem index first
0266     auto volumeIndex = cCache.volumeLinks.find(volume->geometryId());
0267     if (volumeIndex != cCache.volumeLinks.end()) {
0268       std::vector<DetrayMaterialGrid> volumeMaterialGrids = {};
0269       // Now convert the surfaces
0270       for (const auto& surface : selector.surfaces) {
0271         // Find the surfaces to assign
0272         auto vIndex = cCache.volumeIndex(volume);
0273         auto localSurfaceLinks = cCache.localSurfaceLinks.find(vIndex);
0274         if (localSurfaceLinks != cCache.localSurfaceLinks.end()) {
0275           // Find the surface link
0276           auto surfaceIndices =
0277               localSurfaceLinks->second.equal_range(surface->geometryId());
0278
0279           ACTS_VERBOSE(
0280               "DetrayMaterialConverter: assigning to "
0281               << std::distance(surfaceIndices.first, surfaceIndices.second)
0282               << " surfaces with material in volume " << volume->name());
0283           DetrayMaterialGrid materialGrid =
0284               convertGridSurfaceMaterial(*surface->surfaceMaterial(), logger);
0285           // Ignore if an empty payload is returned
0286           if (materialGrid.axes.empty() || materialGrid.bins.empty()) {
0287             continue;
0288           }
0289
0290           // Loop over the equal range and fill one grid each, this is needed
0291           // as the initial portal could be split into multiple surfaces
0292           for (auto itr = surfaceIndices.first; itr != surfaceIndices.second;
0293                ++itr) {
0294             // Fill the surface index
0295             materialGrid.owner_link =
0296                 detray::io::single_link_payload{itr->second};
0297             // Fill the grid
0298             volumeMaterialGrids.push_back(materialGrid);
0299           }
0300         }
0301       }
0302       // Register the grids of this volume
0303       materialGrids.grids.insert({volumeIndex->second, volumeMaterialGrids});
0304     } else {
0305       ACTS_WARNING(
0306           "DetrayMaterialConverter: volume not found in cache, should not "
0307           "happen.");
0308     }
0309   }
0310   // Return the material grids payload
0311   return materialGrids;
0312 }