TrackFinding/src/SeedingAlgorithm.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/SeedingAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/EventData/Seed.hpp"
0013 #include "Acts/Seeding/BinnedGroup.hpp"
0014 #include "Acts/Seeding/SeedFilter.hpp"
0015 #include "Acts/Utilities/Delegate.hpp"
0016 #include "Acts/Utilities/GridBinFinder.hpp"
0017 #include "ActsExamples/EventData/SimSeed.hpp"
0018
0019 #include <cmath>
0020 #include <csignal>
0021 #include <cstddef>
0022 #include <limits>
0023 #include <ostream>
0024 #include <stdexcept>
0025
0026 using namespace Acts::HashedStringLiteral;
0027
0028 namespace ActsExamples {
0029 struct AlgorithmContext;
0030 }  // namespace ActsExamples
0031
0032 ActsExamples::SeedingAlgorithm::SeedingAlgorithm(
0033     ActsExamples::SeedingAlgorithm::Config cfg, Acts::Logging::Level lvl)
0034     : ActsExamples::IAlgorithm("SeedingAlgorithm", lvl), m_cfg(std::move(cfg)) {
0035   using SpacePointProxy_type = typename Acts::SpacePointContainer<
0036       ActsExamples::SpacePointContainer<std::vector<const SimSpacePoint*>>,
0037       Acts::detail::RefHolder>::SpacePointProxyType;
0038
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<SpacePointProxy_type>>(
0044           m_cfg.seedFilterConfig, logger().cloneWithSuffix("SeedFilter"));
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 that the bins required in the custom bin looping
0163     // are contained in the bins defined by the total number of edges
0164
0165     for (std::size_t i : m_cfg.seedFinderConfig.zBinsCustomLooping) {
0166       if (i >= m_cfg.gridConfig.zBinEdges.size()) {
0167         throw std::invalid_argument(
0168             "Inconsistent config zBinsCustomLooping does not contain a subset "
0169             "of bins defined by zBinEdges");
0170       }
0171     }
0172   }
0173
0174   if (m_cfg.useExtraCuts) {
0175     // This function will be applied to select space points during grid filling
0176     m_cfg.seedFinderConfig.spacePointSelector
0177         .connect<itkFastTrackingSPselect>();
0178
0179     // This function will be applied to the doublet compatibility selection
0180     m_cfg.seedFinderConfig.experimentCuts.connect<itkFastTrackingCuts>();
0181   }
0182
0183   using SpacePointProxy_type = typename Acts::SpacePointContainer<
0184       ActsExamples::SpacePointContainer<std::vector<const SimSpacePoint*>>,
0185       Acts::detail::RefHolder>::SpacePointProxyType;
0186
0187   m_bottomBinFinder = std::make_unique<const Acts::GridBinFinder<3ul>>(
0188       m_cfg.numPhiNeighbors, cfg.zBinNeighborsBottom, 0);
0189   m_topBinFinder = std::make_unique<const Acts::GridBinFinder<3ul>>(
0190       m_cfg.numPhiNeighbors, m_cfg.zBinNeighborsTop, 0);
0191
0192   m_seedFinder =
0193       Acts::SeedFinder<SpacePointProxy_type,
0194                        Acts::CylindricalSpacePointGrid<SpacePointProxy_type>>(
0195           m_cfg.seedFinderConfig, logger().cloneWithSuffix("SeedFinder"));
0196 }
0197
0198 ActsExamples::ProcessCode ActsExamples::SeedingAlgorithm::execute(
0199     const AlgorithmContext& ctx) const {
0200   // construct the combined input container of space point pointers from all
0201   // configured input sources.
0202   // pre-compute the total size required so we only need to allocate once
0203   std::size_t nSpacePoints = 0;
0204   for (const auto& isp : m_inputSpacePoints) {
0205     nSpacePoints += (*isp)(ctx).size();
0206   }
0207
0208   std::vector<const SimSpacePoint*> spacePointPtrs;
0209   spacePointPtrs.reserve(nSpacePoints);
0210   for (const auto& isp : m_inputSpacePoints) {
0211     for (const auto& spacePoint : (*isp)(ctx)) {
0212       // since the event store owns the space
0213       // points, their pointers should be stable and
0214       // we do not need to create local copies.
0215       spacePointPtrs.push_back(&spacePoint);
0216     }
0217   }
0218
0219   // Config
0220   Acts::SpacePointContainerConfig spConfig;
0221   spConfig.useDetailedDoubleMeasurementInfo =
0222       m_cfg.seedFinderConfig.useDetailedDoubleMeasurementInfo;
0223   // Options
0224   Acts::SpacePointContainerOptions spOptions;
0225   spOptions.beamPos = {0., 0.};
0226
0227   // Prepare interface SpacePoint backend-ACTS
0228   ActsExamples::SpacePointContainer container(spacePointPtrs);
0229   // Prepare Acts API
0230   Acts::SpacePointContainer<decltype(container), Acts::detail::RefHolder>
0231       spContainer(spConfig, spOptions, container);
0232
0233   using value_type = typename decltype(spContainer)::SpacePointProxyType;
0234   using seed_type = Acts::Seed<value_type>;
0235
0236   Acts::CylindricalSpacePointGrid<value_type> grid =
0237       Acts::CylindricalSpacePointGridCreator::createGrid<value_type>(
0238           m_cfg.gridConfig, m_cfg.gridOptions, logger());
0239
0240   Acts::CylindricalSpacePointGridCreator::fillGrid<value_type>(
0241       m_cfg.seedFinderConfig, m_cfg.seedFinderOptions, grid, spContainer,
0242       logger());
0243
0244   // Compute radius Range
0245   // we rely on the fact the grid is storing the proxies
0246   // with a sorting in the radius
0247   float minRange = std::numeric_limits<float>::max();
0248   float maxRange = std::numeric_limits<float>::lowest();
0249   for (const auto& coll : grid) {
0250     if (coll.empty()) {
0251       continue;
0252     }
0253     const auto* firstEl = coll.front();
0254     const auto* lastEl = coll.back();
0255     minRange = std::min(firstEl->radius(), minRange);
0256     maxRange = std::max(lastEl->radius(), maxRange);
0257   }
0258
0259   std::array<std::vector<std::size_t>, 3ul> navigation;
0260   navigation[1ul] = m_cfg.seedFinderConfig.zBinsCustomLooping;
0261
0262   auto spacePointsGrouping = Acts::CylindricalBinnedGroup<value_type>(
0263       std::move(grid), *m_bottomBinFinder, *m_topBinFinder,
0264       std::move(navigation));
0265
0266   /// variable middle SP radial region of interest
0267   const Acts::Range1D<float> rMiddleSPRange(
0268       std::floor(minRange / 2) * 2 +
0269           m_cfg.seedFinderConfig.deltaRMiddleMinSPRange,
0270       std::floor(maxRange / 2) * 2 -
0271           m_cfg.seedFinderConfig.deltaRMiddleMaxSPRange);
0272
0273   // run the seeding
0274   static thread_local std::vector<seed_type> seeds;
0275   seeds.clear();
0276   static thread_local decltype(m_seedFinder)::SeedingState state;
0277   state.spacePointMutableData.resize(spContainer.size());
0278
0279   for (const auto [bottom, middle, top] : spacePointsGrouping) {
0280     m_seedFinder.createSeedsForGroup(m_cfg.seedFinderOptions, state,
0281                                      spacePointsGrouping.grid(), seeds, bottom,
0282                                      middle, top, rMiddleSPRange);
0283   }
0284
0285   ACTS_DEBUG("Created " << seeds.size() << " track seeds from "
0286                         << spacePointPtrs.size() << " space points");
0287
0288   // we have seeds of proxies
0289   // convert them to seed of external space points
0290   SimSeedContainer SeedContainerForStorage;
0291   SeedContainerForStorage.reserve(seeds.size());
0292   for (const auto& seed : seeds) {
0293     const auto& sps = seed.sp();
0294     SeedContainerForStorage.emplace_back(*sps[0]->externalSpacePoint(),
0295                                          *sps[1]->externalSpacePoint(),
0296                                          *sps[2]->externalSpacePoint());
0297     SeedContainerForStorage.back().setVertexZ(seed.z());
0298     SeedContainerForStorage.back().setQuality(seed.seedQuality());
0299   }
0300
0301   m_outputSeeds(ctx, std::move(SeedContainerForStorage));
0302   return ActsExamples::ProcessCode::SUCCESS;
0303 }