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   return materialPayload;
0112 }
0113
0114 detray::io::grid_payload<detray::io::material_slab_payload,
0115                          detray::io::material_id>
0116 Acts::DetrayMaterialConverter::convertGridSurfaceMaterial(
0117     const ISurfaceMaterial& material, const Logger& logger) {
0118   detray::io::grid_payload<detray::io::material_slab_payload,
0119                            detray::io::material_id>
0120       materialGrid;
0121
0122   // Check the material types
0123   // (1) homogeneous -> skip
0124   auto homogeneousMaterial =
0125       dynamic_cast<const HomogeneousSurfaceMaterial*>(&material);
0126   if (homogeneousMaterial != nullptr) {
0127     ACTS_DEBUG(
0128         "DetrayMaterialConverter: found homogeneous surface material, ignored "
0129         "as this should be handled by the homogeneous material conversion.");
0130     return materialGrid;
0131   }
0132   // (2) - binned material -> convert into grid structure
0133   auto binnedMaterial = dynamic_cast<const BinnedSurfaceMaterial*>(&material);
0134   if (binnedMaterial != nullptr) {
0135     ACTS_VERBOSE("DetrayMaterialConverter: found binned surface material");
0136
0137     // BinUtility modifications
0138     bool swapped = false;
0139     // Get the bin utility (make a copy as we may modify it)
0140     // Detray expects 2-dimensional grid, currently supported are
0141     // x-y, r-phi, phi-z
0142     BinUtility bUtility = binnedMaterial->binUtility();
0143     // Turn the bin value into a 2D grid
0144     if (bUtility.dimensions() == 1u) {
0145       if (bUtility.binningData()[0u].binvalue == AxisDirection::AxisX) {
0146         // Turn to X-Y
0147         bUtility += BinUtility(1u, std::numeric_limits<float>::lowest(),
0148                                std::numeric_limits<float>::max(),
0149                                BinningOption::closed, AxisDirection::AxisY);
0150       } else if (bUtility.binningData()[0u].binvalue == AxisDirection::AxisY) {
0151         // Turn to X-Y
0152         BinUtility nbUtility(1u, std::numeric_limits<float>::lowest(),
0153                              std::numeric_limits<float>::max(),
0154                              BinningOption::closed, AxisDirection::AxisX);
0155         nbUtility += bUtility;
0156         bUtility = std::move(nbUtility);
0157         swapped = true;
0158       } else if (bUtility.binningData()[0u].binvalue == AxisDirection::AxisR) {
0159         // Turn to R-Phi
0160         bUtility += BinUtility(1u, -std::numbers::pi, std::numbers::pi, closed,
0161                                AxisDirection::AxisPhi);
0162       } else if (bUtility.binningData()[0u].binvalue == AxisDirection::AxisZ) {
0163         // Turn to Phi-Z - swap needed
0164         BinUtility nbUtility(1u, -std::numbers::pi, std::numbers::pi, closed,
0165                              AxisDirection::AxisPhi);
0166         nbUtility += bUtility;
0167         bUtility = std::move(nbUtility);
0168         swapped = true;
0169       } else {
0170         std::invalid_argument("Unsupported binning for Detray");
0171       }
0172     } else if (bUtility.dimensions() == 2u &&
0173                bUtility.binningData()[0u].binvalue == AxisDirection::AxisZ &&
0174                bUtility.binningData()[1u].binvalue == AxisDirection::AxisPhi) {
0175       BinUtility nbUtility(bUtility.binningData()[1u]);
0176       nbUtility += BinUtility{bUtility.binningData()[0u]};
0177       bUtility = std::move(nbUtility);
0178       swapped = true;
0179     }
0180
0181     AxisDirection bVal0 = bUtility.binningData()[0u].binvalue;
0182     AxisDirection bVal1 = bUtility.binningData()[1u].binvalue;
0183
0184     // Translate into grid index type
0185     detray::io::material_id gridIndexType = detray::io::material_id::unknown;
0186     if (bVal0 == AxisDirection::AxisR && bVal1 == AxisDirection::AxisPhi) {
0187       gridIndexType = detray::io::material_id::ring2_map;
0188     } else if (bVal0 == AxisDirection::AxisPhi &&
0189                bVal1 == AxisDirection::AxisZ) {
0190       gridIndexType = detray::io::material_id::concentric_cylinder2_map;
0191     } else if (bVal0 == AxisDirection::AxisX && bVal1 == AxisDirection::AxisY) {
0192       gridIndexType = detray::io::material_id::rectangle2_map;
0193     } else {
0194       std::runtime_error(
0195           "DetrayMaterialConverter: Unsupported binning for Detray");
0196     }
0197
0198     detray::io::typed_link_payload<detray::io::material_id> linkPayload{
0199         gridIndexType, 0u};
0200     materialGrid.grid_link = linkPayload;
0201
0202     // Now convert the (modified) bin utility
0203     for (const auto& bData : bUtility.binningData()) {
0204       auto axisPayload = DetrayConversionUtils::convertBinningData(bData);
0205       materialGrid.axes.push_back(axisPayload);
0206     }
0207
0208     // Convert the material slabs from the material matrix
0209     auto materialMatrix = binnedMaterial->fullMaterial();
0210     for (std::size_t ib1 = 0; ib1 < materialMatrix.size(); ++ib1) {
0211       for (std::size_t ib0 = 0; ib0 < materialMatrix[0u].size(); ++ib0) {
0212         // Translate into a local bin
0213         std::size_t lb0 = swapped ? ib1 : ib0;
0214         std::size_t lb1 = swapped ? ib0 : ib1;
0215         detray::io::material_slab_payload slab =
0216             convertMaterialSlab(materialMatrix[ib1][ib0]);
0217         detray::io::grid_bin_payload<detray::io::material_slab_payload> slabBin{
0218             {static_cast<unsigned int>(lb0), static_cast<unsigned int>(lb1)},
0219             {slab}};
0220         // Fill into the grid
0221         materialGrid.bins.push_back(slabBin);
0222       }
0223     }
0224     return materialGrid;
0225   }
0226
0227   if (dynamic_cast<const Acts::ProtoSurfaceMaterial*>(&material) != nullptr ||
0228       dynamic_cast<const Acts::ProtoGridSurfaceMaterial*>(&material) !=
0229           nullptr) {
0230     ACTS_WARNING(
0231         "DetrayMaterialConverter: ProtoSurfaceMaterial and "
0232         "ProtoGridSurfaceMaterial are not being translated, consider to switch "
0233         "material conversion off.");
0234     return materialGrid;
0235   }
0236
0237   throw std::invalid_argument(
0238       "DetrayMaterialConverter: unknown surface material type detected.");
0239 }
0240
0241 detray::io::detector_grids_payload<detray::io::material_slab_payload,
0242                                    detray::io::material_id>
0243 Acts::DetrayMaterialConverter::convertGridSurfaceMaterial(
0244     const DetrayConversionUtils::Cache& cCache,
0245     const Experimental::Detector& detector, const Logger& logger) {
0246   // The material grid payload
0247   detray::io::detector_grids_payload<detray::io::material_slab_payload,
0248                                      detray::io::material_id>
0249       materialGrids;
0250
0251   using DetrayMaterialGrid =
0252       detray::io::grid_payload<detray::io::material_slab_payload,
0253                                detray::io::material_id>;
0254
0255   // Loop over the volumes in order to assign the right volume links
0256   for (const auto& volume : detector.volumes()) {
0257     // Per volume surface selector
0258     MaterialSurfaceSelector selector;
0259     volume->visitSurfaces(selector);
0260     ACTS_VERBOSE("DetrayMaterialConverter: found "
0261                  << selector.surfaces.size()
0262                  << " surfaces/portals with material in volume "
0263                  << volume->name());
0264     // Find the voluem index first
0265     auto volumeIndex = cCache.volumeLinks.find(volume->geometryId());
0266     if (volumeIndex != cCache.volumeLinks.end()) {
0267       std::vector<DetrayMaterialGrid> volumeMaterialGrids = {};
0268       // Now convert the surfaces
0269       for (const auto& surface : selector.surfaces) {
0270         // Find the surfaces to assign
0271         auto vIndex = cCache.volumeIndex(volume);
0272         auto localSurfaceLinks = cCache.localSurfaceLinks.find(vIndex);
0273         if (localSurfaceLinks != cCache.localSurfaceLinks.end()) {
0274           // Find the surface link
0275           auto surfaceIndices =
0276               localSurfaceLinks->second.equal_range(surface->geometryId());
0277
0278           ACTS_VERBOSE(
0279               "DetrayMaterialConverter: assigning to "
0280               << std::distance(surfaceIndices.first, surfaceIndices.second)
0281               << " surfaces with material in volume " << volume->name());
0282           DetrayMaterialGrid materialGrid =
0283               convertGridSurfaceMaterial(*surface->surfaceMaterial(), logger);
0284           // Ignore if an empty payload is returned
0285           if (materialGrid.axes.empty() || materialGrid.bins.empty()) {
0286             continue;
0287           }
0288
0289           // Loop over the equal range and fill one grid each, this is needed
0290           // as the initial portal could be split into multiple surfaces
0291           for (auto itr = surfaceIndices.first; itr != surfaceIndices.second;
0292                ++itr) {
0293             // Fill the surface index
0294             materialGrid.owner_link =
0295                 detray::io::single_link_payload{itr->second};
0296             // Fill the grid
0297             volumeMaterialGrids.push_back(materialGrid);
0298           }
0299         }
0300       }
0301       // Register the grids of this volume
0302       materialGrids.grids.insert({volumeIndex->second, volumeMaterialGrids});
0303     } else {
0304       ACTS_WARNING(
0305           "DetrayMaterialConverter: volume not found in cache, should not "
0306           "happen.");
0307     }
0308   }
0309   // Return the material grids payload
0310   return materialGrids;
0311 }