TrackFinding/src/AdaptiveHoughTransformSeeder.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/AdaptiveHoughTransformSeeder.hpp"
0010
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/EventData/SourceLink.hpp"
0013 #include "Acts/Utilities/ScopedTimer.hpp"
0014 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0015
0016 #include <algorithm>
0017 #include <cmath>
0018 #include <iterator>
0019 #include <map>
0020 #include <numeric>
0021 #include <ostream>
0022 #include <stdexcept>
0023
0024 namespace ActsExamples {
0025
0026 // Helper class describing one section of the accumulator space
0027 AccumulatorSection::AccumulatorSection(float xw, float yw, float xBegin,
0028                                        float yBegin, int div,
0029                                        const std::vector<unsigned> &indices,
0030                                        const std::vector<float> &history)
0031     : m_xSize(xw),
0032       m_ySize(yw),
0033       m_xBegin(xBegin),
0034       m_yBegin(yBegin),
0035       m_divisionLevel(div),
0036       m_indices(indices),
0037       m_history(history) {}
0038
0039 void AccumulatorSection::updateDimensions(float xw, float yw, float xBegin,
0040                                           float yBegin) {
0041   m_xSize = xw;
0042   m_ySize = yw;
0043   m_xBegin = xBegin;
0044   m_yBegin = yBegin;
0045 }
0046
0047 void AccumulatorSection::expand(float xs, float ys) {
0048   m_xBegin = m_xBegin + 0.5f * m_xSize - m_xSize * xs * 0.5f;
0049   m_yBegin = m_yBegin + 0.5f * m_ySize - m_ySize * ys * 0.5f;
0050   m_xSize *= xs;
0051   m_ySize *= ys;
0052 }
0053
0054 AccumulatorSection AccumulatorSection::bottomLeft(float xFraction,
0055                                                   float yFraction) const {
0056   return AccumulatorSection(m_xSize * xFraction, m_ySize * yFraction, m_xBegin,
0057                             m_yBegin, m_divisionLevel + 1, m_indices,
0058                             m_history);
0059 }
0060 AccumulatorSection AccumulatorSection::topLeft(float xFraction,
0061                                                float yFraction) const {
0062   return AccumulatorSection(m_xSize * xFraction, m_ySize * yFraction, m_xBegin,
0063                             m_yBegin + m_ySize - m_ySize * yFraction,
0064                             m_divisionLevel + 1, m_indices, m_history);
0065 }
0066 AccumulatorSection AccumulatorSection::topRight(float xFraction,
0067                                                 float yFraction) const {
0068   return AccumulatorSection(m_xSize * xFraction, m_ySize * yFraction,
0069                             m_xBegin + m_xSize - m_xSize * xFraction,
0070                             m_yBegin + m_ySize - m_ySize * yFraction,
0071                             m_divisionLevel + 1, m_indices, m_history);
0072 }
0073 AccumulatorSection AccumulatorSection::bottomRight(float xFraction,
0074                                                    float yFraction) const {
0075   return AccumulatorSection(m_xSize * xFraction, m_ySize * yFraction,
0076                             m_xBegin + m_xSize - m_xSize * xFraction, m_yBegin,
0077                             m_divisionLevel + 1, m_indices, m_history);
0078 }
0079 AccumulatorSection AccumulatorSection::bottom(float yFraction) const {
0080   return bottomLeft(1.0, yFraction);
0081 }
0082 AccumulatorSection AccumulatorSection::top(float yFraction) const {
0083   return topLeft(1.0, yFraction);
0084 }
0085 AccumulatorSection AccumulatorSection::left(float xFraction) const {
0086   return bottomLeft(xFraction, 1.0);
0087 }
0088 AccumulatorSection AccumulatorSection::right(float xFraction) const {
0089   return bottomRight(xFraction, 1.0);
0090 }
0091
0092 AdaptiveHoughTransformSeeder::AdaptiveHoughTransformSeeder(
0093     ActsExamples::AdaptiveHoughTransformSeeder::Config cfg,
0094     Acts::Logging::Level lvl)
0095     : ActsExamples::IAlgorithm("AdaptiveHoughTransformSeeder", lvl),
0096       m_cfg(std::move(cfg)),
0097       m_logger(Acts::getDefaultLogger("AdaptiveHoughTransformSeeder", lvl)) {
0098   for (const auto &spName : config().inputSpacePoints) {
0099     const auto &handle = m_inputSpacePoints.emplace_back(
0100         std::make_unique<ReadDataHandle<SimSpacePointContainer>>(
0101             this,
0102             "InputSpacePoints#" + std::to_string(m_inputSpacePoints.size())));
0103     handle->initialize(spName);
0104   }
0105   if (config().outputSeeds.empty()) {
0106     throw std::invalid_argument(
0107         "AdaptiveHoughTransformSeeder: Output seeds collection name empty");
0108   }
0109   m_outputSeeds.initialize(config().outputSeeds);
0110 }
0111
0112 // When traversing parameters space in various directions it is useful
0113 // to record the "path".
0114 // AccumulatorSection has simple functionality allowing that in a form of a
0115 // plain vector of floats. These numbers need to be accessed consistently, thus
0116 // this indices.
0117 namespace {
0118 constexpr unsigned int phiSplitMinIndex = 0;
0119 constexpr unsigned int phiSplitWidthIndex = 1;
0120 constexpr unsigned int zSplitMinIndex = 2;
0121 constexpr unsigned int zSplitWidthIndex = 3;
0122 constexpr unsigned int cotThetaSplitMinIndex = 4;
0123 constexpr unsigned int cotThetaSplitWidthIndex = 5;
0124 }  // namespace
0125
0126 ProcessCode AdaptiveHoughTransformSeeder::execute(
0127     const AlgorithmContext &ctx) const {
0128   // get inputs
0129   std::vector<PreprocessedMeasurement> measurements;
0130   preparePreprocessedMeasurements(ctx, measurements);
0131
0132   // prepare initial stack
0133   std::vector<AccumulatorSection> stack1;
0134   fillStackPhiSplit(stack1, measurements);
0135
0136   // split into regions in z_vertex cot theta, there is a lot of duplication and
0137   // it can be optimized better in the future
0138   for (auto &section : stack1) {
0139     section.updateDimensions(2.0f * config().zRange,
0140                              2.0f * config().cotThetaRange, -config().zRange,
0141                              -config().cotThetaRange);
0142   }
0143   std::vector<AccumulatorSection> stack2;
0144   {
0145     Acts::ScopedTimer st("splitInZCotTheta", logger());
0146     processStackZCotThetaSplit(stack1, stack2, measurements);
0147     if (config().doSecondPhase) {
0148       for (AccumulatorSection &section : stack2) {
0149         section.setHistory(zSplitMinIndex, section.xBegin());
0150         section.setHistory(zSplitWidthIndex, section.xSize());
0151         section.setHistory(cotThetaSplitMinIndex, section.yBegin());
0152         section.setHistory(cotThetaSplitWidthIndex, section.ySize());
0153       }
0154     }
0155   }
0156   ACTS_DEBUG("past StackZCotThetaSplit stack1 size "
0157              << stack1.size() << " stack2 size " << stack2.size());
0158   for (auto &section : stack2) {
0159     ACTS_DEBUG("Post z_vertex cot theta regions "
0160                << section.count() << " for region starting from "
0161                << section.xBegin() << " " << section.yBegin());
0162     // now need to change search space into phi - q/pT, section covers phi range
0163     // from initial splitting (therefore needed history)
0164     section.updateDimensions(
0165         section.history(phiSplitWidthIndex), 2.0f * config().qOverPtMin,
0166         section.history(phiSplitMinIndex), -config().qOverPtMin);
0167   }
0168   {
0169     Acts::ScopedTimer st("processQOverPtPhi", logger());
0170     processStackQOverPtPhi(stack2, stack1, measurements);
0171   }
0172   ACTS_DEBUG("past StackQOverPtPhi stack1 size "
0173              << stack1.size() << " stack2 size " << stack2.size());
0174
0175   if (m_cfg.deduplicate) {
0176     Acts::ScopedTimer st("deduplication", logger());
0177     deduplicate(stack1);
0178     ACTS_DEBUG("past deduplication " << stack1.size());
0179   }
0180
0181   // do scan in z_vertex - cot theta
0182   if (config().doSecondPhase) {
0183     Acts::ScopedTimer st("secondPhase", logger());
0184
0185     for (auto &section : stack1) {
0186       section.updateDimensions(section.history(zSplitWidthIndex),
0187                                section.history(cotThetaSplitWidthIndex),
0188                                section.history(zSplitMinIndex),
0189                                section.history(cotThetaSplitMinIndex));
0190     }
0191     processStackZCotTheta(stack1, stack2, measurements);
0192     ACTS_DEBUG("past StackZCotTheta stack1 size "
0193                << stack1.size() << " stack2 size " << stack2.size());
0194
0195     if (m_cfg.deduplicate) {
0196       Acts::ScopedTimer st2("deduplication", logger());
0197       deduplicate(stack2);
0198       ACTS_DEBUG("past second deduplication " << stack2.size());
0199     }
0200   }
0201
0202   std::vector<AccumulatorSection> &solutions =
0203       config().doSecondPhase ? stack2 : stack1;
0204   Acts::ScopedTimer st("seedsMaking", logger());
0205
0206   // post solutions
0207   SimSeedContainer seeds;
0208   seeds.reserve(solutions.size());
0209
0210   ACTS_VERBOSE("Number of solutions " << solutions.size());
0211   makeSeeds(seeds, solutions, measurements);
0212
0213   m_outputSeeds(ctx, std::move(seeds));
0214
0215   return ActsExamples::ProcessCode::SUCCESS;
0216 }
0217
0218 void AdaptiveHoughTransformSeeder::preparePreprocessedMeasurements(
0219     const AlgorithmContext &ctx,
0220     std::vector<PreprocessedMeasurement> &measurements) const {
0221   for (const auto &isp : m_inputSpacePoints) {
0222     const auto &spContainer = (*isp)(ctx);
0223     ACTS_DEBUG("Inserting " << spContainer.size()
0224                             << " space points from collection\"" << isp->key()
0225                             << "\"");
0226
0227     for (const SimSpacePoint &sp : spContainer) {
0228       const double phi = std::atan2(sp.y(), sp.x());
0229       const double invr = 1.0 / sp.r();
0230       measurements.emplace_back(
0231           invr, phi, sp.z(),
0232           Acts::SourceLink(static_cast<const SimSpacePoint *>(&sp)));
0233       // wrap phi by duplicating some seeds
0234       if (phi < -std::numbers::pi + config().phiWrap * std::numbers::pi) {
0235         measurements.emplace_back(
0236             invr, phi + 2 * std::numbers::pi, sp.z(),
0237             Acts::SourceLink(static_cast<const SimSpacePoint *>(&sp)));
0238       }
0239       if (phi > std::numbers::pi - config().phiWrap * std::numbers::pi) {
0240         measurements.emplace_back(
0241             invr, phi - 2 * std::numbers::pi, sp.z(),
0242             Acts::SourceLink(static_cast<const SimSpacePoint *>(&sp)));
0243       }
0244     }
0245   }
0246   ACTS_DEBUG("Total number of " << measurements.size()
0247                                 << " used to account for phi wrapping");
0248 }
0249
0250 void AdaptiveHoughTransformSeeder::fillStackPhiSplit(
0251     std::vector<AccumulatorSection> &stack,
0252     const std::vector<PreprocessedMeasurement> &measurements) const {
0253   Acts::ScopedTimer st("splitInQuadrants", logger());
0254   const int nSplits = 8;
0255   const auto wedgeWidth = static_cast<float>(2.0 * std::numbers::pi / nSplits);
0256
0257   for (int phiIndex = 0; phiIndex < nSplits; phiIndex++) {
0258     const auto startPhi = static_cast<float>(
0259         wedgeWidth * static_cast<float>(phiIndex) - std::numbers::pi);
0260     stack.emplace_back(1.2f * wedgeWidth, 2.0f * config().qOverPtMin, startPhi,
0261                        -config().qOverPtMin);
0262     stack.back().indices().resize(measurements.size());
0263     std::iota(std::begin(stack.back().indices()),
0264               std::end(stack.back().indices()), 0);
0265     stack.back().setHistory(phiSplitMinIndex, startPhi);
0266     stack.back().setHistory(phiSplitWidthIndex, wedgeWidth);
0267
0268     updateSection(stack.back(), measurements, m_qOverPtPhiLineParams);
0269     ACTS_DEBUG("Initial split " << stack.back().count()
0270                                 << " for region starting from " << startPhi);
0271     for (unsigned index : stack.back().indices()) {
0272       const PreprocessedMeasurement &m = measurements[index];
0273       ACTS_DEBUG("phi section "
0274                  << startPhi << " x: " << std::cos(m.phi) / m.invr
0275                  << " y: " << std::sin(m.phi) / m.invr << " z: " << m.z);
0276     }
0277   }
0278 }
0279
0280 void AdaptiveHoughTransformSeeder::processStackQOverPtPhi(
0281     std::vector<AccumulatorSection> &input,
0282     std::vector<AccumulatorSection> &output,
0283     const std::vector<PreprocessedMeasurement> &measurements) const {
0284   struct Stats {
0285     double area{};
0286     int nSections{};
0287     int nLines{};
0288     int discardedByThresholdCut{};
0289     int discardedByCrossingCut{};
0290   };
0291   std::map<int, Stats> sStat;
0292   ExplorationOptions opt;
0293   opt.xMinBinSize = config().phiMinBinSize;
0294   opt.yMinBinSize = config().qOverPtMinBinSize;
0295   opt.lineParamFunctor = m_qOverPtPhiLineParams;
0296   opt.decisionFunctor = [&sStat, &cfg = m_cfg, &opt, this](
0297                             const AccumulatorSection &section,
0298                             const std::vector<PreprocessedMeasurement> &mes) {
0299     using enum ExplorationOptions<PreprocessedMeasurement>::Decision;
0300     if (section.divisionLevel() <= 8) {
0301       return Drill;
0302     }
0303
0304     if (section.count() < cfg.threshold) {
0305       sStat[section.divisionLevel()].discardedByThresholdCut += 1;
0306       return Discard;
0307     }
0308     if (section.count() < 3 * cfg.threshold &&
0309         !passIntersectionsCheck(section, mes, opt.lineParamFunctor,
0310                                 cfg.threshold * (cfg.threshold - 1))) {
0311       sStat[section.divisionLevel()].discardedByCrossingCut += 1;
0312       return Discard;
0313     }
0314
0315     if (section.count() >= cfg.threshold &&
0316         section.count() <= cfg.noiseThreshold &&
0317         section.xSize() <= cfg.phiMinBinSize &&
0318         section.ySize() <= cfg.qOverPtMinBinSize) {
0319       return Accept;
0320     }
0321     sStat[section.divisionLevel()].area += section.xSize() * section.ySize();
0322     sStat[section.divisionLevel()].nSections += 1;
0323     sStat[section.divisionLevel()].nLines += section.count();
0324     return Drill;
0325   };
0326
0327   exploreParametersSpace(input, measurements, opt, output);
0328   const unsigned nl =
0329       std::accumulate(output.begin(), output.end(), 0,
0330                       [](unsigned sum, const AccumulatorSection &s) {
0331                         return sum + s.count();
0332                       });
0333   ACTS_DEBUG("size " << sStat.size());
0334   for (const auto &[div, stats] : sStat) {
0335     ACTS_DEBUG("Area in used by div: "
0336                << div << " area: " << stats.area << " avglines: "
0337                << (stats.nSections > 0 ? stats.nLines / stats.nSections : 0)
0338                << " n sections: " << stats.nSections
0339                << " thresholdCut: " << stats.discardedByThresholdCut
0340                << " crossingCut: " << stats.discardedByCrossingCut);
0341   }
0342   ACTS_DEBUG("Phi - qOverPt scan finds "
0343              << output.size() << " solutions, included measurements " << nl);
0344 }
0345
0346 void AdaptiveHoughTransformSeeder::processStackZCotTheta(
0347     std::vector<AccumulatorSection> &input,
0348     std::vector<AccumulatorSection> &output,
0349     const std::vector<PreprocessedMeasurement> &measurements) const {
0350   ExplorationOptions opt;
0351   opt.xMinBinSize = config().zMinBinSize;
0352   opt.yMinBinSize = config().cotThetaMinBinSize;
0353   opt.lineParamFunctor = m_zCotThetaLineParams;
0354   opt.decisionFunctor = [&cfg = m_cfg](
0355                             const AccumulatorSection &section,
0356                             const std::vector<PreprocessedMeasurement> &) {
0357     using enum ExplorationOptions<PreprocessedMeasurement>::Decision;
0358
0359     if (section.count() < cfg.threshold) {
0360       return Discard;
0361     }
0362     if (section.count() >= cfg.threshold &&
0363         section.xSize() <= cfg.zMinBinSize &&
0364         section.ySize() <= cfg.cotThetaMinBinSize) {
0365       return Accept;
0366     }
0367     return Drill;
0368   };
0369
0370   exploreParametersSpace(input, measurements, opt, output);
0371   const unsigned nl =
0372       std::accumulate(output.begin(), output.end(), 0,
0373                       [](unsigned sum, const AccumulatorSection &s) {
0374                         return sum + s.count();
0375                       });
0376   ACTS_DEBUG("Z - cotTheta scan finds " << output.size()
0377                                         << " included measurements " << nl);
0378 }
0379
0380 void AdaptiveHoughTransformSeeder::processStackZCotThetaSplit(
0381     std::vector<AccumulatorSection> &input,
0382     std::vector<AccumulatorSection> &output,
0383     const std::vector<PreprocessedMeasurement> &measurements) const {
0384   ExplorationOptions opt;
0385   opt.xMinBinSize = 101.0f * Acts::UnitConstants::mm;
0386   opt.yMinBinSize = 10.1f;
0387   opt.lineParamFunctor = m_zCotThetaLineParams;
0388   opt.decisionFunctor = [&cfg = m_cfg, &opt](
0389                             const AccumulatorSection &section,
0390                             const std::vector<PreprocessedMeasurement> &) {
0391     using enum ExplorationOptions<PreprocessedMeasurement>::Decision;
0392     if (section.count() < cfg.threshold) {
0393       return Discard;
0394     }
0395     if (section.count() >= cfg.threshold &&
0396         section.xSize() <= opt.xMinBinSize &&
0397         section.ySize() <= opt.yMinBinSize) {
0398       return Accept;
0399     }
0400     return Drill;
0401   };
0402
0403   exploreParametersSpace(input, measurements, opt, output);
0404   const unsigned nl =
0405       std::accumulate(output.begin(), output.end(), 0,
0406                       [](unsigned sum, const AccumulatorSection &s) {
0407                         return sum + s.count();
0408                       });
0409   ACTS_DEBUG("Z - cotTheta split scan finds "
0410              << output.size() << " included measurements " << nl);
0411 }
0412
0413 void AdaptiveHoughTransformSeeder::makeSeeds(
0414     SimSeedContainer &seeds, const std::vector<AccumulatorSection> &solutions,
0415     const std::vector<PreprocessedMeasurement> &measurements) const {
0416   std::size_t seedIndex = 0;
0417   for (const AccumulatorSection &s : solutions) {
0418     unsigned spIndex = 0;
0419     std::array<const SimSpacePoint *, 3> sp = {nullptr, nullptr, nullptr};
0420     std::vector<unsigned> sortedIndices = s.indices();
0421     std::ranges::sort(sortedIndices, [&measurements](unsigned i1, unsigned i2) {
0422       const auto &m1 = measurements[i1];
0423       const auto &m2 = measurements[i2];
0424       return m1.invr > m2.invr;
0425     });
0426
0427     for (unsigned sidx : sortedIndices) {
0428       const PreprocessedMeasurement &m = measurements[sidx];
0429       sp[spIndex] = m.link.get<const SimSpacePoint *>();
0430       if (spIndex == 0 || std::abs(sp[spIndex]->r() - sp[spIndex - 1]->r()) >
0431                               5. * Acts::UnitConstants::mm) {
0432         spIndex++;
0433       }
0434       if (spIndex >= sp.size()) {
0435         break;
0436       }
0437     }
0438     if (sp[2] == nullptr) {
0439       ACTS_VERBOSE("this solution has less than 3 SP, ignoring");
0440       continue;
0441     }
0442
0443     auto cotThetaEstimate = static_cast<float>((sp[2]->z() - sp[0]->z()) /
0444                                                (sp[2]->r() - sp[0]->r()));
0445     auto cotThetaEstimate01 = static_cast<float>((sp[1]->z() - sp[0]->z()) /
0446                                                  (sp[1]->r() - sp[0]->r()));
0447     auto cotThetaEstimate12 = static_cast<float>((sp[2]->z() - sp[1]->z()) /
0448                                                  (sp[2]->r() - sp[1]->r()));
0449     if (std::abs(cotThetaEstimate01 - cotThetaEstimate12) > 0.1f) {
0450       ACTS_VERBOSE("Large difference in cotTheta " << cotThetaEstimate01 << " "
0451                                                    << cotThetaEstimate12);
0452       continue;
0453     }
0454     seeds.emplace_back(*sp[0], *sp[1], *sp[2]);
0455     auto z = static_cast<float>(sp[1]->z() - sp[1]->r() * cotThetaEstimate);
0456     seeds.back().setVertexZ(z);
0457
0458     ACTS_VERBOSE(seedIndex << ": solution x: " << s.xBegin() << " " << s.xSize()
0459                            << " y: " << s.yBegin() << " " << s.ySize()
0460                            << " nlines: " << s.count() << " vertex z: " << z
0461                            << " r0,1,2: " << sp[0]->r() << " " << sp[1]->r()
0462                            << " " << sp[2]->r() << " cot,01,12 "
0463                            << cotThetaEstimate << " " << cotThetaEstimate01
0464                            << " " << cotThetaEstimate12);
0465
0466     // for the time the quality is fixed
0467     // in the future we can use section count for instance
0468     seeds.back().setQuality(1.0);
0469     seedIndex++;
0470   }
0471 }
0472
0473 bool AdaptiveHoughTransformSeeder::passIntersectionsCheck(
0474     const AccumulatorSection &section,
0475     const std::vector<PreprocessedMeasurement> &measurements,
0476     const LineParamFunctor &lineFunctor, unsigned threshold) const {
0477   unsigned inside = 0;
0478   for (std::size_t idx1 = 0; idx1 < section.count(); ++idx1) {
0479     const auto &m1 = measurements[section.indices()[idx1]];
0480     std::function<float(float)> line1 =
0481         std::bind_front(lineFunctor, std::cref(m1));
0482     for (std::size_t idx2 = idx1 + 1; idx2 < section.count(); ++idx2) {
0483       const auto &m2 = measurements[section.indices()[idx2]];
0484       std::function<float(float)> line2 =
0485           std::bind_front(lineFunctor, std::cref(m2));
0486       if (section.isCrossingInside(line1, line2)) {
0487         inside++;
0488         if (inside >= threshold) {
0489           return true;
0490         }
0491       }
0492     }
0493   }
0494   ACTS_VERBOSE("Number of crossings inside of section " << inside);
0495   return inside >= threshold;
0496 }
0497 void AdaptiveHoughTransformSeeder::deduplicate(
0498     std::vector<AccumulatorSection> &input) const {
0499   std::vector<const AccumulatorSection *> op;
0500   op.reserve(input.size());
0501   std::transform(input.begin(), input.end(), std::back_inserter(op),
0502                  [](const AccumulatorSection &s) { return &s; });
0503
0504   auto binaryPredSort = [](const AccumulatorSection *a,
0505                            const AccumulatorSection *b) {
0506     return a->indices() < b->indices();
0507   };
0508   auto binaryPredUnique = [](const AccumulatorSection *a,
0509                              const AccumulatorSection *b) {
0510     return a->indices() == b->indices();
0511   };
0512
0513   std::ranges::sort(op, binaryPredSort);
0514   auto [rbegin, rend] = std::ranges::unique(op, binaryPredUnique);
0515   std::vector<AccumulatorSection> temp;
0516   for (auto sPtr = std::begin(op); sPtr != rbegin; ++sPtr) {
0517     temp.push_back(**sPtr);
0518   }
0519   input = std::move(temp);
0520 }
0521
0522 }  // namespace ActsExamples