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/Seeding2/BroadTripletSeedFilter.hpp"
0016 #include "Acts/Seeding2/BroadTripletSeedFinder.hpp"
0017 #include "Acts/Seeding2/DoubletSeedFinder.hpp"
0018 #include "Acts/Utilities/Delegate.hpp"
0019 #include "ActsExamples/EventData/SimSeed.hpp"
0020 #include "ActsExamples/EventData/SimSpacePoint.hpp"
0021
0022 #include <cmath>
0023 #include <csignal>
0024 #include <cstddef>
0025 #include <stdexcept>
0026
0027 namespace ActsExamples {
0028
0029 GridTripletSeedingAlgorithm::GridTripletSeedingAlgorithm(
0030     const Config& cfg, Acts::Logging::Level lvl)
0031     : IAlgorithm("GridTripletSeedingAlgorithm", lvl), m_cfg(cfg) {
0032   m_inputSpacePoints.initialize(m_cfg.inputSpacePoints);
0033   m_outputSeeds.initialize(m_cfg.outputSeeds);
0034
0035   // check that the bins required in the custom bin looping
0036   // are contained in the bins defined by the total number of edges
0037   for (std::size_t i : m_cfg.zBinsCustomLooping) {
0038     if (i >= m_cfg.zBinEdges.size()) {
0039       throw std::invalid_argument(
0040           "Inconsistent config zBinsCustomLooping does not contain a subset "
0041           "of bins defined by zBinEdges");
0042     }
0043   }
0044
0045   if (m_cfg.useExtraCuts) {
0046     // This function will be applied to select space points during grid filling
0047     m_spacePointSelector.connect<itkFastTrackingSPselect>();
0048   }
0049
0050   m_gridConfig.minPt = m_cfg.minPt;
0051   // TODO switch to `m_cfg.rMin`
0052   m_gridConfig.rMin = 0;
0053   m_gridConfig.rMax = m_cfg.rMax;
0054   m_gridConfig.zMin = m_cfg.zMin;
0055   m_gridConfig.zMax = m_cfg.zMax;
0056   m_gridConfig.deltaRMax = m_cfg.deltaRMax;
0057   m_gridConfig.cotThetaMax = m_cfg.cotThetaMax;
0058   m_gridConfig.impactMax = m_cfg.impactMax;
0059   m_gridConfig.phiMin = m_cfg.phiMin;
0060   m_gridConfig.phiMax = m_cfg.phiMax;
0061   m_gridConfig.phiBinDeflectionCoverage = m_cfg.phiBinDeflectionCoverage;
0062   m_gridConfig.maxPhiBins = m_cfg.maxPhiBins;
0063   m_gridConfig.rBinEdges = {};
0064   m_gridConfig.zBinEdges = m_cfg.zBinEdges;
0065   m_gridConfig.bFieldInZ = m_cfg.bFieldInZ;
0066   m_gridConfig.bottomBinFinder.emplace(m_cfg.numPhiNeighbors,
0067                                        m_cfg.zBinNeighborsBottom, 0);
0068   m_gridConfig.topBinFinder.emplace(m_cfg.numPhiNeighbors,
0069                                     m_cfg.zBinNeighborsTop, 0);
0070   m_gridConfig.navigation[0ul] = {};
0071   m_gridConfig.navigation[1ul] = m_cfg.zBinsCustomLooping;
0072   m_gridConfig.navigation[2ul] = {};
0073
0074   m_seedFinder = Acts::Experimental::BroadTripletSeedFinder(
0075       logger().cloneWithSuffix("Finder"));
0076
0077   Acts::Experimental::BroadTripletSeedFilter::Config filterConfig;
0078   filterConfig.deltaInvHelixDiameter = m_cfg.deltaInvHelixDiameter;
0079   filterConfig.deltaRMin = m_cfg.deltaRMin;
0080   filterConfig.compatSeedWeight = m_cfg.compatSeedWeight;
0081   filterConfig.impactWeightFactor = m_cfg.impactWeightFactor;
0082   filterConfig.zOriginWeightFactor = m_cfg.zOriginWeightFactor;
0083   filterConfig.maxSeedsPerSpM = m_cfg.maxSeedsPerSpM;
0084   filterConfig.compatSeedLimit = m_cfg.compatSeedLimit;
0085   filterConfig.seedWeightIncrement = m_cfg.seedWeightIncrement;
0086   filterConfig.numSeedIncrement = m_cfg.numSeedIncrement;
0087   filterConfig.seedConfirmation = m_cfg.seedConfirmation;
0088   filterConfig.centralSeedConfirmationRange =
0089       m_cfg.centralSeedConfirmationRange;
0090   filterConfig.forwardSeedConfirmationRange =
0091       m_cfg.forwardSeedConfirmationRange;
0092   filterConfig.maxSeedsPerSpMConf = std::numeric_limits<std::size_t>::max();
0093   filterConfig.maxQualitySeedsPerSpMConf =
0094       std::numeric_limits<std::size_t>::max();
0095   filterConfig.useDeltaRinsteadOfTopRadius = m_cfg.useDeltaRinsteadOfTopRadius;
0096
0097   m_seedFilter = Acts::Experimental::BroadTripletSeedFilter(
0098       filterConfig, logger().cloneWithSuffix("Filter"));
0099 }
0100
0101 ProcessCode GridTripletSeedingAlgorithm::execute(
0102     const AlgorithmContext& ctx) const {
0103   const SimSpacePointContainer& spacePoints = m_inputSpacePoints(ctx);
0104
0105   Acts::Experimental::SpacePointContainer2 coreSpacePoints;
0106   coreSpacePoints.createExtraColumns(
0107       Acts::Experimental::SpacePointKnownExtraColumn::R |
0108       Acts::Experimental::SpacePointKnownExtraColumn::Phi |
0109       Acts::Experimental::SpacePointKnownExtraColumn::VarianceR |
0110       Acts::Experimental::SpacePointKnownExtraColumn::VarianceZ);
0111   coreSpacePoints.reserve(spacePoints.size());
0112   for (const auto& sp : spacePoints) {
0113     // check if the space point passes the selection
0114     if (m_spacePointSelector(sp)) {
0115       auto newSp = coreSpacePoints.createSpacePoint(
0116           std::array<Acts::SourceLink, 1>{Acts::SourceLink(&sp)}, sp.x(),
0117           sp.y(), sp.z());
0118       newSp.r() = sp.r();
0119       newSp.phi() = std::atan2(sp.y(), sp.x());
0120       newSp.varianceR() = sp.varianceR();
0121       newSp.varianceZ() = sp.varianceZ();
0122     }
0123   }
0124
0125   Acts::Experimental::CylindricalSpacePointGrid2 grid(
0126       m_gridConfig, logger().cloneWithSuffix("Grid"));
0127
0128   grid.fill(coreSpacePoints);
0129
0130   // Compute radius Range
0131   // we rely on the fact the grid is storing the proxies
0132   // with a sorting in the radius
0133   const Acts::Range1D<float> rRange = grid.computeRadiusRange(coreSpacePoints);
0134
0135   Acts::Experimental::BroadTripletSeedFinder::Options finderOptions;
0136   finderOptions.bFieldInZ = m_cfg.bFieldInZ;
0137
0138   Acts::Experimental::DoubletSeedFinder::Config bottomDoubletFinderConfig;
0139   bottomDoubletFinderConfig.candidateDirection = Acts::Direction::Backward();
0140   bottomDoubletFinderConfig.deltaRMin = std::isnan(m_cfg.deltaRMaxBottom)
0141                                             ? m_cfg.deltaRMin
0142                                             : m_cfg.deltaRMinBottom;
0143   bottomDoubletFinderConfig.deltaRMax = std::isnan(m_cfg.deltaRMaxBottom)
0144                                             ? m_cfg.deltaRMax
0145                                             : m_cfg.deltaRMaxBottom;
0146   bottomDoubletFinderConfig.deltaZMin = m_cfg.deltaZMin;
0147   bottomDoubletFinderConfig.deltaZMax = m_cfg.deltaZMax;
0148   bottomDoubletFinderConfig.impactMax = m_cfg.impactMax;
0149   bottomDoubletFinderConfig.interactionPointCut = m_cfg.interactionPointCut;
0150   bottomDoubletFinderConfig.collisionRegionMin = m_cfg.collisionRegionMin;
0151   bottomDoubletFinderConfig.collisionRegionMax = m_cfg.collisionRegionMax;
0152   bottomDoubletFinderConfig.cotThetaMax = m_cfg.cotThetaMax;
0153   bottomDoubletFinderConfig.minPt = m_cfg.minPt;
0154   bottomDoubletFinderConfig.helixCutTolerance = m_cfg.helixCutTolerance;
0155   if (m_cfg.useExtraCuts) {
0156     bottomDoubletFinderConfig.experimentCuts.connect<itkFastTrackingCuts>();
0157   }
0158   Acts::Experimental::DoubletSeedFinder bottomDoubletFinder(
0159       Acts::Experimental::DoubletSeedFinder::DerivedConfig(
0160           bottomDoubletFinderConfig, m_cfg.bFieldInZ));
0161
0162   Acts::Experimental::DoubletSeedFinder::Config topDoubletFinderConfig =
0163       bottomDoubletFinderConfig;
0164   topDoubletFinderConfig.candidateDirection = Acts::Direction::Forward();
0165   topDoubletFinderConfig.deltaRMin =
0166       std::isnan(m_cfg.deltaRMaxTop) ? m_cfg.deltaRMin : m_cfg.deltaRMinTop;
0167   topDoubletFinderConfig.deltaRMax =
0168       std::isnan(m_cfg.deltaRMaxTop) ? m_cfg.deltaRMax : m_cfg.deltaRMaxTop;
0169   Acts::Experimental::DoubletSeedFinder topDoubletFinder(
0170       Acts::Experimental::DoubletSeedFinder::DerivedConfig(
0171           topDoubletFinderConfig, m_cfg.bFieldInZ));
0172
0173   Acts::Experimental::BroadTripletSeedFinder::TripletCuts tripletCuts;
0174   tripletCuts.minPt = m_cfg.minPt;
0175   tripletCuts.sigmaScattering = m_cfg.sigmaScattering;
0176   tripletCuts.radLengthPerSeed = m_cfg.radLengthPerSeed;
0177   tripletCuts.maxPtScattering = m_cfg.maxPtScattering;
0178   tripletCuts.impactMax = m_cfg.impactMax;
0179   tripletCuts.helixCutTolerance = m_cfg.helixCutTolerance;
0180   tripletCuts.toleranceParam = m_cfg.toleranceParam;
0181   Acts::Experimental::BroadTripletSeedFinder::DerivedTripletCuts
0182       derivedTripletCuts(tripletCuts, m_cfg.bFieldInZ);
0183
0184   // variable middle SP radial region of interest
0185   Acts::Range1D<float> rMiddleSpRange = {
0186       std::floor(rRange.min() / 2) * 2 + m_cfg.deltaRMiddleMinSPRange,
0187       std::floor(rRange.max() / 2) * 2 - m_cfg.deltaRMiddleMaxSPRange};
0188
0189   // run the seeding
0190   Acts::Experimental::SeedContainer2 seeds;
0191   Acts::Experimental::BroadTripletSeedFinder::State state;
0192   static thread_local Acts::Experimental::BroadTripletSeedFinder::Cache cache;
0193
0194   std::vector<std::span<const Acts::SpacePointIndex2>> bottomSpGroups;
0195   std::span<const Acts::SpacePointIndex2> middleSps;
0196   std::vector<std::span<const Acts::SpacePointIndex2>> topSpGroups;
0197
0198   for (const auto [bottom, middle, top] : grid.binnedGroup()) {
0199     ACTS_VERBOSE("Process middle " << middle);
0200
0201     bottomSpGroups.clear();
0202     for (const auto b : bottom) {
0203       bottomSpGroups.push_back(grid.at(b));
0204     }
0205     middleSps = grid.at(middle);
0206     topSpGroups.clear();
0207     for (const auto t : top) {
0208       topSpGroups.push_back(grid.at(t));
0209     }
0210
0211     // we compute this here since all middle space point candidates belong to
0212     // the same z-bin
0213     auto firstMiddleSp = coreSpacePoints.at(middleSps.front());
0214     auto radiusRangeForMiddle = retrieveRadiusRangeForMiddle(
0215         Acts::Experimental::ConstSpacePointProxy2(firstMiddleSp),
0216         rMiddleSpRange);
0217     ACTS_VERBOSE("Validity range (radius) for the middle space point is ["
0218                  << radiusRangeForMiddle.first << ", "
0219                  << radiusRangeForMiddle.second << "]");
0220
0221     m_seedFinder->createSeedsFromSortedGroups(
0222         finderOptions, state, cache, bottomDoubletFinder, topDoubletFinder,
0223         derivedTripletCuts, *m_seedFilter, coreSpacePoints, bottomSpGroups,
0224         middleSps, topSpGroups, radiusRangeForMiddle, seeds);
0225   }
0226
0227   ACTS_DEBUG("Created " << seeds.size() << " track seeds from "
0228                         << spacePoints.size() << " space points");
0229
0230   // we have seeds of proxies
0231   // convert them to seed of external space points
0232   SimSeedContainer seedContainerForStorage;
0233   seedContainerForStorage.reserve(seeds.size());
0234   for (const auto& seed : seeds) {
0235     auto sps = seed.spacePointIndices();
0236     seedContainerForStorage.emplace_back(*coreSpacePoints.at(sps[0])
0237                                               .sourceLinks()[0]
0238                                               .get<const SimSpacePoint*>(),
0239                                          *coreSpacePoints.at(sps[1])
0240                                               .sourceLinks()[0]
0241                                               .get<const SimSpacePoint*>(),
0242                                          *coreSpacePoints.at(sps[2])
0243                                               .sourceLinks()[0]
0244                                               .get<const SimSpacePoint*>());
0245     seedContainerForStorage.back().setVertexZ(seed.vertexZ());
0246     seedContainerForStorage.back().setQuality(seed.quality());
0247   }
0248
0249   m_outputSeeds(ctx, std::move(seedContainerForStorage));
0250   return ProcessCode::SUCCESS;
0251 }
0252
0253 std::pair<float, float>
0254 GridTripletSeedingAlgorithm::retrieveRadiusRangeForMiddle(
0255     const Acts::Experimental::ConstSpacePointProxy2& spM,
0256     const Acts::Range1D<float>& rMiddleSpRange) const {
0257   if (m_cfg.useVariableMiddleSPRange) {
0258     return {rMiddleSpRange.min(), rMiddleSpRange.max()};
0259   }
0260   if (m_cfg.rRangeMiddleSP.empty()) {
0261     return {m_cfg.rMinMiddle, m_cfg.rMaxMiddle};
0262   }
0263
0264   // get zBin position of the middle SP
0265   auto pVal = std::ranges::lower_bound(m_cfg.zBinEdges, spM.z());
0266   std::size_t zBin = std::distance(m_cfg.zBinEdges.begin(), pVal);
0267   // protects against zM at the limit of zBinEdges
0268   zBin == 0 ? zBin : --zBin;
0269   return {m_cfg.rRangeMiddleSP[zBin][0], m_cfg.rRangeMiddleSP[zBin][1]};
0270 }
0271
0272 }  // namespace ActsExamples