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