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