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