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/VectorMultiTrajectory.hpp"
0019 #include "Acts/EventData/VectorTrackContainer.hpp"
0020 #include "Acts/Geometry/GeometryIdentifier.hpp"
0021 #include "Acts/Propagator/MaterialInteractor.hpp"
0022 #include "Acts/Propagator/Navigator.hpp"
0023 #include "Acts/Propagator/Propagator.hpp"
0024 #include "Acts/Propagator/StandardAborters.hpp"
0025 #include "Acts/Propagator/SympyStepper.hpp"
0026 #include "Acts/Surfaces/PerigeeSurface.hpp"
0027 #include "Acts/Surfaces/Surface.hpp"
0028 #include "Acts/TrackFinding/TrackStateCreator.hpp"
0029 #include "Acts/TrackFitting/GainMatrixUpdater.hpp"
0030 #include "Acts/Utilities/Enumerate.hpp"
0031 #include "Acts/Utilities/Logger.hpp"
0032 #include "Acts/Utilities/TrackHelpers.hpp"
0033 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0034 #include "ActsExamples/EventData/Measurement.hpp"
0035 #include "ActsExamples/EventData/MeasurementCalibration.hpp"
0036 #include "ActsExamples/EventData/Seed.hpp"
0037 #include "ActsExamples/EventData/SpacePoint.hpp"
0038 #include "ActsExamples/EventData/Track.hpp"
0039 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0040 #include "ActsExamples/Framework/ProcessCode.hpp"
0041
0042 #include <functional>
0043 #include <memory>
0044 #include <optional>
0045 #include <ostream>
0046 #include <stdexcept>
0047 #include <unordered_map>
0048 #include <utility>
0049
0050 #include <boost/functional/hash.hpp>
0051
0052 // Specialize std::hash for SeedIdentifier
0053 // This is required to use SeedIdentifier as a key in an `std::unordered_map`.
0054 template <class T, std::size_t N>
0055 struct std::hash<std::array<T, N>> {
0056   std::size_t operator()(const std::array<T, N>& array) const {
0057     std::hash<T> hasher;
0058     std::size_t result = 0;
0059     for (auto&& element : array) {
0060       boost::hash_combine(result, hasher(element));
0061     }
0062     return result;
0063   }
0064 };
0065
0066 namespace ActsExamples {
0067
0068 namespace {
0069
0070 class MeasurementSelector {
0071  public:
0072   using Traj = Acts::VectorMultiTrajectory;
0073
0074   explicit MeasurementSelector(Acts::MeasurementSelector selector)
0075       : m_selector(std::move(selector)) {}
0076
0077   void setSeed(const std::optional<ConstSeedProxy>& seed) { m_seed = seed; }
0078
0079   Acts::Result<std::pair<std::vector<Traj::TrackStateProxy>::iterator,
0080                          std::vector<Traj::TrackStateProxy>::iterator>>
0081   select(std::vector<Traj::TrackStateProxy>& candidates, bool& isOutlier,
0082          const Acts::Logger& logger) const {
0083     if (m_seed.has_value()) {
0084       std::vector<Traj::TrackStateProxy> newCandidates;
0085
0086       for (const auto& candidate : candidates) {
0087         if (isSeedCandidate(candidate)) {
0088           newCandidates.push_back(candidate);
0089         }
0090       }
0091
0092       if (!newCandidates.empty()) {
0093         candidates = std::move(newCandidates);
0094       }
0095     }
0096
0097     return m_selector.select<Acts::VectorMultiTrajectory>(candidates, isOutlier,
0098                                                           logger);
0099   }
0100
0101  private:
0102   Acts::MeasurementSelector m_selector;
0103   std::optional<ConstSeedProxy> m_seed;
0104
0105   bool isSeedCandidate(const Traj::TrackStateProxy& candidate) const {
0106     assert(candidate.hasUncalibratedSourceLink());
0107     assert(m_seed.has_value());
0108
0109     const Acts::SourceLink& sourceLink = candidate.getUncalibratedSourceLink();
0110
0111     for (const ConstSpacePointProxy sp : m_seed->spacePoints()) {
0112       for (const Acts::SourceLink& sl : sp.sourceLinks()) {
0113         if (sourceLink.get<IndexSourceLink>() == sl.get<IndexSourceLink>()) {
0114           return true;
0115         }
0116       }
0117     }
0118
0119     return false;
0120   }
0121 };
0122
0123 /// Source link indices of the bottom, middle, top measurements.
0124 /// In case of strip seeds only the first source link of the pair is used.
0125 using SeedIdentifier = std::array<Index, 3>;
0126
0127 /// Build a seed identifier from a seed.
0128 ///
0129 /// @param seed The seed to build the identifier from.
0130 /// @return The seed identifier.
0131 SeedIdentifier makeSeedIdentifier(const ConstSeedProxy& seed) {
0132   SeedIdentifier result;
0133
0134   for (const auto& [i, spIndex] : Acts::enumerate(seed.spacePointIndices())) {
0135     const ConstSpacePointProxy sp =
0136         seed.container().spacePointContainer().at(spIndex);
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().isHole() ||
0223           trackState.typeFlags().isOutlier()) {
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(
0259     Config config, std::unique_ptr<const Acts::Logger> logger)
0260     : IAlgorithm("TrackFindingAlgorithm", std::move(logger)),
0261       m_cfg(std::move(config)) {
0262   if (m_cfg.inputMeasurements.empty()) {
0263     throw std::invalid_argument("Missing measurements input collection");
0264   }
0265   if (m_cfg.inputInitialTrackParameters.empty()) {
0266     throw std::invalid_argument(
0267         "Missing initial track parameters input collection");
0268   }
0269   if (m_cfg.outputTracks.empty()) {
0270     throw std::invalid_argument("Missing tracks output collection");
0271   }
0272
0273   if (m_cfg.seedDeduplication && m_cfg.inputSeeds.empty()) {
0274     throw std::invalid_argument(
0275         "Missing seeds input collection. This is "
0276         "required for seed deduplication.");
0277   }
0278   if (m_cfg.stayOnSeed && m_cfg.inputSeeds.empty()) {
0279     throw std::invalid_argument(
0280         "Missing seeds input collection. This is "
0281         "required for staying on seed.");
0282   }
0283
0284   if (m_cfg.trackSelectorCfg.has_value()) {
0285     m_trackSelector = std::visit(
0286         [](const auto& cfg) -> std::optional<Acts::TrackSelector> {
0287           return Acts::TrackSelector(cfg);
0288         },
0289         m_cfg.trackSelectorCfg.value());
0290   }
0291
0292   m_inputMeasurements.initialize(m_cfg.inputMeasurements);
0293   m_inputInitialTrackParameters.initialize(m_cfg.inputInitialTrackParameters);
0294   m_inputSeeds.maybeInitialize(m_cfg.inputSeeds);
0295   m_outputTracks.initialize(m_cfg.outputTracks);
0296 }
0297
0298 ProcessCode TrackFindingAlgorithm::execute(const AlgorithmContext& ctx) const {
0299   // Read input data
0300   const auto& measurements = m_inputMeasurements(ctx);
0301   const auto& initialParameters = m_inputInitialTrackParameters(ctx);
0302   const SeedContainer* seeds = nullptr;
0303
0304   if (m_inputSeeds.isInitialized()) {
0305     seeds = &m_inputSeeds(ctx);
0306
0307     if (initialParameters.size() != seeds->size()) {
0308       ACTS_ERROR("Number of initial parameters and seeds do not match. "
0309                  << initialParameters.size() << " != " << seeds->size());
0310     }
0311   }
0312
0313   // Construct a perigee surface as the target surface
0314   auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
0315       Acts::Vector3{0., 0., 0.});
0316
0317   PassThroughCalibrator pcalibrator;
0318   MeasurementCalibratorAdapter calibrator(pcalibrator, measurements);
0319   Acts::GainMatrixUpdater kfUpdater;
0320
0321   using Extensions = Acts::CombinatorialKalmanFilterExtensions<TrackContainer>;
0322
0323   BranchStopper branchStopper(m_cfg);
0324   MeasurementSelector measSel{
0325       Acts::MeasurementSelector(m_cfg.measurementSelectorCfg)};
0326
0327   IndexSourceLinkAccessor slAccessor;
0328   slAccessor.container = &measurements.orderedIndices();
0329
0330   using TrackStateCreatorType =
0331       Acts::TrackStateCreator<IndexSourceLinkAccessor::Iterator,
0332                               TrackContainer>;
0333   TrackStateCreatorType trackStateCreator;
0334   trackStateCreator.sourceLinkAccessor
0335       .template connect<&IndexSourceLinkAccessor::range>(&slAccessor);
0336   trackStateCreator.calibrator
0337       .template connect<&MeasurementCalibratorAdapter::calibrate>(&calibrator);
0338   trackStateCreator.measurementSelector
0339       .template connect<&MeasurementSelector::select>(&measSel);
0340
0341   Extensions extensions;
0342   extensions.updater.connect<&Acts::GainMatrixUpdater::operator()<
0343       typename TrackContainer::TrackStateContainerBackend>>(&kfUpdater);
0344   extensions.branchStopper.connect<&BranchStopper::operator()>(&branchStopper);
0345   extensions.createTrackStates
0346       .template connect<&TrackStateCreatorType ::createTrackStates>(
0347           &trackStateCreator);
0348
0349   Acts::PropagatorPlainOptions firstPropOptions(ctx.geoContext,
0350                                                 ctx.magFieldContext);
0351   firstPropOptions.maxSteps = m_cfg.maxSteps;
0352   firstPropOptions.direction = m_cfg.reverseSearch ? Acts::Direction::Backward()
0353                                                    : Acts::Direction::Forward();
0354   firstPropOptions.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
0355   firstPropOptions.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
0356
0357   Acts::PropagatorPlainOptions secondPropOptions(ctx.geoContext,
0358                                                  ctx.magFieldContext);
0359   secondPropOptions.maxSteps = m_cfg.maxSteps;
0360   secondPropOptions.direction = firstPropOptions.direction.invert();
0361   secondPropOptions.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
0362   secondPropOptions.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
0363
0364   // Set the CombinatorialKalmanFilter options
0365   TrackFinderOptions firstOptions(ctx.geoContext, ctx.magFieldContext,
0366                                   ctx.calibContext, extensions,
0367                                   firstPropOptions);
0368
0369   firstOptions.targetSurface = m_cfg.reverseSearch ? pSurface.get() : nullptr;
0370
0371   TrackFinderOptions secondOptions(ctx.geoContext, ctx.magFieldContext,
0372                                    ctx.calibContext, extensions,
0373                                    secondPropOptions);
0374   secondOptions.targetSurface = m_cfg.reverseSearch ? nullptr : pSurface.get();
0375   secondOptions.skipPrePropagationUpdate = true;
0376
0377   using Extrapolator = Acts::Propagator<Acts::SympyStepper, Acts::Navigator>;
0378   using ExtrapolatorOptions = Extrapolator::template Options<
0379       Acts::ActorList<Acts::MaterialInteractor, Acts::EndOfWorldReached>>;
0380
0381   Extrapolator extrapolator(
0382       Acts::SympyStepper(m_cfg.magneticField),
0383       Acts::Navigator({m_cfg.trackingGeometry},
0384                       logger().cloneWithSuffix("Navigator")),
0385       logger().cloneWithSuffix("Propagator"));
0386
0387   ExtrapolatorOptions extrapolationOptions(ctx.geoContext, ctx.magFieldContext);
0388   extrapolationOptions.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
0389   extrapolationOptions.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
0390
0391   // Perform the track finding for all initial parameters
0392   ACTS_DEBUG("Invoke track finding with " << initialParameters.size()
0393                                           << " seeds.");
0394
0395   auto trackContainer = std::make_shared<Acts::VectorTrackContainer>();
0396   auto trackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
0397
0398   auto trackContainerTemp = std::make_shared<Acts::VectorTrackContainer>();
0399   auto trackStateContainerTemp =
0400       std::make_shared<Acts::VectorMultiTrajectory>();
0401
0402   TrackContainer tracks(trackContainer, trackStateContainer);
0403   TrackContainer tracksTemp(trackContainerTemp, trackStateContainerTemp);
0404
0405   // Note that not all backends support PODs as column types
0406   tracks.addColumn<BranchStopper::BranchState>("MyBranchState");
0407   tracksTemp.addColumn<BranchStopper::BranchState>("MyBranchState");
0408
0409   tracks.addColumn<unsigned int>("trackGroup");
0410   tracksTemp.addColumn<unsigned int>("trackGroup");
0411   Acts::ProxyAccessor<unsigned int> seedNumber("trackGroup");
0412
0413   unsigned int nSeed = 0;
0414
0415   // A map indicating whether a seed has been discovered already
0416   std::unordered_map<SeedIdentifier, bool> discoveredSeeds;
0417
0418   auto addTrack = [&](const TrackProxy& track) {
0419     ++m_nFoundTracks;
0420
0421     // trim the track if requested
0422     if (m_cfg.trimTracks) {
0423       Acts::trimTrack(track, true, true, true, true);
0424     }
0425     Acts::calculateTrackQuantities(track);
0426
0427     if (m_trackSelector.has_value() && !m_trackSelector->isValidTrack(track)) {
0428       return;
0429     }
0430
0431     // flag seeds which are covered by the track
0432     visitSeedIdentifiers(track, [&](const SeedIdentifier& seedIdentifier) {
0433       if (auto it = discoveredSeeds.find(seedIdentifier);
0434           it != discoveredSeeds.end()) {
0435         it->second = true;
0436       }
0437     });
0438
0439     ++m_nSelectedTracks;
0440
0441     auto destProxy = tracks.makeTrack();
0442     // make sure we copy track states!
0443     destProxy.copyFrom(track);
0444   };
0445
0446   if (seeds != nullptr && m_cfg.seedDeduplication) {
0447     // Index the seeds for deduplication
0448     for (const auto& seed : *seeds) {
0449       SeedIdentifier seedIdentifier = makeSeedIdentifier(seed);
0450       discoveredSeeds.emplace(seedIdentifier, false);
0451     }
0452   }
0453
0454   for (std::size_t iSeed = 0; iSeed < initialParameters.size(); ++iSeed) {
0455     m_nTotalSeeds++;
0456
0457     if (seeds != nullptr) {
0458       const ConstSeedProxy seed = seeds->at(iSeed);
0459
0460       if (m_cfg.seedDeduplication) {
0461         SeedIdentifier seedIdentifier = makeSeedIdentifier(seed);
0462         // check if the seed has been discovered already
0463         if (auto it = discoveredSeeds.find(seedIdentifier);
0464             it != discoveredSeeds.end() && it->second) {
0465           m_nDeduplicatedSeeds++;
0466           ACTS_VERBOSE("Skipping seed " << iSeed << " due to deduplication.");
0467           continue;
0468         }
0469       }
0470
0471       if (m_cfg.stayOnSeed) {
0472         measSel.setSeed(seed);
0473       }
0474     }
0475
0476     // Clear trackContainerTemp and trackStateContainerTemp
0477     tracksTemp.clear();
0478
0479     const Acts::BoundTrackParameters& firstInitialParameters =
0480         initialParameters.at(iSeed);
0481     ACTS_VERBOSE("Processing seed " << iSeed << " with initial parameters "
0482                                     << firstInitialParameters);
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);
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           // We are excluding non measurement states and outlier here. Those can
0527           // decrease resolution because only the smoothing corrected the very
0528           // first prediction as filtering is not possible.
0529           if (trackState.typeFlags().isMeasurement()) {
0530             firstMeasurementOpt = trackState;
0531           }
0532         }
0533
0534         if (firstMeasurementOpt.has_value()) {
0535           TrackContainer::TrackStateProxy firstMeasurement{
0536               firstMeasurementOpt.value()};
0537           TrackContainer::ConstTrackStateProxy firstMeasurementConst{
0538               firstMeasurement};
0539
0540           Acts::BoundTrackParameters secondInitialParameters =
0541               trackCandidate.createParametersFromState(firstMeasurementConst);
0542
0543           if (!secondInitialParameters.referenceSurface().insideBounds(
0544                   secondInitialParameters.localPosition())) {
0545             m_nSkippedSecondPass++;
0546             ACTS_DEBUG(
0547                 "Smoothing of first pass fit produced out-of-bounds parameters "
0548                 "relative to the surface. Skipping second pass.");
0549             continue;
0550           }
0551
0552           auto secondRootBranch = tracksTemp.makeTrack();
0553           secondRootBranch.copyFromWithoutStates(trackCandidate);
0554           auto secondResult =
0555               (*m_cfg.findTracks)(secondInitialParameters, secondOptions,
0556                                   tracksTemp, secondRootBranch);
0557
0558           if (!secondResult.ok()) {
0559             ACTS_WARNING("Second track finding failed for seed "
0560                          << iSeed << " with error" << secondResult.error());
0561           } else {
0562             // store the original previous state to restore it later
0563             auto originalFirstMeasurementPrevious = firstMeasurement.previous();
0564
0565             auto& secondTracksForSeed = secondResult.value();
0566             for (auto& secondTrack : secondTracksForSeed) {
0567               // TODO a copy of the track should not be necessary but is the
0568               //      safest way with the current EDM
0569               // TODO a lightweight copy without copying all the track state
0570               //      components might be a solution
0571               auto secondTrackCopy = tracksTemp.makeTrack();
0572               secondTrackCopy.copyFrom(secondTrack);
0573
0574               // Note that this is only valid if there are no branches
0575               // We disallow this by breaking this look after a second track was
0576               // processed
0577               secondTrackCopy.reverseTrackStates(true);
0578
0579               firstMeasurement.previous() =
0580                   secondTrackCopy.outermostTrackState().index();
0581
0582               // Retain tip and stem index of the first track
0583               auto tipIndex = trackCandidate.tipIndex();
0584               auto stemIndex = trackCandidate.stemIndex();
0585               trackCandidate.copyFromWithoutStates(secondTrackCopy);
0586               trackCandidate.tipIndex() = tipIndex;
0587               trackCandidate.stemIndex() = stemIndex;
0588
0589               // finalize the track candidate
0590
0591               bool doExtrapolate = true;
0592
0593               if (!m_cfg.reverseSearch) {
0594                 // these parameters are already extrapolated by the CKF and have
0595                 // the optimal resolution. note that we did not smooth all the
0596                 // states.
0597
0598                 // only extrapolate if we did not do it already
0599                 doExtrapolate = !trackCandidate.hasReferenceSurface();
0600               } else {
0601                 // smooth the full track and extrapolate to the reference
0602
0603                 auto secondSmoothingResult =
0604                     Acts::smoothTrack(ctx.geoContext, trackCandidate, logger());
0605                 if (!secondSmoothingResult.ok()) {
0606                   m_nFailedSmoothing++;
0607                   ACTS_ERROR("Second smoothing for seed "
0608                              << iSeed << " and track " << secondTrack.index()
0609                              << " failed with error "
0610                              << secondSmoothingResult.error());
0611                   continue;
0612                 }
0613
0614                 trackCandidate.reverseTrackStates(true);
0615               }
0616
0617               if (doExtrapolate) {
0618                 auto secondExtrapolationResult =
0619                     Acts::extrapolateTrackToReferenceSurface(
0620                         trackCandidate, *pSurface, extrapolator,
0621                         extrapolationOptions, m_cfg.extrapolationStrategy,
0622                         logger());
0623                 if (!secondExtrapolationResult.ok()) {
0624                   m_nFailedExtrapolation++;
0625                   ACTS_ERROR("Second extrapolation for seed "
0626                              << iSeed << " and track " << secondTrack.index()
0627                              << " failed with error "
0628                              << secondExtrapolationResult.error());
0629                   continue;
0630                 }
0631               }
0632
0633               addTrack(trackCandidate);
0634
0635               ++nSecond;
0636             }
0637
0638             // restore the original previous state
0639             firstMeasurement.previous() = originalFirstMeasurementPrevious;
0640           }
0641         }
0642       }
0643
0644       // if no second track was found, we will use only the first track
0645       if (nSecond == 0) {
0646         // restore the track to the original state
0647         auto tipIndex = trackCandidate.tipIndex();
0648         auto stemIndex = trackCandidate.stemIndex();
0649         trackCandidate.copyFromWithoutStates(firstTrack);
0650         trackCandidate.tipIndex() = tipIndex;
0651         trackCandidate.stemIndex() = stemIndex;
0652
0653         auto firstExtrapolationResult =
0654             Acts::extrapolateTrackToReferenceSurface(
0655                 trackCandidate, *pSurface, extrapolator, extrapolationOptions,
0656                 m_cfg.extrapolationStrategy, logger());
0657         if (!firstExtrapolationResult.ok()) {
0658           m_nFailedExtrapolation++;
0659           ACTS_ERROR("Extrapolation for seed "
0660                      << iSeed << " and track " << firstTrack.index()
0661                      << " failed with error "
0662                      << firstExtrapolationResult.error());
0663           continue;
0664         }
0665
0666         addTrack(trackCandidate);
0667       }
0668     }
0669   }
0670
0671   // Compute shared hits from all the reconstructed tracks
0672   if (m_cfg.computeSharedHits) {
0673     computeSharedHits(tracks, measurements);
0674   }
0675
0676   ACTS_DEBUG("Finalized track finding with " << tracks.size()
0677                                              << " track candidates.");
0678
0679   m_nStoppedBranches += branchStopper.m_nStoppedBranches;
0680
0681   m_memoryStatistics.local().hist +=
0682       tracks.trackStateContainer().statistics().hist;
0683
0684   auto constTrackStateContainer =
0685       std::make_shared<Acts::ConstVectorMultiTrajectory>(
0686           std::move(*trackStateContainer));
0687
0688   auto constTrackContainer = std::make_shared<Acts::ConstVectorTrackContainer>(
0689       std::move(*trackContainer));
0690
0691   ConstTrackContainer constTracks{constTrackContainer,
0692                                   constTrackStateContainer};
0693
0694   m_outputTracks(ctx, std::move(constTracks));
0695   return ProcessCode::SUCCESS;
0696 }
0697
0698 ProcessCode TrackFindingAlgorithm::finalize() {
0699   ACTS_INFO("TrackFindingAlgorithm statistics:");
0700   ACTS_INFO("- total seeds: " << m_nTotalSeeds);
0701   ACTS_INFO("- deduplicated seeds: " << m_nDeduplicatedSeeds);
0702   ACTS_INFO("- failed seeds: " << m_nFailedSeeds);
0703   ACTS_INFO("- failed smoothing: " << m_nFailedSmoothing);
0704   ACTS_INFO("- failed extrapolation: " << m_nFailedExtrapolation);
0705   ACTS_INFO("- failure ratio seeds: " << static_cast<double>(m_nFailedSeeds) /
0706                                              m_nTotalSeeds);
0707   ACTS_INFO("- found tracks: " << m_nFoundTracks);
0708   ACTS_INFO("- selected tracks: " << m_nSelectedTracks);
0709   ACTS_INFO("- stopped branches: " << m_nStoppedBranches);
0710   ACTS_INFO("- skipped second pass: " << m_nSkippedSecondPass);
0711
0712   auto memoryStatistics =
0713       m_memoryStatistics.combine([](const auto& a, const auto& b) {
0714         Acts::VectorMultiTrajectory::Statistics c;
0715         c.hist = a.hist + b.hist;
0716         return c;
0717       });
0718   std::stringstream ss;
0719   memoryStatistics.toStream(ss);
0720   ACTS_DEBUG("Track State memory statistics (averaged):\n" << ss.str());
0721   return ProcessCode::SUCCESS;
0722 }
0723
0724 // TODO this is somewhat duplicated in AmbiguityResolutionAlgorithm.cpp
0725 // TODO we should make a common implementation in the core at some point
0726 void TrackFindingAlgorithm::computeSharedHits(
0727     TrackContainer& tracks, const MeasurementContainer& measurements) const {
0728   // Compute shared hits from all the reconstructed tracks
0729   // Compute nSharedhits and Update ckf results
0730   // hit index -> list of multi traj indexes [traj, meas]
0731
0732   std::vector<std::size_t> firstTrackOnTheHit(
0733       measurements.size(), std::numeric_limits<std::size_t>::max());
0734   std::vector<std::size_t> firstStateOnTheHit(
0735       measurements.size(), std::numeric_limits<std::size_t>::max());
0736
0737   for (auto track : tracks) {
0738     for (auto state : track.trackStatesReversed()) {
0739       if (!state.typeFlags().isMeasurement()) {
0740         continue;
0741       }
0742
0743       std::size_t hitIndex = state.getUncalibratedSourceLink()
0744                                  .template get<IndexSourceLink>()
0745                                  .index();
0746
0747       // Check if hit not already used
0748       if (firstTrackOnTheHit.at(hitIndex) ==
0749           std::numeric_limits<std::size_t>::max()) {
0750         firstTrackOnTheHit.at(hitIndex) = track.index();
0751         firstStateOnTheHit.at(hitIndex) = state.index();
0752         continue;
0753       }
0754
0755       // if already used, control if first track state has been marked
0756       // as shared
0757       std::size_t indexFirstTrack = firstTrackOnTheHit.at(hitIndex);
0758       std::size_t indexFirstState = firstStateOnTheHit.at(hitIndex);
0759
0760       auto firstState = tracks.getTrack(indexFirstTrack)
0761                             .container()
0762                             .trackStateContainer()
0763                             .getTrackState(indexFirstState);
0764       firstState.typeFlags().setIsSharedHit();
0765
0766       // Decorate this track state
0767       state.typeFlags().setIsSharedHit();
0768     }
0769   }
0770 }
0771
0772 }  // namespace ActsExamples