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