src/components/CKFTracking.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck
0003
0004 #include "CKFTracking.h"
0005
0006 // Gaudi
0007 #include "GaudiAlg/GaudiAlgorithm.h"
0008 #include "GaudiKernel/ToolHandle.h"
0009 #include "GaudiAlg/Transformer.h"
0010 #include "GaudiAlg/GaudiTool.h"
0011 #include "GaudiKernel/RndmGenerators.h"
0012 #include "Gaudi/Property.h"
0013
0014 #include "DDRec/CellIDPositionConverter.h"
0015 #include "DDRec/SurfaceManager.h"
0016 #include "DDRec/Surface.h"
0017
0018 #if Acts_VERSION_MAJOR < 36
0019 #include <Acts/EventData/Measurement.hpp>
0020 #endif
0021 #if Acts_VERSION_MAJOR >= 32
0022 #include "Acts/EventData/ProxyAccessor.hpp"
0023 #endif
0024 #include "Acts/Geometry/TrackingGeometry.hpp"
0025 #include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
0026 #include "Acts/Surfaces/PerigeeSurface.hpp"
0027
0028 #include "Acts/TrackFitting/GainMatrixSmoother.hpp"
0029 #include "Acts/TrackFitting/GainMatrixUpdater.hpp"
0030 #include "Acts/Propagator/EigenStepper.hpp"
0031 #if Acts_VERSION_MAJOR >= 34
0032 #if Acts_VERSION_MAJOR >= 37
0033 #include "Acts/Propagator/ActorList.hpp"
0034 #else
0035 #include "Acts/Propagator/AbortList.hpp"
0036 #endif
0037 #include "Acts/Propagator/EigenStepper.hpp"
0038 #include "Acts/Propagator/MaterialInteractor.hpp"
0039 #include "Acts/Propagator/Navigator.hpp"
0040 #endif
0041 #include "Acts/Propagator/Propagator.hpp"
0042 #if Acts_VERSION_MAJOR >= 34
0043 #include "Acts/Propagator/StandardAborters.hpp"
0044 #endif
0045 #include "Acts/Definitions/Common.hpp"
0046 #include "Acts/Utilities/Helpers.hpp"
0047 #include "Acts/Utilities/Logger.hpp"
0048 #if Acts_VERSION_MAJOR >= 34
0049 #include "Acts/Utilities/TrackHelpers.hpp"
0050 #endif
0051 #include "Acts/Definitions/Units.hpp"
0052
0053 #include <k4FWCore/DataHandle.h>
0054 #include <k4Interface/IGeoSvc.h>
0055 #include "JugTrack/IActsGeoSvc.h"
0056 #include "JugTrack/DD4hepBField.h"
0057
0058 #include "ActsExamples/EventData/GeometryContainers.hpp"
0059 #include "ActsExamples/EventData/Measurement.hpp"
0060 #include "ActsExamples/EventData/MeasurementCalibration.hpp"
0061 #include "ActsExamples/EventData/Index.hpp"
0062 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0063 #include "ActsExamples/EventData/Track.hpp"
0064
0065 #include <fmt/core.h>
0066 #include <fmt/ostream.h>
0067
0068 #include "edm4eic/TrackerHitCollection.h"
0069
0070 #include <functional>
0071 #include <stdexcept>
0072 #include <system_error>
0073 #include <type_traits>
0074 #include <vector>
0075 #include <random>
0076 #include <stdexcept>
0077
0078 template<> struct fmt::formatter<std::error_code> : fmt::ostream_formatter {};
0079
0080 static const std::map<int, Acts::Logging::Level> s_msgMap = {
0081     {MSG::DEBUG, Acts::Logging::DEBUG},
0082     {MSG::VERBOSE, Acts::Logging::VERBOSE},
0083     {MSG::INFO, Acts::Logging::INFO},
0084     {MSG::WARNING, Acts::Logging::WARNING},
0085     {MSG::FATAL, Acts::Logging::FATAL},
0086     {MSG::ERROR, Acts::Logging::ERROR},
0087 };
0088
0089 namespace Jug::Reco {
0090
0091   using namespace Acts::UnitLiterals;
0092
0093   CKFTracking::CKFTracking(const std::string& name, ISvcLocator* svcLoc)
0094       : GaudiAlgorithm(name, svcLoc)
0095   {
0096 #if Acts_VERSION_MAJOR < 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR < 1)
0097     declareProperty("inputSourceLinks", m_inputSourceLinks, "");
0098 #endif
0099     declareProperty("inputMeasurements", m_inputMeasurements, "");
0100     declareProperty("inputInitialTrackParameters", m_inputInitialTrackParameters, "");
0101     declareProperty("outputTracks", m_outputTracks, "");
0102     declareProperty("outputTrajectories", m_outputTrajectories, "");
0103   }
0104
0105   StatusCode CKFTracking::initialize()
0106   {
0107     if (GaudiAlgorithm::initialize().isFailure()) {
0108       return StatusCode::FAILURE;
0109     }
0110     m_geoSvc = service("GeoSvc");
0111     if (!m_geoSvc) {
0112       error() << "Unable to locate Geometry Service. "
0113               << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0114       return StatusCode::FAILURE;
0115     }
0116     m_actsGeoSvc = service("ActsGeoSvc");
0117     if (!m_actsGeoSvc) {
0118       error() << "Unable to locate ACTS Geometry Service. "
0119               << "Make sure you have ActsGeoSvc in the right place in the configuration." << endmsg;
0120       return StatusCode::FAILURE;
0121     }
0122
0123     m_BField   = std::dynamic_pointer_cast<const Jug::BField::DD4hepBField>(m_actsGeoSvc->getFieldProvider());
0124     m_fieldctx = Jug::BField::BFieldVariant(m_BField);
0125
0126     // eta bins, chi2 and #sourclinks per surface cutoffs
0127     m_sourcelinkSelectorCfg = {
0128         {Acts::GeometryIdentifier(),
0129             {m_etaBins, m_chi2CutOff,
0130                 {m_numMeasurementsCutOff.begin(), m_numMeasurementsCutOff.end()}
0131             }
0132         },
0133     };
0134     m_trackFinderFunc = CKFTracking::makeCKFTrackingFunction(m_actsGeoSvc->trackingGeometry(), m_BField);
0135     auto im = s_msgMap.find(msgLevel());
0136     if (im != s_msgMap.end()) {
0137         m_actsLoggingLevel = im->second;
0138     }
0139     return StatusCode::SUCCESS;
0140   }
0141
0142   StatusCode CKFTracking::execute()
0143   {
0144     // Read input data
0145 #if Acts_VERSION_MAJOR < 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR < 1)
0146     const auto* const src_links       = m_inputSourceLinks.get();
0147 #endif
0148     const auto* const init_trk_params = m_inputInitialTrackParameters.get();
0149     const auto* const measurements    = m_inputMeasurements.get();
0150
0151     //// Prepare the output data with MultiTrajectory
0152     auto* acts_trajectories = m_outputTrajectories.createAndPut();
0153     acts_trajectories->reserve(init_trk_params->size());
0154
0155     //// Construct a perigee surface as the target surface
0156     auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.});
0157
0158     ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("CKFTracking Logger", m_actsLoggingLevel));
0159
0160 #if Acts_VERSION_MAJOR >= 36
0161     Acts::PropagatorPlainOptions pOptions(m_geoctx, m_fieldctx);
0162 #else
0163     Acts::PropagatorPlainOptions pOptions;
0164 #endif
0165     pOptions.maxSteps = 10000;
0166
0167     ActsExamples::PassThroughCalibrator pcalibrator;
0168     ActsExamples::MeasurementCalibratorAdapter calibrator(pcalibrator, *measurements);
0169     Acts::GainMatrixUpdater kfUpdater;
0170 #if Acts_VERSION_MAJOR < 34
0171     Acts::GainMatrixSmoother kfSmoother;
0172 #endif
0173     Acts::MeasurementSelector measSel{m_sourcelinkSelectorCfg};
0174
0175 #if Acts_VERSION_MAJOR >= 36
0176     Acts::CombinatorialKalmanFilterExtensions<ActsExamples::TrackContainer>
0177         extensions;
0178 #else
0179     Acts::CombinatorialKalmanFilterExtensions<Acts::VectorMultiTrajectory>
0180         extensions;
0181 #endif
0182     extensions.calibrator.connect<
0183         &ActsExamples::MeasurementCalibratorAdapter::calibrate>(
0184         &calibrator);
0185 #if Acts_VERSION_MAJOR >= 36
0186     extensions.updater.connect<
0187         &Acts::GainMatrixUpdater::operator()<
0188         typename ActsExamples::TrackContainer::TrackStateContainerBackend>>(
0189         &kfUpdater);
0190 #else
0191     extensions.updater.connect<
0192         &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>(
0193         &kfUpdater);
0194 #endif
0195 #if Acts_VERSION_MAJOR < 34
0196     extensions.smoother.connect<
0197         &Acts::GainMatrixSmoother::operator()<Acts::VectorMultiTrajectory>>(
0198         &kfSmoother);
0199 #endif
0200     extensions.measurementSelector
0201         .connect<&Acts::MeasurementSelector::select<Acts::VectorMultiTrajectory>>(
0202             &measSel);
0203
0204     ActsExamples::IndexSourceLinkAccessor slAccessor;
0205 #if Acts_VERSION_MAJOR >= 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR >= 1)
0206     slAccessor.container = measurements->orderedIndices();
0207 #else
0208     slAccessor.container = src_links;
0209 #endif
0210     Acts::SourceLinkAccessorDelegate<ActsExamples::IndexSourceLinkAccessor::Iterator>
0211         slAccessorDelegate;
0212     slAccessorDelegate.connect<&ActsExamples::IndexSourceLinkAccessor::range>(&slAccessor);
0213
0214     // Set the CombinatorialKalmanFilter options
0215 #if Acts_VERSION_MAJOR < 34
0216     CKFTracking::TrackFinderOptions options(
0217         m_geoctx, m_fieldctx, m_calibctx, slAccessorDelegate,
0218         extensions, pOptions, &(*pSurface));
0219 #else
0220     CKFTracking::TrackFinderOptions options(
0221         m_geoctx, m_fieldctx, m_calibctx, slAccessorDelegate,
0222         extensions, pOptions);
0223 #endif
0224
0225 #if Acts_VERSION_MAJOR >= 36
0226     using Extrapolator = Acts::Propagator<Acts::EigenStepper<>, Acts::Navigator>;
0227 # if Acts_VERSION_MAJOR >= 37
0228     using ExtrapolatorOptions =
0229         Extrapolator::template Options<Acts::ActorList<Acts::MaterialInteractor,
0230                                                            Acts::EndOfWorldReached>>;
0231 # else
0232     using ExtrapolatorOptions =
0233         Extrapolator::template Options<Acts::ActionList<Acts::MaterialInteractor>,
0234                                            Acts::AbortList<Acts::EndOfWorldReached>>;
0235 # endif
0236     Extrapolator extrapolator(
0237         Acts::EigenStepper<>(m_BField),
0238         Acts::Navigator({m_actsGeoSvc->trackingGeometry()},
0239                         logger().cloneWithSuffix("Navigator")),
0240         logger().cloneWithSuffix("Propagator"));
0241     ExtrapolatorOptions extrapolationOptions(m_geoctx, m_fieldctx);
0242 #elif Acts_VERSION_MAJOR >= 34
0243     Acts::Propagator<Acts::EigenStepper<>, Acts::Navigator> extrapolator(
0244         Acts::EigenStepper<>(m_BField),
0245         Acts::Navigator({m_actsGeoSvc->trackingGeometry()},
0246                         logger().cloneWithSuffix("Navigator")),
0247         logger().cloneWithSuffix("Propagator"));
0248
0249     Acts::PropagatorOptions<Acts::ActionList<Acts::MaterialInteractor>,
0250                             Acts::AbortList<Acts::EndOfWorldReached>>
0251         extrapolationOptions(m_geoctx, m_fieldctx);
0252 #endif
0253
0254     // Create track container
0255     auto trackContainer = std::make_shared<Acts::VectorTrackContainer>();
0256     auto trackStateContainer = std::make_shared<Acts::VectorMultiTrajectory>();
0257     ActsExamples::TrackContainer tracks(trackContainer, trackStateContainer);
0258
0259     // Add seed number column
0260     tracks.addColumn<unsigned int>("seed");
0261 #if Acts_VERSION_MAJOR >= 32
0262     Acts::ProxyAccessor<unsigned int> seedNumber("seed");
0263 #else
0264     Acts::TrackAccessor<unsigned int> seedNumber("seed");
0265 #endif
0266
0267     // Loop over seeds
0268     for (std::size_t iseed = 0; iseed < init_trk_params->size(); ++iseed) {
0269         auto result =
0270             (*m_trackFinderFunc)(init_trk_params->at(iseed), options, tracks);
0271
0272         if (!result.ok()) {
0273             debug() << fmt::format("Track finding failed for seed {} with error {}", iseed, result.error()) << endmsg;
0274             continue;
0275         }
0276
0277         // Set seed number for all found tracks
0278         auto& tracksForSeed = result.value();
0279         for (auto& track : tracksForSeed) {
0280
0281 #if Acts_VERSION_MAJOR >=34
0282             auto smoothingResult = Acts::smoothTrack(m_geoctx, track, logger());
0283             if (!smoothingResult.ok()) {
0284                 ACTS_ERROR("Smoothing for seed "
0285                     << iseed << " and track " << track.index()
0286                     << " failed with error " << smoothingResult.error());
0287                 continue;
0288             }
0289
0290             auto extrapolationResult = Acts::extrapolateTrackToReferenceSurface(
0291                 track, *pSurface, extrapolator, extrapolationOptions,
0292                 Acts::TrackExtrapolationStrategy::firstOrLast, logger());
0293             if (!extrapolationResult.ok()) {
0294                 ACTS_ERROR("Extrapolation for seed "
0295                     << iseed << " and track " << track.index()
0296                     << " failed with error " << extrapolationResult.error());
0297                 continue;
0298             }
0299 #endif
0300
0301             seedNumber(track) = iseed;
0302         }
0303     }
0304
0305     // Move track states and track container to const containers
0306     // NOTE Using the non-const containers leads to references to
0307     // implicitly converted temporaries inside the Trajectories.
0308     auto constTrackStateContainer =
0309         std::make_shared<Acts::ConstVectorMultiTrajectory>(
0310             std::move(*trackStateContainer));
0311
0312     auto constTrackContainer =
0313         std::make_shared<Acts::ConstVectorTrackContainer>(
0314             std::move(*trackContainer));
0315
0316     auto* constTracks = m_outputTracks.put(
0317         std::make_unique<ActsExamples::ConstTrackContainer>(constTrackContainer, constTrackStateContainer)
0318     );
0319
0320     // Seed number column accessor
0321 #if Acts_VERSION_MAJOR >= 32
0322     const Acts::ConstProxyAccessor<unsigned int> constSeedNumber("seed");
0323 #else
0324     const Acts::ConstTrackAccessor<unsigned int> constSeedNumber("seed");
0325 #endif
0326
0327     // Prepare the output data with MultiTrajectory, per seed
0328     ActsExamples::Trajectories::IndexedParameters parameters;
0329     std::vector<Acts::MultiTrajectoryTraits::IndexType> tips;
0330
0331     std::optional<unsigned int> lastSeed;
0332     for (const auto& track : *constTracks) {
0333       if (!lastSeed) {
0334         lastSeed = constSeedNumber(track);
0335       }
0336
0337       if (constSeedNumber(track) != lastSeed.value()) {
0338         // make copies and clear vectors
0339         acts_trajectories->emplace_back(
0340           constTracks->trackStateContainer(),
0341           tips, parameters);
0342
0343         tips.clear();
0344         parameters.clear();
0345       }
0346
0347       lastSeed = constSeedNumber(track);
0348
0349       tips.push_back(track.tipIndex());
0350       parameters.emplace(
0351           std::pair{track.tipIndex(),
0352                     ActsExamples::TrackParameters{track.referenceSurface().getSharedPtr(),
0353                                                   track.parameters(), track.covariance(),
0354                                                   track.particleHypothesis()}});
0355     }
0356
0357     if (tips.empty()) {
0358       info() << fmt::format("Last trajectory is empty") << endmsg;
0359     }
0360
0361     // last entry: move vectors
0362     acts_trajectories->emplace_back(
0363       constTracks->trackStateContainer(),
0364       std::move(tips), std::move(parameters));
0365
0366     return StatusCode::SUCCESS;
0367   }
0368
0369   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0370   DECLARE_COMPONENT(CKFTracking)
0371 } // namespace Jug::Reco