src/Seeding/GbtsNodeStorage.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/Seeding/GbtsNodeStorage.hpp"
0010
0011 #include "Acts/Seeding/GbtsGeometry.hpp"
0012 #include "Acts/SpacePointFormation/detail/StripSpacePointCalibrationImpl.hpp"
0013 #include "Acts/Utilities/MathHelpers.hpp"
0014
0015 #include <algorithm>
0016 #include <cmath>
0017 #include <limits>
0018 #include <numbers>
0019 #include <utility>
0020
0021 namespace Acts::Experimental {
0022
0023 GbtsNodeStorage::GbtsNodeStorage(const Config& config,
0024                                  std::shared_ptr<const GbtsGeometry> geometry,
0025                                  detail::GbtsTauLookupTable tauLut)
0026     : m_cfg(config),
0027       m_geometry(std::move(geometry)),
0028       m_tauLut(std::move(tauLut)),
0029       m_nodes(SpacePointColumns::CopiedFromIndex |
0030               SpacePointColumns::PackedXYZR) {
0031   m_etaBins.resize(m_geometry->numBins());
0032   m_stagedPerBin.resize(m_geometry->numBins());
0033 }
0034
0035 std::optional<std::uint32_t> GbtsNodeStorage::insert(
0036     const SpacePointIndex index, const float x, const float y, const float z,
0037     const std::uint32_t layerIndex, const float clusterWidth,
0038     const float localPositionY) {
0039   const float r = fastHypot(x, y);
0040   const float phi = std::atan2(y, x);
0041   return insert(index, x, y, z, r, phi, layerIndex, clusterWidth,
0042                 localPositionY);
0043 }
0044
0045 std::optional<std::uint32_t> GbtsNodeStorage::insert(
0046     const SpacePointIndex index, const float x, const float y, const float z,
0047     const float r, const float phi, const std::uint32_t layerIndex,
0048     const float clusterWidth, const float localPositionY,
0049     const OuterStripSpacePointCalibrationDetails* strip) {
0050   const detail::GbtsLayer& layer = m_geometry->layerByIndex(layerIndex);
0051   const GbtsLayerDescription& description = layer.layerDescription();
0052
0053   // wide pixel endcap clusters are dropped when the width cuts are on
0054   if (m_cfg.useClusterWidthCuts && description.type == GbtsLayerType::Endcap &&
0055       description.technology == GbtsLayerTechnology::Pixel &&
0056       clusterWidth > m_cfg.maxEndcapClusterWidth) {
0057     return std::nullopt;
0058   }
0059
0060   const std::int32_t binIndex = layer.getEtaBin(z, r);
0061   if (binIndex == -1) {
0062     return std::nullopt;
0063   }
0064
0065   const auto bin = static_cast<std::uint32_t>(binIndex);
0066
0067   std::uint32_t stripIndex = detail::kNoStrip;
0068   // A strip pair on a pixel layer is never read: the seeder takes the strip
0069   // path per bin.
0070   if (strip != nullptr &&
0071       description.technology == GbtsLayerTechnology::Strip) {
0072     stripIndex = static_cast<std::uint32_t>(m_strips.size());
0073     // Derived once here rather than once per pair in the graph: it is six
0074     // cross products and a node takes part in many pairs.
0075     m_strips.push_back(
0076         Acts::detail::deriveOuterStripSpacePointCalibrationDetails(*strip));
0077   }
0078
0079   m_stagedPerBin.at(bin).push_back(static_cast<std::uint32_t>(m_staged.size()));
0080   m_staged.emplace_back(index, x, y, z, r, phi, clusterWidth, localPositionY,
0081                         static_cast<std::uint16_t>(layerIndex), stripIndex);
0082
0083   return bin;
0084 }
0085
0086 void GbtsNodeStorage::extend(
0087     const SpacePointContainer& spacePoints,
0088     const ConstSpacePointColumnProxy<std::uint32_t>& layerColumn,
0089     const ConstSpacePointColumnProxy<float>& clusterWidthColumn,
0090     const ConstSpacePointColumnProxy<float>& localPositionYColumn) {
0091   const bool strips =
0092       spacePoints.hasColumns(SpacePointColumns::StripCalibrationDetails);
0093   m_staged.reserve(m_staged.size() + spacePoints.size());
0094   for (const auto& sp : spacePoints) {
0095     insert(sp, layerColumn, clusterWidthColumn, localPositionYColumn, strips);
0096   }
0097 }
0098
0099 std::vector<std::uint32_t> GbtsNodeStorage::sortBinByPhi(
0100     const std::vector<std::uint32_t>& staged) const {
0101   const std::uint32_t nBuckets = m_cfg.phiSortBuckets;
0102   std::array<std::vector<std::pair<float, std::uint32_t>>,
0103              kMaxPhiSortBuckets + 1>
0104       phiBuckets;
0105
0106   for (const std::uint32_t stagedIdx : staged) {
0107     const float phi = m_staged[stagedIdx].phi;
0108     const auto bIdx = static_cast<std::uint32_t>(
0109         0.5 * nBuckets * (phi / std::numbers::pi_v<float> + 1.0f));
0110     phiBuckets[bIdx].emplace_back(phi, stagedIdx);
0111   }
0112
0113   // Nodes with identical phi are ordered by insertion index.
0114   for (std::uint32_t bucket = 0; bucket <= nBuckets; ++bucket) {
0115     std::ranges::sort(phiBuckets[bucket]);
0116   }
0117
0118   std::vector<std::uint32_t> sorted;
0119   sorted.reserve(staged.size());
0120   for (std::uint32_t bucket = 0; bucket <= nBuckets; ++bucket) {
0121     for (const auto& [phi, stagedIdx] : phiBuckets[bucket]) {
0122       sorted.push_back(stagedIdx);
0123     }
0124   }
0125
0126   return sorted;
0127 }
0128
0129 void GbtsNodeStorage::finalize() {
0130   const auto nNodes = static_cast<std::uint32_t>(m_staged.size());
0131
0132   m_nodes.reserve(nNodes);
0133   m_layers.reserve(nNodes);
0134
0135   // Node order across all bins, used to fill the derived columns below.
0136   std::vector<std::uint32_t> nodeOrder;
0137   nodeOrder.reserve(nNodes);
0138
0139   for (std::uint32_t bin = 0; bin < m_etaBins.size(); ++bin) {
0140     detail::GbtsEtaBinInfo& binInfo = m_etaBins[bin];
0141
0142     binInfo.nodes = {m_nodes.size(), m_nodes.size()};
0143
0144     const std::vector<std::uint32_t>& staged = m_stagedPerBin[bin];
0145     if (staged.empty()) {
0146       continue;
0147     }
0148
0149     const std::vector<std::uint32_t> sorted = sortBinByPhi(staged);
0150
0151     float minRadius = std::numeric_limits<float>::max();
0152     float maxRadius = std::numeric_limits<float>::lowest();
0153
0154     for (const std::uint32_t stagedIdx : sorted) {
0155       const StagedNode& node = m_staged[stagedIdx];
0156
0157       auto newNode = m_nodes.createSpacePoint();
0158       newNode.copiedFromIndex() = node.spacePointIndex;
0159       newNode.xyzr() = std::array<float, 4>{node.x, node.y, node.z, node.r};
0160
0161       m_layers.push_back(node.layer);
0162       nodeOrder.push_back(stagedIdx);
0163
0164       minRadius = std::min(minRadius, node.r);
0165       maxRadius = std::max(maxRadius, node.r);
0166     }
0167
0168     binInfo.nodes.second = m_nodes.size();
0169     binInfo.minRadius = minRadius;
0170     binInfo.maxRadius = maxRadius;
0171     // every node in a bin is on the same layer, so any of them will do
0172     const GbtsLayerDescription& description =
0173         m_geometry->layerDescriptionByIndex(m_staged[staged.front()].layer);
0174     binInfo.layerId = static_cast<std::uint32_t>(description.id);
0175     binInfo.type = description.type;
0176     binInfo.technology = description.technology;
0177   }
0178
0179   // Created now that the container has its final size, so that each column is
0180   // allocated in a single resize.
0181   m_paramsColumn.emplace(
0182       m_nodes.createColumn<detail::GbtsNodeParams>("gbtsNodeParams"));
0183   m_edgeInfoColumn.emplace(
0184       m_nodes.createColumn<detail::GbtsNodeEdgeInfo>("gbtsNodeEdgeInfo"));
0185
0186   // Reordered into node order, so that the pairs a bin reads sit together;
0187   // left empty when nothing carries one, which is what `hasStrips` reports.
0188   std::vector<OuterStripSpacePointCalibrationDetailsDerived> strips;
0189   if (!m_strips.empty()) {
0190     m_stripIndex.assign(nodeOrder.size(), detail::kNoStrip);
0191     strips.reserve(m_strips.size());
0192   }
0193
0194   std::span<detail::GbtsNodeParams> params = m_paramsColumn->data();
0195   for (std::uint32_t node = 0; node < nodeOrder.size(); ++node) {
0196     const StagedNode& staged = m_staged[nodeOrder[node]];
0197
0198     if (staged.strip != detail::kNoStrip) {
0199       m_stripIndex[node] = static_cast<std::uint32_t>(strips.size());
0200       strips.push_back(m_strips[staged.strip]);
0201     }
0202
0203     params[node].phi = staged.phi;
0204     params[node].r = staged.r;
0205     params[node].z = staged.z;
0206
0207     // minTau and maxTau keep their "do not cut" defaults unless the lookup
0208     // table narrows them
0209     if (m_cfg.useClusterWidthCuts) {
0210       applyTauCuts(staged, params[node]);
0211     }
0212   }
0213
0214   m_strips = std::move(strips);
0215
0216   generatePhiIndexing(m_cfg.phiIndexMargin * m_cfg.phiSliceWidth);
0217
0218   m_staged.clear();
0219   m_staged.shrink_to_fit();
0220   m_stagedPerBin.clear();
0221   m_stagedPerBin.shrink_to_fit();
0222 }
0223
0224 void GbtsNodeStorage::applyTauCuts(const StagedNode& staged,
0225                                    detail::GbtsNodeParams& params) const {
0226   const GbtsLayerDescription& description =
0227       m_geometry->layerDescriptionByIndex(staged.layer);
0228
0229   // the table is trained on pixel barrel clusters
0230   if (description.technology != GbtsLayerTechnology::Pixel ||
0231       description.type != GbtsLayerType::Barrel) {
0232     return;
0233   }
0234
0235   // by the reciprocal, not the division: 1/0.05f is exactly 20, the division
0236   // is not, and the difference lands on the bin edges
0237   const auto lutBinIdx =
0238       static_cast<std::int32_t>(
0239           std::floor(staged.clusterWidth * (1.0f / m_cfg.tauLutBinWidth))) -
0240       1;
0241
0242   if (lutBinIdx < 0 ||
0243       lutBinIdx >= static_cast<std::int32_t>(m_tauLut.size())) {
0244     return;
0245   }
0246
0247   const detail::GbtsTauBounds& bounds = m_tauLut[lutBinIdx];
0248
0249   // close to the edge the cluster may be shortened, which the lookup table
0250   // covers with a separate pair of bounds
0251   const float dist2border =
0252       m_cfg.moduleHalfLengthY - std::abs(staged.localPositionY);
0253   const bool nearEdge = dist2border <= m_cfg.moduleEdgeTolerance;
0254
0255   params.minTau = nearEdge ? bounds.minTauNearEdge : bounds.minTau;
0256   params.maxTau = nearEdge ? bounds.maxTauNearEdge : bounds.maxTau;
0257
0258   if (params.maxTau < 0) {
0259     // insufficient training data, do not cut on tau
0260     params.maxTau = std::numeric_limits<float>::infinity();
0261   }
0262 }
0263
0264 void GbtsNodeStorage::generatePhiIndexing(const float dphi) {
0265   const std::span<const detail::GbtsNodeParams> params = m_paramsColumn->data();
0266
0267   for (detail::GbtsEtaBinInfo& bin : m_etaBins) {
0268     if (bin.empty()) {
0269       continue;
0270     }
0271
0272     const SpacePointIndex begin = bin.nodes.first;
0273     const SpacePointIndex end = bin.nodes.second;
0274
0275     for (SpacePointIndex node = begin; node < end; ++node) {
0276       const float phi = params[node].phi;
0277       if (phi <= std::numbers::pi_v<float> - dphi) {
0278         continue;
0279       }
0280       bin.phiNodes.emplace_back(phi - 2 * std::numbers::pi_v<float>, node);
0281     }
0282
0283     for (SpacePointIndex node = begin; node < end; ++node) {
0284       bin.phiNodes.emplace_back(params[node].phi, node);
0285     }
0286
0287     for (SpacePointIndex node = begin; node < end; ++node) {
0288       const float phi = params[node].phi;
0289       if (phi >= -std::numbers::pi_v<float> + dphi) {
0290         break;
0291       }
0292       bin.phiNodes.emplace_back(phi + 2 * std::numbers::pi_v<float>, node);
0293     }
0294   }
0295 }
0296
0297 }  // namespace Acts::Experimental