TrackFitting/src/RefittingAlgorithm.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/TrackFitting/RefittingAlgorithm.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/EventData/BoundTrackParameters.hpp"
0013 #include "Acts/EventData/MultiTrajectory.hpp"
0014 #include "Acts/EventData/SourceLink.hpp"
0015 #include "Acts/EventData/TrackProxy.hpp"
0016 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0017 #include "Acts/EventData/VectorTrackContainer.hpp"
0018 #include "Acts/Surfaces/Surface.hpp"
0019 #include "Acts/Utilities/Result.hpp"
0020 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0021 #include "ActsExamples/TrackFitting/RefittingCalibrator.hpp"
0022 #include "ActsExamples/TrackFitting/TrackFitterFunction.hpp"
0023
0024 #include <algorithm>
0025 #include <optional>
0026 #include <ostream>
0027 #include <stdexcept>
0028 #include <system_error>
0029 #include <utility>
0030 #include <vector>
0031
0032 namespace ActsExamples {
0033
0034 RefittingAlgorithm::RefittingAlgorithm(
0035     Config config, std::unique_ptr<const Acts::Logger> logger)
0036     : IAlgorithm("RefittingAlgorithm", std::move(logger)),
0037       m_cfg(std::move(config)) {
0038   if (m_cfg.inputTracks.empty()) {
0039     throw std::invalid_argument("Missing input tracks collection");
0040   }
0041   if (m_cfg.outputTracks.empty()) {
0042     throw std::invalid_argument("Missing output tracks collection");
0043   }
0044
0045   m_inputTracks.initialize(m_cfg.inputTracks);
0046   m_outputTracks.initialize(m_cfg.outputTracks);
0047 }
0048
0049 ProcessCode RefittingAlgorithm::execute(const AlgorithmContext& ctx) const {
0050   const auto& inputTracks = m_inputTracks(ctx);
0051
0052   auto trackContainer = std::make_shared<Acts::VectorTrackContainer>();
0053   auto trackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
0054   TrackContainer tracks(trackContainer, trackStateContainer);
0055
0056   // Perform the fit for each input track
0057   std::vector<Acts::SourceLink> trackSourceLinks;
0058   std::vector<const Acts::Surface*> surfSequence;
0059   RefittingCalibrator calibrator;
0060
0061   auto itrack = 0ul;
0062   for (const auto& track : inputTracks) {
0063     // Check if you are not in picking mode
0064     ++itrack;
0065     if (m_cfg.pickTrack > -1 &&
0066         static_cast<std::size_t>(m_cfg.pickTrack) != itrack - 1) {
0067       continue;
0068     }
0069
0070     if (!track.hasReferenceSurface()) {
0071       ACTS_VERBOSE("Skip track " << itrack << ": missing ref surface");
0072       continue;
0073     }
0074
0075     TrackFitterFunction::GeneralFitterOptions options{
0076         ctx.geoContext,
0077         ctx.magFieldContext,
0078         ctx.calibContext,
0079         &track.referenceSurface(),
0080         Acts::PropagatorPlainOptions(ctx.geoContext, ctx.magFieldContext),
0081         true};
0082
0083     Acts::BoundTrackParameters initialParams(
0084         track.referenceSurface().getSharedPtr(), track.parameters(),
0085         track.covariance(), track.particleHypothesis());
0086
0087     if (initialParams.covariance()) {
0088       for (auto i = 0ul; i < m_cfg.initialVarInflation.size(); ++i) {
0089         (*initialParams.covariance())(i, i) *= m_cfg.initialVarInflation.at(i);
0090       }
0091     }
0092
0093     trackSourceLinks.clear();
0094     surfSequence.clear();
0095
0096     for (auto state : track.trackStatesReversed()) {
0097       surfSequence.push_back(&state.referenceSurface());
0098
0099       if (!state.hasCalibrated()) {
0100         continue;
0101       }
0102
0103       auto sl = RefittingCalibrator::RefittingSourceLink{state};
0104       trackSourceLinks.push_back(Acts::SourceLink{sl});
0105     }
0106
0107     if (surfSequence.empty()) {
0108       ACTS_WARNING("Empty track " << itrack << " found.");
0109       continue;
0110     }
0111
0112     std::ranges::reverse(surfSequence);
0113
0114     ACTS_VERBOSE("Initial parameters: "
0115                  << initialParams.fourPosition(ctx.geoContext).transpose()
0116                  << " -> " << initialParams.direction().transpose());
0117
0118     ACTS_DEBUG("Invoke direct fitter for track " << itrack);
0119     auto result = (*m_cfg.fit)(trackSourceLinks, initialParams, options,
0120                                calibrator, surfSequence, tracks);
0121
0122     if (result.ok()) {
0123       // Get the fit output object
0124       const auto& refittedTrack = result.value();
0125       if (refittedTrack.hasReferenceSurface()) {
0126         ACTS_VERBOSE("Refitted parameters for track " << itrack);
0127         ACTS_VERBOSE("  " << track.parameters().transpose());
0128         ACTS_VERBOSE("Measurements: " << refittedTrack.nMeasurements());
0129         ACTS_VERBOSE("Outliers: " << refittedTrack.nOutliers());
0130       } else {
0131         ACTS_DEBUG("No refitted parameters for track " << itrack);
0132       }
0133     } else {
0134       ACTS_WARNING("Fit failed for event "
0135                    << ctx.eventNumber << " track " << itrack << " with error: "
0136                    << result.error() << ", " << result.error().message());
0137     }
0138     ++itrack;
0139   }
0140
0141   ACTS_DEBUG("Fitted tracks: " << trackContainer->size());
0142
0143   if (logger().doPrint(Acts::Logging::DEBUG)) {
0144     std::stringstream ss;
0145     trackStateContainer->statistics().toStream(ss);
0146     ACTS_DEBUG(ss.str());
0147   }
0148
0149   ConstTrackContainer constTracks{
0150       std::make_shared<Acts::ConstVectorTrackContainer>(
0151           std::move(*trackContainer)),
0152       std::make_shared<Acts::ConstVectorMultiTrajectory>(
0153           std::move(*trackStateContainer))};
0154
0155   m_outputTracks(ctx, std::move(constTracks));
0156   return ProcessCode::SUCCESS;
0157 }
0158
0159 }  // namespace ActsExamples