TrackFinding/src/SeedingAlgorithmHashing.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/SeedingAlgorithmHashing.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/EventData/Seed.hpp"
0013 #include "ActsPlugins/Hashing/HashingAlgorithm.hpp"
0014 #include "ActsPlugins/Hashing/HashingTraining.hpp"
0015
0016 namespace ActsExamples {
0017
0018 namespace {
0019
0020 // Custom seed comparison function
0021 template <typename external_spacepoint_t>
0022 struct SeedComparison {
0023   bool operator()(const Acts::Seed<external_spacepoint_t>& seed1,
0024                   const Acts::Seed<external_spacepoint_t>& seed2) const {
0025     const auto& sp1 = seed1.sp();
0026     const auto& sp2 = seed2.sp();
0027
0028     for (std::size_t i = 0; i < sp1.size(); ++i) {
0029       if (sp1[i]->z() != sp2[i]->z()) {
0030         return sp1[i]->z() < sp2[i]->z();
0031       }
0032     }
0033
0034     for (std::size_t i = 0; i < sp1.size(); ++i) {
0035       if (sp1[i]->x() != sp2[i]->x()) {
0036         return sp1[i]->x() < sp2[i]->x();
0037       }
0038     }
0039
0040     for (std::size_t i = 0; i < sp1.size(); ++i) {
0041       if (sp1[i]->y() != sp2[i]->y()) {
0042         return sp1[i]->y() < sp2[i]->y();
0043       }
0044     }
0045
0046     return false;
0047   }
0048 };
0049
0050 }  // namespace
0051
0052 SeedingAlgorithmHashing::SeedingAlgorithmHashing(
0053     SeedingAlgorithmHashing::Config cfg, Acts::Logging::Level lvl)
0054     : IAlgorithm("SeedingAlgorithmHashing", lvl), m_cfg(std::move(cfg)) {
0055   using SpacePointProxy_type = typename Acts::SpacePointContainer<
0056       SpacePointContainer<std::vector<const SimSpacePoint*>>,
0057       Acts::detail::RefHolder>::SpacePointProxyType;
0058
0059   // Seed Finder config requires Seed Filter object before conversion to
0060   // internal units
0061   m_cfg.seedFinderConfig.seedFilter =
0062       std::make_unique<Acts::SeedFilter<SpacePointProxy_type>>(
0063           m_cfg.seedFilterConfig);
0064
0065   m_cfg.seedFinderConfig = m_cfg.seedFinderConfig.calculateDerivedQuantities();
0066   m_cfg.seedFinderOptions = m_cfg.seedFinderOptions.calculateDerivedQuantities(
0067       m_cfg.seedFinderConfig);
0068   if (m_cfg.inputSpacePoints.empty()) {
0069     throw std::invalid_argument("Missing space point input collections");
0070   }
0071
0072   for (const auto& spName : m_cfg.inputSpacePoints) {
0073     if (spName.empty()) {
0074       throw std::invalid_argument("Invalid space point input collection");
0075     }
0076
0077     auto& handle = m_inputSpacePoints.emplace_back(
0078         std::make_unique<ReadDataHandle<SimSpacePointContainer>>(
0079             this,
0080             "InputSpacePoints#" + std::to_string(m_inputSpacePoints.size())));
0081     handle->initialize(spName);
0082   }
0083   if (m_cfg.outputSeeds.empty()) {
0084     throw std::invalid_argument("Missing seeds output collection");
0085   }
0086
0087   m_outputSeeds.initialize(m_cfg.outputSeeds);
0088   m_outputBuckets.initialize(m_cfg.outputBuckets);
0089
0090   if (m_cfg.gridConfig.rMax != m_cfg.seedFinderConfig.rMax &&
0091       m_cfg.allowSeparateRMax == false) {
0092     throw std::invalid_argument(
0093         "Inconsistent config rMax: using different values in gridConfig and "
0094         "seedFinderConfig. If values are intentional set allowSeparateRMax to "
0095         "true");
0096   }
0097
0098   if (m_cfg.seedFilterConfig.deltaRMin != m_cfg.seedFinderConfig.deltaRMin) {
0099     throw std::invalid_argument("Inconsistent config deltaRMin");
0100   }
0101
0102   if (m_cfg.gridConfig.deltaRMax != m_cfg.seedFinderConfig.deltaRMax) {
0103     throw std::invalid_argument("Inconsistent config deltaRMax");
0104   }
0105
0106   static_assert(
0107       std::numeric_limits<
0108           decltype(m_cfg.seedFinderConfig.deltaRMaxTopSP)>::has_quiet_NaN,
0109       "Value of deltaRMaxTopSP must support NaN values");
0110
0111   static_assert(
0112       std::numeric_limits<
0113           decltype(m_cfg.seedFinderConfig.deltaRMinTopSP)>::has_quiet_NaN,
0114       "Value of deltaRMinTopSP must support NaN values");
0115
0116   static_assert(
0117       std::numeric_limits<
0118           decltype(m_cfg.seedFinderConfig.deltaRMaxBottomSP)>::has_quiet_NaN,
0119       "Value of deltaRMaxBottomSP must support NaN values");
0120
0121   static_assert(
0122       std::numeric_limits<
0123           decltype(m_cfg.seedFinderConfig.deltaRMinBottomSP)>::has_quiet_NaN,
0124       "Value of deltaRMinBottomSP must support NaN values");
0125
0126   if (std::isnan(m_cfg.seedFinderConfig.deltaRMaxTopSP)) {
0127     m_cfg.seedFinderConfig.deltaRMaxTopSP = m_cfg.seedFinderConfig.deltaRMax;
0128   }
0129
0130   if (std::isnan(m_cfg.seedFinderConfig.deltaRMinTopSP)) {
0131     m_cfg.seedFinderConfig.deltaRMinTopSP = m_cfg.seedFinderConfig.deltaRMin;
0132   }
0133
0134   if (std::isnan(m_cfg.seedFinderConfig.deltaRMaxBottomSP)) {
0135     m_cfg.seedFinderConfig.deltaRMaxBottomSP = m_cfg.seedFinderConfig.deltaRMax;
0136   }
0137
0138   if (std::isnan(m_cfg.seedFinderConfig.deltaRMinBottomSP)) {
0139     m_cfg.seedFinderConfig.deltaRMinBottomSP = m_cfg.seedFinderConfig.deltaRMin;
0140   }
0141
0142   if (m_cfg.gridConfig.zMin != m_cfg.seedFinderConfig.zMin) {
0143     throw std::invalid_argument("Inconsistent config zMin");
0144   }
0145
0146   if (m_cfg.gridConfig.zMax != m_cfg.seedFinderConfig.zMax) {
0147     throw std::invalid_argument("Inconsistent config zMax");
0148   }
0149
0150   if (m_cfg.seedFilterConfig.maxSeedsPerSpM !=
0151       m_cfg.seedFinderConfig.maxSeedsPerSpM) {
0152     throw std::invalid_argument("Inconsistent config maxSeedsPerSpM");
0153   }
0154
0155   if (m_cfg.gridConfig.cotThetaMax != m_cfg.seedFinderConfig.cotThetaMax) {
0156     throw std::invalid_argument("Inconsistent config cotThetaMax");
0157   }
0158
0159   if (m_cfg.gridConfig.minPt != m_cfg.seedFinderConfig.minPt) {
0160     throw std::invalid_argument("Inconsistent config minPt");
0161   }
0162
0163   if (m_cfg.gridOptions.bFieldInZ != m_cfg.seedFinderOptions.bFieldInZ) {
0164     throw std::invalid_argument("Inconsistent config bFieldInZ");
0165   }
0166
0167   if (m_cfg.gridConfig.zBinEdges.size() - 1 != m_cfg.zBinNeighborsTop.size() &&
0168       m_cfg.zBinNeighborsTop.empty() == false) {
0169     throw std::invalid_argument("Inconsistent config zBinNeighborsTop");
0170   }
0171
0172   if (m_cfg.gridConfig.zBinEdges.size() - 1 !=
0173           m_cfg.zBinNeighborsBottom.size() &&
0174       m_cfg.zBinNeighborsBottom.empty() == false) {
0175     throw std::invalid_argument("Inconsistent config zBinNeighborsBottom");
0176   }
0177
0178   m_bottomBinFinder = std::make_unique<const Acts::GridBinFinder<3ul>>(
0179       m_cfg.numPhiNeighbors, m_cfg.zBinNeighborsBottom, 0);
0180   m_topBinFinder = std::make_unique<const Acts::GridBinFinder<3ul>>(
0181       m_cfg.numPhiNeighbors, m_cfg.zBinNeighborsTop, 0);
0182
0183   m_cfg.seedFinderConfig.seedFilter =
0184       std::make_unique<Acts::SeedFilter<SpacePointProxy_type>>(
0185           m_cfg.seedFilterConfig);
0186   m_seedFinder =
0187       Acts::SeedFinder<SpacePointProxy_type,
0188                        Acts::CylindricalSpacePointGrid<SpacePointProxy_type>>(
0189           m_cfg.seedFinderConfig);
0190   m_hashing = ActsPlugins::HashingAlgorithm<const SimSpacePoint*,
0191                                             std::vector<const SimSpacePoint*>>(
0192       m_cfg.hashingConfig);
0193   m_hashingTraining =
0194       ActsPlugins::HashingTrainingAlgorithm<std::vector<const SimSpacePoint*>>(
0195           m_cfg.hashingTrainingConfig);
0196 }
0197
0198 ProcessCode SeedingAlgorithmHashing::execute(
0199     const AlgorithmContext& ctx) const {
0200   ACTS_DEBUG("Start of SeedingAlgorithmHashing execute");
0201   using SpacePointPtrVector = std::vector<const SimSpacePoint*>;
0202
0203   // construct the combined input container of space point pointers from all
0204   // configured input sources.
0205   // pre-compute the total size required so we only need to allocate once
0206   std::size_t nSpacePoints = 0;
0207   for (const auto& isp : m_inputSpacePoints) {
0208     nSpacePoints += (*isp)(ctx).size();
0209   }
0210
0211   SpacePointPtrVector spacePointPtrs;
0212   spacePointPtrs.reserve(nSpacePoints);
0213   for (const auto& isp : m_inputSpacePoints) {
0214     for (const SimSpacePoint& spacePoint : (*isp)(ctx)) {
0215       // since the event store owns the space
0216       // points, their pointers should be stable and
0217       // we do not need to create local copies.
0218       spacePointPtrs.push_back(&spacePoint);
0219     }
0220   }
0221
0222   // Config
0223   Acts::SpacePointContainerConfig spConfig;
0224   spConfig.useDetailedDoubleMeasurementInfo =
0225       m_cfg.seedFinderConfig.useDetailedDoubleMeasurementInfo;
0226   // Options
0227   Acts::SpacePointContainerOptions spOptions;
0228   spOptions.beamPos = {0., 0.};
0229
0230   // Hashing Training
0231   ActsPlugins::AnnoyModel annoyModel =
0232       m_hashingTraining.execute(spacePointPtrs);
0233   // Hashing
0234   static thread_local std::vector<SpacePointPtrVector> bucketsPtrs;
0235   bucketsPtrs.clear();
0236   m_hashing.execute(spacePointPtrs, &annoyModel, bucketsPtrs);
0237
0238   // pre-compute the maximum size required so we only need to allocate once
0239   // doesn't combine the input containers of space point pointers
0240   std::size_t maxNSpacePoints = 0;
0241   std::size_t inSpacePoints = 0;
0242   for (const SpacePointPtrVector& bucket : bucketsPtrs) {
0243     inSpacePoints = bucket.size();
0244     if (inSpacePoints > maxNSpacePoints) {
0245       maxNSpacePoints = inSpacePoints;
0246     }
0247   }
0248
0249   using value_type = typename Acts::SpacePointContainer<
0250       SpacePointContainer<std::vector<const SimSpacePoint*>>,
0251       Acts::detail::RefHolder>::SpacePointProxyType;
0252   using seed_type = Acts::Seed<value_type>;
0253
0254   // Create the set with custom comparison function
0255   static thread_local std::set<SimSeed, SeedComparison<SimSpacePoint>> seedsSet;
0256   seedsSet.clear();
0257   static thread_local decltype(m_seedFinder)::SeedingState state;
0258   state.spacePointMutableData.resize(maxNSpacePoints);
0259
0260   for (SpacePointPtrVector& bucket : bucketsPtrs) {
0261     std::set<seed_type, SeedComparison<value_type>> seedsSetForBucket;
0262     state.spacePointMutableData.clear();
0263     state.spacePointMutableData.resize(maxNSpacePoints);
0264
0265     // Prepare interface SpacePoint backend-ACTS
0266     SpacePointContainer container(bucket);
0267     // Prepare Acts API
0268     Acts::SpacePointContainer<decltype(container), Acts::detail::RefHolder>
0269         spContainer(spConfig, spOptions, container);
0270
0271     // construct the seeding tools
0272     Acts::CylindricalSpacePointGrid<value_type> grid =
0273         Acts::CylindricalSpacePointGridCreator::createGrid<value_type>(
0274             m_cfg.gridConfig, m_cfg.gridOptions);
0275     Acts::CylindricalSpacePointGridCreator::fillGrid(
0276         m_cfg.seedFinderConfig, m_cfg.seedFinderOptions, grid,
0277         spContainer.begin(), spContainer.end());
0278
0279     // Compute radius Range
0280     // we rely on the fact the grid is storing the proxies
0281     // with a sorting in the radius
0282     float minRange = std::numeric_limits<float>::max();
0283     float maxRange = std::numeric_limits<float>::lowest();
0284     for (const auto& coll : grid) {
0285       if (coll.empty()) {
0286         continue;
0287       }
0288       const auto* firstEl = coll.front();
0289       const auto* lastEl = coll.back();
0290       minRange = std::min(firstEl->radius(), minRange);
0291       maxRange = std::max(lastEl->radius(), maxRange);
0292     }
0293
0294     std::array<std::vector<std::size_t>, 3ul> navigation;
0295     navigation[1ul] = m_cfg.seedFinderConfig.zBinsCustomLooping;
0296
0297     // groups spacepoints
0298     auto spacePointsGrouping = Acts::CylindricalBinnedGroup<value_type>(
0299         std::move(grid), *m_bottomBinFinder, *m_topBinFinder,
0300         std::move(navigation));
0301
0302     /// variable middle SP radial region of interest
0303     const Acts::Range1D<float> rMiddleSPRange(
0304         minRange + m_cfg.seedFinderConfig.deltaRMiddleMinSPRange,
0305         maxRange - m_cfg.seedFinderConfig.deltaRMiddleMaxSPRange);
0306
0307     // this creates seeds of proxy, we need to convert it to seed of space
0308     // points
0309     for (const auto [bottom, middle, top] : spacePointsGrouping) {
0310       m_seedFinder.createSeedsForGroup(
0311           m_cfg.seedFinderOptions, state, spacePointsGrouping.grid(),
0312           seedsSetForBucket, bottom, middle, top, rMiddleSPRange);
0313     }
0314
0315     // proxies die when the Acts::SpacePointContainer dies, so we need to
0316     // convert the seeds of space points proxies to seeds of space points
0317     for (const seed_type& seed : seedsSetForBucket) {
0318       const std::array<const value_type*, 3>& sps = seed.sp();
0319       const SimSpacePoint* bottom = sps[0]->externalSpacePoint();
0320       const SimSpacePoint* middle = sps[1]->externalSpacePoint();
0321       const SimSpacePoint* top = sps[2]->externalSpacePoint();
0322
0323       SimSeed toAdd(*bottom, *middle, *top);
0324       toAdd.setVertexZ(seed.z());
0325       toAdd.setQuality(seed.seedQuality());
0326
0327       seedsSet.insert(toAdd);
0328     }
0329   }
0330
0331   // convert the set to a simseed collection
0332   SimSeedContainer seeds;
0333   seeds.reserve(seedsSet.size());
0334   for (const SimSeed& seed : seedsSet) {
0335     seeds.push_back(seed);
0336   }
0337
0338   ACTS_INFO("Created " << seeds.size() << " track seeds from "
0339                        << spacePointPtrs.size() << " space points");
0340
0341   m_outputSeeds(ctx, SimSeedContainer{seeds});
0342   std::vector<SimSpacePointContainer> buckets;
0343   for (const SpacePointPtrVector& bucket : bucketsPtrs) {
0344     SimSpacePointContainer bucketSP;
0345     for (const SimSpacePoint* spacePoint : bucket) {
0346       bucketSP.push_back(*spacePoint);
0347     }
0348     buckets.push_back(bucketSP);
0349   }
0350   m_outputBuckets(ctx, std::vector<SimSpacePointContainer>{buckets});
0351
0352   ACTS_DEBUG("End of SeedingAlgorithmHashing execute");
0353   return ProcessCode::SUCCESS;
0354 }
0355
0356 }  // namespace ActsExamples