TrackFinding/src/TrackFindingAlgorithm.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/TrackFindingAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Direction.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/EventData/MultiTrajectory.hpp"
0015 #include "Acts/EventData/ProxyAccessor.hpp"
0016 #include "Acts/EventData/SourceLink.hpp"
0017 #include "Acts/EventData/TrackContainer.hpp"
0018 #include "Acts/EventData/TrackParameters.hpp"
0019 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0020 #include "Acts/EventData/VectorTrackContainer.hpp"
0021 #include "Acts/Geometry/GeometryIdentifier.hpp"
0022 #include "Acts/Propagator/MaterialInteractor.hpp"
0023 #include "Acts/Propagator/Navigator.hpp"
0024 #include "Acts/Propagator/Propagator.hpp"
0025 #include "Acts/Propagator/StandardAborters.hpp"
0026 #include "Acts/Propagator/SympyStepper.hpp"
0027 #include "Acts/Surfaces/PerigeeSurface.hpp"
0028 #include "Acts/Surfaces/Surface.hpp"
0029 #include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
0030 #include "Acts/TrackFinding/TrackStateCreator.hpp"
0031 #include "Acts/TrackFitting/GainMatrixUpdater.hpp"
0032 #include "Acts/Utilities/Enumerate.hpp"
0033 #include "Acts/Utilities/Logger.hpp"
0034 #include "Acts/Utilities/TrackHelpers.hpp"
0035 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0036 #include "ActsExamples/EventData/Measurement.hpp"
0037 #include "ActsExamples/EventData/MeasurementCalibration.hpp"
0038 #include "ActsExamples/EventData/SimSeed.hpp"
0039 #include "ActsExamples/EventData/Track.hpp"
0040 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0041 #include "ActsExamples/Framework/ProcessCode.hpp"
0042
0043 #include <cmath>
0044 #include <functional>
0045 #include <memory>
0046 #include <optional>
0047 #include <ostream>
0048 #include <stdexcept>
0049 #include <system_error>
0050 #include <unordered_map>
0051 #include <utility>
0052
0053 #include <boost/functional/hash.hpp>
0054
0055 // Specialize std::hash for SeedIdentifier
0056 // This is required to use SeedIdentifier as a key in an `std::unordered_map`.
0057 template <class T, std::size_t N>
0058 struct std::hash<std::array<T, N>> {
0059   std::size_t operator()(const std::array<T, N>& array) const {
0060     std::hash<T> hasher;
0061     std::size_t result = 0;
0062     for (auto&& element : array) {
0063       boost::hash_combine(result, hasher(element));
0064     }
0065     return result;
0066   }
0067 };
0068
0069 namespace ActsExamples {
0070
0071 namespace {
0072
0073 class MeasurementSelector {
0074  public:
0075   using Traj = Acts::VectorMultiTrajectory;
0076
0077   explicit MeasurementSelector(Acts::MeasurementSelector selector)
0078       : m_selector(std::move(selector)) {}
0079
0080   void setSeed(const std::optional<SimSeed>& seed) { m_seed = seed; }
0081
0082   Acts::Result<std::pair<std::vector<Traj::TrackStateProxy>::iterator,
0083                          std::vector<Traj::TrackStateProxy>::iterator>>
0084   select(std::vector<Traj::TrackStateProxy>& candidates, bool& isOutlier,
0085          const Acts::Logger& logger) const {
0086     if (m_seed.has_value()) {
0087       std::vector<Traj::TrackStateProxy> newCandidates;
0088
0089       for (const auto& candidate : candidates) {
0090         if (isSeedCandidate(candidate)) {
0091           newCandidates.push_back(candidate);
0092         }
0093       }
0094
0095       if (!newCandidates.empty()) {
0096         candidates = std::move(newCandidates);
0097       }
0098     }
0099
0100     return m_selector.select<Acts::VectorMultiTrajectory>(candidates, isOutlier,
0101                                                           logger);
0102   }
0103
0104  private:
0105   Acts::MeasurementSelector m_selector;
0106   std::optional<SimSeed> m_seed;
0107
0108   bool isSeedCandidate(const Traj::TrackStateProxy& candidate) const {
0109     assert(candidate.hasUncalibratedSourceLink());
0110
0111     const Acts::SourceLink& sourceLink = candidate.getUncalibratedSourceLink();
0112
0113     for (const auto& sp : m_seed->sp()) {
0114       for (const auto& sl : sp->sourceLinks()) {
0115         if (sourceLink.get<IndexSourceLink>() == sl.get<IndexSourceLink>()) {
0116           return true;
0117         }
0118       }
0119     }
0120
0121     return false;
0122   }
0123 };
0124
0125 /// Source link indices of the bottom, middle, top measurements.
0126 /// In case of strip seeds only the first source link of the pair is used.
0127 using SeedIdentifier = std::array<Index, 3>;
0128
0129 /// Build a seed identifier from a seed.
0130 ///
0131 /// @param seed The seed to build the identifier from.
0132 /// @return The seed identifier.
0133 SeedIdentifier makeSeedIdentifier(const SimSeed& seed) {
0134   SeedIdentifier result;
0135
0136   for (const auto& [i, sp] : Acts::enumerate(seed.sp())) {
0137     const Acts::SourceLink& firstSourceLink = sp->sourceLinks().front();
0138     result.at(i) = firstSourceLink.get<IndexSourceLink>().index();
0139   }
0140
0141   return result;
0142 }
0143
0144 /// Visit all possible seed identifiers of a track.
0145 ///
0146 /// @param track The track to visit the seed identifiers of.
0147 /// @param visitor The visitor to call for each seed identifier.
0148 template <typename Visitor>
0149 void visitSeedIdentifiers(const TrackProxy& track, Visitor visitor) {
0150   // first we collect the source link indices of the track states
0151   std::vector<Index> sourceLinkIndices;
0152   sourceLinkIndices.reserve(track.nMeasurements());
0153   for (const auto& trackState : track.trackStatesReversed()) {
0154     if (!trackState.hasUncalibratedSourceLink()) {
0155       continue;
0156     }
0157     const Acts::SourceLink& sourceLink = trackState.getUncalibratedSourceLink();
0158     sourceLinkIndices.push_back(sourceLink.get<IndexSourceLink>().index());
0159   }
0160
0161   // then we iterate over all possible triplets and form seed identifiers
0162   for (std::size_t i = 0; i < sourceLinkIndices.size(); ++i) {
0163     for (std::size_t j = i + 1; j < sourceLinkIndices.size(); ++j) {
0164       for (std::size_t k = j + 1; k < sourceLinkIndices.size(); ++k) {
0165         // Putting them into reverse order (k, j, i) to compensate for the
0166         // `trackStatesReversed` above.
0167         visitor({sourceLinkIndices.at(k), sourceLinkIndices.at(j),
0168                  sourceLinkIndices.at(i)});
0169       }
0170     }
0171   }
0172 }
0173
0174 class BranchStopper {
0175  public:
0176   using BranchStopperResult =
0177       Acts::CombinatorialKalmanFilterBranchStopperResult;
0178
0179   struct BranchState {
0180     std::size_t nPixelHoles = 0;
0181     std::size_t nStripHoles = 0;
0182   };
0183
0184   static constexpr Acts::ProxyAccessor<BranchState> branchStateAccessor =
0185       Acts::ProxyAccessor<BranchState>(Acts::hashString("MyBranchState"));
0186
0187   mutable std::atomic<std::size_t> m_nStoppedBranches{0};
0188
0189   explicit BranchStopper(const TrackFindingAlgorithm::Config& config)
0190       : m_cfg(config) {}
0191
0192   BranchStopperResult operator()(
0193       const TrackContainer::TrackProxy& track,
0194       const TrackContainer::TrackStateProxy& trackState) const {
0195     if (!m_cfg.trackSelectorCfg.has_value()) {
0196       return BranchStopperResult::Continue;
0197     }
0198
0199     const Acts::TrackSelector::Config* singleConfig = std::visit(
0200         [&](const auto& config) -> const Acts::TrackSelector::Config* {
0201           using T = std::decay_t<decltype(config)>;
0202           if constexpr (std::is_same_v<T, Acts::TrackSelector::Config>) {
0203             return &config;
0204           } else if constexpr (std::is_same_v<
0205                                    T, Acts::TrackSelector::EtaBinnedConfig>) {
0206             double theta = trackState.parameters()[Acts::eBoundTheta];
0207             double eta = Acts::AngleHelpers::etaFromTheta(theta);
0208             return config.hasCuts(eta) ? &config.getCuts(eta) : nullptr;
0209           }
0210         },
0211         *m_cfg.trackSelectorCfg);
0212
0213     if (singleConfig == nullptr) {
0214       ++m_nStoppedBranches;
0215       return BranchStopperResult::StopAndDrop;
0216     }
0217
0218     bool tooManyHolesPS = false;
0219     if (!(m_cfg.pixelVolumeIds.empty() && m_cfg.stripVolumeIds.empty())) {
0220       auto& branchState = branchStateAccessor(track);
0221       // count both holes and outliers as holes for pixel/strip counts
0222       if (trackState.typeFlags().test(Acts::TrackStateFlag::HoleFlag) ||
0223           trackState.typeFlags().test(Acts::TrackStateFlag::OutlierFlag)) {
0224         auto volumeId = trackState.referenceSurface().geometryId().volume();
0225         if (std::find(m_cfg.pixelVolumeIds.begin(), m_cfg.pixelVolumeIds.end(),
0226                       volumeId) != m_cfg.pixelVolumeIds.end()) {
0227           ++branchState.nPixelHoles;
0228         } else if (std::find(m_cfg.stripVolumeIds.begin(),
0229                              m_cfg.stripVolumeIds.end(),
0230                              volumeId) != m_cfg.stripVolumeIds.end()) {
0231           ++branchState.nStripHoles;
0232         }
0233       }
0234       tooManyHolesPS = branchState.nPixelHoles > m_cfg.maxPixelHoles ||
0235                        branchState.nStripHoles > m_cfg.maxStripHoles;
0236     }
0237
0238     bool enoughMeasurements =
0239         track.nMeasurements() >= singleConfig->minMeasurements;
0240     bool tooManyHoles =
0241         track.nHoles() > singleConfig->maxHoles || tooManyHolesPS;
0242     bool tooManyOutliers = track.nOutliers() > singleConfig->maxOutliers;
0243     bool tooManyHolesAndOutliers = (track.nHoles() + track.nOutliers()) >
0244                                    singleConfig->maxHolesAndOutliers;
0245
0246     if (tooManyHoles || tooManyOutliers || tooManyHolesAndOutliers) {
0247       ++m_nStoppedBranches;
0248       return enoughMeasurements ? BranchStopperResult::StopAndKeep
0249                                 : BranchStopperResult::StopAndDrop;
0250     }
0251
0252     return BranchStopperResult::Continue;
0253   }
0254
0255  private:
0256   const TrackFindingAlgorithm::Config& m_cfg;
0257 };
0258
0259 }  // namespace
0260
0261 TrackFindingAlgorithm::TrackFindingAlgorithm(Config config,
0262                                              Acts::Logging::Level level)
0263     : IAlgorithm("TrackFindingAlgorithm", level), m_cfg(std::move(config)) {
0264   if (m_cfg.inputMeasurements.empty()) {
0265     throw std::invalid_argument("Missing measurements input collection");
0266   }
0267   if (m_cfg.inputInitialTrackParameters.empty()) {
0268     throw std::invalid_argument(
0269         "Missing initial track parameters input collection");
0270   }
0271   if (m_cfg.outputTracks.empty()) {
0272     throw std::invalid_argument("Missing tracks output collection");
0273   }
0274
0275   if (m_cfg.seedDeduplication && m_cfg.inputSeeds.empty()) {
0276     throw std::invalid_argument(
0277         "Missing seeds input collection. This is "
0278         "required for seed deduplication.");
0279   }
0280   if (m_cfg.stayOnSeed && m_cfg.inputSeeds.empty()) {
0281     throw std::invalid_argument(
0282         "Missing seeds input collection. This is "
0283         "required for staying on seed.");
0284   }
0285
0286   if (m_cfg.trackSelectorCfg.has_value()) {
0287     m_trackSelector = std::visit(
0288         [](const auto& cfg) -> std::optional<Acts::TrackSelector> {
0289           return Acts::TrackSelector(cfg);
0290         },
0291         m_cfg.trackSelectorCfg.value());
0292   }
0293
0294   m_inputMeasurements.initialize(m_cfg.inputMeasurements);
0295   m_inputInitialTrackParameters.initialize(m_cfg.inputInitialTrackParameters);
0296   m_inputSeeds.maybeInitialize(m_cfg.inputSeeds);
0297   m_outputTracks.initialize(m_cfg.outputTracks);
0298 }
0299
0300 ProcessCode TrackFindingAlgorithm::execute(const AlgorithmContext& ctx) const {
0301   // Read input data
0302   const auto& measurements = m_inputMeasurements(ctx);
0303   const auto& initialParameters = m_inputInitialTrackParameters(ctx);
0304   const SimSeedContainer* seeds = nullptr;
0305
0306   if (m_inputSeeds.isInitialized()) {
0307     seeds = &m_inputSeeds(ctx);
0308
0309     if (initialParameters.size() != seeds->size()) {
0310       ACTS_ERROR("Number of initial parameters and seeds do not match. "
0311                  << initialParameters.size() << " != " << seeds->size());
0312     }
0313   }
0314
0315   // Construct a perigee surface as the target surface
0316   auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
0317       Acts::Vector3{0., 0., 0.});
0318
0319   PassThroughCalibrator pcalibrator;
0320   MeasurementCalibratorAdapter calibrator(pcalibrator, measurements);
0321   Acts::GainMatrixUpdater kfUpdater;
0322
0323   using Extensions = Acts::CombinatorialKalmanFilterExtensions<TrackContainer>;
0324
0325   BranchStopper branchStopper(m_cfg);
0326   MeasurementSelector measSel{
0327       Acts::MeasurementSelector(m_cfg.measurementSelectorCfg)};
0328
0329   IndexSourceLinkAccessor slAccessor;
0330   slAccessor.container = &measurements.orderedIndices();
0331
0332   using TrackStateCreatorType =
0333       Acts::TrackStateCreator<IndexSourceLinkAccessor::Iterator,
0334                               TrackContainer>;
0335   TrackStateCreatorType trackStateCreator;
0336   trackStateCreator.sourceLinkAccessor
0337       .template connect<&IndexSourceLinkAccessor::range>(&slAccessor);
0338   trackStateCreator.calibrator
0339       .template connect<&MeasurementCalibratorAdapter::calibrate>(&calibrator);
0340   trackStateCreator.measurementSelector
0341       .template connect<&MeasurementSelector::select>(&measSel);
0342
0343   Extensions extensions;
0344   extensions.updater.connect<&Acts::GainMatrixUpdater::operator()<
0345       typename TrackContainer::TrackStateContainerBackend>>(&kfUpdater);
0346   extensions.branchStopper.connect<&BranchStopper::operator()>(&branchStopper);
0347   extensions.createTrackStates
0348       .template connect<&TrackStateCreatorType ::createTrackStates>(
0349           &trackStateCreator);
0350
0351   Acts::PropagatorPlainOptions firstPropOptions(ctx.geoContext,
0352                                                 ctx.magFieldContext);
0353   firstPropOptions.maxSteps = m_cfg.maxSteps;
0354   firstPropOptions.direction = m_cfg.reverseSearch ? Acts::Direction::Backward()
0355                                                    : Acts::Direction::Forward();
0356   firstPropOptions.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
0357   firstPropOptions.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
0358
0359   Acts::PropagatorPlainOptions secondPropOptions(ctx.geoContext,
0360                                                  ctx.magFieldContext);
0361   secondPropOptions.maxSteps = m_cfg.maxSteps;
0362   secondPropOptions.direction = firstPropOptions.direction.invert();
0363   secondPropOptions.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
0364   secondPropOptions.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
0365
0366   // Set the CombinatorialKalmanFilter options
0367   TrackFinderOptions firstOptions(ctx.geoContext, ctx.magFieldContext,
0368                                   ctx.calibContext, extensions,
0369                                   firstPropOptions);
0370
0371   firstOptions.targetSurface = m_cfg.reverseSearch ? pSurface.get() : nullptr;
0372
0373   TrackFinderOptions secondOptions(ctx.geoContext, ctx.magFieldContext,
0374                                    ctx.calibContext, extensions,
0375                                    secondPropOptions);
0376   secondOptions.targetSurface = m_cfg.reverseSearch ? nullptr : pSurface.get();
0377   secondOptions.skipPrePropagationUpdate = true;
0378
0379   using Extrapolator = Acts::Propagator<Acts::SympyStepper, Acts::Navigator>;
0380   using ExtrapolatorOptions = Extrapolator::template Options<
0381       Acts::ActorList<Acts::MaterialInteractor, Acts::EndOfWorldReached>>;
0382
0383   Extrapolator extrapolator(
0384       Acts::SympyStepper(m_cfg.magneticField),
0385       Acts::Navigator({m_cfg.trackingGeometry},
0386                       logger().cloneWithSuffix("Navigator")),
0387       logger().cloneWithSuffix("Propagator"));
0388
0389   ExtrapolatorOptions extrapolationOptions(ctx.geoContext, ctx.magFieldContext);
0390   extrapolationOptions.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
0391   extrapolationOptions.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
0392
0393   // Perform the track finding for all initial parameters
0394   ACTS_DEBUG("Invoke track finding with " << initialParameters.size()
0395                                           << " seeds.");
0396
0397   auto trackContainer = std::make_shared<Acts::VectorTrackContainer>();
0398   auto trackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
0399
0400   auto trackContainerTemp = std::make_shared<Acts::VectorTrackContainer>();
0401   auto trackStateContainerTemp =
0402       std::make_shared<Acts::VectorMultiTrajectory>();
0403
0404   TrackContainer tracks(trackContainer, trackStateContainer);
0405   TrackContainer tracksTemp(trackContainerTemp, trackStateContainerTemp);
0406
0407   // Note that not all backends support PODs as column types
0408   tracks.addColumn<BranchStopper::BranchState>("MyBranchState");
0409   tracksTemp.addColumn<BranchStopper::BranchState>("MyBranchState");
0410
0411   tracks.addColumn<unsigned int>("trackGroup");
0412   tracksTemp.addColumn<unsigned int>("trackGroup");
0413   Acts::ProxyAccessor<unsigned int> seedNumber("trackGroup");
0414
0415   unsigned int nSeed = 0;
0416
0417   // A map indicating whether a seed has been discovered already
0418   std::unordered_map<SeedIdentifier, bool> discoveredSeeds;
0419
0420   auto addTrack = [&](const TrackProxy& track) {
0421     ++m_nFoundTracks;
0422
0423     // trim the track if requested
0424     if (m_cfg.trimTracks) {
0425       Acts::trimTrack(track, true, true, true, true);
0426     }
0427     Acts::calculateTrackQuantities(track);
0428
0429     if (m_trackSelector.has_value() && !m_trackSelector->isValidTrack(track)) {
0430       return;
0431     }
0432
0433     // flag seeds which are covered by the track
0434     visitSeedIdentifiers(track, [&](const SeedIdentifier& seedIdentifier) {
0435       if (auto it = discoveredSeeds.find(seedIdentifier);
0436           it != discoveredSeeds.end()) {
0437         it->second = true;
0438       }
0439     });
0440
0441     ++m_nSelectedTracks;
0442
0443     auto destProxy = tracks.makeTrack();
0444     // make sure we copy track states!
0445     destProxy.copyFrom(track, true);
0446   };
0447
0448   if (seeds != nullptr && m_cfg.seedDeduplication) {
0449     // Index the seeds for deduplication
0450     for (const auto& seed : *seeds) {
0451       SeedIdentifier seedIdentifier = makeSeedIdentifier(seed);
0452       discoveredSeeds.emplace(seedIdentifier, false);
0453     }
0454   }
0455
0456   for (std::size_t iSeed = 0; iSeed < initialParameters.size(); ++iSeed) {
0457     m_nTotalSeeds++;
0458
0459     if (seeds != nullptr) {
0460       const SimSeed& seed = seeds->at(iSeed);
0461
0462       if (m_cfg.seedDeduplication) {
0463         SeedIdentifier seedIdentifier = makeSeedIdentifier(seed);
0464         // check if the seed has been discovered already
0465         if (auto it = discoveredSeeds.find(seedIdentifier);
0466             it != discoveredSeeds.end() && it->second) {
0467           m_nDeduplicatedSeeds++;
0468           ACTS_VERBOSE("Skipping seed " << iSeed << " due to deduplication.");
0469           continue;
0470         }
0471       }
0472
0473       if (m_cfg.stayOnSeed) {
0474         measSel.setSeed(seed);
0475       }
0476     }
0477
0478     // Clear trackContainerTemp and trackStateContainerTemp
0479     tracksTemp.clear();
0480
0481     const Acts::BoundTrackParameters& firstInitialParameters =
0482         initialParameters.at(iSeed);
0483
0484     auto firstRootBranch = tracksTemp.makeTrack();
0485     auto firstResult = (*m_cfg.findTracks)(firstInitialParameters, firstOptions,
0486                                            tracksTemp, firstRootBranch);
0487     nSeed++;
0488
0489     if (!firstResult.ok()) {
0490       m_nFailedSeeds++;
0491       ACTS_WARNING("Track finding failed for seed " << iSeed << " with error"
0492                                                     << firstResult.error());
0493       continue;
0494     }
0495
0496     auto& firstTracksForSeed = firstResult.value();
0497     for (auto& firstTrack : firstTracksForSeed) {
0498       // TODO a copy of the track should not be necessary but is the safest way
0499       //      with the current EDM
0500       // TODO a lightweight copy without copying all the track state components
0501       //      might be a solution
0502       auto trackCandidate = tracksTemp.makeTrack();
0503       trackCandidate.copyFrom(firstTrack, true);
0504
0505       Acts::Result<void> firstSmoothingResult{
0506           Acts::smoothTrack(ctx.geoContext, trackCandidate, logger())};
0507       if (!firstSmoothingResult.ok()) {
0508         m_nFailedSmoothing++;
0509         ACTS_ERROR("First smoothing for seed "
0510                    << iSeed << " and track " << firstTrack.index()
0511                    << " failed with error " << firstSmoothingResult.error());
0512         continue;
0513       }
0514
0515       // number of second tracks found
0516       std::size_t nSecond = 0;
0517
0518       // Set the seed number, this number decrease by 1 since the seed number
0519       // has already been updated
0520       seedNumber(trackCandidate) = nSeed - 1;
0521
0522       if (m_cfg.twoWay) {
0523         std::optional<Acts::VectorMultiTrajectory::TrackStateProxy>
0524             firstMeasurementOpt;
0525         for (auto trackState : trackCandidate.trackStatesReversed()) {
0526           bool isMeasurement = trackState.typeFlags().test(
0527               Acts::TrackStateFlag::MeasurementFlag);
0528           bool isOutlier =
0529               trackState.typeFlags().test(Acts::TrackStateFlag::OutlierFlag);
0530           // We are excluding non measurement states and outlier here. Those can
0531           // decrease resolution because only the smoothing corrected the very
0532           // first prediction as filtering is not possible.
0533           if (isMeasurement && !isOutlier) {
0534             firstMeasurementOpt = trackState;
0535           }
0536         }
0537
0538         if (firstMeasurementOpt.has_value()) {
0539           TrackContainer::TrackStateProxy firstMeasurement{
0540               firstMeasurementOpt.value()};
0541           TrackContainer::ConstTrackStateProxy firstMeasurementConst{
0542               firstMeasurement};
0543
0544           Acts::BoundTrackParameters secondInitialParameters =
0545               trackCandidate.createParametersFromState(firstMeasurementConst);
0546
0547           if (!secondInitialParameters.referenceSurface().insideBounds(
0548                   secondInitialParameters.localPosition())) {
0549             m_nSkippedSecondPass++;
0550             ACTS_DEBUG(
0551                 "Smoothing of first pass fit produced out-of-bounds parameters "
0552                 "relative to the surface. Skipping second pass.");
0553             continue;
0554           }
0555
0556           auto secondRootBranch = tracksTemp.makeTrack();
0557           secondRootBranch.copyFrom(trackCandidate, false);
0558           auto secondResult =
0559               (*m_cfg.findTracks)(secondInitialParameters, secondOptions,
0560                                   tracksTemp, secondRootBranch);
0561
0562           if (!secondResult.ok()) {
0563             ACTS_WARNING("Second track finding failed for seed "
0564                          << iSeed << " with error" << secondResult.error());
0565           } else {
0566             // store the original previous state to restore it later
0567             auto originalFirstMeasurementPrevious = firstMeasurement.previous();
0568
0569             auto& secondTracksForSeed = secondResult.value();
0570             for (auto& secondTrack : secondTracksForSeed) {
0571               // TODO a copy of the track should not be necessary but is the
0572               //      safest way with the current EDM
0573               // TODO a lightweight copy without copying all the track state
0574               //      components might be a solution
0575               auto secondTrackCopy = tracksTemp.makeTrack();
0576               secondTrackCopy.copyFrom(secondTrack, true);
0577
0578               // Note that this is only valid if there are no branches
0579               // We disallow this by breaking this look after a second track was
0580               // processed
0581               secondTrackCopy.reverseTrackStates(true);
0582
0583               firstMeasurement.previous() =
0584                   secondTrackCopy.outermostTrackState().index();
0585
0586               trackCandidate.copyFrom(secondTrackCopy, false);
0587
0588               // finalize the track candidate
0589
0590               bool doExtrapolate = true;
0591
0592               if (!m_cfg.reverseSearch) {
0593                 // these parameters are already extrapolated by the CKF and have
0594                 // the optimal resolution. note that we did not smooth all the
0595                 // states.
0596
0597                 // only extrapolate if we did not do it already
0598                 doExtrapolate = !trackCandidate.hasReferenceSurface();
0599               } else {
0600                 // smooth the full track and extrapolate to the reference
0601
0602                 auto secondSmoothingResult =
0603                     Acts::smoothTrack(ctx.geoContext, trackCandidate, logger());
0604                 if (!secondSmoothingResult.ok()) {
0605                   m_nFailedSmoothing++;
0606                   ACTS_ERROR("Second smoothing for seed "
0607                              << iSeed << " and track " << secondTrack.index()
0608                              << " failed with error "
0609                              << secondSmoothingResult.error());
0610                   continue;
0611                 }
0612
0613                 trackCandidate.reverseTrackStates(true);
0614               }
0615
0616               if (doExtrapolate) {
0617                 auto secondExtrapolationResult =
0618                     Acts::extrapolateTrackToReferenceSurface(
0619                         trackCandidate, *pSurface, extrapolator,
0620                         extrapolationOptions, m_cfg.extrapolationStrategy,
0621                         logger());
0622                 if (!secondExtrapolationResult.ok()) {
0623                   m_nFailedExtrapolation++;
0624                   ACTS_ERROR("Second extrapolation for seed "
0625                              << iSeed << " and track " << secondTrack.index()
0626                              << " failed with error "
0627                              << secondExtrapolationResult.error());
0628                   continue;
0629                 }
0630               }
0631
0632               addTrack(trackCandidate);
0633
0634               ++nSecond;
0635             }
0636
0637             // restore the original previous state
0638             firstMeasurement.previous() = originalFirstMeasurementPrevious;
0639           }
0640         }
0641       }
0642
0643       // if no second track was found, we will use only the first track
0644       if (nSecond == 0) {
0645         // restore the track to the original state
0646         trackCandidate.copyFrom(firstTrack, false);
0647
0648         auto firstExtrapolationResult =
0649             Acts::extrapolateTrackToReferenceSurface(
0650                 trackCandidate, *pSurface, extrapolator, extrapolationOptions,
0651                 m_cfg.extrapolationStrategy, logger());
0652         if (!firstExtrapolationResult.ok()) {
0653           m_nFailedExtrapolation++;
0654           ACTS_ERROR("Extrapolation for seed "
0655                      << iSeed << " and track " << firstTrack.index()
0656                      << " failed with error "
0657                      << firstExtrapolationResult.error());
0658           continue;
0659         }
0660
0661         addTrack(trackCandidate);
0662       }
0663     }
0664   }
0665
0666   // Compute shared hits from all the reconstructed tracks
0667   if (m_cfg.computeSharedHits) {
0668     computeSharedHits(tracks, measurements);
0669   }
0670
0671   ACTS_DEBUG("Finalized track finding with " << tracks.size()
0672                                              << " track candidates.");
0673
0674   m_nStoppedBranches += branchStopper.m_nStoppedBranches;
0675
0676   m_memoryStatistics.local().hist +=
0677       tracks.trackStateContainer().statistics().hist;
0678
0679   auto constTrackStateContainer =
0680       std::make_shared<Acts::ConstVectorMultiTrajectory>(
0681           std::move(*trackStateContainer));
0682
0683   auto constTrackContainer = std::make_shared<Acts::ConstVectorTrackContainer>(
0684       std::move(*trackContainer));
0685
0686   ConstTrackContainer constTracks{constTrackContainer,
0687                                   constTrackStateContainer};
0688
0689   m_outputTracks(ctx, std::move(constTracks));
0690   return ProcessCode::SUCCESS;
0691 }
0692
0693 ProcessCode TrackFindingAlgorithm::finalize() {
0694   ACTS_INFO("TrackFindingAlgorithm statistics:");
0695   ACTS_INFO("- total seeds: " << m_nTotalSeeds);
0696   ACTS_INFO("- deduplicated seeds: " << m_nDeduplicatedSeeds);
0697   ACTS_INFO("- failed seeds: " << m_nFailedSeeds);
0698   ACTS_INFO("- failed smoothing: " << m_nFailedSmoothing);
0699   ACTS_INFO("- failed extrapolation: " << m_nFailedExtrapolation);
0700   ACTS_INFO("- failure ratio seeds: " << static_cast<double>(m_nFailedSeeds) /
0701                                              m_nTotalSeeds);
0702   ACTS_INFO("- found tracks: " << m_nFoundTracks);
0703   ACTS_INFO("- selected tracks: " << m_nSelectedTracks);
0704   ACTS_INFO("- stopped branches: " << m_nStoppedBranches);
0705   ACTS_INFO("- skipped second pass: " << m_nSkippedSecondPass);
0706
0707   auto memoryStatistics =
0708       m_memoryStatistics.combine([](const auto& a, const auto& b) {
0709         Acts::VectorMultiTrajectory::Statistics c;
0710         c.hist = a.hist + b.hist;
0711         return c;
0712       });
0713   std::stringstream ss;
0714   memoryStatistics.toStream(ss);
0715   ACTS_DEBUG("Track State memory statistics (averaged):\n" << ss.str());
0716   return ProcessCode::SUCCESS;
0717 }
0718
0719 // TODO this is somewhat duplicated in AmbiguityResolutionAlgorithm.cpp
0720 // TODO we should make a common implementation in the core at some point
0721 void TrackFindingAlgorithm::computeSharedHits(
0722     TrackContainer& tracks, const MeasurementContainer& measurements) const {
0723   // Compute shared hits from all the reconstructed tracks
0724   // Compute nSharedhits and Update ckf results
0725   // hit index -> list of multi traj indexes [traj, meas]
0726
0727   std::vector<std::size_t> firstTrackOnTheHit(
0728       measurements.size(), std::numeric_limits<std::size_t>::max());
0729   std::vector<std::size_t> firstStateOnTheHit(
0730       measurements.size(), std::numeric_limits<std::size_t>::max());
0731
0732   for (auto track : tracks) {
0733     for (auto state : track.trackStatesReversed()) {
0734       if (!state.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) {
0735         continue;
0736       }
0737
0738       std::size_t hitIndex = state.getUncalibratedSourceLink()
0739                                  .template get<IndexSourceLink>()
0740                                  .index();
0741
0742       // Check if hit not already used
0743       if (firstTrackOnTheHit.at(hitIndex) ==
0744           std::numeric_limits<std::size_t>::max()) {
0745         firstTrackOnTheHit.at(hitIndex) = track.index();
0746         firstStateOnTheHit.at(hitIndex) = state.index();
0747         continue;
0748       }
0749
0750       // if already used, control if first track state has been marked
0751       // as shared
0752       int indexFirstTrack = firstTrackOnTheHit.at(hitIndex);
0753       int indexFirstState = firstStateOnTheHit.at(hitIndex);
0754
0755       auto firstState = tracks.getTrack(indexFirstTrack)
0756                             .container()
0757                             .trackStateContainer()
0758                             .getTrackState(indexFirstState);
0759       if (!firstState.typeFlags().test(Acts::TrackStateFlag::SharedHitFlag)) {
0760         firstState.typeFlags().set(Acts::TrackStateFlag::SharedHitFlag);
0761       }
0762
0763       // Decorate this track state
0764       state.typeFlags().set(Acts::TrackStateFlag::SharedHitFlag);
0765     }
0766   }
0767 }
0768
0769 }  // namespace ActsExamples