acts/src/TrackFindingAlgorithm.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck
0003
0004 #include "TrackFindingAlgorithm.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   TrackFindingAlgorithm::TrackFindingAlgorithm(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 TrackFindingAlgorithm::initialize()
0075   {
0076     warning() << "Deprecated algorithm, use CKFTracking instead" << endmsg;
0077     if (GaudiAlgorithm::initialize().isFailure()) {
0078       return StatusCode::FAILURE;
0079     }
0080     m_geoSvc = service("GeoSvc");
0081     if (!m_geoSvc) {
0082       error() << "Unable to locate Geometry Service. "
0083               << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0084       return StatusCode::FAILURE;
0085     }
0086
0087     m_BField   = std::dynamic_pointer_cast<const Jug::BField::DD4hepBField>(m_geoSvc->getFieldProvider());
0088     m_fieldctx = Jug::BField::BFieldVariant(m_BField);
0089
0090     // eta bins, chi2 and #sourclinks per surface cutoffs
0091     m_sourcelinkSelectorCfg = {
0092         {Acts::GeometryIdentifier(),
0093             {m_etaBins, m_chi2CutOff,
0094                 {m_numMeasurementsCutOff.begin(), m_numMeasurementsCutOff.end()}
0095             }
0096         },
0097     };
0098     m_trackFinderFunc = TrackFindingAlgorithm::makeTrackFinderFunction(m_geoSvc->trackingGeometry(), m_BField);
0099     auto im = s_msgMap.find(msgLevel());
0100     if (im != s_msgMap.end()) {
0101         m_actsLoggingLevel = im->second;
0102     }
0103     return StatusCode::SUCCESS;
0104   }
0105
0106   StatusCode TrackFindingAlgorithm::execute()
0107   {
0108     // Read input data
0109     const auto* const src_links       = m_inputSourceLinks.get();
0110     const auto* const init_trk_params = m_inputInitialTrackParameters.get();
0111     const auto* const measurements    = m_inputMeasurements.get();
0112
0113     //// Prepare the output data with MultiTrajectory
0114     // TrajectoryContainer trajectories;
0115     auto* trajectories = m_outputTrajectories.createAndPut();
0116     trajectories->reserve(init_trk_params->size());
0117
0118     //// Construct a perigee surface as the target surface
0119     auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.});
0120
0121     ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("TrackFindingAlgorithm Logger", m_actsLoggingLevel));
0122
0123     Acts::PropagatorPlainOptions pOptions;
0124     pOptions.maxSteps = 10000;
0125
0126     MeasurementCalibrator calibrator{*measurements};
0127     Acts::GainMatrixUpdater kfUpdater;
0128     Acts::GainMatrixSmoother kfSmoother;
0129     Acts::MeasurementSelector measSel{m_sourcelinkSelectorCfg};
0130
0131     Acts::CombinatorialKalmanFilterExtensions<Acts::VectorMultiTrajectory> 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     TrackFindingAlgorithm::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(TrackFindingAlgorithm)
0179
0180 } // namespace Jug::Reco