TrackFinding/src/GridTripletSeedingAlgorithm.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 "ActsExamples/TrackFinding/GridTripletSeedingAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Direction.hpp"
0012 #include "Acts/EventData/SeedContainer2.hpp"
0013 #include "Acts/EventData/SourceLink.hpp"
0014 #include "Acts/EventData/SpacePointContainer2.hpp"
0015 #include "Acts/EventData/Types.hpp"
0016 #include "Acts/Seeding2/BroadTripletSeedFilter.hpp"
0017 #include "Acts/Seeding2/DoubletSeedFinder.hpp"
0018 #include "Acts/Seeding2/TripletSeedFinder.hpp"
0019 #include "Acts/Utilities/Delegate.hpp"
0020 #include "ActsExamples/EventData/SimSeed.hpp"
0021 #include "ActsExamples/EventData/SimSpacePoint.hpp"
0022
0023 #include <cmath>
0024 #include <csignal>
0025 #include <cstddef>
0026 #include <stdexcept>
0027
0028 namespace ActsExamples {
0029
0030 namespace {
0031
0032 static inline bool itkFastTrackingCuts(
0033     const Acts::ConstSpacePointProxy2& /*middle*/,
0034     const Acts::ConstSpacePointProxy2& other, float cotTheta,
0035     bool isBottomCandidate) {
0036   static float rMin = 45;
0037   static float cotThetaMax = 1.5;
0038
0039   if (isBottomCandidate && other.zr()[1] < rMin &&
0040       (cotTheta > cotThetaMax || cotTheta < -cotThetaMax)) {
0041     return false;
0042   }
0043   return true;
0044 }
0045
0046 static inline bool itkFastTrackingSPselect(const SimSpacePoint& sp) {
0047   // At small r we remove points beyond |z| > 200.
0048   float r = sp.r();
0049   float zabs = std::abs(sp.z());
0050   if (zabs > 200. && r < 45.) {
0051     return false;
0052   }
0053
0054   // Remove space points beyond eta=4 if their z is larger than the max seed
0055   // z0 (150.)
0056   float cotTheta = 27.2899;  // corresponds to eta=4
0057   if ((zabs - 150.) > cotTheta * r) {
0058     return false;
0059   }
0060   return true;
0061 }
0062
0063 }  // namespace
0064
0065 GridTripletSeedingAlgorithm::GridTripletSeedingAlgorithm(
0066     const Config& cfg, Acts::Logging::Level lvl)
0067     : IAlgorithm("GridTripletSeedingAlgorithm", lvl), m_cfg(cfg) {
0068   m_inputSpacePoints.initialize(m_cfg.inputSpacePoints);
0069   m_outputSeeds.initialize(m_cfg.outputSeeds);
0070
0071   // check that the bins required in the custom bin looping
0072   // are contained in the bins defined by the total number of edges
0073   for (std::size_t i : m_cfg.zBinsCustomLooping) {
0074     if (i >= m_cfg.zBinEdges.size()) {
0075       throw std::invalid_argument(
0076           "Inconsistent config zBinsCustomLooping does not contain a subset "
0077           "of bins defined by zBinEdges");
0078     }
0079   }
0080
0081   if (m_cfg.useExtraCuts) {
0082     // This function will be applied to select space points during grid filling
0083     m_spacePointSelector.connect<itkFastTrackingSPselect>();
0084   }
0085
0086   m_gridConfig.minPt = m_cfg.minPt;
0087   // TODO switch to `m_cfg.rMin`
0088   m_gridConfig.rMin = 0;
0089   m_gridConfig.rMax = m_cfg.rMax;
0090   m_gridConfig.zMin = m_cfg.zMin;
0091   m_gridConfig.zMax = m_cfg.zMax;
0092   m_gridConfig.deltaRMax = m_cfg.deltaRMax;
0093   m_gridConfig.cotThetaMax = m_cfg.cotThetaMax;
0094   m_gridConfig.impactMax = m_cfg.impactMax;
0095   m_gridConfig.phiMin = m_cfg.phiMin;
0096   m_gridConfig.phiMax = m_cfg.phiMax;
0097   m_gridConfig.phiBinDeflectionCoverage = m_cfg.phiBinDeflectionCoverage;
0098   m_gridConfig.maxPhiBins = m_cfg.maxPhiBins;
0099   m_gridConfig.rBinEdges = {};
0100   m_gridConfig.zBinEdges = m_cfg.zBinEdges;
0101   m_gridConfig.bFieldInZ = m_cfg.bFieldInZ;
0102   m_gridConfig.bottomBinFinder.emplace(m_cfg.numPhiNeighbors,
0103                                        m_cfg.zBinNeighborsBottom, 0);
0104   m_gridConfig.topBinFinder.emplace(m_cfg.numPhiNeighbors,
0105                                     m_cfg.zBinNeighborsTop, 0);
0106   m_gridConfig.navigation[0ul] = {};
0107   m_gridConfig.navigation[1ul] = m_cfg.zBinsCustomLooping;
0108   m_gridConfig.navigation[2ul] = {};
0109
0110   m_filterConfig.deltaInvHelixDiameter = m_cfg.deltaInvHelixDiameter;
0111   m_filterConfig.deltaRMin = m_cfg.deltaRMin;
0112   m_filterConfig.compatSeedWeight = m_cfg.compatSeedWeight;
0113   m_filterConfig.impactWeightFactor = m_cfg.impactWeightFactor;
0114   m_filterConfig.zOriginWeightFactor = m_cfg.zOriginWeightFactor;
0115   m_filterConfig.maxSeedsPerSpM = m_cfg.maxSeedsPerSpM;
0116   m_filterConfig.compatSeedLimit = m_cfg.compatSeedLimit;
0117   m_filterConfig.seedWeightIncrement = m_cfg.seedWeightIncrement;
0118   m_filterConfig.numSeedIncrement = m_cfg.numSeedIncrement;
0119   m_filterConfig.seedConfirmation = m_cfg.seedConfirmation;
0120   m_filterConfig.centralSeedConfirmationRange =
0121       m_cfg.centralSeedConfirmationRange;
0122   m_filterConfig.forwardSeedConfirmationRange =
0123       m_cfg.forwardSeedConfirmationRange;
0124   m_filterConfig.maxSeedsPerSpMConf = m_cfg.maxSeedsPerSpMConf;
0125   m_filterConfig.maxQualitySeedsPerSpMConf = m_cfg.maxQualitySeedsPerSpMConf;
0126   m_filterConfig.useDeltaRinsteadOfTopRadius =
0127       m_cfg.useDeltaRinsteadOfTopRadius;
0128
0129   m_filterLogger = logger().cloneWithSuffix("Filter");
0130
0131   m_seedFinder = Acts::TripletSeeder(logger().cloneWithSuffix("Finder"));
0132 }
0133
0134 ProcessCode GridTripletSeedingAlgorithm::execute(
0135     const AlgorithmContext& ctx) const {
0136   const SimSpacePointContainer& spacePoints = m_inputSpacePoints(ctx);
0137
0138   Acts::CylindricalSpacePointGrid2 grid(m_gridConfig,
0139                                         logger().cloneWithSuffix("Grid"));
0140
0141   for (std::size_t i = 0; i < spacePoints.size(); ++i) {
0142     const auto& sp = spacePoints[i];
0143
0144     // check if the space point passes the selection
0145     if (m_spacePointSelector.connected() && !m_spacePointSelector(sp)) {
0146       continue;
0147     }
0148
0149     float phi = std::atan2(sp.y(), sp.x());
0150     grid.insert(i, phi, sp.z(), sp.r());
0151   }
0152
0153   for (std::size_t i = 0; i < grid.numberOfBins(); ++i) {
0154     std::ranges::sort(grid.at(i), [&](const Acts::SpacePointIndex2& a,
0155                                       const Acts::SpacePointIndex2& b) {
0156       return spacePoints[a].r() < spacePoints[b].r();
0157     });
0158   }
0159
0160   Acts::SpacePointContainer2 coreSpacePoints(
0161       Acts::SpacePointColumns::SourceLinks | Acts::SpacePointColumns::XY |
0162       Acts::SpacePointColumns::ZR | Acts::SpacePointColumns::VarianceZ |
0163       Acts::SpacePointColumns::VarianceR);
0164   coreSpacePoints.reserve(grid.numberOfSpacePoints());
0165   std::vector<Acts::SpacePointIndexRange2> gridSpacePointRanges;
0166   gridSpacePointRanges.reserve(grid.numberOfBins());
0167   for (std::size_t i = 0; i < grid.numberOfBins(); ++i) {
0168     std::uint32_t begin = coreSpacePoints.size();
0169     for (Acts::SpacePointIndex2 spIndex : grid.at(i)) {
0170       const SimSpacePoint& sp = spacePoints[spIndex];
0171
0172       auto newSp = coreSpacePoints.createSpacePoint();
0173       newSp.assignSourceLinks(
0174           std::array<Acts::SourceLink, 1>{Acts::SourceLink(&sp)});
0175       newSp.xy() = std::array<float, 2>{static_cast<float>(sp.x()),
0176                                         static_cast<float>(sp.y())};
0177       newSp.zr() = std::array<float, 2>{static_cast<float>(sp.z()),
0178                                         static_cast<float>(sp.r())};
0179       newSp.varianceZ() = static_cast<float>(sp.varianceZ());
0180       newSp.varianceR() = static_cast<float>(sp.varianceR());
0181     }
0182     std::uint32_t end = coreSpacePoints.size();
0183     gridSpacePointRanges.emplace_back(begin, end);
0184   }
0185
0186   // Compute radius range. We rely on the fact the grid is storing the proxies
0187   // with a sorting in the radius
0188   const Acts::Range1D<float> rRange = [&]() -> Acts::Range1D<float> {
0189     float minRange = std::numeric_limits<float>::max();
0190     float maxRange = std::numeric_limits<float>::lowest();
0191     for (const Acts::SpacePointIndexRange2& range : gridSpacePointRanges) {
0192       if (range.first == range.second) {
0193         continue;
0194       }
0195       auto first = coreSpacePoints[range.first];
0196       auto last = coreSpacePoints[range.second - 1];
0197       minRange = std::min(first.zr()[1], minRange);
0198       maxRange = std::max(last.zr()[1], maxRange);
0199     }
0200     return {minRange, maxRange};
0201   }();
0202
0203   Acts::DoubletSeedFinder::Config bottomDoubletFinderConfig;
0204   bottomDoubletFinderConfig.spacePointsSortedByRadius = true;
0205   bottomDoubletFinderConfig.candidateDirection = Acts::Direction::Backward();
0206   bottomDoubletFinderConfig.deltaRMin = std::isnan(m_cfg.deltaRMaxBottom)
0207                                             ? m_cfg.deltaRMin
0208                                             : m_cfg.deltaRMinBottom;
0209   bottomDoubletFinderConfig.deltaRMax = std::isnan(m_cfg.deltaRMaxBottom)
0210                                             ? m_cfg.deltaRMax
0211                                             : m_cfg.deltaRMaxBottom;
0212   bottomDoubletFinderConfig.deltaZMin = m_cfg.deltaZMin;
0213   bottomDoubletFinderConfig.deltaZMax = m_cfg.deltaZMax;
0214   bottomDoubletFinderConfig.impactMax = m_cfg.impactMax;
0215   bottomDoubletFinderConfig.interactionPointCut = m_cfg.interactionPointCut;
0216   bottomDoubletFinderConfig.collisionRegionMin = m_cfg.collisionRegionMin;
0217   bottomDoubletFinderConfig.collisionRegionMax = m_cfg.collisionRegionMax;
0218   bottomDoubletFinderConfig.cotThetaMax = m_cfg.cotThetaMax;
0219   bottomDoubletFinderConfig.minPt = m_cfg.minPt;
0220   bottomDoubletFinderConfig.helixCutTolerance = m_cfg.helixCutTolerance;
0221   if (m_cfg.useExtraCuts) {
0222     bottomDoubletFinderConfig.experimentCuts.connect<itkFastTrackingCuts>();
0223   }
0224   auto bottomDoubletFinder =
0225       Acts::DoubletSeedFinder::create(Acts::DoubletSeedFinder::DerivedConfig(
0226           bottomDoubletFinderConfig, m_cfg.bFieldInZ));
0227
0228   Acts::DoubletSeedFinder::Config topDoubletFinderConfig =
0229       bottomDoubletFinderConfig;
0230   topDoubletFinderConfig.candidateDirection = Acts::Direction::Forward();
0231   topDoubletFinderConfig.deltaRMin =
0232       std::isnan(m_cfg.deltaRMaxTop) ? m_cfg.deltaRMin : m_cfg.deltaRMinTop;
0233   topDoubletFinderConfig.deltaRMax =
0234       std::isnan(m_cfg.deltaRMaxTop) ? m_cfg.deltaRMax : m_cfg.deltaRMaxTop;
0235   auto topDoubletFinder =
0236       Acts::DoubletSeedFinder::create(Acts::DoubletSeedFinder::DerivedConfig(
0237           topDoubletFinderConfig, m_cfg.bFieldInZ));
0238
0239   Acts::TripletSeedFinder::Config tripletFinderConfig;
0240   tripletFinderConfig.useStripInfo = false;
0241   tripletFinderConfig.sortedByCotTheta = true;
0242   tripletFinderConfig.minPt = m_cfg.minPt;
0243   tripletFinderConfig.sigmaScattering = m_cfg.sigmaScattering;
0244   tripletFinderConfig.radLengthPerSeed = m_cfg.radLengthPerSeed;
0245   tripletFinderConfig.impactMax = m_cfg.impactMax;
0246   tripletFinderConfig.helixCutTolerance = m_cfg.helixCutTolerance;
0247   tripletFinderConfig.toleranceParam = m_cfg.toleranceParam;
0248   auto tripletFinder =
0249       Acts::TripletSeedFinder::create(Acts::TripletSeedFinder::DerivedConfig(
0250           tripletFinderConfig, m_cfg.bFieldInZ));
0251
0252   // variable middle SP radial region of interest
0253   Acts::Range1D<float> rMiddleSpRange = {
0254       std::floor(rRange.min() / 2) * 2 + m_cfg.deltaRMiddleMinSPRange,
0255       std::floor(rRange.max() / 2) * 2 - m_cfg.deltaRMiddleMaxSPRange};
0256
0257   // run the seeding
0258   Acts::SeedContainer2 seeds;
0259   Acts::BroadTripletSeedFilter::State filterState;
0260   Acts::BroadTripletSeedFilter::Cache filterCache;
0261   Acts::BroadTripletSeedFilter seedFilter(m_filterConfig, filterState,
0262                                           filterCache, *m_filterLogger);
0263   static thread_local Acts::TripletSeeder::Cache cache;
0264
0265   std::vector<Acts::SpacePointContainer2::ConstRange> bottomSpRanges;
0266   std::optional<Acts::SpacePointContainer2::ConstRange> middleSpRange;
0267   std::vector<Acts::SpacePointContainer2::ConstRange> topSpRanges;
0268
0269   for (const auto [bottom, middle, top] : grid.binnedGroup()) {
0270     ACTS_VERBOSE("Process middle " << middle);
0271
0272     bottomSpRanges.clear();
0273     for (const auto b : bottom) {
0274       bottomSpRanges.push_back(
0275           coreSpacePoints.range(gridSpacePointRanges.at(b)).asConst());
0276     }
0277     middleSpRange =
0278         coreSpacePoints.range(gridSpacePointRanges.at(middle)).asConst();
0279     topSpRanges.clear();
0280     for (const auto t : top) {
0281       topSpRanges.push_back(
0282           coreSpacePoints.range(gridSpacePointRanges.at(t)).asConst());
0283     }
0284
0285     if (middleSpRange->empty()) {
0286       ACTS_DEBUG("No middle space points in this group, skipping");
0287       continue;
0288     }
0289
0290     // we compute this here since all middle space point candidates belong to
0291     // the same z-bin
0292     Acts::ConstSpacePointProxy2 firstMiddleSp = middleSpRange->front();
0293     std::pair<float, float> radiusRangeForMiddle =
0294         retrieveRadiusRangeForMiddle(firstMiddleSp, rMiddleSpRange);
0295     ACTS_VERBOSE("Validity range (radius) for the middle space point is ["
0296                  << radiusRangeForMiddle.first << ", "
0297                  << radiusRangeForMiddle.second << "]");
0298
0299     m_seedFinder->createSeedsFromGroups(
0300         cache, *bottomDoubletFinder, *topDoubletFinder, *tripletFinder,
0301         seedFilter, coreSpacePoints, bottomSpRanges, *middleSpRange,
0302         topSpRanges, radiusRangeForMiddle, seeds);
0303   }
0304
0305   ACTS_DEBUG("Created " << seeds.size() << " track seeds from "
0306                         << spacePoints.size() << " space points");
0307
0308   // we have seeds of proxies
0309   // convert them to seed of external space points
0310   SimSeedContainer seedContainerForStorage;
0311   seedContainerForStorage.reserve(seeds.size());
0312   for (const auto& seed : seeds) {
0313     auto sps = seed.spacePointIndices();
0314     seedContainerForStorage.emplace_back(*coreSpacePoints.at(sps[0])
0315                                               .sourceLinks()[0]
0316                                               .get<const SimSpacePoint*>(),
0317                                          *coreSpacePoints.at(sps[1])
0318                                               .sourceLinks()[0]
0319                                               .get<const SimSpacePoint*>(),
0320                                          *coreSpacePoints.at(sps[2])
0321                                               .sourceLinks()[0]
0322                                               .get<const SimSpacePoint*>());
0323     seedContainerForStorage.back().setVertexZ(seed.vertexZ());
0324     seedContainerForStorage.back().setQuality(seed.quality());
0325   }
0326
0327   m_outputSeeds(ctx, std::move(seedContainerForStorage));
0328   return ProcessCode::SUCCESS;
0329 }
0330
0331 std::pair<float, float>
0332 GridTripletSeedingAlgorithm::retrieveRadiusRangeForMiddle(
0333     const Acts::ConstSpacePointProxy2& spM,
0334     const Acts::Range1D<float>& rMiddleSpRange) const {
0335   if (m_cfg.useVariableMiddleSPRange) {
0336     return {rMiddleSpRange.min(), rMiddleSpRange.max()};
0337   }
0338   if (m_cfg.rRangeMiddleSP.empty()) {
0339     return {m_cfg.rMinMiddle, m_cfg.rMaxMiddle};
0340   }
0341
0342   // get zBin position of the middle SP
0343   auto pVal = std::ranges::lower_bound(m_cfg.zBinEdges, spM.zr()[0]);
0344   std::size_t zBin = std::distance(m_cfg.zBinEdges.begin(), pVal);
0345   // protects against zM at the limit of zBinEdges
0346   zBin == 0 ? zBin : --zBin;
0347   return {m_cfg.rRangeMiddleSP[zBin][0], m_cfg.rRangeMiddleSP[zBin][1]};
0348 }
0349
0350 }  // namespace ActsExamples