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/Definitions/Units.hpp"
0013 #include "Acts/EventData/BoundTrackParameters.hpp"
0014 #include "Acts/EventData/MultiTrajectory.hpp"
0015 #include "Acts/EventData/SourceLink.hpp"
0016 #include "Acts/EventData/TrackProxy.hpp"
0017 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0018 #include "Acts/EventData/VectorTrackContainer.hpp"
0019 #include "Acts/Surfaces/PerigeeSurface.hpp"
0020 #include "Acts/Surfaces/Surface.hpp"
0021 #include "Acts/Utilities/Result.hpp"
0022 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0023 #include "ActsExamples/TrackFitting/RefittingCalibrator.hpp"
0024 #include "ActsExamples/TrackFitting/TrackFitterFunction.hpp"
0025
0026 #include <algorithm>
0027 #include <optional>
0028 #include <ostream>
0029 #include <stdexcept>
0030 #include <system_error>
0031 #include <utility>
0032 #include <vector>
0033
0034 namespace ActsExamples {
0035
0036 RefittingAlgorithm::RefittingAlgorithm(
0037     Config config, std::unique_ptr<const Acts::Logger> logger)
0038     : IAlgorithm("RefittingAlgorithm", std::move(logger)),
0039       m_cfg(std::move(config)) {
0040   if (m_cfg.inputTracks.empty()) {
0041     throw std::invalid_argument("Missing input tracks collection");
0042   }
0043   if (m_cfg.outputTracks.empty()) {
0044     throw std::invalid_argument("Missing output tracks collection");
0045   }
0046
0047   m_inputTracks.initialize(m_cfg.inputTracks);
0048   m_outputTracks.initialize(m_cfg.outputTracks);
0049 }
0050
0051 ProcessCode RefittingAlgorithm::execute(const AlgorithmContext& ctx) const {
0052   const auto& inputTracks = m_inputTracks(ctx);
0053
0054   auto trackContainer = std::make_shared<Acts::VectorTrackContainer>();
0055   auto trackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
0056   TrackContainer tracks(trackContainer, trackStateContainer);
0057
0058   auto perigeeSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
0059       Acts::Vector3{0., 0., 0.});
0060
0061   // The following code is only necessary if a beamspot constraint is in use but
0062   // unguarded for lifetime and simplicity reasons. The Core KF does not support
0063   // a beamspot constraint by itself but it can be achieved by injecting a fake
0064   // measurement on the perigee surface and processing it first in the fitter.
0065   // This can lead to problems if the beamspot constraint moves the initial
0066   // parameters too far and we start to miss following surfaces and
0067   // measurements. But in a refitting scenario with increased bounds and direct
0068   // navigation it *should* work OK.
0069
0070   auto beamSpotVectorTrackStateContainer =
0071       std::make_shared<Acts::VectorMultiTrajectory>();
0072   auto beamSpotTrackState = beamSpotVectorTrackStateContainer->makeTrackState();
0073
0074   const Acts::Vector2 beamSpotMeasValue{0., 0.};
0075
0076   beamSpotTrackState.setReferenceSurface(perigeeSurface);
0077
0078   if (m_cfg.beamSpotConstraint.has_value()) {
0079     ACTS_DEBUG("Using provided beam spot constraint matrix");
0080     beamSpotTrackState.allocateCalibrated(beamSpotMeasValue,
0081                                           m_cfg.beamSpotConstraint.value());
0082   } else {
0083     ACTS_DEBUG("No beam spot constraint provided, using zero matrix");
0084     beamSpotTrackState.allocateCalibrated(beamSpotMeasValue,
0085                                           Acts::SquareMatrix2::Zero());
0086   }
0087
0088   Acts::SourceLink testSL{42};
0089   beamSpotTrackState.setUncalibratedSourceLink(std::move(testSL));
0090
0091   auto beamSpotConstVectorTrackStateContainer =
0092       std::make_shared<Acts::ConstVectorMultiTrajectory>(
0093           std::move(*beamSpotVectorTrackStateContainer));
0094
0095   auto beamSpotConstTrackState =
0096       beamSpotConstVectorTrackStateContainer->getTrackState(
0097           beamSpotTrackState.index());
0098
0099   // Perform the fit for each input track
0100   std::vector<Acts::SourceLink> trackSourceLinks;
0101   std::vector<const Acts::Surface*> surfSequence;
0102   RefittingCalibrator calibrator;
0103
0104   auto itrack = 0ul;
0105   for (const auto& track : inputTracks) {
0106     // Check if you are not in picking mode
0107     ++itrack;
0108     if (m_cfg.pickTrack > -1 &&
0109         static_cast<std::size_t>(m_cfg.pickTrack) != itrack - 1) {
0110       continue;
0111     }
0112
0113     if (!track.hasReferenceSurface()) {
0114       ACTS_VERBOSE("Skip track " << itrack << ": missing ref surface");
0115       continue;
0116     }
0117
0118     TrackFitterFunction::GeneralFitterOptions options{
0119         ctx.geoContext,
0120         ctx.magFieldContext,
0121         ctx.calibContext,
0122         perigeeSurface.get(),
0123         Acts::PropagatorPlainOptions(ctx.geoContext, ctx.magFieldContext),
0124         true};
0125
0126     Acts::BoundTrackParameters initialParams(
0127         track.referenceSurface().getSharedPtr(), track.parameters(),
0128         track.covariance(), track.particleHypothesis());
0129
0130     if (initialParams.covariance()) {
0131       for (auto i = 0ul; i < m_cfg.initialVarInflation.size(); ++i) {
0132         (*initialParams.covariance())(i, i) *= m_cfg.initialVarInflation.at(i);
0133       }
0134     }
0135
0136     trackSourceLinks.clear();
0137     surfSequence.clear();
0138
0139     for (auto state : track.trackStatesReversed()) {
0140       surfSequence.push_back(&state.referenceSurface());
0141
0142       if (!state.hasCalibrated()) {
0143         continue;
0144       }
0145
0146       auto sl = RefittingCalibrator::RefittingSourceLink{state};
0147       trackSourceLinks.push_back(Acts::SourceLink{sl});
0148     }
0149
0150     if (surfSequence.empty()) {
0151       ACTS_DEBUG("Empty track " << itrack << " found.");
0152       continue;
0153     }
0154
0155     if (m_cfg.beamSpotConstraint.has_value()) {
0156       RefittingCalibrator::RefittingSourceLink beamSpotSL{
0157           beamSpotConstTrackState};
0158       trackSourceLinks.emplace_back(Acts::SourceLink{beamSpotSL});
0159       surfSequence.push_back(perigeeSurface.get());
0160     }
0161
0162     std::ranges::reverse(surfSequence);
0163
0164     ACTS_VERBOSE("Initial parameters: "
0165                  << initialParams.fourPosition(ctx.geoContext).transpose()
0166                  << " -> " << initialParams.direction().transpose());
0167
0168     ACTS_DEBUG("Invoke direct fitter for track " << itrack);
0169     auto result = (*m_cfg.fit)(trackSourceLinks, initialParams, options,
0170                                calibrator, surfSequence, tracks);
0171
0172     if (result.ok()) {
0173       // Get the fit output object
0174       const auto& refittedTrack = result.value();
0175       if (refittedTrack.hasReferenceSurface()) {
0176         ACTS_VERBOSE("Refitted parameters for track " << itrack);
0177         ACTS_VERBOSE("  " << track.parameters().transpose());
0178         ACTS_VERBOSE("Measurements: " << refittedTrack.nMeasurements());
0179         ACTS_VERBOSE("Outliers: " << refittedTrack.nOutliers());
0180       } else {
0181         ACTS_DEBUG("No refitted parameters for track " << itrack);
0182       }
0183     } else {
0184       ACTS_DEBUG("Fit failed for event "
0185                  << ctx.eventNumber << " track " << itrack << " with error: "
0186                  << result.error() << ", " << result.error().message());
0187     }
0188     ++itrack;
0189   }
0190
0191   ACTS_DEBUG("Fitted tracks: " << trackContainer->size());
0192
0193   if (logger().doPrint(Acts::Logging::DEBUG)) {
0194     std::stringstream ss;
0195     trackStateContainer->statistics().toStream(ss);
0196     ACTS_DEBUG(ss.str());
0197   }
0198
0199   ConstTrackContainer constTracks{
0200       std::make_shared<Acts::ConstVectorTrackContainer>(
0201           std::move(*trackContainer)),
0202       std::make_shared<Acts::ConstVectorMultiTrajectory>(
0203           std::move(*trackStateContainer))};
0204
0205   m_outputTracks(ctx, std::move(constTracks));
0206   return ProcessCode::SUCCESS;
0207 }
0208
0209 }  // namespace ActsExamples