TrackFinding/src/SeedingAlgorithm.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2024 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/TrackFinding/SeedingAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/EventData/SpacePointData.hpp"
0013 #include "Acts/Geometry/Extent.hpp"
0014 #include "Acts/Seeding/BinnedGroup.hpp"
0015 #include "Acts/Seeding/InternalSpacePoint.hpp"
0016 #include "Acts/Seeding/SeedFilter.hpp"
0017 #include "Acts/Utilities/BinningType.hpp"
0018 #include "Acts/Utilities/Delegate.hpp"
0019 #include "Acts/Utilities/Grid.hpp"
0020 #include "Acts/Utilities/GridBinFinder.hpp"
0021 #include "Acts/Utilities/Helpers.hpp"
0022 #include "ActsExamples/EventData/SimSeed.hpp"
0023
0024 #include <cmath>
0025 #include <csignal>
0026 #include <cstddef>
0027 #include <iterator>
0028 #include <limits>
0029 #include <ostream>
0030 #include <stdexcept>
0031
0032 namespace ActsExamples {
0033 struct AlgorithmContext;
0034 }  // namespace ActsExamples
0035
0036 ActsExamples::SeedingAlgorithm::SeedingAlgorithm(
0037     ActsExamples::SeedingAlgorithm::Config cfg, Acts::Logging::Level lvl)
0038     : ActsExamples::IAlgorithm("SeedingAlgorithm", lvl), m_cfg(std::move(cfg)) {
0039   // Seed Finder config requires Seed Filter object before conversion to
0040   // internal units
0041   m_cfg.seedFilterConfig = m_cfg.seedFilterConfig.toInternalUnits();
0042   m_cfg.seedFinderConfig.seedFilter =
0043       std::make_unique<Acts::SeedFilter<SimSpacePoint>>(m_cfg.seedFilterConfig);
0044
0045   m_cfg.seedFinderConfig =
0046       m_cfg.seedFinderConfig.toInternalUnits().calculateDerivedQuantities();
0047   m_cfg.seedFinderOptions =
0048       m_cfg.seedFinderOptions.toInternalUnits().calculateDerivedQuantities(
0049           m_cfg.seedFinderConfig);
0050   m_cfg.gridConfig = m_cfg.gridConfig.toInternalUnits();
0051   m_cfg.gridOptions = m_cfg.gridOptions.toInternalUnits();
0052   if (m_cfg.inputSpacePoints.empty()) {
0053     throw std::invalid_argument("Missing space point input collections");
0054   }
0055
0056   for (const auto& spName : m_cfg.inputSpacePoints) {
0057     if (spName.empty()) {
0058       throw std::invalid_argument("Invalid space point input collection");
0059     }
0060
0061     auto& handle = m_inputSpacePoints.emplace_back(
0062         std::make_unique<ReadDataHandle<SimSpacePointContainer>>(
0063             this,
0064             "InputSpacePoints#" + std::to_string(m_inputSpacePoints.size())));
0065     handle->initialize(spName);
0066   }
0067   if (m_cfg.outputSeeds.empty()) {
0068     throw std::invalid_argument("Missing seeds output collection");
0069   }
0070
0071   m_outputSeeds.initialize(m_cfg.outputSeeds);
0072
0073   if (m_cfg.gridConfig.rMax != m_cfg.seedFinderConfig.rMax &&
0074       m_cfg.allowSeparateRMax == false) {
0075     throw std::invalid_argument(
0076         "Inconsistent config rMax: using different values in gridConfig and "
0077         "seedFinderConfig. If values are intentional set allowSeparateRMax to "
0078         "true");
0079   }
0080
0081   if (m_cfg.seedFilterConfig.deltaRMin != m_cfg.seedFinderConfig.deltaRMin) {
0082     throw std::invalid_argument("Inconsistent config deltaRMin");
0083   }
0084
0085   if (m_cfg.gridConfig.deltaRMax != m_cfg.seedFinderConfig.deltaRMax) {
0086     throw std::invalid_argument("Inconsistent config deltaRMax");
0087   }
0088
0089   static_assert(
0090       std::numeric_limits<
0091           decltype(m_cfg.seedFinderConfig.deltaRMaxTopSP)>::has_quiet_NaN,
0092       "Value of deltaRMaxTopSP must support NaN values");
0093
0094   static_assert(
0095       std::numeric_limits<
0096           decltype(m_cfg.seedFinderConfig.deltaRMinTopSP)>::has_quiet_NaN,
0097       "Value of deltaRMinTopSP must support NaN values");
0098
0099   static_assert(
0100       std::numeric_limits<
0101           decltype(m_cfg.seedFinderConfig.deltaRMaxBottomSP)>::has_quiet_NaN,
0102       "Value of deltaRMaxBottomSP must support NaN values");
0103
0104   static_assert(
0105       std::numeric_limits<
0106           decltype(m_cfg.seedFinderConfig.deltaRMinBottomSP)>::has_quiet_NaN,
0107       "Value of deltaRMinBottomSP must support NaN values");
0108
0109   if (std::isnan(m_cfg.seedFinderConfig.deltaRMaxTopSP)) {
0110     m_cfg.seedFinderConfig.deltaRMaxTopSP = m_cfg.seedFinderConfig.deltaRMax;
0111   }
0112
0113   if (std::isnan(m_cfg.seedFinderConfig.deltaRMinTopSP)) {
0114     m_cfg.seedFinderConfig.deltaRMinTopSP = m_cfg.seedFinderConfig.deltaRMin;
0115   }
0116
0117   if (std::isnan(m_cfg.seedFinderConfig.deltaRMaxBottomSP)) {
0118     m_cfg.seedFinderConfig.deltaRMaxBottomSP = m_cfg.seedFinderConfig.deltaRMax;
0119   }
0120
0121   if (std::isnan(m_cfg.seedFinderConfig.deltaRMinBottomSP)) {
0122     m_cfg.seedFinderConfig.deltaRMinBottomSP = m_cfg.seedFinderConfig.deltaRMin;
0123   }
0124
0125   if (m_cfg.gridConfig.zMin != m_cfg.seedFinderConfig.zMin) {
0126     throw std::invalid_argument("Inconsistent config zMin");
0127   }
0128
0129   if (m_cfg.gridConfig.zMax != m_cfg.seedFinderConfig.zMax) {
0130     throw std::invalid_argument("Inconsistent config zMax");
0131   }
0132
0133   if (m_cfg.seedFilterConfig.maxSeedsPerSpM !=
0134       m_cfg.seedFinderConfig.maxSeedsPerSpM) {
0135     throw std::invalid_argument("Inconsistent config maxSeedsPerSpM");
0136   }
0137
0138   if (m_cfg.gridConfig.cotThetaMax != m_cfg.seedFinderConfig.cotThetaMax) {
0139     throw std::invalid_argument("Inconsistent config cotThetaMax");
0140   }
0141
0142   if (m_cfg.gridConfig.minPt != m_cfg.seedFinderConfig.minPt) {
0143     throw std::invalid_argument("Inconsistent config minPt");
0144   }
0145
0146   if (m_cfg.gridOptions.bFieldInZ != m_cfg.seedFinderOptions.bFieldInZ) {
0147     throw std::invalid_argument("Inconsistent config bFieldInZ");
0148   }
0149
0150   if (m_cfg.gridConfig.zBinEdges.size() - 1 != m_cfg.zBinNeighborsTop.size() &&
0151       m_cfg.zBinNeighborsTop.empty() == false) {
0152     throw std::invalid_argument("Inconsistent config zBinNeighborsTop");
0153   }
0154
0155   if (m_cfg.gridConfig.zBinEdges.size() - 1 !=
0156           m_cfg.zBinNeighborsBottom.size() &&
0157       m_cfg.zBinNeighborsBottom.empty() == false) {
0158     throw std::invalid_argument("Inconsistent config zBinNeighborsBottom");
0159   }
0160
0161   if (!m_cfg.seedFinderConfig.zBinsCustomLooping.empty()) {
0162     // check if zBinsCustomLooping contains numbers from 1 to the total number
0163     // of bin in zBinEdges
0164     for (std::size_t i = 1; i != m_cfg.gridConfig.zBinEdges.size(); i++) {
0165       if (std::find(m_cfg.seedFinderConfig.zBinsCustomLooping.begin(),
0166                     m_cfg.seedFinderConfig.zBinsCustomLooping.end(),
0167                     i) == m_cfg.seedFinderConfig.zBinsCustomLooping.end()) {
0168         throw std::invalid_argument(
0169             "Inconsistent config zBinsCustomLooping does not contain the same "
0170             "bins as zBinEdges");
0171       }
0172     }
0173   }
0174
0175   if (m_cfg.seedFinderConfig.useDetailedDoubleMeasurementInfo) {
0176     m_cfg.seedFinderConfig.getTopHalfStripLength.connect(
0177         [](const void*, const SimSpacePoint& sp) -> float {
0178           return sp.topHalfStripLength();
0179         });
0180
0181     m_cfg.seedFinderConfig.getBottomHalfStripLength.connect(
0182         [](const void*, const SimSpacePoint& sp) -> float {
0183           return sp.bottomHalfStripLength();
0184         });
0185
0186     m_cfg.seedFinderConfig.getTopStripDirection.connect(
0187         [](const void*, const SimSpacePoint& sp) -> Acts::Vector3 {
0188           return sp.topStripDirection();
0189         });
0190
0191     m_cfg.seedFinderConfig.getBottomStripDirection.connect(
0192         [](const void*, const SimSpacePoint& sp) -> Acts::Vector3 {
0193           return sp.bottomStripDirection();
0194         });
0195
0196     m_cfg.seedFinderConfig.getStripCenterDistance.connect(
0197         [](const void*, const SimSpacePoint& sp) -> Acts::Vector3 {
0198           return sp.stripCenterDistance();
0199         });
0200
0201     m_cfg.seedFinderConfig.getTopStripCenterPosition.connect(
0202         [](const void*, const SimSpacePoint& sp) -> Acts::Vector3 {
0203           return sp.topStripCenterPosition();
0204         });
0205   }
0206
0207   m_bottomBinFinder = std::make_unique<const Acts::GridBinFinder<2ul>>(
0208       m_cfg.numPhiNeighbors, m_cfg.zBinNeighborsBottom);
0209   m_topBinFinder = std::make_unique<const Acts::GridBinFinder<2ul>>(
0210       m_cfg.numPhiNeighbors, m_cfg.zBinNeighborsTop);
0211
0212   m_cfg.seedFinderConfig.seedFilter =
0213       std::make_unique<Acts::SeedFilter<SimSpacePoint>>(m_cfg.seedFilterConfig);
0214   m_seedFinder =
0215       Acts::SeedFinder<SimSpacePoint,
0216                        Acts::CylindricalSpacePointGrid<SimSpacePoint>>(
0217           m_cfg.seedFinderConfig);
0218 }
0219
0220 ActsExamples::ProcessCode ActsExamples::SeedingAlgorithm::execute(
0221     const AlgorithmContext& ctx) const {
0222   // construct the combined input container of space point pointers from all
0223   // configured input sources.
0224   // pre-compute the total size required so we only need to allocate once
0225   std::size_t nSpacePoints = 0;
0226   for (const auto& isp : m_inputSpacePoints) {
0227     nSpacePoints += (*isp)(ctx).size();
0228   }
0229
0230   std::vector<const SimSpacePoint*> spacePointPtrs;
0231   spacePointPtrs.reserve(nSpacePoints);
0232   for (const auto& isp : m_inputSpacePoints) {
0233     for (const auto& spacePoint : (*isp)(ctx)) {
0234       // since the event store owns the space
0235       // points, their pointers should be stable and
0236       // we do not need to create local copies.
0237       spacePointPtrs.push_back(&spacePoint);
0238     }
0239   }
0240
0241   // construct the seeding tools
0242   // covariance tool, extracts covariances per spacepoint as required
0243   auto extractGlobalQuantities = [=](const SimSpacePoint& sp, float, float,
0244                                      float) {
0245     Acts::Vector3 position{sp.x(), sp.y(), sp.z()};
0246     Acts::Vector2 covariance{sp.varianceR(), sp.varianceZ()};
0247     return std::make_tuple(position, covariance, sp.t());
0248   };
0249
0250   // extent used to store r range for middle spacepoint
0251   Acts::Extent rRangeSPExtent;
0252
0253   Acts::CylindricalSpacePointGrid<SimSpacePoint> grid =
0254       Acts::CylindricalSpacePointGridCreator::createGrid<SimSpacePoint>(
0255           m_cfg.gridConfig, m_cfg.gridOptions);
0256   Acts::CylindricalSpacePointGridCreator::fillGrid(
0257       m_cfg.seedFinderConfig, m_cfg.seedFinderOptions, grid,
0258       spacePointPtrs.begin(), spacePointPtrs.end(), extractGlobalQuantities,
0259       rRangeSPExtent);
0260
0261   std::array<std::vector<std::size_t>, 2ul> navigation;
0262   navigation[1ul] = m_cfg.seedFinderConfig.zBinsCustomLooping;
0263
0264   auto spacePointsGrouping = Acts::CylindricalBinnedGroup<SimSpacePoint>(
0265       std::move(grid), *m_bottomBinFinder, *m_topBinFinder,
0266       std::move(navigation));
0267
0268   // safely clamp double to float
0269   float up = Acts::clampValue<float>(
0270       std::floor(rRangeSPExtent.max(Acts::binR) / 2) * 2);
0271
0272   /// variable middle SP radial region of interest
0273   const Acts::Range1D<float> rMiddleSPRange(
0274       std::floor(rRangeSPExtent.min(Acts::binR) / 2) * 2 +
0275           m_cfg.seedFinderConfig.deltaRMiddleMinSPRange,
0276       up - m_cfg.seedFinderConfig.deltaRMiddleMaxSPRange);
0277
0278   // run the seeding
0279   static thread_local SimSeedContainer seeds;
0280   seeds.clear();
0281   static thread_local decltype(m_seedFinder)::SeedingState state;
0282   state.spacePointData.resize(
0283       spacePointPtrs.size(),
0284       m_cfg.seedFinderConfig.useDetailedDoubleMeasurementInfo);
0285
0286   if (m_cfg.seedFinderConfig.useDetailedDoubleMeasurementInfo) {
0287     for (std::size_t grid_glob_bin(0);
0288          grid_glob_bin < spacePointsGrouping.grid().size(); ++grid_glob_bin) {
0289       const auto& collection = spacePointsGrouping.grid().at(grid_glob_bin);
0290       for (const auto& sp : collection) {
0291         std::size_t index = sp->index();
0292
0293         const float topHalfStripLength =
0294             m_cfg.seedFinderConfig.getTopHalfStripLength(sp->sp());
0295         const float bottomHalfStripLength =
0296             m_cfg.seedFinderConfig.getBottomHalfStripLength(sp->sp());
0297         const Acts::Vector3 topStripDirection =
0298             m_cfg.seedFinderConfig.getTopStripDirection(sp->sp());
0299         const Acts::Vector3 bottomStripDirection =
0300             m_cfg.seedFinderConfig.getBottomStripDirection(sp->sp());
0301
0302         state.spacePointData.setTopStripVector(
0303             index, topHalfStripLength * topStripDirection);
0304         state.spacePointData.setBottomStripVector(
0305             index, bottomHalfStripLength * bottomStripDirection);
0306         state.spacePointData.setStripCenterDistance(
0307             index, m_cfg.seedFinderConfig.getStripCenterDistance(sp->sp()));
0308         state.spacePointData.setTopStripCenterPosition(
0309             index, m_cfg.seedFinderConfig.getTopStripCenterPosition(sp->sp()));
0310       }
0311     }
0312   }
0313
0314   for (const auto [bottom, middle, top] : spacePointsGrouping) {
0315     m_seedFinder.createSeedsForGroup(
0316         m_cfg.seedFinderOptions, state, spacePointsGrouping.grid(),
0317         std::back_inserter(seeds), bottom, middle, top, rMiddleSPRange);
0318   }
0319
0320   ACTS_DEBUG("Created " << seeds.size() << " track seeds from "
0321                         << spacePointPtrs.size() << " space points");
0322
0323   m_outputSeeds(ctx, SimSeedContainer{seeds});
0324   return ActsExamples::ProcessCode::SUCCESS;
0325 }