src/Geometry/RzGridVolumeLookup.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/Geometry/RzGridVolumeLookup.hpp"
0010
0011 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
0012 #include "Acts/Geometry/TrackingVolume.hpp"
0013 #include "Acts/Utilities/Axis.hpp"
0014 #include "Acts/Utilities/MultiAxis.hpp"
0015 #include "Acts/Utilities/VectorHelpers.hpp"
0016
0017 #include <algorithm>
0018 #include <array>
0019 #include <cmath>
0020 #include <limits>
0021 #include <numbers>
0022 #include <numeric>
0023 #include <optional>
0024 #include <span>
0025 #include <stdexcept>
0026 #include <vector>
0027
0028 #include <Eigen/Geometry>
0029
0030 #include "VolumeLookupImpl.hpp"
0031
0032 namespace Acts {
0033 namespace {
0034 constexpr std::size_t noRzNode = std::numeric_limits<std::size_t>::max();
0035
0036 // Unlike Acts::AxisAlignedBoundingBox::intersect, Eigen's AlignedBox includes
0037 // the upper edges, as required for conservative RZ candidate selection.
0038 // Conservative envelope for one subtree candidate.
0039 struct RzVolumeCandidate {
0040   const TrackingVolume* m_volume;
0041   Eigen::AlignedBox2d m_box;
0042   bool m_indexable = false;
0043   // Stop indexing here and dispatch to this volume's explicitly set lookup.
0044   bool m_delegate = false;
0045
0046   RzVolumeCandidate(const GeometryContext& gctx, const TrackingVolume& v,
0047                     double maxTolerance)
0048       : m_volume(&v) {
0049     const auto* cylinder =
0050         dynamic_cast<const CylinderVolumeBounds*>(&v.volumeBounds());
0051     const auto& trf = v.localToGlobalTransform(gctx);
0052     // Only fixed full cylinders with near-identity rotations are indexed.
0053     if (!v.isAlignable() && cylinder != nullptr &&
0054         trf.linear().isIdentity(1e-12) &&
0055         cylinder->get(CylinderVolumeBounds::eHalfPhiSector) ==
0056             std::numbers::pi) {
0057       const double offset = VectorHelpers::perp(trf.translation());
0058       const double hz = cylinder->get(CylinderVolumeBounds::eHalfLengthZ);
0059       // Bound displacement from the accepted rotation, including tolerance.
0060       // This keeps the envelope conservative even for tiny nonzero tilts.
0061       const double extent =
0062           std::hypot(cylinder->get(CylinderVolumeBounds::eMaxR) + maxTolerance,
0063                      hz + maxTolerance);
0064       const double margin =
0065           maxTolerance + 1e-9 +
0066           (trf.linear() - RotationMatrix3::Identity()).norm() * extent;
0067       m_box.min() =
0068           Vector2(cylinder->get(CylinderVolumeBounds::eMinR) - offset - margin,
0069                   trf.translation().z() - hz - margin);
0070       m_box.max() =
0071           Vector2(cylinder->get(CylinderVolumeBounds::eMaxR) + offset + margin,
0072                   trf.translation().z() + hz + margin);
0073       m_indexable = true;
0074     }
0075   }
0076 };
0077
0078 // Grid edges are the sorted unique rectangle edges. Each rectangle occupies
0079 // every intersected cell, including its upper-edge bin. Candidate lists retain
0080 // hierarchy order and use CSR storage: offsets delimit slices in one array.
0081 struct RzGrid {
0082   using Config = RzGridVolumeLookup::Config;
0083   using RzAxis = Axis<AxisType::Variable, AxisBoundaryType::Bound>;
0084   using RzAxes = MultiAxis<RzAxis, RzAxis>;
0085   std::optional<RzAxes> m_axes;
0086   std::vector<std::size_t> m_offsets;
0087   std::vector<std::size_t> m_candidates;
0088
0089   bool build(const std::vector<RzVolumeCandidate>& children,
0090              const Config& config) {
0091     std::array<std::vector<double>, 2> edges;
0092     for (int axis = 0; axis < 2; ++axis) {
0093       auto& e = edges[axis];
0094       for (const auto& child : children) {
0095         if (!child.m_indexable) {
0096           continue;
0097         }
0098         e.push_back(child.m_box.min()[axis]);
0099         e.push_back(child.m_box.max()[axis]);
0100       }
0101       std::ranges::sort(e);
0102       const auto duplicates = std::ranges::unique(e);
0103       e.erase(duplicates.begin(), duplicates.end());
0104     }
0105     // N edges give N-1 regular bins plus two reserved under-/overflow bins.
0106     const auto nr = edges[0].size() + 1;
0107     const auto nz = edges[1].size() + 1;
0108     // Check nr*nz against the limit without overflowing the product.
0109     if (nr > config.maxGridCells / nz) {
0110       return false;
0111     }
0112     m_axes.emplace(RzAxis(std::move(edges[0])), RzAxis(std::move(edges[1])));
0113     // One extra offset marks the end of the last cell's candidate slice.
0114     m_offsets.resize(m_axes->getNTotalBins(true) + 1);
0115     auto cells = [&](const Eigen::AlignedBox2d& b, auto&& visit) {
0116       const auto lower = m_axes->getLocalBinsFromPoint(b.min());
0117       const auto upper = m_axes->getLocalBinsFromPoint(b.max());
0118       // Include the bin at the upper edge. Bound axes clamp the final edge to
0119       // the last regular bin, matching the inclusive rectangle containment.
0120       for (auto r = lower[0]; r <= upper[0]; ++r) {
0121         for (auto z = lower[1]; z <= upper[1]; ++z) {
0122           visit(m_axes->getGlobalBinFromLocalBins({r, z}));
0123         }
0124       }
0125     };
0126     // Shift counts by one, prefix-sum to slice offsets, then fill with separate
0127     // cursors so the offsets remain unchanged.
0128     for (const auto& child : children) {
0129       if (!child.m_indexable) {
0130         continue;
0131       }
0132       cells(child.m_box, [&](auto cell) { ++m_offsets[cell + 1]; });
0133     }
0134     std::partial_sum(m_offsets.begin(), m_offsets.end(), m_offsets.begin());
0135     m_candidates.resize(m_offsets.back());
0136     auto cursor = m_offsets;
0137     for (std::size_t i = 0; i < children.size(); ++i) {
0138       if (!children[i].m_indexable) {
0139         continue;
0140       }
0141       cells(children[i].m_box,
0142             [&](auto cell) { m_candidates[cursor[cell]++] = i; });
0143     }
0144     return true;
0145   }
0146
0147   std::span<const std::size_t> candidates(const Vector2& key) const {
0148     // The caller checks the inclusive domain; MultiAxis::isInside() would
0149     // exclude the upper edge, while Bound axes alone also clamp outside points.
0150     const auto cell = m_axes->getGlobalBinFromPoint(key);
0151     return std::span{m_candidates}.subspan(
0152         m_offsets[cell], m_offsets[cell + 1] - m_offsets[cell]);
0153   }
0154 };
0155 }  // namespace
0156 // One spatial index covers the subtree. Candidates are ordered in postorder:
0157 // descendants before their parent fallback, and an entire earlier sibling's
0158 // subtree before later siblings. Geometry construction must ensure each child
0159 // is contained within its parent, so candidates need no ancestor checks.
0160 struct RzGridVolumeLookup::Impl {
0161   const TrackingVolume& m_root;
0162   Config m_config;
0163   bool m_gridLimit = false;
0164   bool m_hasIndex = false;
0165   std::size_t m_alignable = 0;
0166   std::size_t m_unsupported = 0;
0167   Eigen::AlignedBox2d m_domain;
0168   std::vector<RzVolumeCandidate> m_children;
0169   std::vector<std::size_t> m_fallbackCandidates;
0170   // One grid for this lookup's entire collected subtree, not one per volume.
0171   RzGrid m_grid;
0172
0173   Impl(const GeometryContext& gctx, const TrackingVolume& root,
0174        const Config& cfg)
0175       : m_root(root), m_config(cfg) {
0176     // A bare TrackingVolume has no generation flag. Reject Gen1-only storage
0177     // wherever it is encountered in this index's scope, including the root.
0178     auto requireGen3 = [](const TrackingVolume& volume) {
0179       if (volume.confinedVolumes() || volume.confinedLayers() != nullptr ||
0180           !volume.denseVolumes().empty()) {
0181         throw std::invalid_argument(
0182             "RZ volume lookup does not support Gen1 storage in volume '" +
0183             volume.volumeName() + "'");
0184       }
0185     };
0186     requireGen3(root);
0187     // Collect geometry independently of the spatial grid. Postorder IDs
0188     // encode recursive-scan priority; build() below indexes their RZ envelopes.
0189     auto collect = [&](auto&& self, const TrackingVolume& volume) -> void {
0190       requireGen3(volume);
0191       RzVolumeCandidate entry(gctx, volume, cfg.maxTolerance);
0192       entry.m_delegate = volume.volumeLookup() != nullptr;
0193       if (!entry.m_delegate) {
0194         for (const auto& child : volume.volumes()) {
0195           self(self, child);
0196         }
0197       }
0198       if (entry.m_indexable) {
0199         m_domain.extend(entry.m_box);
0200         m_hasIndex = true;
0201       } else {
0202         m_fallbackCandidates.push_back(m_children.size());
0203         if (volume.isAlignable()) {
0204           ++m_alignable;
0205         } else {
0206           ++m_unsupported;
0207         }
0208       }
0209       m_children.push_back(std::move(entry));
0210     };
0211     for (const auto& child : root.volumes()) {
0212       collect(collect, child);
0213     }
0214     // Build the single grid after collecting the whole indexed scope.
0215     if (m_hasIndex) {
0216       m_gridLimit = !m_grid.build(m_children, cfg);
0217     }
0218   }
0219
0220   void warnAboutFallback(const TrackingVolume& volume, double fraction,
0221                          const Logger& logger) const {
0222     // A threshold of 1 suppresses this warning.
0223     if (fraction == 1 || !logger.doPrint(Logging::WARNING)) {
0224       return;
0225     }
0226     const auto total = m_children.size();
0227     const auto affected = m_gridLimit ? total : m_fallbackCandidates.size();
0228     // Avoid division by zero; warn only above the configured fraction.
0229     if (total == 0 ||
0230         static_cast<double>(affected) / static_cast<double>(total) <=
0231             fraction) {
0232       return;
0233     }
0234     ACTS_WARNING("RZ volume lookup for '"
0235                  << volume.volumeName() << "': " << affected << "/" << total
0236                  << " volumes in the indexed scope require linear fallback ("
0237                  << 100. * static_cast<double>(affected) /
0238                         static_cast<double>(total)
0239                  << "%, warning threshold " << 100. * fraction
0240                  << "%). Reasons: alignable=" << m_alignable
0241                  << ", unsupported shape/orientation=" << m_unsupported
0242                  << ", grid cell limit=" << m_gridLimit
0243                  << ". Explicit descendant lookup boundaries remain active.");
0244   }
0245
0246   const TrackingVolume* find(const GeometryContext& gctx, const Vector3& p,
0247                              const Vector2& key, double tolerance) const {
0248     if (m_gridLimit) {
0249       return detail::linearVolumeScan(m_root, gctx, p, tolerance);
0250     }
0251     if (!m_root.inside(gctx, p, tolerance)) {
0252       return nullptr;
0253     }
0254     // IDs encode hierarchy priority. Each ordered candidate list can stop
0255     // at its first exact match; choose the earlier match across both lists.
0256     std::size_t best = noRzNode;
0257     // Return true to keep scanning this ordered list, false to stop it.
0258     auto consider = [&](std::size_t id) {
0259       // Later IDs cannot beat a match already found in either list.
0260       if (id >= best) {
0261         return false;
0262       }
0263       const auto& entry = m_children[id];
0264       // Check the conservative box when available, then exact containment.
0265       // A miss continues the search; a match stops this ordered list.
0266       if ((entry.m_indexable && !entry.m_box.contains(key)) ||
0267           !entry.m_volume->inside(gctx, p, tolerance)) {
0268         return true;
0269       }
0270       best = id;
0271       // The other list may still contain a higher-priority match.
0272       return false;
0273     };
0274     for (auto id : m_fallbackCandidates) {
0275       if (!consider(id)) {
0276         break;
0277       }
0278     }
0279     if (m_hasIndex && m_domain.contains(key)) {
0280       for (auto id : m_grid.candidates(key)) {
0281         if (!consider(id)) {
0282           break;
0283         }
0284       }
0285     }
0286     if (best == noRzNode) {
0287       return &m_root;
0288     }
0289     const auto& entry = m_children[best];
0290     return entry.m_delegate
0291                ? entry.m_volume->lowestTrackingVolume(gctx, p, tolerance)
0292                : entry.m_volume;
0293   }
0294 };
0295
0296 RzGridVolumeLookup::RzGridVolumeLookup(const GeometryContext& gctx,
0297                                        const TrackingVolume& world,
0298                                        const Config& config,
0299                                        const Logger& logger) {
0300   if (!std::isfinite(config.fallbackWarningFraction) ||
0301       config.fallbackWarningFraction < 0 ||
0302       config.fallbackWarningFraction > 1 ||
0303       !std::isfinite(config.maxTolerance) || config.maxTolerance < 0 ||
0304       config.maxGridCells == 0) {
0305     throw std::invalid_argument("Invalid RZ volume lookup configuration");
0306   }
0307   m_impl = std::make_unique<Impl>(gctx, world, config);
0308   m_impl->warnAboutFallback(world, config.fallbackWarningFraction, logger);
0309 }
0310
0311 RzGridVolumeLookup::~RzGridVolumeLookup() = default;
0312
0313 const RzGridVolumeLookup::Config& RzGridVolumeLookup::config() const {
0314   return m_impl->m_config;
0315 }
0316
0317 const TrackingVolume* RzGridVolumeLookup::findLowest(
0318     const GeometryContext& gctx, const Vector3& position,
0319     double tolerance) const {
0320   if (!std::isfinite(tolerance) || tolerance < 0 ||
0321       tolerance > m_impl->m_config.maxTolerance) {
0322     return detail::linearVolumeScan(m_impl->m_root, gctx, position, tolerance);
0323   }
0324   return m_impl->find(gctx, position,
0325                       Vector2(VectorHelpers::perp(position), position.z()),
0326                       tolerance);
0327 }
0328 }  // namespace Acts