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/SeedFilter.hpp"
0014 #include "Acts/Utilities/Delegate.hpp"
0015 #include "Acts/Utilities/GridBinFinder.hpp"
0016 #include "ActsExamples/EventData/SimSeed.hpp"
0017
0018 #include <cmath>
0019 #include <csignal>
0020 #include <cstddef>
0021 #include <limits>
0022 #include <ostream>
0023 #include <stdexcept>
0024
0025 using namespace Acts::HashedStringLiteral;
0026
0027 namespace ActsExamples {
0028
0029 SeedingAlgorithm::SeedingAlgorithm(SeedingAlgorithm::Config cfg,
0030                                    Acts::Logging::Level lvl)
0031     : IAlgorithm("SeedingAlgorithm", lvl), m_cfg(std::move(cfg)) {
0032   using SpacePointProxy_type = typename Acts::SpacePointContainer<
0033       SpacePointContainer<std::vector<const SimSpacePoint*>>,
0034       Acts::detail::RefHolder>::SpacePointProxyType;
0035
0036   // Seed Finder config requires Seed Filter object before conversion to
0037   // internal units
0038   m_cfg.seedFinderConfig.seedFilter =
0039       std::make_unique<Acts::SeedFilter<SpacePointProxy_type>>(
0040           m_cfg.seedFilterConfig, logger().cloneWithSuffix("SeedFilter"));
0041   m_cfg.seedFinderConfig = m_cfg.seedFinderConfig.calculateDerivedQuantities();
0042   m_cfg.seedFinderOptions = m_cfg.seedFinderOptions.calculateDerivedQuantities(
0043       m_cfg.seedFinderConfig);
0044   if (m_cfg.inputSpacePoints.empty()) {
0045     throw std::invalid_argument("Missing space point input collections");
0046   }
0047
0048   for (const auto& spName : m_cfg.inputSpacePoints) {
0049     if (spName.empty()) {
0050       throw std::invalid_argument("Invalid space point input collection");
0051     }
0052
0053     auto& handle = m_inputSpacePoints.emplace_back(
0054         std::make_unique<ReadDataHandle<SimSpacePointContainer>>(
0055             this,
0056             "InputSpacePoints#" + std::to_string(m_inputSpacePoints.size())));
0057     handle->initialize(spName);
0058   }
0059   if (m_cfg.outputSeeds.empty()) {
0060     throw std::invalid_argument("Missing seeds output collection");
0061   }
0062
0063   m_outputSeeds.initialize(m_cfg.outputSeeds);
0064
0065   if (m_cfg.gridConfig.rMax != m_cfg.seedFinderConfig.rMax &&
0066       m_cfg.allowSeparateRMax == false) {
0067     throw std::invalid_argument(
0068         "Inconsistent config rMax: using different values in gridConfig and "
0069         "seedFinderConfig. If values are intentional set allowSeparateRMax to "
0070         "true");
0071   }
0072
0073   if (m_cfg.seedFilterConfig.deltaRMin != m_cfg.seedFinderConfig.deltaRMin) {
0074     throw std::invalid_argument("Inconsistent config deltaRMin");
0075   }
0076
0077   if (m_cfg.gridConfig.deltaRMax != m_cfg.seedFinderConfig.deltaRMax) {
0078     throw std::invalid_argument("Inconsistent config deltaRMax");
0079   }
0080
0081   static_assert(
0082       std::numeric_limits<
0083           decltype(m_cfg.seedFinderConfig.deltaRMaxTopSP)>::has_quiet_NaN,
0084       "Value of deltaRMaxTopSP must support NaN values");
0085
0086   static_assert(
0087       std::numeric_limits<
0088           decltype(m_cfg.seedFinderConfig.deltaRMinTopSP)>::has_quiet_NaN,
0089       "Value of deltaRMinTopSP must support NaN values");
0090
0091   static_assert(
0092       std::numeric_limits<
0093           decltype(m_cfg.seedFinderConfig.deltaRMaxBottomSP)>::has_quiet_NaN,
0094       "Value of deltaRMaxBottomSP must support NaN values");
0095
0096   static_assert(
0097       std::numeric_limits<
0098           decltype(m_cfg.seedFinderConfig.deltaRMinBottomSP)>::has_quiet_NaN,
0099       "Value of deltaRMinBottomSP must support NaN values");
0100
0101   if (std::isnan(m_cfg.seedFinderConfig.deltaRMaxTopSP)) {
0102     m_cfg.seedFinderConfig.deltaRMaxTopSP = m_cfg.seedFinderConfig.deltaRMax;
0103   }
0104
0105   if (std::isnan(m_cfg.seedFinderConfig.deltaRMinTopSP)) {
0106     m_cfg.seedFinderConfig.deltaRMinTopSP = m_cfg.seedFinderConfig.deltaRMin;
0107   }
0108
0109   if (std::isnan(m_cfg.seedFinderConfig.deltaRMaxBottomSP)) {
0110     m_cfg.seedFinderConfig.deltaRMaxBottomSP = m_cfg.seedFinderConfig.deltaRMax;
0111   }
0112
0113   if (std::isnan(m_cfg.seedFinderConfig.deltaRMinBottomSP)) {
0114     m_cfg.seedFinderConfig.deltaRMinBottomSP = m_cfg.seedFinderConfig.deltaRMin;
0115   }
0116
0117   if (m_cfg.gridConfig.zMin != m_cfg.seedFinderConfig.zMin) {
0118     throw std::invalid_argument("Inconsistent config zMin");
0119   }
0120
0121   if (m_cfg.gridConfig.zMax != m_cfg.seedFinderConfig.zMax) {
0122     throw std::invalid_argument("Inconsistent config zMax");
0123   }
0124
0125   if (m_cfg.seedFilterConfig.maxSeedsPerSpM !=
0126       m_cfg.seedFinderConfig.maxSeedsPerSpM) {
0127     throw std::invalid_argument("Inconsistent config maxSeedsPerSpM");
0128   }
0129
0130   if (m_cfg.gridConfig.cotThetaMax != m_cfg.seedFinderConfig.cotThetaMax) {
0131     throw std::invalid_argument("Inconsistent config cotThetaMax");
0132   }
0133
0134   if (m_cfg.gridConfig.minPt != m_cfg.seedFinderConfig.minPt) {
0135     throw std::invalid_argument("Inconsistent config minPt");
0136   }
0137
0138   if (m_cfg.gridOptions.bFieldInZ != m_cfg.seedFinderOptions.bFieldInZ) {
0139     throw std::invalid_argument("Inconsistent config bFieldInZ");
0140   }
0141
0142   if (m_cfg.gridConfig.zBinEdges.size() - 1 != m_cfg.zBinNeighborsTop.size() &&
0143       m_cfg.zBinNeighborsTop.empty() == false) {
0144     throw std::invalid_argument("Inconsistent config zBinNeighborsTop");
0145   }
0146
0147   if (m_cfg.gridConfig.zBinEdges.size() - 1 !=
0148           m_cfg.zBinNeighborsBottom.size() &&
0149       m_cfg.zBinNeighborsBottom.empty() == false) {
0150     throw std::invalid_argument("Inconsistent config zBinNeighborsBottom");
0151   }
0152
0153   if (!m_cfg.seedFinderConfig.zBinsCustomLooping.empty()) {
0154     // check that the bins required in the custom bin looping
0155     // are contained in the bins defined by the total number of edges
0156
0157     for (std::size_t i : m_cfg.seedFinderConfig.zBinsCustomLooping) {
0158       if (i >= m_cfg.gridConfig.zBinEdges.size()) {
0159         throw std::invalid_argument(
0160             "Inconsistent config zBinsCustomLooping does not contain a subset "
0161             "of bins defined by zBinEdges");
0162       }
0163     }
0164   }
0165
0166   if (m_cfg.useExtraCuts) {
0167     // This function will be applied to select space points during grid filling
0168     m_cfg.seedFinderConfig.spacePointSelector
0169         .connect<itkFastTrackingSPselect>();
0170
0171     // This function will be applied to the doublet compatibility selection
0172     m_cfg.seedFinderConfig.experimentCuts.connect<itkFastTrackingCuts>();
0173   }
0174
0175   using SpacePointProxy_type = typename Acts::SpacePointContainer<
0176       SpacePointContainer<std::vector<const SimSpacePoint*>>,
0177       Acts::detail::RefHolder>::SpacePointProxyType;
0178
0179   m_bottomBinFinder = std::make_unique<const Acts::GridBinFinder<3ul>>(
0180       m_cfg.numPhiNeighbors, cfg.zBinNeighborsBottom, 0);
0181   m_topBinFinder = std::make_unique<const Acts::GridBinFinder<3ul>>(
0182       m_cfg.numPhiNeighbors, m_cfg.zBinNeighborsTop, 0);
0183
0184   m_seedFinder =
0185       Acts::SeedFinder<SpacePointProxy_type,
0186                        Acts::CylindricalSpacePointGrid<SpacePointProxy_type>>(
0187           m_cfg.seedFinderConfig, logger().cloneWithSuffix("SeedFinder"));
0188 }
0189
0190 ProcessCode SeedingAlgorithm::execute(const AlgorithmContext& ctx) const {
0191   // construct the combined input container of space point pointers from all
0192   // configured input sources.
0193   // pre-compute the total size required so we only need to allocate once
0194   std::size_t nSpacePoints = 0;
0195   for (const auto& isp : m_inputSpacePoints) {
0196     nSpacePoints += (*isp)(ctx).size();
0197   }
0198
0199   std::vector<const SimSpacePoint*> spacePointPtrs;
0200   spacePointPtrs.reserve(nSpacePoints);
0201   for (const auto& isp : m_inputSpacePoints) {
0202     for (const auto& spacePoint : (*isp)(ctx)) {
0203       // since the event store owns the space
0204       // points, their pointers should be stable and
0205       // we do not need to create local copies.
0206       spacePointPtrs.push_back(&spacePoint);
0207     }
0208   }
0209
0210   // Config
0211   Acts::SpacePointContainerConfig spConfig;
0212   spConfig.useDetailedDoubleMeasurementInfo =
0213       m_cfg.seedFinderConfig.useDetailedDoubleMeasurementInfo;
0214   // Options
0215   Acts::SpacePointContainerOptions spOptions;
0216   spOptions.beamPos = {0., 0.};
0217
0218   // Prepare interface SpacePoint backend-ACTS
0219   SpacePointContainer container(spacePointPtrs);
0220   // Prepare Acts API
0221   Acts::SpacePointContainer<decltype(container), Acts::detail::RefHolder>
0222       spContainer(spConfig, spOptions, container);
0223
0224   using value_type = typename decltype(spContainer)::SpacePointProxyType;
0225   using seed_type = Acts::Seed<value_type>;
0226
0227   Acts::CylindricalSpacePointGrid<value_type> grid =
0228       Acts::CylindricalSpacePointGridCreator::createGrid<value_type>(
0229           m_cfg.gridConfig, m_cfg.gridOptions, logger());
0230
0231   Acts::CylindricalSpacePointGridCreator::fillGrid<value_type>(
0232       m_cfg.seedFinderConfig, m_cfg.seedFinderOptions, grid, spContainer,
0233       logger());
0234
0235   // Compute radius Range
0236   // we rely on the fact the grid is storing the proxies
0237   // with a sorting in the radius
0238   float minRange = std::numeric_limits<float>::max();
0239   float maxRange = std::numeric_limits<float>::lowest();
0240   for (const auto& coll : grid) {
0241     if (coll.empty()) {
0242       continue;
0243     }
0244     const auto* firstEl = coll.front();
0245     const auto* lastEl = coll.back();
0246     minRange = std::min(firstEl->radius(), minRange);
0247     maxRange = std::max(lastEl->radius(), maxRange);
0248   }
0249
0250   std::array<std::vector<std::size_t>, 3ul> navigation;
0251   navigation[1ul] = m_cfg.seedFinderConfig.zBinsCustomLooping;
0252
0253   auto spacePointsGrouping = Acts::CylindricalBinnedGroup<value_type>(
0254       std::move(grid), *m_bottomBinFinder, *m_topBinFinder,
0255       std::move(navigation));
0256
0257   /// variable middle SP radial region of interest
0258   const Acts::Range1D<float> rMiddleSPRange(
0259       std::floor(minRange / 2) * 2 +
0260           m_cfg.seedFinderConfig.deltaRMiddleMinSPRange,
0261       std::floor(maxRange / 2) * 2 -
0262           m_cfg.seedFinderConfig.deltaRMiddleMaxSPRange);
0263
0264   // run the seeding
0265   static thread_local std::vector<seed_type> seeds;
0266   seeds.clear();
0267   static thread_local decltype(m_seedFinder)::SeedingState state;
0268   state.spacePointMutableData.resize(spContainer.size());
0269
0270   for (const auto [bottom, middle, top] : spacePointsGrouping) {
0271     m_seedFinder.createSeedsForGroup(m_cfg.seedFinderOptions, state,
0272                                      spacePointsGrouping.grid(), seeds, bottom,
0273                                      middle, top, rMiddleSPRange);
0274   }
0275
0276   ACTS_DEBUG("Created " << seeds.size() << " track seeds from "
0277                         << spacePointPtrs.size() << " space points");
0278
0279   // we have seeds of proxies
0280   // convert them to seed of external space points
0281   SimSeedContainer SeedContainerForStorage;
0282   SeedContainerForStorage.reserve(seeds.size());
0283   for (const auto& seed : seeds) {
0284     const auto& sps = seed.sp();
0285     SeedContainerForStorage.emplace_back(*sps[0]->externalSpacePoint(),
0286                                          *sps[1]->externalSpacePoint(),
0287                                          *sps[2]->externalSpacePoint());
0288     SeedContainerForStorage.back().setVertexZ(seed.z());
0289     SeedContainerForStorage.back().setQuality(seed.seedQuality());
0290   }
0291
0292   m_outputSeeds(ctx, std::move(SeedContainerForStorage));
0293   return ProcessCode::SUCCESS;
0294 }
0295
0296 }  // namespace ActsExamples