Detray/src/DetrayGeometryConverter.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 "ActsPlugins/Detray/DetrayGeometryConverter.hpp"
0010
0011 #include "Acts/Detector/Detector.hpp"
0012 #include "Acts/Detector/DetectorVolume.hpp"
0013 #include "Acts/Detector/Portal.hpp"
0014 #include "Acts/Navigation/PortalNavigation.hpp"
0015 #include "Acts/Surfaces/CylinderBounds.hpp"
0016 #include "Acts/Surfaces/CylinderSurface.hpp"
0017 #include "Acts/Surfaces/DiscSurface.hpp"
0018 #include "Acts/Surfaces/RadialBounds.hpp"
0019 #include "Acts/Surfaces/RegularSurface.hpp"
0020 #include "Acts/Surfaces/Surface.hpp"
0021 #include "Acts/Surfaces/SurfaceBounds.hpp"
0022 #include "ActsPlugins/Detray/DetrayConversionUtils.hpp"
0023 #include "ActsPlugins/Json/DetrayJsonHelper.hpp"
0024
0025 #include <algorithm>
0026
0027 #include <detray/io/frontend/detector_writer.hpp>
0028
0029 using namespace Acts;
0030 using namespace detray;
0031
0032 detray::io::mask_payload ActsPlugins::DetrayGeometryConverter::convertMask(
0033     const SurfaceBounds& bounds, bool portal) {
0034   detray::io::mask_payload maskPayload;
0035   auto [shape, boundaries] = DetrayJsonHelper::maskFromBounds(bounds, portal);
0036   maskPayload.shape = static_cast<io::mask_payload::mask_shape>(shape);
0037   maskPayload.boundaries = static_cast<std::vector<real_io>>(boundaries);
0038   // default maskPayload.volume_link
0039
0040   return maskPayload;
0041 }
0042
0043 detray::io::surface_payload
0044 ActsPlugins::DetrayGeometryConverter::convertSurface(
0045     const GeometryContext& gctx, const Surface& surface, bool portal) {
0046   detray::io::surface_payload surfacePayload;
0047
0048   surfacePayload.transform =
0049       DetrayConversionUtils::convertTransform(surface.transform(gctx));
0050   surfacePayload.source = surface.geometryId().value();
0051   surfacePayload.barcode = std::nullopt;
0052   surfacePayload.type = static_cast<detray::surface_id>(
0053       portal ? surface_id::e_portal
0054              : (surface.geometryId().sensitive() > 0
0055                     ? detray::surface_id::e_sensitive
0056                     : detray::surface_id::e_passive));
0057   surfacePayload.masks = {convertMask(surface.bounds())};
0058   return surfacePayload;
0059 }
0060
0061 std::vector<detray::io::surface_payload>
0062 ActsPlugins::DetrayGeometryConverter::convertPortal(
0063     DetrayConversionUtils::Cache& cCache, const GeometryContext& gctx,
0064     const Experimental::Portal& portal, std::size_t ip,
0065     const Experimental::DetectorVolume& volume,
0066     const std::vector<OrientedSurface>& orientedSurfaces) {
0067   std::vector<detray::io::surface_payload> portals{};
0068
0069   const RegularSurface& surface = portal.surface();
0070   const auto& volumeLinks = portal.portalNavigation();
0071
0072   // First assumption for outside link (along direction)
0073   std::size_t outside = 1u;
0074
0075   // Find out if you need to take the outside or inside volume
0076   // for planar surfaces that's easy
0077   if (surface.type() != Surface::SurfaceType::Cylinder) {
0078     // Get the two volume center
0079     const auto volumeCenter = volume.transform(gctx).translation();
0080     const auto surfaceCenter = surface.center(gctx);
0081     const auto surfaceNormal = surface.normal(gctx, surfaceCenter);
0082     // Get the direction from the volume to the surface, correct link
0083     const auto volumeToSurface = surfaceCenter - volumeCenter;
0084     if (volumeToSurface.dot(surfaceNormal) < 0.) {
0085       outside = 0u;
0086     }
0087   } else {
0088     // This is a cylinder portal, inner cover reverses the normal
0089     if (ip == 3u) {
0090       outside = 0u;
0091     }
0092   }
0093
0094   const auto& outsideLink = volumeLinks[outside];
0095   // Grab the corresponding volume link
0096   // If it is a single link, we are done
0097   const auto* instance = outsideLink.instance();
0098   // Single link cast
0099   auto singleLink =
0100       dynamic_cast<const Experimental::SingleDetectorVolumeNavigation*>(
0101           instance);
0102
0103   auto [surfaceAdjusted, insidePointer] = orientedSurfaces[ip];
0104   // Assign the geometry id to the surface
0105   surfaceAdjusted->assignGeometryId(surface.geometryId());
0106
0107   // Single link detected - just write it out, we use the oriented surface
0108   // in order to make sure the size is adjusted
0109   if (singleLink != nullptr) {
0110     // Single link can be written out
0111     std::size_t vLink = cCache.volumeIndex(singleLink->object());
0112     auto portalPayload = convertSurface(gctx, *surfaceAdjusted, true);
0113     portalPayload.masks.at(0).volume_link.link = vLink;
0114     portals.push_back(portalPayload);
0115   } else {
0116     // Multi link detected - 1D
0117     auto multiLink1D =
0118         dynamic_cast<const Experimental::BoundVolumesGrid1Navigation*>(
0119             instance);
0120     if (multiLink1D != nullptr) {
0121       // Resolve the multi link 1D
0122       auto boundaries =
0123           multiLink1D->indexedUpdater.grid.axes()[0u]->getBinEdges();
0124       const auto& cast = multiLink1D->indexedUpdater.casts[0u];
0125       const auto& transform = multiLink1D->indexedUpdater.transform;
0126       const auto& volumes = multiLink1D->indexedUpdater.extractor.dVolumes;
0127
0128       // Apply the correction from local to global boundaries
0129       double gCorr = VectorHelpers::cast(transform.translation(), cast);
0130       std::ranges::for_each(boundaries, [&gCorr](double& b) { b -= gCorr; });
0131
0132       // Get the volume indices
0133       auto surfaceType = surfaceAdjusted->type();
0134       std::vector<unsigned int> vIndices = {};
0135       for (const auto& v : volumes) {
0136         vIndices.push_back(cCache.volumeIndex(v));
0137       }
0138
0139       // Pick the surface dimension
0140       std::array<double, 2u> clipRange = {0., 0.};
0141       std::vector<double> boundValues = surfaceAdjusted->bounds().values();
0142       if (surfaceType == Surface::SurfaceType::Cylinder &&
0143           cast == AxisDirection::AxisZ) {
0144         double zPosition = surfaceAdjusted->center(gctx).z();
0145         clipRange = {
0146             zPosition - boundValues[CylinderBounds::BoundValues::eHalfLengthZ],
0147             zPosition + boundValues[CylinderBounds::BoundValues::eHalfLengthZ]};
0148       } else if (surfaceType == Surface::SurfaceType::Disc &&
0149                  cast == AxisDirection::AxisR) {
0150         clipRange = {boundValues[RadialBounds::BoundValues::eMinR],
0151                      boundValues[RadialBounds::BoundValues::eMaxR]};
0152       } else {
0153         throw std::runtime_error(
0154             "PortalDetrayGeometryConverter: surface type not (yet) supported "
0155             "for "
0156             "detray "
0157             "conversion, only cylinder and disc are currently supported.");
0158       }
0159
0160       // Need to clip the parameter space to the surface dimension
0161       std::vector<unsigned int> clippedIndices = {};
0162       std::vector<double> clippedBoundaries = {};
0163       clippedBoundaries.push_back(clipRange[0u]);
0164       for (const auto [ib, b] : enumerate(boundaries)) {
0165         if (ib > 0) {
0166           unsigned int vI = vIndices[ib - 1u];
0167           double highEdge = boundaries[ib];
0168           if (boundaries[ib - 1] >= clipRange[1u]) {
0169             break;
0170           }
0171           if (highEdge <= clipRange[0u] ||
0172               std::abs(highEdge - clipRange[0u]) < 1e-5) {
0173             continue;
0174           }
0175           if (highEdge > clipRange[1u]) {
0176             highEdge = clipRange[1u];
0177           }
0178           clippedIndices.push_back(vI);
0179           clippedBoundaries.push_back(highEdge);
0180         }
0181       }
0182       // Interpret the clipped information
0183       //
0184       // Clipped cylinder case
0185       if (surfaceType == Surface::SurfaceType::Cylinder) {
0186         for (auto [ib, b] : enumerate(clippedBoundaries)) {
0187           if (ib > 0) {
0188             // Create sub surfaces
0189             std::array<double, CylinderBounds::BoundValues::eSize>
0190                 subBoundValues = {};
0191             for (auto [ibv, bv] : enumerate(boundValues)) {
0192               subBoundValues[ibv] = bv;
0193             }
0194             subBoundValues[CylinderBounds::BoundValues::eHalfLengthZ] =
0195                 0.5 * (b - clippedBoundaries[ib - 1u]);
0196             auto subBounds = std::make_shared<CylinderBounds>(subBoundValues);
0197             auto subTransform = Transform3::Identity();
0198             subTransform.pretranslate(Vector3(
0199                 0., 0.,
0200                 clippedBoundaries[ib - 1u] +
0201                     subBoundValues[CylinderBounds::BoundValues::eHalfLengthZ]));
0202
0203             auto subSurface =
0204                 Surface::makeShared<CylinderSurface>(subTransform, subBounds);
0205             subSurface->assignGeometryId(surface.geometryId());
0206
0207             auto portalPayload = convertSurface(gctx, *subSurface, true);
0208             portalPayload.masks.at(0).volume_link.link =
0209                 clippedIndices[ib - 1u];
0210             portals.push_back(portalPayload);
0211           }
0212         }
0213       } else {
0214         for (auto [ib, b] : enumerate(clippedBoundaries)) {
0215           if (ib > 0) {
0216             // Create sub surfaces
0217             std::array<double, RadialBounds::BoundValues::eSize>
0218                 subBoundValues = {};
0219             for (auto [ibv, bv] : enumerate(boundValues)) {
0220               subBoundValues[ibv] = bv;
0221             }
0222             subBoundValues[RadialBounds::BoundValues::eMinR] =
0223                 clippedBoundaries[ib - 1u];
0224             subBoundValues[RadialBounds::BoundValues::eMaxR] = b;
0225             auto subBounds = std::make_shared<RadialBounds>(subBoundValues);
0226             auto subSurface = Surface::makeShared<DiscSurface>(
0227                 portal.surface().transform(gctx), subBounds);
0228
0229             subSurface->assignGeometryId(surface.geometryId());
0230             auto portalPayload = convertSurface(gctx, *subSurface, true);
0231             portalPayload.masks.at(0).volume_link.link =
0232                 clippedIndices[ib - 1u];
0233             portals.push_back(portalPayload);
0234           }
0235         }
0236       }
0237
0238     } else {
0239       // Write surface with invalid link
0240       auto portalPayload = convertSurface(gctx, *surfaceAdjusted, true);
0241       using NavigationLink =
0242           typename DetrayHostDetector::surface_type::navigation_link;
0243       portalPayload.masks.at(0).volume_link.link =
0244           std::numeric_limits<NavigationLink>::max();
0245
0246       portals.push_back(portalPayload);
0247     }
0248   }
0249   return portals;
0250 }
0251
0252 detray::io::volume_payload ActsPlugins::DetrayGeometryConverter::convertVolume(
0253     DetrayConversionUtils::Cache& cCache, const GeometryContext& gctx,
0254     const Experimental::DetectorVolume& volume, const Logger& logger) {
0255   ACTS_DEBUG("DetrayGeometryConverter: converting volume "
0256              << volume.name() << " with " << volume.surfaces().size()
0257              << " surfaces and " << volume.portals().size() << " portals");
0258
0259   detray::io::volume_payload volumePayload;
0260   std::size_t volumeIndex = cCache.volumeIndex(&volume);
0261   volumePayload.name = volume.name();
0262   volumePayload.index.link = volumeIndex;
0263   volumePayload.transform =
0264       DetrayConversionUtils::convertTransform(volume.transform(gctx));
0265
0266   // Remember the link
0267   cCache.volumeLinks[volume.geometryId()] = volumePayload.index.link;
0268
0269   std::multimap<GeometryIdentifier, unsigned long> localSurfaceLinks;
0270
0271   // iterate over surfaces and portals keeping the same surf_pd.index_in_coll
0272   std::size_t sIndex = 0;
0273   for (const auto surface : volume.surfaces()) {
0274     io::surface_payload surfacePayload = convertSurface(gctx, *surface, false);
0275     // Set the index in the collection & remember it in the cache
0276     surfacePayload.index_in_coll = sIndex++;
0277     localSurfaceLinks.insert(
0278         {surface->geometryId(), surfacePayload.index_in_coll.value()});
0279     // Set mask to volume link
0280     surfacePayload.masks.at(0).volume_link.link =
0281         volumePayload.index.link;  // link surface' mask to volume
0282     volumePayload.surfaces.push_back(surfacePayload);
0283   }
0284
0285   auto orientedSurfaces =
0286       volume.volumeBounds().orientedSurfaces(volume.transform(gctx));
0287
0288   // Iterate over portals
0289   int portalCounter = 0;
0290   for (const auto& [ip, p] : enumerate(volume.portals())) {
0291     auto portals =
0292         convertPortal(cCache, gctx, *p, ip, volume, orientedSurfaces);
0293     ACTS_VERBOSE(" > portal " << ip << " split into " << portals.size()
0294                               << " surfaces");
0295     GeometryIdentifier geoID = p->surface().geometryId();
0296     std::ranges::for_each(portals, [&](auto& portalPayload) {
0297       // Set the index in the collection & remember it in the cache
0298       portalPayload.index_in_coll = sIndex++;
0299       localSurfaceLinks.insert({geoID, portalPayload.index_in_coll.value()});
0300       // Add it to the surfaces
0301       volumePayload.surfaces.push_back(portalPayload);
0302       portalCounter++;
0303     });
0304   }
0305   cCache.localSurfaceLinks[volumeIndex] = localSurfaceLinks;
0306   ACTS_DEBUG(" > " << volume.portals().size()
0307                    << " initial ACTS portals split into final " << portalCounter
0308                    << " detray portals");
0309   ACTS_VERBOSE(" > Local surface link cache has " << localSurfaceLinks.size()
0310                                                   << " entries");
0311
0312   return volumePayload;
0313 }
0314
0315 detray::io::detector_payload
0316 ActsPlugins::DetrayGeometryConverter::convertDetector(
0317     DetrayConversionUtils::Cache& cCache, const GeometryContext& gctx,
0318     const Experimental::Detector& detector, const Logger& logger) {
0319   ACTS_DEBUG("DetrayGeometryConverter: converting detector"
0320              << detector.name() << " with " << detector.volumes().size()
0321              << " volumes.");
0322
0323   // The detector payload which will be handed back
0324   detray::io::detector_payload detectorPayload;
0325
0326   for (const auto volume : detector.volumes()) {
0327     detectorPayload.volumes.push_back(
0328         convertVolume(cCache, gctx, *volume, logger));
0329   }
0330
0331   return detectorPayload;
0332 }