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