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