Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:24:17 

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 "ActsPlugins/Detray/DetrayPayloadConverter.hpp"
0010 
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/CompositePortalLink.hpp"
0013 #include "Acts/Geometry/GeometryContext.hpp"
0014 #include "Acts/Geometry/GeometryIdentifier.hpp"
0015 #include "Acts/Geometry/GridPortalLink.hpp"
0016 #include "Acts/Geometry/TrackingGeometry.hpp"
0017 #include "Acts/Geometry/TrivialPortalLink.hpp"
0018 #include "Acts/Geometry/VolumeBounds.hpp"
0019 #include "Acts/Material/ISurfaceMaterial.hpp"
0020 #include "Acts/Navigation/INavigationPolicy.hpp"
0021 #include "Acts/Surfaces/AnnulusBounds.hpp"
0022 #include "Acts/Surfaces/CylinderBounds.hpp"
0023 #include "Acts/Surfaces/RadialBounds.hpp"
0024 #include "Acts/Surfaces/RectangleBounds.hpp"
0025 #include "Acts/Surfaces/Surface.hpp"
0026 #include "Acts/Surfaces/SurfaceBounds.hpp"
0027 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0028 #include "Acts/Utilities/Helpers.hpp"
0029 #include "Acts/Utilities/Logger.hpp"
0030 #include "ActsPlugins/Detray/DetrayConversionUtils.hpp"
0031 
0032 #include <optional>
0033 
0034 #include <detray/geometry/shapes/annulus2D.hpp>
0035 #include <detray/geometry/shapes/concentric_cylinder2D.hpp>
0036 #include <detray/geometry/shapes/cylinder2D.hpp>
0037 #include <detray/geometry/shapes/rectangle2D.hpp>
0038 #include <detray/geometry/shapes/ring2D.hpp>
0039 #include <detray/geometry/shapes/trapezoid2D.hpp>
0040 #include <detray/io/frontend/definitions.hpp>
0041 #include <detray/io/frontend/payloads.hpp>
0042 
0043 using namespace Acts;
0044 
0045 namespace ActsPlugins {
0046 
0047 DetrayPayloadConverter::DetrayPayloadConverter(
0048     const Config& config, std::unique_ptr<const Logger> logger)
0049     : m_cfg(config), m_logger(std::move(logger)) {}
0050 
0051 namespace {
0052 using enum detray::io::shape_id;
0053 
0054 detray::io::mask_payload convertBounds(const AnnulusBounds& annulus) {
0055   using enum detray::annulus2D::boundaries;
0056   using enum AnnulusBounds::BoundValues;
0057 
0058   detray::io::mask_payload payload;
0059   payload.shape = annulus2;
0060   payload.boundaries.resize(e_size);
0061   payload.boundaries.at(e_min_r) = annulus.get(eMinR);
0062   payload.boundaries.at(e_max_r) = annulus.get(eMaxR);
0063   payload.boundaries.at(e_min_phi_rel) = annulus.get(eMinPhiRel);
0064   payload.boundaries.at(e_max_phi_rel) = annulus.get(eMaxPhiRel);
0065   payload.boundaries.at(e_average_phi) = annulus.get(eAveragePhi);
0066   payload.boundaries.at(e_shift_x) = annulus.get(eOriginX);
0067   payload.boundaries.at(e_shift_y) = annulus.get(eOriginY);
0068 
0069   return payload;
0070 }
0071 
0072 detray::io::mask_payload convertBounds(const RectangleBounds& rectangle) {
0073   using enum RectangleBounds::BoundValues;
0074   using enum detray::rectangle2D::boundaries;
0075 
0076   detray::io::mask_payload payload;
0077   payload.shape = rectangle2;
0078   payload.boundaries.resize(e_size);
0079 
0080   double minX = rectangle.get(eMinX);
0081   double maxX = rectangle.get(eMaxX);
0082   double minY = rectangle.get(eMinY);
0083   double maxY = rectangle.get(eMaxY);
0084 
0085   if (minX != -maxX || minY != -maxY) {
0086     throw std::runtime_error(
0087         "Rectangle bounds are not symmetric, detray cannot handle this");
0088   }
0089 
0090   payload.boundaries.at(e_half_x) = maxX;
0091   payload.boundaries.at(e_half_y) = maxY;
0092 
0093   return payload;
0094 }
0095 
0096 detray::io::mask_payload convertBounds(const CylinderBounds& cylinder,
0097                                        detray::io::shape_id shape) {
0098   using enum CylinderBounds::BoundValues;
0099 
0100   detray::io::mask_payload payload;
0101   payload.shape = shape;
0102   if (shape == portal_cylinder2) {
0103     using enum detray::concentric_cylinder2D::boundaries;
0104     payload.boundaries.resize(e_size);
0105     payload.boundaries.at(e_r) = cylinder.get(eR);
0106     double hlZ = cylinder.get(eHalfLengthZ);
0107     payload.boundaries.at(e_lower_z) = -hlZ;
0108     payload.boundaries.at(e_upper_z) = hlZ;
0109   } else {
0110     using enum detray::cylinder2D::boundaries;
0111     payload.boundaries.resize(e_size);
0112     payload.boundaries.at(e_r) = cylinder.get(eR);
0113     double hlZ = cylinder.get(eHalfLengthZ);
0114     payload.boundaries.at(e_lower_z) = -hlZ;
0115     payload.boundaries.at(e_upper_z) = hlZ;
0116   }
0117   return payload;
0118 }
0119 
0120 detray::io::mask_payload convertBounds(const TrapezoidBounds& trapezoid) {
0121   using enum TrapezoidBounds::BoundValues;
0122   using enum detray::trapezoid2D::boundaries;
0123 
0124   detray::io::mask_payload payload;
0125   payload.shape = trapezoid2;
0126 
0127   payload.boundaries.resize(e_size);
0128   payload.boundaries.at(e_half_length_0) = trapezoid.get(eHalfLengthXnegY);
0129   payload.boundaries.at(e_half_length_1) = trapezoid.get(eHalfLengthXposY);
0130   payload.boundaries.at(e_half_length_2) = trapezoid.get(eHalfLengthY);
0131   payload.boundaries.at(e_divisor) = 1 / (2 * trapezoid.get(eHalfLengthY));
0132 
0133   return payload;
0134 }
0135 
0136 detray::io::mask_payload convertBounds(const RadialBounds& radial) {
0137   using enum RadialBounds::BoundValues;
0138   using enum detray::ring2D::boundaries;
0139 
0140   if (!radial.coversFullAzimuth()) {
0141     throw std::runtime_error(
0142         "Radial bounds do not cover full azimuth, detray cannot handle this");
0143   }
0144 
0145   if (radial.get(eAveragePhi) != 0.) {
0146     throw std::runtime_error(
0147         "Radial bounds have an average phi, detray cannot handle this");
0148   }
0149 
0150   detray::io::mask_payload payload;
0151   payload.shape = ring2;
0152 
0153   payload.boundaries.resize(e_size);
0154   payload.boundaries.at(e_inner_r) = radial.get(eMinR);
0155   payload.boundaries.at(e_outer_r) = radial.get(eMaxR);
0156 
0157   return payload;
0158 }
0159 
0160 }  // namespace
0161 
0162 detray::io::mask_payload DetrayPayloadConverter::convertMask(
0163     const SurfaceBounds& bounds, bool forPortal) {
0164   detray::io::mask_payload payload;
0165 
0166   switch (bounds.type()) {
0167     using enum SurfaceBounds::BoundsType;
0168     using enum detray::io::shape_id;
0169     case eAnnulus:
0170       payload = convertBounds(dynamic_cast<const AnnulusBounds&>(bounds));
0171       break;
0172     case eRectangle:
0173       payload = convertBounds(dynamic_cast<const RectangleBounds&>(bounds));
0174       break;
0175     case eCylinder:
0176       payload = convertBounds(dynamic_cast<const CylinderBounds&>(bounds),
0177                               forPortal ? portal_cylinder2 : cylinder2);
0178       break;
0179     case eTrapezoid:
0180       payload = convertBounds(dynamic_cast<const TrapezoidBounds&>(bounds));
0181       break;
0182     case eDisc:
0183       if (auto* radial = dynamic_cast<const RadialBounds*>(&bounds);
0184           radial != nullptr) {
0185         payload = convertBounds(*radial);
0186       } else {
0187         throw std::runtime_error(
0188             "Disc bounds type but not radial bounds currently unsupported");
0189       }
0190       break;
0191     default:
0192       payload.shape = unknown;
0193       break;
0194   }
0195 
0196   return payload;
0197 }
0198 
0199 detray::io::surface_payload DetrayPayloadConverter::convertSurface(
0200     const GeometryContext& gctx, const Surface& surface, bool portal) const {
0201   detray::io::surface_payload payload;
0202 
0203   payload.transform =
0204       DetrayConversionUtils::convertTransform(surface.transform(gctx));
0205   payload.source = surface.geometryId().value();
0206   payload.barcode = std::nullopt;
0207 
0208   bool isSensitive = false;
0209   if (m_cfg.sensitiveStrategy ==
0210       DetrayPayloadConverter::Config::SensitiveStrategy::Identifier) {
0211     isSensitive = surface.geometryId().sensitive() > 0;
0212   } else {
0213     isSensitive = surface.associatedDetectorElement() != nullptr;
0214   }
0215 
0216   if (portal) {
0217     payload.type = detray::surface_id::e_portal;
0218   } else {
0219     payload.type = isSensitive ? detray::surface_id::e_sensitive
0220                                : detray::surface_id::e_passive;
0221   }
0222   payload.masks = {convertMask(surface.bounds(), portal)};
0223   return payload;
0224 }
0225 
0226 detray::io::volume_payload DetrayPayloadConverter::convertVolume(
0227     const TrackingVolume& volume) const {
0228   detray::io::volume_payload payload;
0229   payload.transform =
0230       DetrayConversionUtils::convertTransform(volume.transform());
0231   payload.name = volume.volumeName();
0232   switch (volume.volumeBounds().type()) {
0233     using enum VolumeBounds::BoundsType;
0234     using enum detray::volume_id;
0235     case eCylinder:
0236       payload.type = e_cylinder;
0237       break;
0238     case eCuboid:
0239       payload.type = e_cuboid;
0240       break;
0241     case eTrapezoid:
0242       payload.type = e_trapezoid;
0243       break;
0244     case eCone:
0245       payload.type = e_cone;
0246       break;
0247     default:
0248       payload.type = e_unknown;
0249       break;
0250   }
0251   return payload;
0252 }
0253 
0254 namespace {
0255 std::vector<const TrivialPortalLink*> decomposeToTrivials(
0256     const PortalLinkBase& link, const Logger& logger) {
0257   std::vector<const TrivialPortalLink*> trivials;
0258 
0259   if (auto* trivial = dynamic_cast<const TrivialPortalLink*>(&link);
0260       trivial != nullptr) {
0261     trivials.push_back(trivial);
0262   } else if (auto* composite = dynamic_cast<const CompositePortalLink*>(&link);
0263              composite != nullptr) {
0264     ACTS_VERBOSE("Converting composite portal link with "
0265                  << composite->links().size() << " sub-links");
0266     for (const auto& subLink : composite->links()) {
0267       const auto* subTrivial = dynamic_cast<const TrivialPortalLink*>(&subLink);
0268 
0269       if (subTrivial == nullptr) {
0270         throw std::runtime_error(
0271             "Composite portal link contains non-trivial portal links");
0272       } else {
0273         trivials.push_back(subTrivial);
0274       }
0275     }
0276   } else if (auto* grid = dynamic_cast<const GridPortalLink*>(&link);
0277              grid != nullptr) {
0278     ACTS_VERBOSE("Converting grid portal link with "
0279                  << grid->artifactPortalLinks().size() << " link artifacts");
0280     for (const auto& artifact : grid->artifactPortalLinks()) {
0281       trivials.push_back(&artifact);
0282     }
0283   } else {
0284     throw std::runtime_error(
0285         "Unknown portal link type, detray cannot handle this");
0286   }
0287 
0288   return trivials;
0289 }
0290 /// Check compatibility between ACTS surface type and detray material_id
0291 /// @returns pair<is_compatible, expected_material_id>
0292 std::pair<bool, detray::io::material_id> isGridMaterialCompatible(
0293     Surface::SurfaceType surfaceType, detray::io::material_id materialId) {
0294   using enum Surface::SurfaceType;
0295   using enum detray::io::material_id;
0296 
0297   switch (surfaceType) {
0298     case Cylinder:
0299       return {materialId == concentric_cylinder2_map, concentric_cylinder2_map};
0300     case Disc:
0301       return {materialId == ring2_map, ring2_map};
0302     case Plane:
0303       return {materialId == rectangle2_map, rectangle2_map};
0304     default:
0305       // For other surface types, we don't have specific checks yet
0306       return {true, unknown};
0307   }
0308 }
0309 
0310 }  // namespace
0311 
0312 void DetrayPayloadConverter::handlePortalLink(
0313     const GeometryContext& gctx, const TrackingVolume& volume,
0314     detray::io::volume_payload& volPayload,
0315     const std::function<std::size_t(const TrackingVolume*)>& volumeLookup,
0316     std::unordered_map<const Surface*, std::size_t>& surfaceIndices,
0317     const PortalLinkBase& link) const {
0318   std::vector<const TrivialPortalLink*> trivials =
0319       decomposeToTrivials(link, logger());
0320 
0321   // If ANY of the trivials point at the current volume, we don't handle this
0322   // portal link at all, otherwise we would get a self-referencing volume link
0323 
0324   if (std::ranges::any_of(
0325           trivials, [&](const auto* t) { return &t->volume() == &volume; })) {
0326     ACTS_VERBOSE("At least one trivial link points at this volume ("
0327                  << volume.volumeName() << ") => skipping");
0328     return;
0329   }
0330 
0331   for (const auto* trivial : trivials) {
0332     ACTS_VERBOSE("Converting trivial portal link registered to volume "
0333                  << volume.volumeName());
0334     ACTS_VERBOSE(
0335         "Portal link surface is: " << trivial->surface().toStream(gctx));
0336     if (&trivial->volume() == &volume) {
0337       ACTS_VERBOSE("~> points at this volume (" << volume.volumeName()
0338                                                 << ") => skipping");
0339       return;
0340     }
0341 
0342     ACTS_VERBOSE("~> points at different volume ("
0343                  << trivial->volume().volumeName()
0344                  << ") => adding link to this volume (" << volume.volumeName()
0345                  << ")");
0346 
0347     // add the surface (including mask first)
0348     auto& srfPayload = volPayload.surfaces.emplace_back(
0349         convertSurface(gctx, trivial->surface(), true));
0350     srfPayload.index_in_coll = volPayload.surfaces.size() - 1;
0351 
0352     // lookup the target volume index (we already converted this)
0353     ACTS_VERBOSE("Target volume index for "
0354                  << trivial->volume().volumeName() << ": "
0355                  << volumeLookup(&trivial->volume()));
0356     std::size_t targetVolumeIndex = volumeLookup(&trivial->volume());
0357     srfPayload.masks.at(0).volume_link.link = targetVolumeIndex;
0358     surfaceIndices[&trivial->surface()] = srfPayload.index_in_coll.value();
0359   }
0360 }
0361 
0362 void DetrayPayloadConverter::makeEndOfWorld(
0363     const GeometryContext& gctx, detray::io::volume_payload& volPayload,
0364     std::unordered_map<const Surface*, std::size_t>& surfaceIndices,
0365     const Surface& surface) const {
0366   ACTS_VERBOSE("Adding end of world surface");
0367   auto& srfPayload =
0368       volPayload.surfaces.emplace_back(convertSurface(gctx, surface, true));
0369   srfPayload.index_in_coll = volPayload.surfaces.size() - 1;
0370 
0371   // Marker for end of world is MAX
0372   srfPayload.masks.at(0).volume_link.link =
0373       std::numeric_limits<std::size_t>::max();
0374 
0375   surfaceIndices[&surface] = srfPayload.index_in_coll.value();
0376 }
0377 
0378 void DetrayPayloadConverter::handlePortal(
0379     const GeometryContext& gctx, const TrackingVolume& volume,
0380     detray::io::volume_payload& volPayload,
0381     const std::function<std::size_t(const TrackingVolume*)>& volumeLookup,
0382     std::unordered_map<const Surface*, std::size_t>& surfaceIndices,
0383     const Portal& portal) const {
0384   auto* lAlong = portal.getLink(Direction::AlongNormal());
0385   auto* lOpposite = portal.getLink(Direction::OppositeNormal());
0386 
0387   if (lAlong == nullptr && lOpposite == nullptr) {
0388     // Sanity check: this shouldn't happen
0389     throw std::runtime_error("Portal link is not symmetric");
0390   }
0391 
0392   if (lAlong != nullptr) {
0393     handlePortalLink(gctx, volume, volPayload, volumeLookup, surfaceIndices,
0394                      *lAlong);
0395   } else {
0396     // can't both be nullptr
0397     assert(lOpposite != nullptr);
0398     makeEndOfWorld(gctx, volPayload, surfaceIndices, lOpposite->surface());
0399   }
0400 
0401   if (lOpposite != nullptr) {
0402     handlePortalLink(gctx, volume, volPayload, volumeLookup, surfaceIndices,
0403                      *lOpposite);
0404   } else {
0405     // can't both be nullptr
0406     assert(lAlong != nullptr);
0407     makeEndOfWorld(gctx, volPayload, surfaceIndices, lAlong->surface());
0408   }
0409 }
0410 
0411 namespace {
0412 
0413 constexpr static detray::io::material_slab_payload s_dummyMaterialSlab{
0414     .type = detray::io::material_id::slab,
0415     .index_in_coll = std::numeric_limits<std::size_t>::max(),
0416     .thickness = 42,
0417     .mat = {42, 42, 42, 42, 42, 42, 42},
0418 };
0419 
0420 }  // namespace
0421 
0422 std::pair<std::vector<detray::io::grid_payload<
0423               detray::io::material_slab_payload, detray::io::material_id>>,
0424           detray::io::material_volume_payload>
0425 DetrayPayloadConverter::convertMaterial(
0426     const TrackingVolume& volume,
0427     const std::unordered_map<const Surface*, std::size_t>& surfaceIndices,
0428     detray::io::volume_payload& volPayload) const {
0429   ACTS_DEBUG("Converting material for volume " << volume.volumeName());
0430   std::vector<detray::io::grid_payload<detray::io::material_slab_payload,
0431                                        detray::io::material_id>>
0432       grids;
0433   detray::io::material_volume_payload homogeneous;
0434   homogeneous.volume_link.link = volPayload.index.link;
0435 
0436   // @HACK: Detray does not like homoegeneous material only on SOME surfaces.
0437   ACTS_INFO("Adding dummy material slabs to homogeneous collection for "
0438             << volPayload.surfaces.size()
0439             << " surfaces (detray "
0440                "hack)");
0441   for (const auto& surface : volPayload.surfaces) {
0442     auto& slabPayload = homogeneous.mat_slabs.emplace_back(s_dummyMaterialSlab);
0443     slabPayload.index_in_coll = surface.index_in_coll.value();
0444     slabPayload.surface.link = surface.index_in_coll.value();
0445   }
0446 
0447   std::map<std::size_t, const ISurfaceMaterial*> srfIdxToMaterial;
0448 
0449   auto assignMaterial = [&](const ISurfaceMaterial* material,
0450                             DetraySurfaceMaterial& detrayMaterial,
0451                             std::size_t srfIdx, const Surface& surface) {
0452     auto handleHomogeneous =
0453         [&](const detray::io::material_slab_payload& slab) {
0454           // Given the pseudo slabs from before, we need to either update an
0455           // existing slab or create a new one.
0456 
0457           auto it = std::ranges::find_if(
0458               homogeneous.mat_slabs, [srfIdx](const auto& matslab) {
0459                 return matslab.surface.link == srfIdx;
0460               });
0461 
0462           if (it != homogeneous.mat_slabs.end()) {
0463             ACTS_VERBOSE("Adding slab to homogeneous material for surface "
0464                          << srfIdx);
0465             auto& targetSlab = *it;
0466             targetSlab = slab;
0467             targetSlab.index_in_coll = srfIdx;
0468             targetSlab.surface.link = srfIdx;
0469           } else {
0470             ACTS_VERBOSE("Updating slab in homogeneous material for surface "
0471                          << srfIdx);
0472             homogeneous.mat_slabs.emplace_back(slab);
0473             homogeneous.mat_slabs.back().index_in_coll = srfIdx;
0474             homogeneous.mat_slabs.back().surface.link = srfIdx;
0475           }
0476 
0477           auto sit = srfIdxToMaterial.find(srfIdx);
0478           if (sit != srfIdxToMaterial.end() && sit->second != material) {
0479             ACTS_ERROR("Surface "
0480                        << srfIdx
0481                        << " already has a different material assigned");
0482             throw std::runtime_error("Material mismatch for surface");
0483           }
0484         };
0485 
0486     auto handleGrid = [&](const detray::io::grid_payload<
0487                           detray::io::material_slab_payload,
0488                           detray::io::material_id>& grid) {
0489       ACTS_DEBUG("Assigning grid material to surface " << srfIdx);
0490       auto it = srfIdxToMaterial.find(srfIdx);
0491 
0492       if (it != srfIdxToMaterial.end()) {
0493         // Surface already has some material assigned
0494         if (it->second != material) {
0495           ACTS_ERROR("Surface "
0496                      << srfIdx << " already has a different material assigned");
0497           throw std::runtime_error("Material mismatch for surface");
0498         }
0499 
0500         // It's the same material again, we just skip adding a duplicate
0501         return;
0502       }
0503 
0504       // @TODO: Add a consistency check between the surface type and the
0505       //        axis types: detray's consistency check will find this but
0506       //        the error message is not trivial. In this location, we can
0507       //        print out the surface id and guide debugging!
0508 
0509       // Check compatibility between surface type and grid material_id
0510       auto [isCompatible, expectedId] =
0511           isGridMaterialCompatible(surface.type(), grid.grid_link.type);
0512 
0513       if (!isCompatible) {
0514         ACTS_ERROR("Grid material compatibility error for surface "
0515                    << surface.geometryId() << " (index " << srfIdx << "): "
0516                    << "Surface type " << surface.type()
0517                    << " is incompatible with material_id '"
0518                    << toUnderlying(grid.grid_link.type) << "'. Expected '"
0519                    << toUnderlying(expectedId)
0520                    << "'. This indicates the BinUtility axis definition "
0521                    << "does not match the surface geometry.");
0522         throw std::runtime_error(
0523             "Grid material has incompatible axis definition for surface type");
0524       }
0525 
0526       grids.emplace_back(grid);
0527       grids.back().owner_link.link = srfIdx;
0528     };
0529 
0530     std::visit(overloaded{handleHomogeneous, handleGrid}, detrayMaterial);
0531 
0532     srfIdxToMaterial[srfIdx] = material;
0533   };
0534 
0535   auto printSurfaceInfo = [&](DetraySurfaceMaterial& detrayMaterial,
0536                               const Surface& surface) {
0537     auto handleHomogeneous = [&](const detray::io::material_slab_payload&) {
0538       ACTS_VERBOSE("Surface " << surface.geometryId()
0539                               << " has homogeneous material");
0540     };
0541 
0542     auto handleGrid =
0543         [&](const detray::io::grid_payload<detray::io::material_slab_payload,
0544                                            detray::io::material_id>&) {
0545           ACTS_VERBOSE("Surface " << surface.geometryId()
0546                                   << " has grid material");
0547         };
0548     std::visit(overloaded{handleHomogeneous, handleGrid}, detrayMaterial);
0549   };
0550 
0551   ACTS_VERBOSE("Looping over " << volume.surfaces().size()
0552                                << " surfaces in volume " << volPayload.name);
0553   for (const auto& surface : volume.surfaces()) {
0554     auto srfIt = surfaceIndices.find(&surface);
0555 
0556     if (srfIt == surfaceIndices.end()) {
0557       ACTS_ERROR("Surface " << surface.geometryId().value()
0558                             << " not found in volume " << volPayload.name
0559                             << ". This is a bug in the conversion.");
0560       throw std::runtime_error("Surface not found in volume");
0561     }
0562 
0563     std::size_t srfIdx = srfIt->second;
0564 
0565     if (surface.surfaceMaterial() == nullptr) {
0566       continue;
0567     }
0568 
0569     std::optional detrayMaterial =
0570         m_cfg.convertSurfaceMaterial(*surface.surfaceMaterial());
0571 
0572     if (!detrayMaterial.has_value()) {
0573       continue;
0574     }
0575 
0576     printSurfaceInfo(*detrayMaterial, surface);
0577 
0578     assignMaterial(surface.surfaceMaterial(), *detrayMaterial, srfIdx, surface);
0579   }
0580   ACTS_VERBOSE("Looping over " << volume.portals().size()
0581                                << " portals in volume " << volPayload.name);
0582 
0583   // Portals need special treatment: we have decomposed them to their trivial
0584   // portal links, and only registered a subset to them as well. We again need
0585   // to decompose here, and only look for the portals that are actually found in
0586   // the volume payload
0587   for (const auto& portal : volume.portals()) {
0588     // First check, if the combined portal surface has material assigned at all,
0589     // if not there's nothing to do
0590     const auto* surfaceMaterial = portal.surface().surfaceMaterial();
0591     if (surfaceMaterial == nullptr) {
0592       continue;
0593     }
0594 
0595     std::optional detrayMaterial =
0596         m_cfg.convertSurfaceMaterial(*surfaceMaterial);
0597 
0598     // Portal surface material reports it does not apply to detray, skip
0599     if (!detrayMaterial.has_value()) {
0600       continue;
0601     }
0602 
0603     printSurfaceInfo(*detrayMaterial, portal.surface());
0604 
0605     // Have valid detray material, now we need to find the surfaces that are
0606     // actually there in detray
0607 
0608     for (auto dir : {Direction::AlongNormal(), Direction::OppositeNormal()}) {
0609       const auto* link = portal.getLink(dir);
0610 
0611       if (link == nullptr) {
0612         continue;
0613       }
0614       ACTS_VERBOSE("Processing dir=" << dir << " for portal on surface "
0615                                      << portal.surface().geometryId());
0616 
0617       std::vector<const TrivialPortalLink*> trivials =
0618           decomposeToTrivials(*link, logger());
0619 
0620       ACTS_VERBOSE("Portal link in volume "
0621                    << volPayload.name << " has produced " << trivials.size()
0622                    << " trivials");
0623       for (const auto* trivial : trivials) {
0624         auto srfIt = surfaceIndices.find(&trivial->surface());
0625 
0626         if (srfIt == surfaceIndices.end()) {
0627           // This trivial was not converted, skip
0628           continue;
0629         }
0630 
0631         std::size_t srfIdx = srfIt->second;
0632 
0633         ACTS_VERBOSE("Trivial portal link in volume "
0634                      << volPayload.name << " has surface "
0635                      << trivial->surface().geometryId() << " detray idx "
0636                      << srfIdx);
0637 
0638         // Assign (a copy of) the detray material to the surface payload
0639         // associated with this trivial
0640         assignMaterial(surfaceMaterial, *detrayMaterial, srfIdx,
0641                        trivial->surface());
0642       }
0643     }
0644   }
0645 
0646   return {grids, homogeneous};
0647 }
0648 
0649 DetrayPayloadConverter::Payloads
0650 DetrayPayloadConverter::convertTrackingGeometry(
0651     const GeometryContext& gctx, const TrackingGeometry& geometry) const {
0652   ACTS_INFO("Converting tracking geometry to detray format");
0653 
0654   if (geometry.geometryVersion() != TrackingGeometry::GeometryVersion::Gen3) {
0655     ACTS_WARNING(
0656         "Only Gen3 tracking geometries are supported. Gen1 geometries will "
0657         "give wrong results");
0658   }
0659 
0660   if (m_cfg.beampipeVolume == nullptr) {
0661     throw std::runtime_error(
0662         "Beampipe volume not set. This is needed to ensure detray receives the "
0663         "beampip volume where it expects it");
0664   }
0665 
0666   Payloads payloads;
0667   payloads.detector = std::make_unique<detray::io::detector_payload>();
0668   detray::io::detector_payload& detPayload = *payloads.detector;
0669 
0670   payloads.homogeneousMaterial =
0671       std::make_unique<detray::io::detector_homogeneous_material_payload>();
0672   detray::io::detector_homogeneous_material_payload& dthmPayload =
0673       *payloads.homogeneousMaterial;
0674 
0675   payloads.materialGrids = std::make_unique<detray::io::detector_grids_payload<
0676       detray::io::material_slab_payload, detray::io::material_id>>();
0677 
0678   detray::io::detector_grids_payload<detray::io::material_slab_payload,
0679                                      detray::io::material_id>& materialGrids =
0680       *payloads.materialGrids;
0681 
0682   payloads.surfaceGrids = std::make_unique<
0683       detray::io::detector_grids_payload<std::size_t, detray::io::accel_id>>();
0684 
0685   detray::io::detector_grids_payload<std::size_t, detray::io::accel_id>&
0686       surfaceGrids = *payloads.surfaceGrids;
0687 
0688   std::unordered_map<const TrackingVolume*, std::size_t> volumeIds;
0689 
0690   auto lookup = [&volumeIds](const TrackingVolume* v) {
0691     return volumeIds.at(v);
0692   };
0693 
0694   std::unordered_map<const TrackingVolume*,
0695                      std::unordered_map<const Surface*, std::size_t>>
0696       volumeSurfaceIndices;
0697 
0698   geometry.apply([&](const TrackingVolume& volume) {
0699     auto& volPayload = detPayload.volumes.emplace_back(convertVolume(volume));
0700     volPayload.index.link = detPayload.volumes.size() - 1;
0701     volumeIds[&volume] = volPayload.index.link;
0702 
0703     ACTS_DEBUG("Volume " << volume.volumeName() << " has index "
0704                          << volPayload.index.link);
0705 
0706     auto& surfaceIndices = volumeSurfaceIndices[&volume];
0707 
0708     for (auto& surface : volume.surfaces()) {
0709       auto& srfPayload =
0710           volPayload.surfaces.emplace_back(convertSurface(gctx, surface));
0711       srfPayload.index_in_coll = volPayload.surfaces.size() - 1;
0712       srfPayload.masks.at(0).volume_link.link = volPayload.index.link;
0713       surfaceIndices[&surface] = srfPayload.index_in_coll.value();
0714     }
0715   });
0716 
0717   // Run again over volumes, can lookup volume index from pointer now
0718   geometry.apply([&](const TrackingVolume& volume) {
0719     auto& volPayload = detPayload.volumes.at(volumeIds.at(&volume));
0720     auto& surfaceIndices = volumeSurfaceIndices[&volume];
0721 
0722     for (const auto& portal : volume.portals()) {
0723       handlePortal(gctx, volume, volPayload, lookup, surfaceIndices, portal);
0724     }
0725 
0726     ACTS_DEBUG("Volume " << volume.volumeName() << " (detray idx: "
0727                          << volPayload.index.link << ") has "
0728                          << volPayload.surfaces.size() << " total surfaces");
0729 
0730     std::size_t nPortals =
0731         std::ranges::count_if(volPayload.surfaces, [](const auto& srfPayload) {
0732           return srfPayload.type == detray::surface_id::e_portal;
0733         });
0734     ACTS_DEBUG("-> portals:        " << nPortals);
0735     std::size_t nSensitives =
0736         std::ranges::count_if(volPayload.surfaces, [](const auto& srfPayload) {
0737           return srfPayload.type == detray::surface_id::e_sensitive;
0738         });
0739     ACTS_DEBUG("-> sensitives:     " << nSensitives);
0740     ACTS_DEBUG("-> other surfaces: " << volPayload.surfaces.size() - nPortals -
0741                                             nSensitives);
0742 
0743     for (const auto& [surface, idx] : surfaceIndices) {
0744       ACTS_VERBOSE("Surface " << surface->geometryId() << " (&: " << surface
0745                               << ", type: " << surface->type()
0746                               << ") has detray index " << idx);
0747     }
0748 
0749     // Portals have produced surfaces and are added in volume payload, handle
0750     // material now
0751 
0752     auto [grids, homogeneous] =
0753         convertMaterial(volume, surfaceIndices, volPayload);
0754 
0755     ACTS_DEBUG("Volume " << volume.volumeName() << " (detray idx: "
0756                          << volPayload.index.link << ") has "
0757                          << homogeneous.mat_slabs.size() << " material slabs");
0758 
0759     if (!homogeneous.mat_slabs.empty()) {
0760       // Only add if it's not empty (it might be)
0761       // NOTE: Currently, it'll always be populated by at least the homogeneous
0762       // NOTE: Volume association is internal to
0763       // `detray::io::material_volume_payload`
0764       dthmPayload.volumes.emplace_back(std::move(homogeneous));
0765     }
0766 
0767     ACTS_DEBUG("Volume " << volume.volumeName()
0768                          << " (detray idx: " << volPayload.index.link
0769                          << ") has " << grids.size() << " material grids");
0770     if (!grids.empty()) {
0771       // Only add if we have grids
0772       // NOTE: Volume association is EXTERNAL, i.e. we need to fill a map keyed
0773       // by the volume index
0774       materialGrids.grids[volPayload.index.link] = std::move(grids);
0775     }
0776 
0777     // Look for navigation policies that we need to convert!
0778     const auto* navPolicy = volume.navigationPolicy();
0779     if (navPolicy != nullptr) {
0780       // Create surface lookup function for this volume
0781       auto surfaceLookupFn =
0782           [&surfaceIndices](const Surface* surface) -> std::size_t {
0783         auto it = surfaceIndices.find(surface);
0784         if (it == surfaceIndices.end()) {
0785           throw std::runtime_error("Surface not found in surface indices map");
0786         }
0787         return it->second;
0788       };
0789 
0790       std::optional<DetraySurfaceGrid> detrayGrid = std::nullopt;
0791 
0792       navPolicy->visit([&](const INavigationPolicy& policy) {
0793         if (detrayGrid.has_value()) {
0794           ACTS_ERROR("Volume "
0795                      << volume.volumeName()
0796                      << " has more than one detray-convertible navigation "
0797                         "policy. This cannot currently be handled.");
0798           throw std::runtime_error{
0799               "Multiple detray-compatible navigation policies"};
0800         }
0801 
0802         detrayGrid = m_cfg.convertNavigationPolicy(policy, gctx,
0803                                                    surfaceLookupFn, logger());
0804       });
0805 
0806       if (detrayGrid.has_value()) {
0807         ACTS_DEBUG("Volume " << volume.volumeName()
0808                              << " (detray idx: " << volPayload.index.link
0809                              << ") has navigation policy which produced "
0810                              << detrayGrid->bins.size() << " populated bins");
0811 
0812         detrayGrid->owner_link.link = volPayload.index.link;
0813 
0814         // Add the surface grid to the payload
0815         surfaceGrids.grids[volPayload.index.link].push_back(*detrayGrid);
0816       }
0817       // per volume, we have a VECTOR of grids: what are they? are they always
0818       // tied to a surface? which one?
0819     }
0820   });
0821 
0822   // HACK: Beampipe MUST have index 0
0823   std::size_t beampipeIdx = volumeIds.at(m_cfg.beampipeVolume);
0824   ACTS_DEBUG("Beampipe volume (" << m_cfg.beampipeVolume->volumeName()
0825                                  << ") index: " << beampipeIdx);
0826   ACTS_DEBUG("Volume at index 0 is " << detPayload.volumes.at(0).name);
0827 
0828   // Swap beampipe and world volumes self-index
0829   std::swap(detPayload.volumes.at(0).index.link,
0830             detPayload.volumes.at(beampipeIdx).index.link);
0831 
0832   // Swap beampipe and world volumes location in vector
0833   std::swap(detPayload.volumes.at(0), detPayload.volumes.at(beampipeIdx));
0834 
0835   // Adjust volume indices in surfaces after swapping
0836   for (auto& vol : detPayload.volumes) {
0837     for (auto& srf : vol.surfaces) {
0838       auto& mask = srf.masks.at(0);
0839       if (mask.volume_link.link == beampipeIdx) {
0840         mask.volume_link.link = 0;
0841       } else if (mask.volume_link.link == 0) {
0842         mask.volume_link.link = beampipeIdx;
0843       }
0844     }
0845   }
0846 
0847   {
0848     // Possibly swap homogeneous material entries in vector if they both exist
0849     auto find = [](std::size_t id) {
0850       return [id](const auto& vol) { return vol.volume_link.link == id; };
0851     };
0852 
0853     auto beampipeIt =
0854         std::ranges::find_if(dthmPayload.volumes, find(beampipeIdx));
0855     auto worldIt = std::ranges::find_if(dthmPayload.volumes, find(0));
0856 
0857     if (beampipeIt != dthmPayload.volumes.end() &&
0858         worldIt != dthmPayload.volumes.end()) {
0859       // BOTH world and beampipe have homogoenous material: swap them
0860       ACTS_DEBUG("Swapping beampipe and world homogoenous material entries");
0861       std::swap(*beampipeIt, *worldIt);
0862     }
0863 
0864     // Retarget the entries, regardless of whether there is an entry for only
0865     // one of them
0866     for (auto& mat : dthmPayload.volumes) {
0867       if (mat.volume_link.link == beampipeIdx) {
0868         ACTS_DEBUG("Reassigning beampipe homogoenous material to index 0");
0869         mat.volume_link.link = 0;
0870       } else if (mat.volume_link.link == 0) {
0871         ACTS_DEBUG("Reassigning world homogoenous material to beampipe index "
0872                    << beampipeIdx);
0873         mat.volume_link.link = beampipeIdx;
0874       }
0875     }
0876   }
0877 
0878   {
0879     // Adjust material grids after swapping
0880     auto beampipeGridIt = materialGrids.grids.find(beampipeIdx);
0881     auto worldGridIt = materialGrids.grids.find(0);
0882 
0883     if (beampipeGridIt != materialGrids.grids.end() &&
0884         worldGridIt != materialGrids.grids.end()) {
0885       // BOTH world and beampipe have grid specifiers: swap them
0886       ACTS_DEBUG("Swapping beampipe and world material grid specifiers");
0887       std::swap(beampipeGridIt->second, worldGridIt->second);
0888     } else if (beampipeGridIt != materialGrids.grids.end()) {
0889       // ONLY beampipe has grid specifier: move it to world
0890       ACTS_DEBUG("Moving beampipe material grid specifier to world");
0891       materialGrids.grids[0] = std::move(beampipeGridIt->second);
0892       materialGrids.grids.erase(beampipeGridIt);
0893     } else if (worldGridIt != materialGrids.grids.end()) {
0894       // ONLY world has grid specifier: move it to beampipe
0895       ACTS_DEBUG("Moving world material grid specifier to beampipe");
0896       materialGrids.grids[beampipeIdx] = std::move(worldGridIt->second);
0897       materialGrids.grids.erase(worldGridIt);
0898     }
0899   }
0900 
0901   {
0902     // Adjust surface grids after swapping
0903     // @NOTE: The beampipe should **generally** not have a surface grid, but
0904     //        let's be safe and swap them regardless
0905 
0906     auto beampipeGridIt = surfaceGrids.grids.find(beampipeIdx);
0907     auto worldGridIt = surfaceGrids.grids.find(0);
0908 
0909     if (beampipeGridIt != surfaceGrids.grids.end() &&
0910         worldGridIt != surfaceGrids.grids.end()) {
0911       // BOTH world and beampipe have grid specifiers: swap them
0912       ACTS_DEBUG("Swapping beampipe and world surface grid specifiers");
0913       std::swap(beampipeGridIt->second, worldGridIt->second);
0914     } else if (beampipeGridIt != surfaceGrids.grids.end()) {
0915       // ONLY beampipe has grid specifier: move it to world
0916       ACTS_DEBUG("Moving beampipe surface grid specifier to world");
0917       surfaceGrids.grids[0] = std::move(beampipeGridIt->second);
0918       surfaceGrids.grids.erase(beampipeGridIt);
0919     } else if (worldGridIt != surfaceGrids.grids.end()) {
0920       // ONLY world has grid specifier: move it to beampipe
0921       ACTS_DEBUG("Moving world surface grid specifier to beampipe");
0922       surfaceGrids.grids[beampipeIdx] = std::move(worldGridIt->second);
0923       surfaceGrids.grids.erase(worldGridIt);
0924     }
0925   }
0926 
0927   // This needs to happen after swapping so that the indices are correct
0928   payloads.names = {{0, "Detector"}};
0929   for (const auto& volume : detPayload.volumes) {
0930     payloads.names.emplace(volume.index.link + 1, volume.name);
0931   }
0932 
0933   ACTS_DEBUG("Collected " << detPayload.volumes.size() << " volumes");
0934 
0935   return payloads;
0936 }
0937 
0938 }  // namespace ActsPlugins
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!