acts/src/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 #include "Acts/Geometry/TrackingGeometry.hpp"
0019 #include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
0020 #include "Acts/Surfaces/PerigeeSurface.hpp"
0021
0022 #include "Acts/TrackFitting/GainMatrixSmoother.hpp"
0023 #include "Acts/TrackFitting/GainMatrixUpdater.hpp"
0024 #include "Acts/Propagator/EigenStepper.hpp"
0025 #include "Acts/Propagator/Navigator.hpp"
0026 #include "Acts/Propagator/Propagator.hpp"
0027 #include "Acts/Definitions/Common.hpp"
0028 #include "Acts/Utilities/Helpers.hpp"
0029 #include "Acts/Utilities/Logger.hpp"
0030 #include "Acts/Definitions/Units.hpp"
0031
0032 #include "JugBase/DataHandle.h"
0033 #include "JugBase/IGeoSvc.h"
0034 #include "JugBase/BField/DD4hepBField.h"
0035
0036 #include "JugTrack/GeometryContainers.hpp"
0037 #include "JugTrack/Measurement.hpp"
0038 #include "JugTrack/Index.hpp"
0039 #include "JugTrack/IndexSourceLink.hpp"
0040 #include "JugTrack/Track.hpp"
0041
0042
0043 #include "edm4eic/TrackerHitCollection.h"
0044
0045 #include <functional>
0046 #include <stdexcept>
0047 #include <vector>
0048 #include <random>
0049 #include <stdexcept>
0050
0051
0052 static const std::map<int, Acts::Logging::Level> s_msgMap = {
0053     {MSG::DEBUG, Acts::Logging::DEBUG},
0054     {MSG::VERBOSE, Acts::Logging::VERBOSE},
0055     {MSG::INFO, Acts::Logging::INFO},
0056     {MSG::WARNING, Acts::Logging::WARNING},
0057     {MSG::FATAL, Acts::Logging::FATAL},
0058     {MSG::ERROR, Acts::Logging::ERROR},
0059 };
0060
0061 namespace Jug::Reco {
0062
0063   using namespace Acts::UnitLiterals;
0064
0065   CKFTracking::CKFTracking(const std::string& name, ISvcLocator* svcLoc)
0066       : GaudiAlgorithm(name, svcLoc)
0067   {
0068     declareProperty("inputSourceLinks", m_inputSourceLinks, "");
0069     declareProperty("inputMeasurements", m_inputMeasurements, "");
0070     declareProperty("inputInitialTrackParameters", m_inputInitialTrackParameters, "");
0071     declareProperty("outputTrajectories", m_outputTrajectories, "");
0072   }
0073
0074   StatusCode CKFTracking::initialize()
0075   {
0076     if (GaudiAlgorithm::initialize().isFailure()) {
0077       return StatusCode::FAILURE;
0078     }
0079     m_geoSvc = service("GeoSvc");
0080     if (!m_geoSvc) {
0081       error() << "Unable to locate Geometry Service. "
0082               << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0083       return StatusCode::FAILURE;
0084     }
0085
0086     m_BField   = std::dynamic_pointer_cast<const Jug::BField::DD4hepBField>(m_geoSvc->getFieldProvider());
0087     m_fieldctx = Jug::BField::BFieldVariant(m_BField);
0088
0089     // eta bins, chi2 and #sourclinks per surface cutoffs
0090     m_sourcelinkSelectorCfg = {
0091         {Acts::GeometryIdentifier(),
0092             {m_etaBins, m_chi2CutOff,
0093                 {m_numMeasurementsCutOff.begin(), m_numMeasurementsCutOff.end()}
0094             }
0095         },
0096     };
0097     m_trackFinderFunc = CKFTracking::makeCKFTrackingFunction(m_geoSvc->trackingGeometry(), m_BField);
0098     auto im = s_msgMap.find(msgLevel());
0099     if (im != s_msgMap.end()) {
0100         m_actsLoggingLevel = im->second;
0101     }
0102     return StatusCode::SUCCESS;
0103   }
0104
0105   StatusCode CKFTracking::execute()
0106   {
0107     // Read input data
0108     const auto* const src_links       = m_inputSourceLinks.get();
0109     const auto* const init_trk_params = m_inputInitialTrackParameters.get();
0110     const auto* const measurements    = m_inputMeasurements.get();
0111
0112     //// Prepare the output data with MultiTrajectory
0113     // TrajectoryContainer trajectories;
0114     auto* trajectories = m_outputTrajectories.createAndPut();
0115     trajectories->reserve(init_trk_params->size());
0116
0117     //// Construct a perigee surface as the target surface
0118     auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.});
0119
0120     ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("CKFTracking Logger", m_actsLoggingLevel));
0121
0122     Acts::PropagatorPlainOptions pOptions;
0123     pOptions.maxSteps = 10000;
0124
0125     MeasurementCalibrator calibrator{*measurements};
0126     Acts::GainMatrixUpdater kfUpdater;
0127     Acts::GainMatrixSmoother kfSmoother;
0128     Acts::MeasurementSelector measSel{m_sourcelinkSelectorCfg};
0129
0130     Acts::CombinatorialKalmanFilterExtensions<Acts::VectorMultiTrajectory>
0131         extensions;
0132     extensions.calibrator.connect<&MeasurementCalibrator::calibrate>(&calibrator);
0133     extensions.updater.connect<
0134         &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>(
0135         &kfUpdater);
0136     extensions.smoother.connect<
0137         &Acts::GainMatrixSmoother::operator()<Acts::VectorMultiTrajectory>>(
0138         &kfSmoother);
0139     extensions.measurementSelector
0140         .connect<&Acts::MeasurementSelector::select<Acts::VectorMultiTrajectory>>(
0141             &measSel);
0142
0143     IndexSourceLinkAccessor slAccessor;
0144     slAccessor.container = src_links;
0145     Acts::SourceLinkAccessorDelegate<IndexSourceLinkAccessor::Iterator>
0146         slAccessorDelegate;
0147     slAccessorDelegate.connect<&IndexSourceLinkAccessor::range>(&slAccessor);
0148
0149     // Set the CombinatorialKalmanFilter options
0150     CKFTracking::TrackFinderOptions options(
0151         m_geoctx, m_fieldctx, m_calibctx, slAccessorDelegate,
0152         extensions, Acts::LoggerWrapper{logger()}, pOptions, &(*pSurface));
0153
0154     auto results = (*m_trackFinderFunc)(*init_trk_params, options);
0155
0156     for (std::size_t iseed = 0; iseed < init_trk_params->size(); ++iseed) {
0157
0158       auto& result = results[iseed];
0159
0160       if (result.ok()) {
0161         // Get the track finding output object
0162         auto& trackFindingOutput = result.value();
0163         // Create a SimMultiTrajectory
0164         trajectories->emplace_back(std::move(trackFindingOutput.fittedStates),
0165                                    std::move(trackFindingOutput.lastMeasurementIndices),
0166                                    std::move(trackFindingOutput.fittedParameters));
0167       } else {
0168         if (msgLevel(MSG::DEBUG)) {
0169           debug() << "Track finding failed for truth seed " << iseed << "with error: " << result.error() << endmsg;
0170         }
0171       }
0172     }
0173
0174     return StatusCode::SUCCESS;
0175   }
0176
0177   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0178   DECLARE_COMPONENT(CKFTracking)
0179 } // namespace Jug::Reco