TrackFindingGnn/src/TrackFindingFromProtoTracksAlgorithm.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/TrackFindingGnn/TrackFindingFromProtoTracksAlgorithm.hpp"
0010
0011 #include "Acts/EventData/ProxyAccessor.hpp"
0012 #include "Acts/TrackFinding/TrackStateCreator.hpp"
0013 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0014 #include "ActsExamples/EventData/MeasurementCalibration.hpp"
0015 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0016
0017 #include <algorithm>
0018 #include <ranges>
0019
0020 #include <boost/accumulators/accumulators.hpp>
0021 #include <boost/accumulators/statistics.hpp>
0022
0023 using namespace Acts;
0024
0025 namespace {
0026
0027 using namespace ActsExamples;
0028
0029 struct ProtoTrackSourceLinkAccessor
0030     : GeometryIdMultisetAccessor<IndexSourceLink> {
0031   using BaseIterator = GeometryIdMultisetAccessor<IndexSourceLink>::Iterator;
0032   using Iterator = SourceLinkAdapterIterator<BaseIterator>;
0033
0034   std::unique_ptr<const Logger> loggerPtr;
0035   Container protoTrackSourceLinks;
0036
0037   // get the range of elements with requested geoId
0038   std::pair<Iterator, Iterator> range(const Surface& surface) const {
0039     const auto& logger = *loggerPtr;
0040
0041     if (protoTrackSourceLinks.contains(surface.geometryId())) {
0042       auto [begin, end] =
0043           protoTrackSourceLinks.equal_range(surface.geometryId());
0044       ACTS_VERBOSE("Select " << std::distance(begin, end)
0045                              << " source-links from proto track on "
0046                              << surface.geometryId());
0047       return {Iterator{begin}, Iterator{end}};
0048     }
0049
0050     assert(container != nullptr);
0051     auto [begin, end] = container->equal_range(surface.geometryId());
0052     ACTS_VERBOSE("Select " << std::distance(begin, end)
0053                            << " source-links from collection on "
0054                            << surface.geometryId());
0055     return {Iterator{begin}, Iterator{end}};
0056   }
0057 };
0058
0059 }  // namespace
0060
0061 namespace ActsExamples {
0062
0063 TrackFindingFromProtoTracksAlgorithm::TrackFindingFromProtoTracksAlgorithm(
0064     Config cfg, std::unique_ptr<const Acts::Logger> logger)
0065     : IAlgorithm(cfg.tag + "CkfFromProtoTracks", std::move(logger)),
0066       m_cfg(cfg) {
0067   m_inputInitialTrackParameters.initialize(m_cfg.inputInitialTrackParameters);
0068   m_inputMeasurements.initialize(m_cfg.inputMeasurements);
0069   m_inputProtoTracks.initialize(m_cfg.inputProtoTracks);
0070   m_outputTracks.initialize(m_cfg.outputTracks);
0071 }
0072
0073 ProcessCode TrackFindingFromProtoTracksAlgorithm::execute(
0074     const AlgorithmContext& ctx) const {
0075   const auto& measurements = m_inputMeasurements(ctx);
0076   const auto& protoTracks = m_inputProtoTracks(ctx);
0077   const auto& initialParameters = m_inputInitialTrackParameters(ctx);
0078
0079   if (initialParameters.size() != protoTracks.size()) {
0080     ACTS_FATAL("Inconsistent number of parameters and proto tracks");
0081     return ProcessCode::ABORT;
0082   }
0083
0084   // Construct a perigee surface as the target surface
0085   auto pSurface = Surface::makeShared<PerigeeSurface>(Vector3{0., 0., 0.});
0086
0087   PropagatorPlainOptions pOptions(ctx.geoContext, ctx.magFieldContext);
0088   pOptions.maxSteps = 10000;
0089
0090   PassThroughCalibrator pcalibrator;
0091   MeasurementCalibratorAdapter calibrator(pcalibrator, measurements);
0092   GainMatrixUpdater kfUpdater;
0093   GainMatrixSmoother kfSmoother;
0094   MeasurementSelector measSel{m_cfg.measurementSelectorCfg};
0095
0096   // The source link accessor
0097   ProtoTrackSourceLinkAccessor sourceLinkAccessor;
0098   sourceLinkAccessor.loggerPtr = logger().clone("SourceLinkAccessor");
0099   sourceLinkAccessor.container = &measurements.orderedIndices();
0100
0101   using TrackStateCreatorType =
0102       TrackStateCreator<IndexSourceLinkAccessor::Iterator, TrackContainer>;
0103   TrackStateCreatorType trackStateCreator;
0104   trackStateCreator.sourceLinkAccessor
0105       .template connect<&ProtoTrackSourceLinkAccessor::range>(
0106           &sourceLinkAccessor);
0107   trackStateCreator.calibrator
0108       .connect<&MeasurementCalibratorAdapter::calibrate>(&calibrator);
0109   trackStateCreator.measurementSelector.connect<&MeasurementSelector::select<
0110       typename TrackContainer::TrackStateContainerBackend>>(&measSel);
0111
0112   CombinatorialKalmanFilterExtensions<TrackContainer> extensions;
0113   extensions.updater.connect<&GainMatrixUpdater::operator()<
0114       typename TrackContainer::TrackStateContainerBackend>>(&kfUpdater);
0115   extensions.createTrackStates
0116       .template connect<&TrackStateCreatorType ::createTrackStates>(
0117           &trackStateCreator);
0118
0119   // Set the CombinatorialKalmanFilter options
0120   TrackFindingAlgorithm::TrackFinderOptions options(
0121       ctx.geoContext, ctx.magFieldContext, ctx.calibContext, extensions,
0122       pOptions, &(*pSurface));
0123
0124   // Perform the track finding for all initial parameters
0125   ACTS_DEBUG("Invoke track finding with " << initialParameters.size()
0126                                           << " seeds.");
0127
0128   auto trackContainer = std::make_shared<VectorTrackContainer>();
0129   auto trackStateContainer = std::make_shared<VectorMultiTrajectory>();
0130
0131   TrackContainer tracks(trackContainer, trackStateContainer);
0132
0133   tracks.addColumn<unsigned int>("trackGroup");
0134   ProxyAccessor<unsigned int> seedNumber("trackGroup");
0135
0136   std::size_t nSeed = 0;
0137   std::size_t nFailed = 0;
0138
0139   std::vector<std::size_t> nTracksPerSeeds;
0140   nTracksPerSeeds.reserve(initialParameters.size());
0141
0142   for (auto i = 0ul; i < initialParameters.size(); ++i) {
0143     sourceLinkAccessor.protoTrackSourceLinks.clear();
0144
0145     // Fill the source links via their indices from the container
0146     for (const auto hitIndex : protoTracks.at(i)) {
0147       if (auto it = measurements.orderedIndices().nth(hitIndex);
0148           it != measurements.orderedIndices().end()) {
0149         sourceLinkAccessor.protoTrackSourceLinks.insert(*it);
0150       } else {
0151         ACTS_FATAL("Proto track " << i << " contains invalid hit index"
0152                                   << hitIndex);
0153         return ProcessCode::ABORT;
0154       }
0155     }
0156
0157     auto rootBranch = tracks.makeTrack();
0158     auto result = (*m_cfg.findTracks)(initialParameters.at(i), options, tracks,
0159                                       rootBranch);
0160     nSeed++;
0161
0162     if (!result.ok()) {
0163       nFailed++;
0164       ACTS_WARNING("Track finding failed for proto track " << i << " with error"
0165                                                            << result.error());
0166       continue;
0167     }
0168
0169     auto& tracksForSeed = result.value();
0170
0171     nTracksPerSeeds.push_back(tracksForSeed.size());
0172
0173     for (auto& track : tracksForSeed) {
0174       // Set the seed number, this number decrease by 1 since the seed number
0175       // has already been updated
0176       seedNumber(track) = nSeed - 1;
0177     }
0178   }
0179
0180   {
0181     std::lock_guard<std::mutex> guard(m_mutex);
0182
0183     std::copy(nTracksPerSeeds.begin(), nTracksPerSeeds.end(),
0184               std::back_inserter(m_nTracksPerSeeds));
0185   }
0186
0187   // TODO The computeSharedHits function is still a member function of
0188   // TrackFindingAlgorithm, but could also be a free function. Uncomment this
0189   // once this is done.
0190   // Compute shared hits from all the reconstructed tracks if
0191   // (m_cfg.computeSharedHits) {
0192   //   computeSharedHits(measurements, tracks);
0193   // }
0194
0195   ACTS_INFO("Event " << ctx.eventNumber << ": " << nFailed << " / " << nSeed
0196                      << " failed (" << ((100.f * nFailed) / nSeed) << "%)");
0197   ACTS_DEBUG("Finalized track finding with " << tracks.size()
0198                                              << " track candidates.");
0199   auto constTrackStateContainer = std::make_shared<ConstVectorMultiTrajectory>(
0200       std::move(*trackStateContainer));
0201
0202   auto constTrackContainer =
0203       std::make_shared<ConstVectorTrackContainer>(std::move(*trackContainer));
0204
0205   ConstTrackContainer constTracks{constTrackContainer,
0206                                   constTrackStateContainer};
0207
0208   m_outputTracks(ctx, std::move(constTracks));
0209   return ProcessCode::SUCCESS;
0210 }
0211
0212 ProcessCode TrackFindingFromProtoTracksAlgorithm::finalize() {
0213   assert(std::distance(m_nTracksPerSeeds.begin(), m_nTracksPerSeeds.end()) > 0);
0214
0215   ACTS_INFO("TrackFindingFromProtoTracksAlgorithm statistics:");
0216   namespace ba = boost::accumulators;
0217   using Accumulator = ba::accumulator_set<
0218       float, ba::features<ba::tag::sum, ba::tag::mean, ba::tag::variance>>;
0219
0220   Accumulator totalAcc;
0221   std::ranges::for_each(m_nTracksPerSeeds,
0222                         [&](auto v) { totalAcc(static_cast<float>(v)); });
0223   ACTS_INFO("- total number tracks: " << ba::sum(totalAcc));
0224   ACTS_INFO("- avg tracks per seed: " << ba::mean(totalAcc) << " +- "
0225                                       << std::sqrt(ba::variance(totalAcc)));
0226
0227   return {};
0228 }
0229
0230 }  // namespace ActsExamples