Utilities/src/TrackExtrapolationAlgorithm.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/Utilities/TrackExtrapolationAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Direction.hpp"
0012 #include "Acts/EventData/BoundTrackParameters.hpp"
0013 #include "Acts/EventData/TrackContainer.hpp"
0014 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0015 #include "Acts/EventData/VectorTrackContainer.hpp"
0016 #include "Acts/Geometry/TrackingGeometry.hpp"
0017 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0018 #include "Acts/Propagator/ActorList.hpp"
0019 #include "Acts/Propagator/MaterialInteractor.hpp"
0020 #include "Acts/Propagator/Navigator.hpp"
0021 #include "Acts/Propagator/Propagator.hpp"
0022 #include "Acts/Propagator/StandardAborters.hpp"
0023 #include "Acts/Propagator/SympyStepper.hpp"
0024 #include "Acts/Surfaces/Surface.hpp"
0025 #include "Acts/Utilities/Result.hpp"
0026
0027 #include <memory>
0028 #include <stdexcept>
0029 #include <utility>
0030
0031 namespace ActsExamples {
0032
0033 TrackExtrapolationAlgorithm::TrackExtrapolationAlgorithm(
0034     Config config, std::unique_ptr<const Acts::Logger> logger)
0035     : IAlgorithm("TrackExtrapolationAlgorithm", std::move(logger)),
0036       m_cfg(std::move(config)) {
0037   if (m_cfg.inputTracks.empty()) {
0038     throw std::invalid_argument("Missing input track collection");
0039   }
0040   if (m_cfg.outputTracks.empty()) {
0041     throw std::invalid_argument("Missing output track collection");
0042   }
0043   if (m_cfg.targetSurface == nullptr) {
0044     throw std::invalid_argument("Missing target surface");
0045   }
0046   if (m_cfg.trackingGeometry == nullptr) {
0047     throw std::invalid_argument("Missing tracking geometry");
0048   }
0049   if (m_cfg.magneticField == nullptr) {
0050     throw std::invalid_argument("Missing magnetic field");
0051   }
0052
0053   m_inputTracks.initialize(m_cfg.inputTracks);
0054   m_outputTracks.initialize(m_cfg.outputTracks);
0055 }
0056
0057 ProcessCode TrackExtrapolationAlgorithm::execute(
0058     const AlgorithmContext& ctx) const {
0059   const ConstTrackContainer& inputTracks = m_inputTracks(ctx);
0060
0061   using Propagator = Acts::Propagator<Acts::SympyStepper, Acts::Navigator>;
0062   using Options = Propagator::Options<
0063       Acts::ActorList<Acts::MaterialInteractor, Acts::EndOfWorldReached>>;
0064
0065   const Propagator propagator(
0066       Acts::SympyStepper(m_cfg.magneticField),
0067       Acts::Navigator({m_cfg.trackingGeometry},
0068                       logger().cloneWithSuffix("Navigator")),
0069       logger().cloneWithSuffix("Propagator"));
0070
0071   Options options(ctx.recoGeoContext, ctx.magFieldContext);
0072   options.constrainToVolumeIds = m_cfg.constrainToVolumeIds;
0073   options.endOfWorldVolumeIds = m_cfg.endOfWorldVolumeIds;
0074
0075   // `Acts::extrapolateTrackToReferenceSurface` split up, so the states are read
0076   // off the input track and only the parameters are written to the output one
0077   auto extrapolate = [&](const ConstTrackProxy& track)
0078       -> Acts::Result<Acts::BoundTrackParameters> {
0079     auto findResult = Acts::findTrackStateForExtrapolation(
0080         ctx.recoGeoContext, track, *m_cfg.targetSurface, m_cfg.strategy,
0081         logger());
0082     if (!findResult.ok()) {
0083       return findResult.error();
0084     }
0085     const auto& [trackState, distance] = *findResult;
0086
0087     Options trackOptions = options;
0088     trackOptions.direction =
0089         Acts::Direction::fromScalarZeroAsPositive(distance);
0090
0091     auto propagateResult =
0092         propagator.propagate<Options, Acts::ForcedSurfaceReached>(
0093             track.createParametersFromState(trackState), *m_cfg.targetSurface,
0094             trackOptions);
0095     if (!propagateResult.ok()) {
0096       return propagateResult.error();
0097     }
0098
0099     return propagateResult->endParameters.value();
0100   };
0101
0102   // The tip and stem indices point into the input state backend, which the
0103   // output container takes over below. The empty backend here is only because
0104   // `Acts::TrackContainer` cannot pair a mutable track backend with a
0105   // read-only state one.
0106   auto trackBackend = std::make_shared<Acts::VectorTrackContainer>();
0107   TrackContainer extrapolated{trackBackend,
0108                               std::make_shared<Acts::VectorMultiTrajectory>()};
0109   extrapolated.ensureDynamicColumns(inputTracks);
0110
0111   std::size_t nFailed = 0;
0112
0113   for (const auto& track : inputTracks) {
0114     // one output track per input track, so the indices stay the same and any
0115     // truth matching of the input remains valid
0116     auto destination = extrapolated.makeTrack();
0117     destination.copyFromShallow(track);
0118
0119     const auto result = extrapolate(track);
0120     if (!result.ok()) {
0121       ACTS_DEBUG("Extrapolation of track " << track.index() << " failed with "
0122                                            << result.error());
0123       // no parameters on the target surface
0124       destination.setReferenceSurface(nullptr);
0125       ++nFailed;
0126       continue;
0127     }
0128
0129     destination.setReferenceSurface(m_cfg.targetSurface);
0130     destination.parameters() = result->parameters();
0131     destination.covariance() = result->covariance().value();
0132   }
0133
0134   m_nTotalTracks += inputTracks.size();
0135   m_nFailedTracks += nFailed;
0136
0137   if (nFailed > 0) {
0138     ACTS_DEBUG(nFailed << " tracks could not be extrapolated and are left "
0139                           "without a reference surface");
0140   }
0141
0142   ConstTrackContainer outputTracks{
0143       std::make_shared<Acts::ConstVectorTrackContainer>(
0144           std::move(*trackBackend)),
0145       inputTracks.trackStateContainerHolder()};
0146
0147   ACTS_DEBUG("Extrapolated " << (outputTracks.size() - nFailed) << " of "
0148                              << outputTracks.size() << " tracks");
0149
0150   m_outputTracks(ctx, std::move(outputTracks));
0151
0152   return ProcessCode::SUCCESS;
0153 }
0154
0155 ProcessCode TrackExtrapolationAlgorithm::finalize() {
0156   ACTS_INFO("TrackExtrapolationAlgorithm statistics:");
0157   ACTS_INFO("- total tracks: " << m_nTotalTracks);
0158   ACTS_INFO("- failed tracks: " << m_nFailedTracks);
0159   if (m_nTotalTracks > 0) {
0160     ACTS_INFO("- failure ratio: " << static_cast<double>(m_nFailedTracks) /
0161                                          m_nTotalTracks);
0162   }
0163
0164   return ProcessCode::SUCCESS;
0165 }
0166
0167 }  // namespace ActsExamples