EIC code displayed by LXR master/​acts/​Core/​src/​Seeding/​GbtsGraphBuilder.cpp	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

Warning, file /acts/Core/src/Seeding/GbtsGraphBuilder.cpp was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

 // This file is part of the ACTS project.
 //
 // Copyright (C) 2016 CERN for the benefit of the ACTS project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
 #include "Acts/Seeding/GbtsGraphBuilder.hpp"
 
 #include "Acts/SpacePointFormation/detail/StripSpacePointCalibrationImpl.hpp"
 #include "Acts/Utilities/MathHelpers.hpp"
 
 #include <algorithm>
 #include <array>
 #include <cmath>
 #include <cstdint>
 #include <numbers>
 #include <span>
 #include <stdexcept>
 #include <utility>
 #include <vector>
 
 namespace Acts::Experimental {
 
 namespace {
 
 /// Sliding window in phi used to define range used for edge creation.
 ///
 /// Covers one non-empty source eta bin, whose phi-ordered node list it holds
 /// directly so that the innermost loop does not reach through the bin.
 struct SlidingWindow {
   /// phi-ordered nodes of the bin, including the wrap-around duplicates
   const std::pair<float, SpacePointIndex>* phiNodes{};
   /// number of entries in @c phiNodes
   std::uint32_t numPhiNodes{};
   /// sliding window position
   std::uint32_t firstIt{};
   /// window half-width;
   float deltaPhi{};
   /// Inside-out pixel barrel ordinal of the bin's layer, -1 for the rest.
   std::int32_t barrelOrder{-1};
   /// Type of the bin's layer.
   GbtsLayerType type{};
   /// Technology of the bin's layer.
   GbtsLayerTechnology technology{};
 };
 
 /// Check a triplet against the pT and d0 cuts.
 /// @param cfg Graph builder configuration
 /// @param nodeView View of the node positions and layers
 /// @param candidateTriplet The three graph nodes
 /// @param tripletMinPt Minimum transverse momentum
 /// @param tauRatio Tau ratio of the triplet
 /// @param tauRatioCut Tau ratio cut threshold
 /// @param bFieldInZ Magnetic field in z, in GeV/(e*mm)
 /// @return Whether the triplet is accepted
 bool validateTriplet(const GbtsGraphBuilder::Config& cfg,
                      const detail::GbtsNodeView& nodeView,
                      const std::array<SpacePointIndex, 3>& candidateTriplet,
                      const float tripletMinPt, const float tauRatio,
                      const float tauRatioCut, const float bFieldInZ) {
   // conformal mapping with the center at the middle spacepoint
 
   std::array<float, 2> u{};
   std::array<float, 2> v{};
 
   const detail::GbtsNodeProxy n0 = nodeView[candidateTriplet[1]];
 
   const float x0 = n0.x();
   const float y0 = n0.y();
 
   const float r0 = n0.r();
 
   const float cosA = x0 / r0;
 
   const float sinA = y0 / r0;
 
   for (std::uint32_t k = 0; k < 2; k++) {
     const std::uint32_t spIdx = (k == 1) ? 2 : k;
 
     const detail::GbtsNodeProxy nk = nodeView[candidateTriplet[spIdx]];
 
     const float dx = nk.x() - x0;
 
     const float dy = nk.y() - y0;
 
     const float r2Inv = 1.0f / (dx * dx + dy * dy);
 
     const float xn = dx * cosA + dy * sinA;
 
     const float yn = -dx * sinA + dy * cosA;
 
     u[k] = xn * r2Inv;
     v[k] = yn * r2Inv;
   }
 
   const float du = u[0] - u[1];
 
   if (du == 0.0) {
     return false;
   }
 
   const float A = (v[0] - v[1]) / du;
 
   const float B = v[1] - A * u[1];
 
   const float d0 = r0 * (B * r0 - A);
 
   if (std::abs(d0) > cfg.d0Max) {
     return false;
   }
 
   if (B != 0.0) {  // straight-line track is OK
 
     const float R = std::sqrt(1 + A * A) / B;
 
     // 1T magnetic field used
     const float pT = std::abs(bFieldInZ * R / 2);
 
     if (pT < tripletMinPt) {
       return false;
     }
 
     if (pT > 5 * tripletMinPt) {  // relatively high-pT track
 
       if (tauRatio > 0.9f * tauRatioCut) {
         return false;
       }
     }
   }
 
   return true;
 }
 
 /// Check to see if z0 of segment is within the expected z range of the
 /// beamspot
 /// @param cfg Graph builder configuration
 /// @param z0BitMask Sets allowed bins of allowed z value
 /// @param z0 Estimated z0 of segments z value at beamspot
 /// @param z0HistoCoeff Scalfactor that converts z coodindate into bin index
 /// @return Whether segment is within beamspot range
 bool checkZ0BitMask(const GbtsGraphBuilder::Config& cfg,
                     const std::uint16_t z0BitMask, const float z0,
                     const float z0HistoCoeff) {
   if (z0BitMask == 0) {
     return true;
   }
 
   // z0 is not yet range checked here -- doubletFilterRZ runs after this -- so
   // a bin outside the histogram is reachable and must not reach the shift.
   const auto isSet = [z0BitMask](const std::int32_t bin) {
     return bin >= 0 && bin < detail::kGbtsZ0HistogramBins &&
            ((z0BitMask >> bin) & 1) != 0;
   };
 
   const float dz = z0 - cfg.minZ0;
   const std::int32_t z0BinIndex = static_cast<std::int32_t>(z0HistoCoeff * dz);
 
   if (isSet(z0BinIndex)) {
     return true;
   }
 
   // check adjacent bins as well
 
   const float dzm = dz - cfg.z0Resolution;
 
   std::int32_t nextBin = static_cast<std::int32_t>(z0HistoCoeff * dzm);
 
   if (nextBin != z0BinIndex && isSet(nextBin)) {
     return true;
   }
 
   const float dzp = dz + cfg.z0Resolution;
 
   nextBin = static_cast<std::int32_t>(z0HistoCoeff * dzp);
 
   if (nextBin != z0BinIndex && isSet(nextBin)) {
     return true;
   }
 
   return false;
 }
 
 }  // namespace
 
 GbtsGraphBuilder::GbtsGraphBuilder(const Config& config,
                                    std::shared_ptr<const GbtsGeometry> geometry,
                                    std::unique_ptr<const Acts::Logger> logger)
     : m_cfg(config),
       m_geometry(std::move(geometry)),
       m_logger(std::move(logger)) {
   // buildTheGraph pre-computes the loosest tau ratio threshold it can apply,
   // which assumes the corrections only ever widen the cut.
   if (m_cfg.tauRatioCorr < 0) {
     throw std::invalid_argument(
         "GbtsGraphBuilder: tauRatioCorr must not be negative");
   }
 
   if (m_cfg.tauRatioCorrStrip < 0) {
     throw std::invalid_argument(
         "GbtsGraphBuilder: tauRatioCorrStrip must not be negative");
   }
 }
 
 GbtsGraph GbtsGraphBuilder::buildTheGraph(const GbtsRoiDescriptor& roi,
                                           GbtsNodeStorage& nodeStorage,
                                           const float bFieldInZ) const {
   GbtsGraph graph;
   std::vector<detail::GbtsEdge>& edgeStorage = graph.edgeStorage;
 
   // used to calculate Z cut on doublets
   const float cutZMinU =
       m_cfg.minZ0 + m_cfg.maxOuterRadius * static_cast<float>(roi.dzdrMin());
   const float cutZMaxU =
       m_cfg.maxZ0 + m_cfg.maxOuterRadius * static_cast<float>(roi.dzdrMax());
 
   const float tripletPtMin = m_cfg.tripletPtFraction * m_cfg.minPt;
 
   const float ptScale = m_cfg.tuningPt / m_cfg.minPt;
 
   const float curvatureCutHighEta = m_cfg.maxCurvatureHighEta * ptScale;
   const float curvatureCutLowEta = m_cfg.maxCurvatureLowEta * ptScale;
 
   // the loosest tau ratio threshold the triplet matching can apply
   const float maxTauRatioCut =
       m_cfg.tauRatioCut + (m_cfg.useAdaptiveCuts ? m_cfg.tauRatioCorr : 0.0f) +
       (nodeStorage.hasStrips() ? m_cfg.tauRatioCorrStrip : 0.0f);
 
   // the default sliding window along phi. Taken from the node storage so that
   // the windows and the phi indexing they slide over cannot disagree.
   const float deltaPhi0 = 0.5f * nodeStorage.config().phiSliceWidth;
 
   std::uint32_t nConnections = 0;
 
   // Scales the maximum number of edges by the number of space points using the
   // factor maxEdgesPerSP
   const auto maxEdges = static_cast<std::size_t>(
       m_cfg.maxEdgesPerSP * static_cast<float>(nodeStorage.numberOfNodes()));
 
   edgeStorage.reserve(maxEdges);
 
   // scale factor to get indexes of binned beamspot
   const float z0HistoCoeff =
       detail::kGbtsZ0HistogramBins / (m_cfg.maxZ0 - m_cfg.minZ0 + 1e-6f);
 
   const detail::GbtsNodeView nodeView = nodeStorage.nodeView();
   const std::span<const detail::GbtsNodeParams> params =
       nodeStorage.nodeParams();
   const std::span<detail::GbtsNodeEdgeInfo> edgeInfo =
       nodeStorage.nodeEdgeInfo();
 
   // reused across bin groups so that the windows are allocated once
   std::vector<SlidingWindow> phiSlidingWindow;
 
   const bool calibrate = m_cfg.calibrateStrips && nodeStorage.hasStrips();
 
   // Put a strip node back where a direction says it crossed. Both ends, since
   // the nominal position sits on the inner strip and the calibrated one on the
   // outer. A direction that misses either strip is no crossing at all.
   const auto calibrateNode = [&](const SpacePointIndex node,
                                  const std::array<float, 3>& direction,
                                  float& r, float& z) {
     std::array<float, 3> point{};
     if (!Acts::detail::calibrateOuterStripSpacePoint(
             direction, nodeStorage.strip(node), point,
             m_cfg.maxStripLengthFraction)) {
       return false;
     }
     r = fastHypot(point[0], point[1]);
     z = point[2];
     return true;
   };
 
   // loop over bin groups
   for (const auto& bg : m_geometry->binGroups()) {
     const detail::GbtsEtaBinInfo& B1 = nodeStorage.etaBin(bg.bin);
 
     if (B1.empty()) {
       continue;
     }
 
     const float rb1 = B1.minRadius;
 
     const std::int32_t barrelOrder1 = B1.barrelOrder;
 
     const bool isPixel1 = B1.technology == GbtsLayerTechnology::Pixel;
     const bool isPixelBarrel1 = barrelOrder1 >= 0;
 
     const bool useZ0Histogram =
         barrelOrder1 >= 0 && barrelOrder1 <= m_cfg.z0HistogramMaxBarrelOrder;
     const bool useMatchBeforeCreate =
         m_cfg.matchBeforeCreate && barrelOrder1 >= 0 &&
         barrelOrder1 <= m_cfg.matchBeforeCreateMaxBarrelOrder;
 
     // prepare a sliding window for each non-empty bin2 in the group
 
     phiSlidingWindow.clear();
 
     // loop over n2 eta-bins in L2 layers
     for (const std::uint32_t b2Idx : bg.links) {
       const detail::GbtsEtaBinInfo& B2 = nodeStorage.etaBin(b2Idx);
 
       if (B2.empty()) {
         continue;
       }
 
       const float rb2 = B2.maxRadius;
 
       float deltaPhi = deltaPhi0;  // the default
 
       // override the default window width
       if (m_cfg.useEtaBinning) {
         const float absDr = std::fabs(rb2 - rb1);
         if (absDr < m_cfg.phiWindowSplitDeltaRadius) {
           deltaPhi = m_cfg.phiWindowNearOffset +
                      m_cfg.phiWindowNearSlope * ptScale * absDr;
         } else {
           deltaPhi = m_cfg.phiWindowFarOffset +
                      m_cfg.phiWindowFarSlope * ptScale * absDr;
         }
       }
 
       SlidingWindow& window = phiSlidingWindow.emplace_back();
       window.phiNodes = B2.phiNodes.data();
       window.numPhiNodes = static_cast<std::uint32_t>(B2.phiNodes.size());
       window.deltaPhi = deltaPhi;
       window.barrelOrder = B2.barrelOrder;
       window.type = B2.type;
       window.technology = B2.technology;
     }
 
     // in GBTSv3 the outer loop goes over n1 nodes in the Layer 1 bin
     for (SpacePointIndex n1Idx = B1.nodes.first; n1Idx < B1.nodes.second;
          ++n1Idx) {
       // initialization using the top watermark of the edge storage
       edgeInfo[n1Idx].firstEdge =
           static_cast<std::uint32_t>(edgeStorage.size());
 
       // the counter for the incoming graph edges created for n1
       std::uint16_t numCreatedEdges = 0;
 
       bool isConnected = false;
 
       std::array<std::uint8_t, detail::kGbtsZ0HistogramBins> z0Histo = {};
 
       const detail::GbtsNodeParams& n1pars = params[n1Idx];
 
       const float phi1 = n1pars.phi;
       const float r1 = n1pars.r;
       const float z1 = n1pars.z;
 
       // the chord of a pair, which only the strip path below reads
       float x1 = 0.f;
       float y1 = 0.f;
       if (calibrate) {
         const std::array<float, 4>& position = nodeView.positions[n1Idx];
         x1 = position[0];
         y1 = position[1];
       }
 
       // the intermediate loop over sliding windows
       for (auto& slw : phiSlidingWindow) {
         const std::int32_t barrelOrder2 = slw.barrelOrder;
 
         const bool isPixel2 = slw.technology == GbtsLayerTechnology::Pixel;
         const bool isPixelBarrel2 = barrelOrder2 >= 0;
 
         const bool stripPair = calibrate && (!isPixel1 || !isPixel2);
 
         const float deltaPhi = slw.deltaPhi;
 
         // sliding window phi1 +/- deltaPhi
 
         const float minPhi = phi1 - deltaPhi;
         const float maxPhi = phi1 + deltaPhi;
 
         // the inner loop over n2 nodes using sliding window
         for (std::uint32_t n2PhiIdx = slw.firstIt; n2PhiIdx < slw.numPhiNodes;
              ++n2PhiIdx) {
           const float phi2 = slw.phiNodes[n2PhiIdx].first;
 
           if (phi2 < minPhi) {
             // update the window position
             slw.firstIt = n2PhiIdx;
             continue;
           }
           if (phi2 > maxPhi) {
             // break and go to the next window
             break;
           }
 
           const SpacePointIndex n2Idx = slw.phiNodes[n2PhiIdx].second;
 
           const detail::GbtsNodeEdgeInfo& n2Info = edgeInfo[n2Idx];
 
           const std::uint16_t nodeInfo = n2Info.isConnected;
 
           // skip isolated nodes as their incoming edges lead to nowhere
           if (useZ0Histogram && (nodeInfo == 0)) {
             continue;
           }
 
           const std::uint32_t n2FirstEdge = n2Info.firstEdge;
           const std::uint16_t n2NumEdges = n2Info.numEdges;
           const std::uint32_t n2LastEdge = n2FirstEdge + n2NumEdges;
 
           const detail::GbtsNodeParams& n2pars = params[n2Idx];
 
           const float r2 = n2pars.r;
 
           float dr = r2 - r1;
 
           // On the nominal radii, so an endcap pair the slide would have
           // opened up is lost here. The cheap reject is worth more.
           if (dr < m_cfg.minDeltaRadius) {
             continue;
           }
 
           const float z2 = n2pars.z;
 
           // the ends as the pair puts them: nominal, or slid along the strip
           // when resolved. Azimuth is kept as it was -- exactly right in the
           // barrel, and to the stereo angle in the endcap, where the strip
           // slid along is a hair off radial.
           float r1c = r1;
           float z1c = z1;
           float r2c = r2;
           float z2c = z2;
 
           if (stripPair) {
             // A circle's chord bisects the tangent directions, so turning the
             // chord by the pair's own azimuth difference gives each tangent.
             const std::array<float, 4>& position = nodeView.positions[n2Idx];
             const float dx = position[0] - x1;
             const float dy = position[1] - y1;
             const float dzc = z2 - z1;
             const float turn = phi2 - phi1;
 
             bool resolved = true;
             if (!isPixel1) {
               resolved = calibrateNode(
                   n1Idx, {dx + dy * turn, dy - dx * turn, dzc}, r1c, z1c);
             }
             if (resolved && !isPixel2) {
               resolved = calibrateNode(
                   n2Idx, {dx - dy * turn, dy + dx * turn, dzc}, r2c, z2c);
             }
             if (!resolved) {
               continue;
             }
 
             dr = r2c - r1c;
             if (dr < m_cfg.minDeltaRadius) {
               continue;
             }
           }
 
           const float dz = z2c - z1c;
           const float tau = dz / dr;
           const float ftau = std::fabs(tau);
           if (ftau > m_cfg.maxAbsTau) {
             continue;
           }
 
           if (ftau < n1pars.minTau) {
             continue;
           }
           if (ftau > n1pars.maxTau) {
             continue;
           }
 
           if (ftau < n2pars.minTau) {
             continue;
           }
           if (ftau > n2pars.maxTau) {
             continue;
           }
 
           const float z0 = z1c - r1c * tau;
 
           // check against a non-empty z0 histogram
           if (useZ0Histogram &&
               !checkZ0BitMask(m_cfg, nodeInfo, z0, z0HistoCoeff)) {
             continue;
           }
 
           if (m_cfg.doubletFilterRZ) {
             if (z0 < m_cfg.minZ0 || z0 > m_cfg.maxZ0) {
               continue;
             }
 
             const float zOuter = z0 + m_cfg.maxOuterRadius * tau;
 
             if (zOuter < cutZMinU || zOuter > cutZMaxU) {
               continue;
             }
           }
 
           const float curv = (phi2 - phi1) / dr;
 
           if (std::abs(curv) > (ftau < m_cfg.curvatureSplitAbsTau
                                     ? curvatureCutLowEta
                                     : curvatureCutHighEta)) {
             continue;
           }
 
           const float hypotTau = fastHypot(1, tau);
           const float expEta = hypotTau - tau;
           // 1 / expEta, since (hypotTau - tau) * (hypotTau + tau) == 1. The
           // sum is also the better conditioned form for large tau.
           const float invExpEta = hypotTau + tau;
 
           // match edge candidate against edges incoming to n2
           if (useMatchBeforeCreate) {
             // we must have enough incoming edges to decide
             bool isGood = n2NumEdges <= m_cfg.matchBeforeCreateMaxEdges;
 
             if (!isGood) {
               const float uat1 = invExpEta;
 
               for (std::uint32_t n2InIdx = n2FirstEdge; n2InIdx < n2LastEdge;
                    ++n2InIdx) {
                 const float tau2 = edgeStorage[n2InIdx].p[0];
                 const float tauRatio = tau2 * uat1 - 1.0f;
 
                 if (std::abs(tauRatio) > m_cfg.tauRatioPrecut) {  // bad match
                   continue;
                 }
                 isGood = true;  // good match found
                 break;
               }
             }
 
             if (!isGood) {  // no match found, skip creating [n1 <- n2] edge
               continue;
             }
           }
 
           const float dPhi2 = curv * r2c;
           const float dPhi1 = curv * r1c;
 
           if (edgeStorage.size() < maxEdges) {
             const auto outEdgeIdx =
                 static_cast<std::uint32_t>(edgeStorage.size());
 
             edgeStorage.emplace_back(n1Idx, n2Idx, barrelOrder2, expEta, curv,
                                      phi1 + dPhi1);
 
             ++numCreatedEdges;
 
             const float uat2 = invExpEta;
             const float phi2u = phi2 + dPhi2;
             const float curv2 = curv;
 
             // looking for neighbours of the new edge
             for (std::uint32_t inEdgeIdx = n2FirstEdge; inEdgeIdx < n2LastEdge;
                  ++inEdgeIdx) {
               detail::GbtsEdge* pS = &edgeStorage[inEdgeIdx];
 
               const float absTauRatio = std::abs(pS->p[0] * uat2 - 1.0f);
 
               // rejects most candidates before the layer bookkeeping below
               if (absTauRatio > maxTauRatioCut) {
                 continue;
               }
 
               if (pS->nNei >= detail::kGbtsMaxEdgeNeighbours) {
                 continue;
               }
 
               const std::int32_t barrelOrder3 = pS->n2BarrelOrder;
 
               const bool isPixelBarrel3 = barrelOrder3 >= 0;
 
               float addTauRatioCorr = 0;
 
               if (m_cfg.useAdaptiveCuts) {
                 if (isPixelBarrel1 && isPixelBarrel2 && isPixelBarrel3) {
                   // three radially consecutive layers, none skipped
                   const bool noGap = (barrelOrder2 - barrelOrder1) == 1 &&
                                      (barrelOrder3 - barrelOrder2) == 1;
 
                   // assume more scattering due to the layer in between
                   if (!noGap) {
                     addTauRatioCorr = m_cfg.tauRatioCorr;
                   }
                 } else {
                   bool mixedTriplet =
                       isPixelBarrel1 && isPixelBarrel2 && !isPixelBarrel3;
                   if (mixedTriplet) {
                     addTauRatioCorr = m_cfg.tauRatioCorr;
                   }
                 }
               }
               // The two doublets sharing a strip node resolved it separately,
               // so a triplet through a strip may disagree on tau by more. Any
               // of the three: the outer two carry their end's error into tau.
               if (m_cfg.tauRatioCorrStrip > 0.f &&
                   (!isPixel1 || !isPixel2 ||
                    nodeView.strip(pS->n2) != nullptr)) {
                 addTauRatioCorr += m_cfg.tauRatioCorrStrip;
               }
 
               // bad match
               if (absTauRatio > m_cfg.tauRatioCut + addTauRatioCorr) {
                 continue;
               }
 
               float dPhi = phi2u - pS->p[2];
 
               if (dPhi < -std::numbers::pi_v<float>) {
                 dPhi += 2 * std::numbers::pi_v<float>;
               } else if (dPhi > std::numbers::pi_v<float>) {
                 dPhi -= 2 * std::numbers::pi_v<float>;
               }
 
               if (std::abs(dPhi) > m_cfg.cutDPhiMax) {
                 continue;
               }
 
               const float dcurv = curv2 - pS->p[1];
 
               if (dcurv < -m_cfg.cutDCurvMax || dcurv > m_cfg.cutDCurvMax) {
                 continue;
               }
 
               // final check: cuts on pT and d0
               if (m_cfg.validateTriplets && isPixelBarrel1 && isPixelBarrel2 &&
                   isPixelBarrel3) {
                 const std::array<SpacePointIndex, 3> candidateTriplet = {
                     n1Idx, n2Idx, pS->n2};
 
                 if (!validateTriplet(m_cfg, nodeView, candidateTriplet,
                                      tripletPtMin, absTauRatio,
                                      m_cfg.tauRatioCut, bFieldInZ)) {
                   continue;
                 }
               }
 
               pS->vNei[pS->nNei] = outEdgeIdx;
               ++pS->nNei;
 
               isConnected = true;  // there is at least one good match
 
               // edge confirmed - update z0 histogram. Only doubletFilterRZ
               // holds z0 to the histogram range, and it is optional.
               if (const auto z0BinIndex = static_cast<std::int32_t>(
                       z0HistoCoeff * (z0 - m_cfg.minZ0));
                   z0BinIndex >= 0 &&
                   z0BinIndex < detail::kGbtsZ0HistogramBins) {
                 ++z0Histo[static_cast<std::size_t>(z0BinIndex)];
               }
 
               nConnections++;
             }
           }
         }  // loop over n2 (outer) nodes inside a sliding window on n2 bin
       }  // loop over sliding windows associated with n2 bins
 
       // updating the n1 node attributes
 
       edgeInfo[n1Idx].numEdges = numCreatedEdges;
       if (isConnected) {
         std::uint16_t z0BitMask = 0x0;
 
         for (std::int32_t bIdx = 0; bIdx < detail::kGbtsZ0HistogramBins;
              ++bIdx) {
           if (z0Histo[static_cast<std::size_t>(bIdx)] == 0) {
             continue;
           }
 
           z0BitMask |= (1 << bIdx);
         }
 
         // non-zero mask indicates that there is at least one connected edge
         edgeInfo[n1Idx].isConnected = z0BitMask;
       }
 
     }  // loop over n1 (inner) nodes
   }  // loop over bin groups: a single n1 bin and multiple n2 bins
 
   if (edgeStorage.size() >= maxEdges) {
     ACTS_WARNING(
         "Maximum number of graph edges exceeded - possible efficiency loss "
         << edgeStorage.size());
   }
   graph.nConnections = nConnections;
   return graph;
 }
 
 std::uint32_t GbtsGraph::runCCA(const std::uint32_t maxIterations) {
   const auto nEdges = static_cast<std::uint32_t>(edgeStorage.size());
 
   std::uint32_t maxLevel = 0;
 
   std::uint32_t iter = 0;
 
   std::vector<detail::GbtsEdge*> vOld;
 
   for (std::uint32_t edgeIndex = 0; edgeIndex < nEdges; ++edgeIndex) {
     detail::GbtsEdge* pS = &(edgeStorage[edgeIndex]);
     if (pS->nNei == 0) {
       continue;
     }
 
     // TODO: increment level for segments as they already have at least one
     // neighbour
     vOld.push_back(pS);
   }
 
   std::vector<detail::GbtsEdge*> vNew;
   vNew.reserve(vOld.size());
 
   // generate proposals
   for (; iter < maxIterations; iter++) {
     vNew.clear();
 
     for (detail::GbtsEdge* pS : vOld) {
       std::int32_t nextLevel = pS->level;
 
       for (std::uint32_t nIdx = 0; nIdx < pS->nNei; ++nIdx) {
         const std::uint32_t nextEdgeIdx = pS->vNei[nIdx];
 
         const detail::GbtsEdge* pN = &(edgeStorage[nextEdgeIdx]);
 
         if (pS->level == pN->level) {
           nextLevel = pS->level + 1;
           vNew.push_back(pS);
           break;
         }
       }
 
       // proposal
       pS->next = static_cast<std::int8_t>(nextLevel);
     }
 
     // update
 
     std::uint32_t nChanges = 0;
 
     for (auto pS : vNew) {
       if (pS->next != pS->level) {
         nChanges++;
         pS->level = pS->next;
         // levels only grow from zero here, so the cast is safe
         maxLevel = std::max(maxLevel, static_cast<std::uint32_t>(pS->level));
       }
     }
 
     if (nChanges == 0) {
       break;
     }
 
     vOld.swap(vNew);
     vNew.clear();
   }
 
   return maxLevel;
 }
 
 std::vector<detail::GbtsEdge*> GbtsGraph::extractChainHeads(
     const std::uint8_t minSeedLevel, const bool addTriplets,
     const float maxAbsEtaAddTriplets) {
   const auto nEdges = static_cast<std::uint32_t>(edgeStorage.size());
 
   // `addTriplets` accepts a chain one level short. Signed: an uncollected
   // edge sits at level -1 and `minSeedLevel` may be configured to 0.
   const std::int8_t minLevelAddTriplets =
       static_cast<std::int8_t>(minSeedLevel - 1);
   std::vector<detail::GbtsEdge*> vChainHeads;
 
   vChainHeads.reserve(nEdges / 2);
 
   for (std::uint32_t edgeIndex = 0; edgeIndex < nEdges; ++edgeIndex) {
     detail::GbtsEdge* pS = &edgeStorage[edgeIndex];
 
     if (!addTriplets) {
       if (pS->level < minSeedLevel) {
         continue;
       }
     } else {  // eta-dependent cut
       const float edgeAbsEta = std::abs(-std::log(pS->p[0]));
 
       if (edgeAbsEta > maxAbsEtaAddTriplets) {
         if (pS->level < minSeedLevel) {
           continue;
         }
       } else {
         if (pS->level < minLevelAddTriplets) {
           continue;
         }
       }
     }
 
     vChainHeads.push_back(pS);
   }
 
   std::ranges::sort(vChainHeads, std::ranges::greater{},
                     [](const detail::GbtsEdge* e) { return e->level; });
 
   return vChainHeads;
 }
 
 }  // namespace Acts::Experimental

This page was automatically generated by the 2.3.7 LXR engine. The LXR team