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 (Acts::rangeContainsValue(m_cfg.pixelVolumeIds, volumeId)) {
0226           ++branchState.nPixelHoles;
0227         } else if (Acts::rangeContainsValue(m_cfg.stripVolumeIds, volumeId)) {
0228           ++branchState.nStripHoles;
0229         }
0230       }
0231       tooManyHolesPS = branchState.nPixelHoles > m_cfg.maxPixelHoles ||
0232                        branchState.nStripHoles > m_cfg.maxStripHoles;
0233     }
0234
0235     bool enoughMeasurements =
0236         track.nMeasurements() >= singleConfig->minMeasurements;
0237     bool tooManyHoles =
0238         track.nHoles() > singleConfig->maxHoles || tooManyHolesPS;
0239     bool tooManyOutliers = track.nOutliers() > singleConfig->maxOutliers;
0240     bool tooManyHolesAndOutliers = (track.nHoles() + track.nOutliers()) >
0241                                    singleConfig->maxHolesAndOutliers;
0242
0243     if (tooManyHoles || tooManyOutliers || tooManyHolesAndOutliers) {
0244       ++m_nStoppedBranches;
0245       return enoughMeasurements ? BranchStopperResult::StopAndKeep
0246                                 : BranchStopperResult::StopAndDrop;
0247     }
0248
0249     return BranchStopperResult::Continue;
0250   }
0251
0252  private:
0253   const TrackFindingAlgorithm::Config& m_cfg;
0254 };
0255
0256 }  // namespace
0257
0258 TrackFindingAlgorithm::TrackFindingAlgorithm(Config config,
0259                                              Acts::Logging::Level level)
0260     : IAlgorithm("TrackFindingAlgorithm", level), m_cfg(std::move(config)) {
0261   if (m_cfg.inputMeasurements.empty()) {
0262     throw std::invalid_argument("Missing measurements input collection");
0263   }
0264   if (m_cfg.inputInitialTrackParameters.empty()) {
0265     throw std::invalid_argument(
0266         "Missing initial track parameters input collection");
0267   }
0268   if (m_cfg.outputTracks.empty()) {
0269     throw std::invalid_argument("Missing tracks output collection");
0270   }
0271
0272   if (m_cfg.seedDeduplication && m_cfg.inputSeeds.empty()) {
0273     throw std::invalid_argument(
0274         "Missing seeds input collection. This is "
0275         "required for seed deduplication.");
0276   }
0277   if (m_cfg.stayOnSeed && m_cfg.inputSeeds.empty()) {
0278     throw std::invalid_argument(
0279         "Missing seeds input collection. This is "
0280         "required for staying on seed.");
0281   }
0282
0283   if (m_cfg.trackSelectorCfg.has_value()) {
0284     m_trackSelector = std::visit(
0285         [](const auto& cfg) -> std::optional<Acts::TrackSelector> {
0286           return Acts::TrackSelector(cfg);
0287         },
0288         m_cfg.trackSelectorCfg.value());
0289   }
0290
0291   m_inputMeasurements.initialize(m_cfg.inputMeasurements);
0292   m_inputInitialTrackParameters.initialize(m_cfg.inputInitialTrackParameters);
0293   m_inputSeeds.maybeInitialize(m_cfg.inputSeeds);
0294   m_outputTracks.initialize(m_cfg.outputTracks);
0295 }
0296
0297 ProcessCode TrackFindingAlgorithm::execute(const AlgorithmContext& ctx) const {
0298   // Read input data
0299   const auto& measurements = m_inputMeasurements(ctx);
0300   const auto& initialParameters = m_inputInitialTrackParameters(ctx);
0301   const SimSeedContainer* seeds = nullptr;
0302
0303   if (m_inputSeeds.isInitialized()) {
0304     seeds = &m_inputSeeds(ctx);
0305
0306     if (initialParameters.size() != seeds->size()) {
0307       ACTS_ERROR("Number of initial parameters and seeds do not match. "
0308                  << initialParameters.size() << " != " << seeds->size());
0309     }
0310   }
0311
0312   // Construct a perigee surface as the target surface
0313   auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
0314       Acts::Vector3{0., 0., 0.});
0315
0316   PassThroughCalibrator pcalibrator;
0317   MeasurementCalibratorAdapter calibrator(pcalibrator, measurements);
0318   Acts::GainMatrixUpdater kfUpdater;
0319
0320   using Extensions = Acts::CombinatorialKalmanFilterExtensions<TrackContainer>;
0321
0322   BranchStopper branchStopper(m_cfg);
0323   MeasurementSelector measSel{
0324       Acts::MeasurementSelector(m_cfg.measurementSelectorCfg)};
0325
0326   IndexSourceLinkAccessor slAccessor;
0327   slAccessor.container = &measurements.orderedIndices();
0328
0329   using TrackStateCreatorType =
0330       Acts::TrackStateCreator<IndexSourceLinkAccessor::Iterator,
0331                               TrackContainer>;
0332   TrackStateCreatorType trackStateCreator;
0333   trackStateCreator.sourceLinkAccessor
0334       .template connect<&IndexSourceLinkAccessor::range>(&slAccessor);
0335   trackStateCreator.calibrator
0336       .template connect<&MeasurementCalibratorAdapter::calibrate>(&calibrator);
0337   trackStateCreator.measurementSelector
0338       .template connect<&MeasurementSelector::select>(&measSel);
0339
0340   Extensions extensions;
0341   extensions.updater.connect<&Acts::GainMatrixUpdater::operator()<
0342       typename TrackContainer::TrackStateContainerBackend>>(&kfUpdater);
0343   extensions.branchStopper.connect<&BranchStopper::operator()>(&branchStopper);
0344   extensions.createTrackStates
0345       .template connect<&TrackStateCreatorType ::createTrackStates>(
0346           &trackStateCreator);
0347
0348   Acts::PropagatorPlainOptions firstPropOptions(ctx.geoContext,
0349                                                 ctx.magFieldContext);
0350   firstPropOptions.maxSteps = m_cfg.maxSteps;
0351   firstPropOptions.direction = m_cfg.reverseSearch ? Acts::Direction::Backward()
0352                                                    : Acts::Direction::Forward();
0353   firstPropOptions.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
0354   firstPropOptions.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
0355
0356   Acts::PropagatorPlainOptions secondPropOptions(ctx.geoContext,
0357                                                  ctx.magFieldContext);
0358   secondPropOptions.maxSteps = m_cfg.maxSteps;
0359   secondPropOptions.direction = firstPropOptions.direction.invert();
0360   secondPropOptions.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
0361   secondPropOptions.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
0362
0363   // Set the CombinatorialKalmanFilter options
0364   TrackFinderOptions firstOptions(ctx.geoContext, ctx.magFieldContext,
0365                                   ctx.calibContext, extensions,
0366                                   firstPropOptions);
0367
0368   firstOptions.targetSurface = m_cfg.reverseSearch ? pSurface.get() : nullptr;
0369
0370   TrackFinderOptions secondOptions(ctx.geoContext, ctx.magFieldContext,
0371                                    ctx.calibContext, extensions,
0372                                    secondPropOptions);
0373   secondOptions.targetSurface = m_cfg.reverseSearch ? nullptr : pSurface.get();
0374   secondOptions.skipPrePropagationUpdate = true;
0375
0376   using Extrapolator = Acts::Propagator<Acts::SympyStepper, Acts::Navigator>;
0377   using ExtrapolatorOptions = Extrapolator::template Options<
0378       Acts::ActorList<Acts::MaterialInteractor, Acts::EndOfWorldReached>>;
0379
0380   Extrapolator extrapolator(
0381       Acts::SympyStepper(m_cfg.magneticField),
0382       Acts::Navigator({m_cfg.trackingGeometry},
0383                       logger().cloneWithSuffix("Navigator")),
0384       logger().cloneWithSuffix("Propagator"));
0385
0386   ExtrapolatorOptions extrapolationOptions(ctx.geoContext, ctx.magFieldContext);
0387   extrapolationOptions.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
0388   extrapolationOptions.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
0389
0390   // Perform the track finding for all initial parameters
0391   ACTS_DEBUG("Invoke track finding with " << initialParameters.size()
0392                                           << " seeds.");
0393
0394   auto trackContainer = std::make_shared<Acts::VectorTrackContainer>();
0395   auto trackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
0396
0397   auto trackContainerTemp = std::make_shared<Acts::VectorTrackContainer>();
0398   auto trackStateContainerTemp =
0399       std::make_shared<Acts::VectorMultiTrajectory>();
0400
0401   TrackContainer tracks(trackContainer, trackStateContainer);
0402   TrackContainer tracksTemp(trackContainerTemp, trackStateContainerTemp);
0403
0404   // Note that not all backends support PODs as column types
0405   tracks.addColumn<BranchStopper::BranchState>("MyBranchState");
0406   tracksTemp.addColumn<BranchStopper::BranchState>("MyBranchState");
0407
0408   tracks.addColumn<unsigned int>("trackGroup");
0409   tracksTemp.addColumn<unsigned int>("trackGroup");
0410   Acts::ProxyAccessor<unsigned int> seedNumber("trackGroup");
0411
0412   unsigned int nSeed = 0;
0413
0414   // A map indicating whether a seed has been discovered already
0415   std::unordered_map<SeedIdentifier, bool> discoveredSeeds;
0416
0417   auto addTrack = [&](const TrackProxy& track) {
0418     ++m_nFoundTracks;
0419
0420     // trim the track if requested
0421     if (m_cfg.trimTracks) {
0422       Acts::trimTrack(track, true, true, true, true);
0423     }
0424     Acts::calculateTrackQuantities(track);
0425
0426     if (m_trackSelector.has_value() && !m_trackSelector->isValidTrack(track)) {
0427       return;
0428     }
0429
0430     // flag seeds which are covered by the track
0431     visitSeedIdentifiers(track, [&](const SeedIdentifier& seedIdentifier) {
0432       if (auto it = discoveredSeeds.find(seedIdentifier);
0433           it != discoveredSeeds.end()) {
0434         it->second = true;
0435       }
0436     });
0437
0438     ++m_nSelectedTracks;
0439
0440     auto destProxy = tracks.makeTrack();
0441     // make sure we copy track states!
0442     destProxy.copyFrom(track);
0443   };
0444
0445   if (seeds != nullptr && m_cfg.seedDeduplication) {
0446     // Index the seeds for deduplication
0447     for (const auto& seed : *seeds) {
0448       SeedIdentifier seedIdentifier = makeSeedIdentifier(seed);
0449       discoveredSeeds.emplace(seedIdentifier, false);
0450     }
0451   }
0452
0453   for (std::size_t iSeed = 0; iSeed < initialParameters.size(); ++iSeed) {
0454     m_nTotalSeeds++;
0455
0456     if (seeds != nullptr) {
0457       const SimSeed& seed = seeds->at(iSeed);
0458
0459       if (m_cfg.seedDeduplication) {
0460         SeedIdentifier seedIdentifier = makeSeedIdentifier(seed);
0461         // check if the seed has been discovered already
0462         if (auto it = discoveredSeeds.find(seedIdentifier);
0463             it != discoveredSeeds.end() && it->second) {
0464           m_nDeduplicatedSeeds++;
0465           ACTS_VERBOSE("Skipping seed " << iSeed << " due to deduplication.");
0466           continue;
0467         }
0468       }
0469
0470       if (m_cfg.stayOnSeed) {
0471         measSel.setSeed(seed);
0472       }
0473     }
0474
0475     // Clear trackContainerTemp and trackStateContainerTemp
0476     tracksTemp.clear();
0477
0478     const Acts::BoundTrackParameters& firstInitialParameters =
0479         initialParameters.at(iSeed);
0480     ACTS_VERBOSE("Processing seed " << iSeed << " with initial parameters "
0481                                     << firstInitialParameters);
0482
0483     auto firstRootBranch = tracksTemp.makeTrack();
0484     auto firstResult = (*m_cfg.findTracks)(firstInitialParameters, firstOptions,
0485                                            tracksTemp, firstRootBranch);
0486     nSeed++;
0487
0488     if (!firstResult.ok()) {
0489       m_nFailedSeeds++;
0490       ACTS_WARNING("Track finding failed for seed " << iSeed << " with error"
0491                                                     << firstResult.error());
0492       continue;
0493     }
0494
0495     auto& firstTracksForSeed = firstResult.value();
0496     for (auto& firstTrack : firstTracksForSeed) {
0497       // TODO a copy of the track should not be necessary but is the safest way
0498       //      with the current EDM
0499       // TODO a lightweight copy without copying all the track state components
0500       //      might be a solution
0501       auto trackCandidate = tracksTemp.makeTrack();
0502       trackCandidate.copyFrom(firstTrack);
0503
0504       Acts::Result<void> firstSmoothingResult{
0505           Acts::smoothTrack(ctx.geoContext, trackCandidate, logger())};
0506       if (!firstSmoothingResult.ok()) {
0507         m_nFailedSmoothing++;
0508         ACTS_ERROR("First smoothing for seed "
0509                    << iSeed << " and track " << firstTrack.index()
0510                    << " failed with error " << firstSmoothingResult.error());
0511         continue;
0512       }
0513
0514       // number of second tracks found
0515       std::size_t nSecond = 0;
0516
0517       // Set the seed number, this number decrease by 1 since the seed number
0518       // has already been updated
0519       seedNumber(trackCandidate) = nSeed - 1;
0520
0521       if (m_cfg.twoWay) {
0522         std::optional<Acts::VectorMultiTrajectory::TrackStateProxy>
0523             firstMeasurementOpt;
0524         for (auto trackState : trackCandidate.trackStatesReversed()) {
0525           bool isMeasurement = trackState.typeFlags().test(
0526               Acts::TrackStateFlag::MeasurementFlag);
0527           bool isOutlier =
0528               trackState.typeFlags().test(Acts::TrackStateFlag::OutlierFlag);
0529           // We are excluding non measurement states and outlier here. Those can
0530           // decrease resolution because only the smoothing corrected the very
0531           // first prediction as filtering is not possible.
0532           if (isMeasurement && !isOutlier) {
0533             firstMeasurementOpt = trackState;
0534           }
0535         }
0536
0537         if (firstMeasurementOpt.has_value()) {
0538           TrackContainer::TrackStateProxy firstMeasurement{
0539               firstMeasurementOpt.value()};
0540           TrackContainer::ConstTrackStateProxy firstMeasurementConst{
0541               firstMeasurement};
0542
0543           Acts::BoundTrackParameters secondInitialParameters =
0544               trackCandidate.createParametersFromState(firstMeasurementConst);
0545
0546           if (!secondInitialParameters.referenceSurface().insideBounds(
0547                   secondInitialParameters.localPosition())) {
0548             m_nSkippedSecondPass++;
0549             ACTS_DEBUG(
0550                 "Smoothing of first pass fit produced out-of-bounds parameters "
0551                 "relative to the surface. Skipping second pass.");
0552             continue;
0553           }
0554
0555           auto secondRootBranch = tracksTemp.makeTrack();
0556           secondRootBranch.copyFromWithoutStates(trackCandidate);
0557           auto secondResult =
0558               (*m_cfg.findTracks)(secondInitialParameters, secondOptions,
0559                                   tracksTemp, secondRootBranch);
0560
0561           if (!secondResult.ok()) {
0562             ACTS_WARNING("Second track finding failed for seed "
0563                          << iSeed << " with error" << secondResult.error());
0564           } else {
0565             // store the original previous state to restore it later
0566             auto originalFirstMeasurementPrevious = firstMeasurement.previous();
0567
0568             auto& secondTracksForSeed = secondResult.value();
0569             for (auto& secondTrack : secondTracksForSeed) {
0570               // TODO a copy of the track should not be necessary but is the
0571               //      safest way with the current EDM
0572               // TODO a lightweight copy without copying all the track state
0573               //      components might be a solution
0574               auto secondTrackCopy = tracksTemp.makeTrack();
0575               secondTrackCopy.copyFrom(secondTrack);
0576
0577               // Note that this is only valid if there are no branches
0578               // We disallow this by breaking this look after a second track was
0579               // processed
0580               secondTrackCopy.reverseTrackStates(true);
0581
0582               firstMeasurement.previous() =
0583                   secondTrackCopy.outermostTrackState().index();
0584
0585               // Retain tip and stem index of the first track
0586               auto tipIndex = trackCandidate.tipIndex();
0587               auto stemIndex = trackCandidate.stemIndex();
0588               trackCandidate.copyFromWithoutStates(secondTrackCopy);
0589               trackCandidate.tipIndex() = tipIndex;
0590               trackCandidate.stemIndex() = stemIndex;
0591
0592               // finalize the track candidate
0593
0594               bool doExtrapolate = true;
0595
0596               if (!m_cfg.reverseSearch) {
0597                 // these parameters are already extrapolated by the CKF and have
0598                 // the optimal resolution. note that we did not smooth all the
0599                 // states.
0600
0601                 // only extrapolate if we did not do it already
0602                 doExtrapolate = !trackCandidate.hasReferenceSurface();
0603               } else {
0604                 // smooth the full track and extrapolate to the reference
0605
0606                 auto secondSmoothingResult =
0607                     Acts::smoothTrack(ctx.geoContext, trackCandidate, logger());
0608                 if (!secondSmoothingResult.ok()) {
0609                   m_nFailedSmoothing++;
0610                   ACTS_ERROR("Second smoothing for seed "
0611                              << iSeed << " and track " << secondTrack.index()
0612                              << " failed with error "
0613                              << secondSmoothingResult.error());
0614                   continue;
0615                 }
0616
0617                 trackCandidate.reverseTrackStates(true);
0618               }
0619
0620               if (doExtrapolate) {
0621                 auto secondExtrapolationResult =
0622                     Acts::extrapolateTrackToReferenceSurface(
0623                         trackCandidate, *pSurface, extrapolator,
0624                         extrapolationOptions, m_cfg.extrapolationStrategy,
0625                         logger());
0626                 if (!secondExtrapolationResult.ok()) {
0627                   m_nFailedExtrapolation++;
0628                   ACTS_ERROR("Second extrapolation for seed "
0629                              << iSeed << " and track " << secondTrack.index()
0630                              << " failed with error "
0631                              << secondExtrapolationResult.error());
0632                   continue;
0633                 }
0634               }
0635
0636               addTrack(trackCandidate);
0637
0638               ++nSecond;
0639             }
0640
0641             // restore the original previous state
0642             firstMeasurement.previous() = originalFirstMeasurementPrevious;
0643           }
0644         }
0645       }
0646
0647       // if no second track was found, we will use only the first track
0648       if (nSecond == 0) {
0649         // restore the track to the original state
0650         auto tipIndex = trackCandidate.tipIndex();
0651         auto stemIndex = trackCandidate.stemIndex();
0652         trackCandidate.copyFromWithoutStates(firstTrack);
0653         trackCandidate.tipIndex() = tipIndex;
0654         trackCandidate.stemIndex() = stemIndex;
0655
0656         auto firstExtrapolationResult =
0657             Acts::extrapolateTrackToReferenceSurface(
0658                 trackCandidate, *pSurface, extrapolator, extrapolationOptions,
0659                 m_cfg.extrapolationStrategy, logger());
0660         if (!firstExtrapolationResult.ok()) {
0661           m_nFailedExtrapolation++;
0662           ACTS_ERROR("Extrapolation for seed "
0663                      << iSeed << " and track " << firstTrack.index()
0664                      << " failed with error "
0665                      << firstExtrapolationResult.error());
0666           continue;
0667         }
0668
0669         addTrack(trackCandidate);
0670       }
0671     }
0672   }
0673
0674   // Compute shared hits from all the reconstructed tracks
0675   if (m_cfg.computeSharedHits) {
0676     computeSharedHits(tracks, measurements);
0677   }
0678
0679   ACTS_DEBUG("Finalized track finding with " << tracks.size()
0680                                              << " track candidates.");
0681
0682   m_nStoppedBranches += branchStopper.m_nStoppedBranches;
0683
0684   m_memoryStatistics.local().hist +=
0685       tracks.trackStateContainer().statistics().hist;
0686
0687   auto constTrackStateContainer =
0688       std::make_shared<Acts::ConstVectorMultiTrajectory>(
0689           std::move(*trackStateContainer));
0690
0691   auto constTrackContainer = std::make_shared<Acts::ConstVectorTrackContainer>(
0692       std::move(*trackContainer));
0693
0694   ConstTrackContainer constTracks{constTrackContainer,
0695                                   constTrackStateContainer};
0696
0697   m_outputTracks(ctx, std::move(constTracks));
0698   return ProcessCode::SUCCESS;
0699 }
0700
0701 ProcessCode TrackFindingAlgorithm::finalize() {
0702   ACTS_INFO("TrackFindingAlgorithm statistics:");
0703   ACTS_INFO("- total seeds: " << m_nTotalSeeds);
0704   ACTS_INFO("- deduplicated seeds: " << m_nDeduplicatedSeeds);
0705   ACTS_INFO("- failed seeds: " << m_nFailedSeeds);
0706   ACTS_INFO("- failed smoothing: " << m_nFailedSmoothing);
0707   ACTS_INFO("- failed extrapolation: " << m_nFailedExtrapolation);
0708   ACTS_INFO("- failure ratio seeds: " << static_cast<double>(m_nFailedSeeds) /
0709                                              m_nTotalSeeds);
0710   ACTS_INFO("- found tracks: " << m_nFoundTracks);
0711   ACTS_INFO("- selected tracks: " << m_nSelectedTracks);
0712   ACTS_INFO("- stopped branches: " << m_nStoppedBranches);
0713   ACTS_INFO("- skipped second pass: " << m_nSkippedSecondPass);
0714
0715   auto memoryStatistics =
0716       m_memoryStatistics.combine([](const auto& a, const auto& b) {
0717         Acts::VectorMultiTrajectory::Statistics c;
0718         c.hist = a.hist + b.hist;
0719         return c;
0720       });
0721   std::stringstream ss;
0722   memoryStatistics.toStream(ss);
0723   ACTS_DEBUG("Track State memory statistics (averaged):\n" << ss.str());
0724   return ProcessCode::SUCCESS;
0725 }
0726
0727 // TODO this is somewhat duplicated in AmbiguityResolutionAlgorithm.cpp
0728 // TODO we should make a common implementation in the core at some point
0729 void TrackFindingAlgorithm::computeSharedHits(
0730     TrackContainer& tracks, const MeasurementContainer& measurements) const {
0731   // Compute shared hits from all the reconstructed tracks
0732   // Compute nSharedhits and Update ckf results
0733   // hit index -> list of multi traj indexes [traj, meas]
0734
0735   std::vector<std::size_t> firstTrackOnTheHit(
0736       measurements.size(), std::numeric_limits<std::size_t>::max());
0737   std::vector<std::size_t> firstStateOnTheHit(
0738       measurements.size(), std::numeric_limits<std::size_t>::max());
0739
0740   for (auto track : tracks) {
0741     for (auto state : track.trackStatesReversed()) {
0742       if (!state.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) {
0743         continue;
0744       }
0745
0746       std::size_t hitIndex = state.getUncalibratedSourceLink()
0747                                  .template get<IndexSourceLink>()
0748                                  .index();
0749
0750       // Check if hit not already used
0751       if (firstTrackOnTheHit.at(hitIndex) ==
0752           std::numeric_limits<std::size_t>::max()) {
0753         firstTrackOnTheHit.at(hitIndex) = track.index();
0754         firstStateOnTheHit.at(hitIndex) = state.index();
0755         continue;
0756       }
0757
0758       // if already used, control if first track state has been marked
0759       // as shared
0760       int indexFirstTrack = firstTrackOnTheHit.at(hitIndex);
0761       int indexFirstState = firstStateOnTheHit.at(hitIndex);
0762
0763       auto firstState = tracks.getTrack(indexFirstTrack)
0764                             .container()
0765                             .trackStateContainer()
0766                             .getTrackState(indexFirstState);
0767       if (!firstState.typeFlags().test(Acts::TrackStateFlag::SharedHitFlag)) {
0768         firstState.typeFlags().set(Acts::TrackStateFlag::SharedHitFlag);
0769       }
0770
0771       // Decorate this track state
0772       state.typeFlags().set(Acts::TrackStateFlag::SharedHitFlag);
0773     }
0774   }
0775 }
0776
0777 }  // namespace ActsExamples