Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-10-14 08:00:50

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/OrthogonalTripletSeedingAlgorithm.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/CylindricalSpacePointKDTree.hpp"
0018 #include "Acts/Seeding2/DoubletSeedFinder.hpp"
0019 #include "Acts/Seeding2/TripletSeedFinder.hpp"
0020 #include "Acts/Utilities/Delegate.hpp"
0021 #include "ActsExamples/EventData/SimSeed.hpp"
0022 #include "ActsExamples/EventData/SimSpacePoint.hpp"
0023 
0024 #include <cmath>
0025 #include <csignal>
0026 #include <cstddef>
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 OrthogonalTripletSeedingAlgorithm::OrthogonalTripletSeedingAlgorithm(
0066     const Config &cfg, Acts::Logging::Level lvl)
0067     : IAlgorithm("OrthogonalTripletSeedingAlgorithm", lvl), m_cfg(cfg) {
0068   m_inputSpacePoints.initialize(m_cfg.inputSpacePoints);
0069   m_outputSeeds.initialize(m_cfg.outputSeeds);
0070 
0071   if (m_cfg.useExtraCuts) {
0072     // This function will be applied to select space points during grid filling
0073     m_spacePointSelector.connect<itkFastTrackingSPselect>();
0074   }
0075 
0076   m_filterConfig.deltaInvHelixDiameter = m_cfg.deltaInvHelixDiameter;
0077   m_filterConfig.deltaRMin = m_cfg.deltaRMin;
0078   m_filterConfig.compatSeedWeight = m_cfg.compatSeedWeight;
0079   m_filterConfig.impactWeightFactor = m_cfg.impactWeightFactor;
0080   m_filterConfig.zOriginWeightFactor = m_cfg.zOriginWeightFactor;
0081   m_filterConfig.maxSeedsPerSpM = m_cfg.maxSeedsPerSpM;
0082   m_filterConfig.compatSeedLimit = m_cfg.compatSeedLimit;
0083   m_filterConfig.seedWeightIncrement = m_cfg.seedWeightIncrement;
0084   m_filterConfig.numSeedIncrement = m_cfg.numSeedIncrement;
0085   m_filterConfig.seedConfirmation = m_cfg.seedConfirmation;
0086   m_filterConfig.centralSeedConfirmationRange =
0087       m_cfg.centralSeedConfirmationRange;
0088   m_filterConfig.forwardSeedConfirmationRange =
0089       m_cfg.forwardSeedConfirmationRange;
0090   m_filterConfig.maxSeedsPerSpMConf = m_cfg.maxSeedsPerSpMConf;
0091   m_filterConfig.maxQualitySeedsPerSpMConf = m_cfg.maxQualitySeedsPerSpMConf;
0092   m_filterConfig.useDeltaRinsteadOfTopRadius =
0093       m_cfg.useDeltaRinsteadOfTopRadius;
0094 
0095   m_filterLogger = logger().cloneWithSuffix("Filter");
0096 
0097   m_seedFinder = Acts::TripletSeeder(logger().cloneWithSuffix("Finder"));
0098 }
0099 
0100 ProcessCode OrthogonalTripletSeedingAlgorithm::execute(
0101     const AlgorithmContext &ctx) const {
0102   const SimSpacePointContainer &spacePoints = m_inputSpacePoints(ctx);
0103 
0104   Acts::SpacePointContainer2 coreSpacePoints(
0105       Acts::SpacePointColumns::SourceLinks | Acts::SpacePointColumns::XY |
0106       Acts::SpacePointColumns::ZR | Acts::SpacePointColumns::Phi |
0107       Acts::SpacePointColumns::VarianceZ | Acts::SpacePointColumns::VarianceR);
0108   coreSpacePoints.reserve(spacePoints.size());
0109 
0110   Acts::Experimental::CylindricalSpacePointKDTreeBuilder kdTreeBuilder;
0111   kdTreeBuilder.reserve(spacePoints.size());
0112 
0113   Acts::Extent rRangeSPExtent;
0114 
0115   for (std::size_t i = 0; i < spacePoints.size(); ++i) {
0116     const auto &sp = spacePoints[i];
0117 
0118     // check if the space point passes the selection
0119     if (m_spacePointSelector.connected() && !m_spacePointSelector(sp)) {
0120       continue;
0121     }
0122 
0123     Acts::SpacePointIndex2 newSpIndex = coreSpacePoints.size();
0124     auto newSp = coreSpacePoints.createSpacePoint();
0125     newSp.assignSourceLinks(
0126         std::array<Acts::SourceLink, 1>{Acts::SourceLink(&sp)});
0127     newSp.xy() = std::array<float, 2>{static_cast<float>(sp.x()),
0128                                       static_cast<float>(sp.y())};
0129     newSp.zr() = std::array<float, 2>{static_cast<float>(sp.z()),
0130                                       static_cast<float>(sp.r())};
0131     newSp.phi() = static_cast<float>(std::atan2(sp.y(), sp.x()));
0132     newSp.varianceZ() = static_cast<float>(sp.varianceZ());
0133     newSp.varianceR() = static_cast<float>(sp.varianceR());
0134 
0135     kdTreeBuilder.insert(newSpIndex, newSp.phi(), newSp.zr()[1], newSp.zr()[0]);
0136 
0137     rRangeSPExtent.extend({newSp.xy()[0], newSp.xy()[1], newSp.zr()[0]});
0138   }
0139 
0140   ACTS_VERBOSE("Created k-d tree populated with " << kdTreeBuilder.size()
0141                                                   << " space points");
0142 
0143   Acts::Experimental::CylindricalSpacePointKDTree kdTree =
0144       kdTreeBuilder.build();
0145 
0146   Acts::Experimental::CylindricalSpacePointKDTree::Options lhOptions;
0147   lhOptions.rMax = m_cfg.rMax;
0148   lhOptions.zMin = m_cfg.zMin;
0149   lhOptions.zMax = m_cfg.zMax;
0150   lhOptions.phiMin = m_cfg.phiMin;
0151   lhOptions.phiMax = m_cfg.phiMax;
0152   lhOptions.deltaRMin = std::isnan(m_cfg.deltaRMinBottom)
0153                             ? m_cfg.deltaRMin
0154                             : m_cfg.deltaRMinBottom;
0155   lhOptions.deltaRMax = std::isnan(m_cfg.deltaRMaxBottom)
0156                             ? m_cfg.deltaRMax
0157                             : m_cfg.deltaRMaxBottom;
0158   lhOptions.collisionRegionMin = m_cfg.collisionRegionMin;
0159   lhOptions.collisionRegionMax = m_cfg.collisionRegionMax;
0160   lhOptions.cotThetaMax = m_cfg.cotThetaMax;
0161   lhOptions.deltaPhiMax = m_cfg.deltaPhiMax;
0162   Acts::Experimental::CylindricalSpacePointKDTree::Options hlOptions =
0163       lhOptions;
0164   hlOptions.deltaRMin =
0165       std::isnan(m_cfg.deltaRMinTop) ? m_cfg.deltaRMin : m_cfg.deltaRMinTop;
0166   hlOptions.deltaRMax =
0167       std::isnan(m_cfg.deltaRMaxTop) ? m_cfg.deltaRMax : m_cfg.deltaRMaxTop;
0168 
0169   Acts::DoubletSeedFinder::Config bottomDoubletFinderConfig;
0170   bottomDoubletFinderConfig.spacePointsSortedByRadius = false;
0171   bottomDoubletFinderConfig.candidateDirection = Acts::Direction::Backward();
0172   bottomDoubletFinderConfig.deltaRMin = std::isnan(m_cfg.deltaRMaxBottom)
0173                                             ? m_cfg.deltaRMin
0174                                             : m_cfg.deltaRMinBottom;
0175   bottomDoubletFinderConfig.deltaRMax = std::isnan(m_cfg.deltaRMaxBottom)
0176                                             ? m_cfg.deltaRMax
0177                                             : m_cfg.deltaRMaxBottom;
0178   bottomDoubletFinderConfig.deltaZMin = m_cfg.deltaZMin;
0179   bottomDoubletFinderConfig.deltaZMax = m_cfg.deltaZMax;
0180   bottomDoubletFinderConfig.impactMax = m_cfg.impactMax;
0181   bottomDoubletFinderConfig.interactionPointCut = m_cfg.interactionPointCut;
0182   bottomDoubletFinderConfig.collisionRegionMin = m_cfg.collisionRegionMin;
0183   bottomDoubletFinderConfig.collisionRegionMax = m_cfg.collisionRegionMax;
0184   bottomDoubletFinderConfig.cotThetaMax = m_cfg.cotThetaMax;
0185   bottomDoubletFinderConfig.minPt = m_cfg.minPt;
0186   bottomDoubletFinderConfig.helixCutTolerance = m_cfg.helixCutTolerance;
0187   if (m_cfg.useExtraCuts) {
0188     bottomDoubletFinderConfig.experimentCuts.connect<itkFastTrackingCuts>();
0189   }
0190   auto bottomDoubletFinder =
0191       Acts::DoubletSeedFinder::create(Acts::DoubletSeedFinder::DerivedConfig(
0192           bottomDoubletFinderConfig, m_cfg.bFieldInZ));
0193 
0194   Acts::DoubletSeedFinder::Config topDoubletFinderConfig =
0195       bottomDoubletFinderConfig;
0196   topDoubletFinderConfig.candidateDirection = Acts::Direction::Forward();
0197   topDoubletFinderConfig.deltaRMin =
0198       std::isnan(m_cfg.deltaRMaxTop) ? m_cfg.deltaRMin : m_cfg.deltaRMinTop;
0199   topDoubletFinderConfig.deltaRMax =
0200       std::isnan(m_cfg.deltaRMaxTop) ? m_cfg.deltaRMax : m_cfg.deltaRMaxTop;
0201   auto topDoubletFinder =
0202       Acts::DoubletSeedFinder::create(Acts::DoubletSeedFinder::DerivedConfig(
0203           topDoubletFinderConfig, m_cfg.bFieldInZ));
0204 
0205   Acts::TripletSeedFinder::Config tripletFinderConfig;
0206   tripletFinderConfig.useStripInfo = false;
0207   tripletFinderConfig.sortedByCotTheta = true;
0208   tripletFinderConfig.minPt = m_cfg.minPt;
0209   tripletFinderConfig.sigmaScattering = m_cfg.sigmaScattering;
0210   tripletFinderConfig.radLengthPerSeed = m_cfg.radLengthPerSeed;
0211   tripletFinderConfig.impactMax = m_cfg.impactMax;
0212   tripletFinderConfig.helixCutTolerance = m_cfg.helixCutTolerance;
0213   tripletFinderConfig.toleranceParam = m_cfg.toleranceParam;
0214   auto tripletFinder =
0215       Acts::TripletSeedFinder::create(Acts::TripletSeedFinder::DerivedConfig(
0216           tripletFinderConfig, m_cfg.bFieldInZ));
0217 
0218   // variable middle SP radial region of interest
0219   const Acts::Range1D<float> rMiddleSPRange(
0220       std::floor(rRangeSPExtent.min(Acts::AxisDirection::AxisR) / 2) * 2 +
0221           m_cfg.deltaRMiddleMinSPRange,
0222       std::floor(rRangeSPExtent.max(Acts::AxisDirection::AxisR) / 2) * 2 -
0223           m_cfg.deltaRMiddleMaxSPRange);
0224 
0225   // run the seeding
0226   Acts::SeedContainer2 seeds;
0227   Acts::BroadTripletSeedFilter::State filterState;
0228   Acts::BroadTripletSeedFilter::Cache filterCache;
0229   Acts::BroadTripletSeedFilter seedFilter(m_filterConfig, filterState,
0230                                           filterCache, *m_filterLogger);
0231 
0232   static thread_local Acts::TripletSeeder::Cache cache;
0233   static thread_local Acts::Experimental::CylindricalSpacePointKDTree::
0234       Candidates candidates;
0235 
0236   // Run the seeding algorithm by iterating over all the points in the tree
0237   // and seeing what happens if we take them to be our middle spacepoint.
0238   for (const auto &middle : kdTree) {
0239     ACTS_VERBOSE("Process middle " << middle.second);
0240 
0241     const auto spM = coreSpacePoints.at(middle.second).asConst();
0242 
0243     // Cut: Ensure that the middle spacepoint lies within a valid r-region for
0244     // middle points.
0245     const float rM = spM.zr()[1];
0246     if (m_cfg.useVariableMiddleSPRange) {
0247       if (rM < rMiddleSPRange.min() || rM > rMiddleSPRange.max()) {
0248         continue;
0249       }
0250     } else {
0251       if (rM > m_cfg.rMaxMiddle || rM < m_cfg.rMinMiddle) {
0252         continue;
0253       }
0254     }
0255 
0256     // remove all middle SPs outside phi and z region of interest
0257     const float zM = spM.zr()[0];
0258     if (zM < m_cfg.zOutermostLayers.first ||
0259         zM > m_cfg.zOutermostLayers.second) {
0260       continue;
0261     }
0262     if (const float phiM = spM.phi();
0263         phiM > m_cfg.phiMax || phiM < m_cfg.phiMin) {
0264       continue;
0265     }
0266 
0267     std::size_t nTopSeedConf = 0;
0268     if (m_cfg.seedConfirmation) {
0269       // check if middle SP is in the central or forward region
0270       Acts::SeedConfirmationRangeConfig seedConfRange =
0271           (zM > m_cfg.centralSeedConfirmationRange.zMaxSeedConf ||
0272            zM < m_cfg.centralSeedConfirmationRange.zMinSeedConf)
0273               ? m_cfg.forwardSeedConfirmationRange
0274               : m_cfg.centralSeedConfirmationRange;
0275       // set the minimum number of top SP depending on whether the middle SP is
0276       // in the central or forward region
0277       nTopSeedConf = rM > seedConfRange.rMaxSeedConf
0278                          ? seedConfRange.nTopForLargeR
0279                          : seedConfRange.nTopForSmallR;
0280     }
0281 
0282     candidates.clear();
0283     kdTree.validTuples(lhOptions, hlOptions, spM, nTopSeedConf, candidates);
0284 
0285     Acts::SpacePointContainer2::ConstSubset bottomSps =
0286         coreSpacePoints.subset(candidates.bottom_lh_v).asConst();
0287     Acts::SpacePointContainer2::ConstSubset topSps =
0288         coreSpacePoints.subset(candidates.top_lh_v).asConst();
0289     m_seedFinder->createSeedsFromGroup(
0290         cache, *bottomDoubletFinder, *topDoubletFinder, *tripletFinder,
0291         seedFilter, coreSpacePoints, bottomSps, spM, topSps, seeds);
0292 
0293     bottomSps = coreSpacePoints.subset(candidates.bottom_hl_v).asConst();
0294     topSps = coreSpacePoints.subset(candidates.top_hl_v).asConst();
0295     m_seedFinder->createSeedsFromGroup(
0296         cache, *bottomDoubletFinder, *topDoubletFinder, *tripletFinder,
0297         seedFilter, coreSpacePoints, bottomSps, spM, topSps, seeds);
0298   }
0299 
0300   ACTS_DEBUG("Created " << seeds.size() << " track seeds from "
0301                         << spacePoints.size() << " space points");
0302 
0303   // we have seeds of proxies
0304   // convert them to seed of external space points
0305   SimSeedContainer seedContainerForStorage;
0306   seedContainerForStorage.reserve(seeds.size());
0307   for (const auto &seed : seeds) {
0308     auto sps = seed.spacePointIndices();
0309     seedContainerForStorage.emplace_back(*coreSpacePoints.at(sps[0])
0310                                               .sourceLinks()[0]
0311                                               .get<const SimSpacePoint *>(),
0312                                          *coreSpacePoints.at(sps[1])
0313                                               .sourceLinks()[0]
0314                                               .get<const SimSpacePoint *>(),
0315                                          *coreSpacePoints.at(sps[2])
0316                                               .sourceLinks()[0]
0317                                               .get<const SimSpacePoint *>());
0318     seedContainerForStorage.back().setVertexZ(seed.vertexZ());
0319     seedContainerForStorage.back().setQuality(seed.quality());
0320   }
0321 
0322   m_outputSeeds(ctx, std::move(seedContainerForStorage));
0323   return ProcessCode::SUCCESS;
0324 }
0325 
0326 }  // namespace ActsExamples