File indexing completed on 2025-10-14 08:00:50
0001
0002
0003
0004
0005
0006
0007
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 & ,
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
0048 float r = sp.r();
0049 float zabs = std::abs(sp.z());
0050 if (zabs > 200. && r < 45.) {
0051 return false;
0052 }
0053
0054
0055
0056 float cotTheta = 27.2899;
0057 if ((zabs - 150.) > cotTheta * r) {
0058 return false;
0059 }
0060 return true;
0061 }
0062
0063 }
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
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
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
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
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
0237
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
0244
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
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
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
0276
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
0304
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 }