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 "Acts/EventData/SpacePointData.hpp"
0014 #include "Acts/Plugins/Hashing/HashingAlgorithm.hpp"
0015 #include "Acts/Plugins/Hashing/HashingTraining.hpp"
0016 #include "Acts/Seeding/BinnedGroup.hpp"
0017 #include "Acts/Seeding/SeedFilter.hpp"
0018 #include "Acts/Seeding/SeedFinder.hpp"
0019 #include "Acts/Surfaces/Surface.hpp"
0020 #include "Acts/Utilities/BinningType.hpp"
0021 #include "Acts/Utilities/Delegate.hpp"
0022 #include "Acts/Utilities/Grid.hpp"
0023 #include "Acts/Utilities/GridBinFinder.hpp"
0024 #include "Acts/Utilities/Helpers.hpp"
0025 #include "ActsExamples/EventData/ProtoTrack.hpp"
0026 #include "ActsExamples/EventData/SimSeed.hpp"
0027 #include "ActsExamples/Framework/WhiteBoard.hpp"
0028
0029 #include <cmath>
0030 #include <csignal>
0031
0032 namespace ActsExamples {
0033 struct AlgorithmContext;
0034 }  // namespace ActsExamples
0035
0036 ActsExamples::SeedingAlgorithmHashing::SeedingAlgorithmHashing(
0037     ActsExamples::SeedingAlgorithmHashing::Config cfg, Acts::Logging::Level lvl)
0038     : ActsExamples::IAlgorithm("SeedingAlgorithmHashing", lvl),
0039       m_cfg(std::move(cfg)) {
0040   using SpacePointProxy_type = typename Acts::SpacePointContainer<
0041       ActsExamples::SpacePointContainer<std::vector<const SimSpacePoint*>>,
0042       Acts::detail::RefHolder>::SpacePointProxyType;
0043
0044   // Seed Finder config requires Seed Filter object before conversion to
0045   // internal units
0046   m_cfg.seedFilterConfig = m_cfg.seedFilterConfig.toInternalUnits();
0047   m_cfg.seedFinderConfig.seedFilter =
0048       std::make_unique<Acts::SeedFilter<SpacePointProxy_type>>(
0049           m_cfg.seedFilterConfig);
0050
0051   m_cfg.seedFinderConfig =
0052       m_cfg.seedFinderConfig.toInternalUnits().calculateDerivedQuantities();
0053   m_cfg.seedFinderOptions =
0054       m_cfg.seedFinderOptions.toInternalUnits().calculateDerivedQuantities(
0055           m_cfg.seedFinderConfig);
0056   m_cfg.gridConfig = m_cfg.gridConfig.toInternalUnits();
0057   m_cfg.gridOptions = m_cfg.gridOptions.toInternalUnits();
0058   if (m_cfg.inputSpacePoints.empty()) {
0059     throw std::invalid_argument("Missing space point input collections");
0060   }
0061
0062   for (const auto& spName : m_cfg.inputSpacePoints) {
0063     if (spName.empty()) {
0064       throw std::invalid_argument("Invalid space point input collection");
0065     }
0066
0067     auto& handle = m_inputSpacePoints.emplace_back(
0068         std::make_unique<ReadDataHandle<SimSpacePointContainer>>(
0069             this,
0070             "InputSpacePoints#" + std::to_string(m_inputSpacePoints.size())));
0071     handle->initialize(spName);
0072   }
0073   if (m_cfg.outputSeeds.empty()) {
0074     throw std::invalid_argument("Missing seeds output collection");
0075   }
0076
0077   m_outputSeeds.initialize(m_cfg.outputSeeds);
0078   m_outputBuckets.initialize(m_cfg.outputBuckets);
0079
0080   if (m_cfg.gridConfig.rMax != m_cfg.seedFinderConfig.rMax &&
0081       m_cfg.allowSeparateRMax == false) {
0082     throw std::invalid_argument(
0083         "Inconsistent config rMax: using different values in gridConfig and "
0084         "seedFinderConfig. If values are intentional set allowSeparateRMax to "
0085         "true");
0086   }
0087
0088   if (m_cfg.seedFilterConfig.deltaRMin != m_cfg.seedFinderConfig.deltaRMin) {
0089     throw std::invalid_argument("Inconsistent config deltaRMin");
0090   }
0091
0092   if (m_cfg.gridConfig.deltaRMax != m_cfg.seedFinderConfig.deltaRMax) {
0093     throw std::invalid_argument("Inconsistent config deltaRMax");
0094   }
0095
0096   static_assert(
0097       std::numeric_limits<
0098           decltype(m_cfg.seedFinderConfig.deltaRMaxTopSP)>::has_quiet_NaN,
0099       "Value of deltaRMaxTopSP must support NaN values");
0100
0101   static_assert(
0102       std::numeric_limits<
0103           decltype(m_cfg.seedFinderConfig.deltaRMinTopSP)>::has_quiet_NaN,
0104       "Value of deltaRMinTopSP must support NaN values");
0105
0106   static_assert(
0107       std::numeric_limits<
0108           decltype(m_cfg.seedFinderConfig.deltaRMaxBottomSP)>::has_quiet_NaN,
0109       "Value of deltaRMaxBottomSP must support NaN values");
0110
0111   static_assert(
0112       std::numeric_limits<
0113           decltype(m_cfg.seedFinderConfig.deltaRMinBottomSP)>::has_quiet_NaN,
0114       "Value of deltaRMinBottomSP must support NaN values");
0115
0116   if (std::isnan(m_cfg.seedFinderConfig.deltaRMaxTopSP)) {
0117     m_cfg.seedFinderConfig.deltaRMaxTopSP = m_cfg.seedFinderConfig.deltaRMax;
0118   }
0119
0120   if (std::isnan(m_cfg.seedFinderConfig.deltaRMinTopSP)) {
0121     m_cfg.seedFinderConfig.deltaRMinTopSP = m_cfg.seedFinderConfig.deltaRMin;
0122   }
0123
0124   if (std::isnan(m_cfg.seedFinderConfig.deltaRMaxBottomSP)) {
0125     m_cfg.seedFinderConfig.deltaRMaxBottomSP = m_cfg.seedFinderConfig.deltaRMax;
0126   }
0127
0128   if (std::isnan(m_cfg.seedFinderConfig.deltaRMinBottomSP)) {
0129     m_cfg.seedFinderConfig.deltaRMinBottomSP = m_cfg.seedFinderConfig.deltaRMin;
0130   }
0131
0132   if (m_cfg.gridConfig.zMin != m_cfg.seedFinderConfig.zMin) {
0133     throw std::invalid_argument("Inconsistent config zMin");
0134   }
0135
0136   if (m_cfg.gridConfig.zMax != m_cfg.seedFinderConfig.zMax) {
0137     throw std::invalid_argument("Inconsistent config zMax");
0138   }
0139
0140   if (m_cfg.seedFilterConfig.maxSeedsPerSpM !=
0141       m_cfg.seedFinderConfig.maxSeedsPerSpM) {
0142     throw std::invalid_argument("Inconsistent config maxSeedsPerSpM");
0143   }
0144
0145   if (m_cfg.gridConfig.cotThetaMax != m_cfg.seedFinderConfig.cotThetaMax) {
0146     throw std::invalid_argument("Inconsistent config cotThetaMax");
0147   }
0148
0149   if (m_cfg.gridConfig.minPt != m_cfg.seedFinderConfig.minPt) {
0150     throw std::invalid_argument("Inconsistent config minPt");
0151   }
0152
0153   if (m_cfg.gridOptions.bFieldInZ != m_cfg.seedFinderOptions.bFieldInZ) {
0154     throw std::invalid_argument("Inconsistent config bFieldInZ");
0155   }
0156
0157   if (m_cfg.gridConfig.zBinEdges.size() - 1 != m_cfg.zBinNeighborsTop.size() &&
0158       m_cfg.zBinNeighborsTop.empty() == false) {
0159     throw std::invalid_argument("Inconsistent config zBinNeighborsTop");
0160   }
0161
0162   if (m_cfg.gridConfig.zBinEdges.size() - 1 !=
0163           m_cfg.zBinNeighborsBottom.size() &&
0164       m_cfg.zBinNeighborsBottom.empty() == false) {
0165     throw std::invalid_argument("Inconsistent config zBinNeighborsBottom");
0166   }
0167
0168   if (m_cfg.useExtraCuts) {
0169     // This function will be applied to select space points during grid filling
0170     m_cfg.seedFinderConfig.spacePointSelector
0171         .connect<itkFastTrackingSPselect>();
0172
0173     // This function will be applied to the doublet compatibility selection
0174     m_cfg.seedFinderConfig.experimentCuts.connect<itkFastTrackingCuts>();
0175   }
0176
0177   m_bottomBinFinder = std::make_unique<const Acts::GridBinFinder<3ul>>(
0178       m_cfg.numPhiNeighbors, m_cfg.zBinNeighborsBottom, 0);
0179   m_topBinFinder = std::make_unique<const Acts::GridBinFinder<3ul>>(
0180       m_cfg.numPhiNeighbors, m_cfg.zBinNeighborsTop, 0);
0181
0182   m_cfg.seedFinderConfig.seedFilter =
0183       std::make_unique<Acts::SeedFilter<SpacePointProxy_type>>(
0184           m_cfg.seedFilterConfig);
0185   m_seedFinder =
0186       Acts::SeedFinder<SpacePointProxy_type,
0187                        Acts::CylindricalSpacePointGrid<SpacePointProxy_type>>(
0188           m_cfg.seedFinderConfig);
0189   m_hashing = Acts::HashingAlgorithm<const SimSpacePoint*,
0190                                      std::vector<const SimSpacePoint*>>(
0191       m_cfg.hashingConfig);
0192   m_hashingTraining =
0193       Acts::HashingTrainingAlgorithm<std::vector<const SimSpacePoint*>>(
0194           m_cfg.hashingTrainingConfig);
0195 }
0196
0197 ActsExamples::ProcessCode ActsExamples::SeedingAlgorithmHashing::execute(
0198     const AlgorithmContext& ctx) const {
0199   ACTS_DEBUG("Start of SeedingAlgorithmHashing execute");
0200   using SpacePointPtrVector = std::vector<const SimSpacePoint*>;
0201
0202   // construct the combined input container of space point pointers from all
0203   // configured input sources.
0204   // pre-compute the total size required so we only need to allocate once
0205   std::size_t nSpacePoints = 0;
0206   for (const auto& isp : m_inputSpacePoints) {
0207     nSpacePoints += (*isp)(ctx).size();
0208   }
0209
0210   SpacePointPtrVector spacePointPtrs;
0211   spacePointPtrs.reserve(nSpacePoints);
0212   for (const auto& isp : m_inputSpacePoints) {
0213     for (const SimSpacePoint& spacePoint : (*isp)(ctx)) {
0214       // since the event store owns the space
0215       // points, their pointers should be stable and
0216       // we do not need to create local copies.
0217       spacePointPtrs.push_back(&spacePoint);
0218     }
0219   }
0220
0221   // Config
0222   Acts::SpacePointContainerConfig spConfig;
0223   spConfig.useDetailedDoubleMeasurementInfo =
0224       m_cfg.seedFinderConfig.useDetailedDoubleMeasurementInfo;
0225   // Options
0226   Acts::SpacePointContainerOptions spOptions;
0227   spOptions.beamPos = {0., 0.};
0228
0229   // Hashing Training
0230   Acts::AnnoyModel annoyModel = m_hashingTraining.execute(spacePointPtrs);
0231   // Hashing
0232   static thread_local std::vector<SpacePointPtrVector> bucketsPtrs;
0233   bucketsPtrs.clear();
0234   m_hashing.execute(spacePointPtrs, &annoyModel, bucketsPtrs);
0235
0236   // pre-compute the maximum size required so we only need to allocate once
0237   // doesn't combine the input containers of space point pointers
0238   std::size_t maxNSpacePoints = 0;
0239   std::size_t inSpacePoints = 0;
0240   for (const SpacePointPtrVector& bucket : bucketsPtrs) {
0241     inSpacePoints = bucket.size();
0242     if (inSpacePoints > maxNSpacePoints) {
0243       maxNSpacePoints = inSpacePoints;
0244     }
0245   }
0246
0247   using value_type = typename Acts::SpacePointContainer<
0248       ActsExamples::SpacePointContainer<std::vector<const SimSpacePoint*>>,
0249       Acts::detail::RefHolder>::SpacePointProxyType;
0250   using seed_type = Acts::Seed<value_type>;
0251
0252   // Create the set with custom comparison function
0253   static thread_local std::set<ActsExamples::SimSeed,
0254                                SeedComparison<SimSpacePoint>>
0255       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     ActsExamples::SpacePointContainer<SpacePointPtrVector> 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       ActsExamples::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   ActsExamples::SimSeedContainer seeds;
0333   seeds.reserve(seedsSet.size());
0334   for (const ActsExamples::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 ActsExamples::ProcessCode::SUCCESS;
0354 }