TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp

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