components/disabled/TrackFittingAlgorithm.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck
0003
0004 //
0005 #include "TrackFittingAlgorithm.h"
0006
0007 // Gaudi
0008 #include "Gaudi/Algorithm.h"
0009 #include "GaudiKernel/ToolHandle.h"
0010 #include "GaudiKernel/RndmGenerators.h"
0011 #include "Gaudi/Property.h"
0012
0013 #include "DDRec/CellIDPositionConverter.h"
0014 #include "DDRec/SurfaceManager.h"
0015 #include "DDRec/Surface.h"
0016
0017 #include "Acts/Geometry/TrackingGeometry.hpp"
0018 #include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
0019 #include "Acts/Surfaces/PerigeeSurface.hpp"
0020 #include "Acts/Propagator/EigenStepper.hpp"
0021 #include "Acts/Propagator/Navigator.hpp"
0022 #include "Acts/Propagator/Propagator.hpp"
0023 #include "Acts/Definitions/Common.hpp"
0024 #include "Acts/Utilities/Helpers.hpp"
0025 #include "Acts/Utilities/Logger.hpp"
0026 #include "Acts/Definitions/Units.hpp"
0027
0028 #include <k4FWCore/DataHandle.h>
0029 #include <k4Interface/IGeoSvc.h>
0030 #include "JugTrack/IActsGeoSvc.h"
0031 #include "JugTrack/DD4hepBField.h"
0032
0033 #include "ActsExamples/EventData/GeometryContainers.hpp"
0034 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0035 #include "ActsExamples/EventData/Track.hpp"
0036 #include "ActsExamples/EventData/Measurement.hpp"
0037 #include "ActsExamples/EventData/MeasurementCalibration.hpp"
0038
0039 #include "edm4eic/TrackerHitCollection.h"
0040
0041 #include <functional>
0042 #include <stdexcept>
0043 #include <vector>
0044 #include <random>
0045 #include <stdexcept>
0046
0047 namespace Jug::Reco {
0048
0049   using namespace Acts::UnitLiterals;
0050
0051   TrackFittingAlgorithm::TrackFittingAlgorithm(const std::string& name, ISvcLocator* svcLoc)
0052       : Gaudi::Algorithm(name, svcLoc)
0053   {
0054     declareProperty("inputSourceLinks", m_inputSourceLinks, "");
0055     declareProperty("initialTrackParameters", m_initialTrackParameters, "");
0056     declareProperty("inputMeasurements", m_inputMeasurements, "");
0057     declareProperty("inputProtoTracks", m_inputProtoTracks, "");
0058     declareProperty("foundTracks", m_foundTracks, "");
0059     declareProperty("outputTrajectories", m_outputTrajectories, "");
0060   }
0061
0062   StatusCode TrackFittingAlgorithm::initialize()
0063   {
0064     if (Gaudi::Algorithm::initialize().isFailure()) {
0065       return StatusCode::FAILURE;
0066     }
0067     m_geoSvc = service("GeoSvc");
0068     if (!m_geoSvc) {
0069       error() << "Unable to locate Geometry Service. "
0070               << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0071       return StatusCode::FAILURE;
0072     }
0073     m_actsGeoSvc = service("ActsGeoSvc");
0074     if (!m_actsGeoSvc) {
0075       error() << "Unable to locate ACTS Geometry Service. "
0076               << "Make sure you have ActsGeoSvc in the right place in the configuration." << endmsg;
0077       return StatusCode::FAILURE;
0078     }
0079
0080     m_BField   = std::dynamic_pointer_cast<const Jug::BField::DD4hepBField>(m_actsGeoSvc->getFieldProvider());
0081     m_fieldctx = Jug::BField::BFieldVariant(m_BField);
0082
0083     m_trackFittingFunc = makeTrackFittingFunction(m_actsGeoSvc->trackingGeometry(), m_BField);
0084     return StatusCode::SUCCESS;
0085   }
0086
0087   StatusCode TrackFittingAlgorithm::execute(const EventContext&) const
0088   {
0089     // Read input data
0090     const auto* const sourceLinks       = m_inputSourceLinks.get();
0091     const auto* const initialParameters = m_initialTrackParameters.get();
0092     const auto* const measurements      = m_inputMeasurements.get();
0093     const auto* const protoTracks       = m_inputProtoTracks.get();
0094     ACTS_LOCAL_LOGGER(Acts::getDefaultLogger("TrackFittingAlgorithm Logger", Acts::Logging::INFO));
0095
0096     // Consistency cross checks
0097     if (protoTracks->size() != initialParameters->size()) {
0098       ACTS_FATAL("Inconsistent number of proto tracks and initial parameters");
0099       return StatusCode::FAILURE;
0100     }
0101
0102     // TrajectoryContainer trajectories;
0103     auto* trajectories = m_outputTrajectories.createAndPut();
0104     trajectories->reserve(initialParameters->size());
0105
0106     // Construct a perigee surface as the target surface
0107     auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0., 0., 0.});
0108
0109     Acts::PropagatorPlainOptions pOptions;
0110     pOptions.maxSteps = 10000;
0111     // kfOptions.multipleScattering = m_cfg.multipleScattering;
0112     // kfOptions.energyLoss         = m_cfg.energyLoss;
0113
0114     Acts::KalmanFitterExtensions<Acts::VectorMultiTrajectory> extensions;
0115     ActsExamples::PassThroughCalibrator pcalibrator;
0116     ActsExamples::MeasurementCalibratorAdapter calibrator(pcalibrator, *measurements);
0117     extensions.calibrator.connect<&ActsExamples::MeasurementCalibratorAdapter::calibrate>(
0118         &calibrator);
0119     Acts::GainMatrixUpdater kfUpdater;
0120     Acts::GainMatrixSmoother kfSmoother;
0121     extensions.updater.connect<&Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>(
0122         &kfUpdater);
0123     extensions.smoother.connect<&Acts::GainMatrixSmoother::operator()<Acts::VectorMultiTrajectory>>(
0124         &kfSmoother);
0125
0126     Acts::KalmanFitterOptions kfOptions(
0127         m_geoctx, m_fieldctx, m_calibctx, extensions,
0128         Acts::PropagatorPlainOptions(),
0129         &(*pSurface));
0130
0131     // used for processing the data
0132     std::vector<ActsExamples::IndexSourceLink> trackSourceLinks;
0133     std::vector<const Acts::Surface*> surfSequence;
0134
0135     if (msgLevel(MSG::DEBUG)) {
0136       debug() << "initialParams size:  " << initialParameters->size() << endmsg;
0137       debug() << "measurements size:  " << measurements->size() << endmsg;
0138     }
0139
0140     // Perform the track finding for each starting parameter
0141     // @TODO: use seeds from track seeding algorithm as starting parameter
0142     // initial track params and proto tracks might likely have the same size.
0143     //for (std::size_t iseed = 0; iseed < init_trk_params->size(); ++iseed) {
0144     for (std::size_t itrack = 0; itrack < (*protoTracks).size(); ++itrack) {
0145
0146       const auto& protoTrack    = (*protoTracks)[itrack];
0147       const auto& initialParams = (*initialParameters)[itrack];
0148
0149       if (msgLevel(MSG::DEBUG)) {
0150         debug() << "protoTrack size:  " << protoTrack.size() << endmsg;
0151         debug() << "sourceLinks size:  " << sourceLinks->size() << endmsg;
0152       }
0153
0154       trackSourceLinks.clear();
0155       trackSourceLinks.reserve(protoTrack.size());
0156
0157       for (auto hitIndex : protoTrack) {
0158         if (auto it = sourceLinks->nth(hitIndex); it != sourceLinks->end()) {
0159           const ActsExamples::IndexSourceLink& sourceLink = *it;
0160           trackSourceLinks.push_back(std::cref(sourceLink));
0161         } else {
0162           ACTS_FATAL("Proto track " << itrack << " contains invalid hit index"
0163                                     << hitIndex);
0164           return StatusCode::FAILURE;
0165         }
0166       }
0167
0168       if (msgLevel(MSG::DEBUG)) {
0169         debug() << "Invoke track fitting ...  " << itrack << endmsg;
0170       }
0171       auto result = fitTrack(trackSourceLinks, initialParams, kfOptions);
0172       if (msgLevel(MSG::DEBUG)) {
0173         debug() << "fitting done." << endmsg;
0174       }
0175       if (result.ok())
0176       {
0177         // Get the fit output object
0178         auto& fitOutput = result.value();
0179         // The track entry indices container. One element here.
0180         std::vector<Acts::MultiTrajectoryTraits::IndexType> trackTips;
0181         trackTips.reserve(1);
0182         trackTips.emplace_back(fitOutput.lastMeasurementIndex);
0183         // The fitted parameters container. One element (at most) here.
0184         ActsExamples::Trajectories::IndexedParameters indexedParams;
0185         if (fitOutput.fittedParameters) {
0186           const auto& params = fitOutput.fittedParameters.value();
0187           ACTS_VERBOSE("Fitted paramemeters for track " << itrack);
0188           ACTS_VERBOSE("  " << params.parameters().transpose());
0189           // Push the fitted parameters to the container
0190           indexedParams.emplace(fitOutput.lastMeasurementIndex, params);
0191         } else {
0192           ACTS_DEBUG("No fitted paramemeters for track " << itrack);
0193         }
0194         // store the result
0195         trajectories->emplace_back(std::move(fitOutput.fittedStates), std::move(trackTips),
0196                                   std::move(indexedParams));
0197       } else {
0198         ACTS_WARNING("Fit failed for track " << itrack << " with error" << result.error());
0199         // Fit failed. Add an empty result so the output container has
0200         // the same number of entries as the input.
0201         trajectories->push_back(ActsExamples::Trajectories());
0202       }
0203     }
0204
0205     return StatusCode::SUCCESS;
0206   }
0207
0208   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0209   DECLARE_COMPONENT(TrackFittingAlgorithm)
0210
0211 } // namespace Jug::Reco