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