From aee1808efa1bcf43c1873763f28ecdda169ddeba Mon Sep 17 00:00:00 2001 From: Dmitry Romanov Date: Tue, 19 Nov 2024 10:28:33 -0500 Subject: [PATCH] CKF Tracking --- source/tdis/CMakeLists.txt | 3 +- source/tdis/io/PodioWriteProcessor.hpp | 4 +- source/tdis/tdis_main.cpp | 10 +- source/tdis/tracking/CKFTracking.cc | 375 ++++++++++++++++++ source/tdis/tracking/CKFTracking.h | 91 +++++ .../tracking/TruthTrackParameterFactory.h | 163 +++++--- 6 files changed, 589 insertions(+), 57 deletions(-) create mode 100644 source/tdis/tracking/CKFTracking.cc create mode 100644 source/tdis/tracking/CKFTracking.h diff --git a/source/tdis/CMakeLists.txt b/source/tdis/CMakeLists.txt index 2c17dec..6b2323e 100644 --- a/source/tdis/CMakeLists.txt +++ b/source/tdis/CMakeLists.txt @@ -50,7 +50,8 @@ add_executable(tdis # tracking/Measurement2DFactory.h tracking/TruthTrackParameterFactory.h tracking/KalmanFittingFactory.h - + tracking/CKFTracking.h + tracking/CKFTracking.cc ) # ---------- FIND REQUIRED PACKAGES ------------- diff --git a/source/tdis/io/PodioWriteProcessor.hpp b/source/tdis/io/PodioWriteProcessor.hpp index 1278d51..ba11c75 100644 --- a/source/tdis/io/PodioWriteProcessor.hpp +++ b/source/tdis/io/PodioWriteProcessor.hpp @@ -75,7 +75,9 @@ class PodioWriteProcessor : public JEventProcessor { "DigitizedMtpcMcHit", // Digitized hits - "TrackerHit" + "TrackerHit", + "Measurement2D", + "TruthTrackInitParameters", }; PodioWriteProcessor(JApplication * app); diff --git a/source/tdis/tdis_main.cpp b/source/tdis/tdis_main.cpp index 0937994..33de276 100644 --- a/source/tdis/tdis_main.cpp +++ b/source/tdis/tdis_main.cpp @@ -9,9 +9,9 @@ // #include // #include +#include #include #include -#include #include @@ -21,6 +21,7 @@ #include "services/LogService.hpp" #include "tracking/ActsGeometryService.h" #include "tracking/ReconstructedHitFactory.h" +#include "tracking/TruthTrackParameterFactory.h" // #include "tracking/Measurement2DFactory.h" struct ProgramArguments { @@ -123,6 +124,13 @@ int main(int argc, char* argv[]) { {"TrackerHit", "Measurement2D"}); app.Add(reco_hit_generator); + auto truth_track_init_generator = new JOmniFactoryGeneratorT(); + truth_track_init_generator->AddWiring( + "TruthTrackParameterGenerator", + {"DigitizedMtpcMcTrack"}, + {"TruthTrackInitParameters"}); + app.Add(truth_track_init_generator); + // auto measurement_2d_generator = new JOmniFactoryGeneratorT(); // measurement_2d_generator->AddWiring("TrackerHitGenerator", {"TrackerHit"}, {"Measurement2D"}); // app.Add(measurement_2d_generator); diff --git a/source/tdis/tracking/CKFTracking.cc b/source/tdis/tracking/CKFTracking.cc new file mode 100644 index 0000000..93d33fe --- /dev/null +++ b/source/tdis/tracking/CKFTracking.cc @@ -0,0 +1,375 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck, Dmitry Romanov, Shujie Li + +#include "CKFTracking.h" + +#include +#include +#include +#include +#if Acts_VERSION_MAJOR < 36 +#include +#endif +#include +#include +#if Acts_VERSION_MAJOR >= 32 +#include "Acts/EventData/ProxyAccessor.hpp" +#endif +#include +#include +#include +#include +#include +#include +#if Acts_VERSION_MAJOR >= 34 +#include "Acts/Propagator/AbortList.hpp" +#include "Acts/Propagator/EigenStepper.hpp" +#include "Acts/Propagator/MaterialInteractor.hpp" +#include "Acts/Propagator/Navigator.hpp" +#endif +#include +#if Acts_VERSION_MAJOR >= 34 +#include "Acts/Propagator/StandardAborters.hpp" +#endif +#include +#include +#include +#include +#include +#if Acts_VERSION_MAJOR >= 34 +#include "Acts/Utilities/TrackHelpers.hpp" +#endif +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "ActsGeometryProvider.h" +#include "DD4hepBField.h" +#include "extensions/spdlog/SpdlogFormatters.h" // IWYU pragma: keep +#include "extensions/spdlog/SpdlogToActs.h" + +namespace eicrecon { + + using namespace Acts::UnitLiterals; + + // This array relates the Acts and EDM4eic covariance matrices, including + // the unit conversion to get from Acts units into EDM4eic units. + // + // Note: std::map is not constexpr, so we use a constexpr std::array + // std::array initialization need double braces since arrays are aggregates + // ref: https://en.cppreference.com/w/cpp/language/aggregate_initialization + static constexpr std::array, 6> edm4eic_indexed_units{{ + {Acts::eBoundLoc0, Acts::UnitConstants::mm}, + {Acts::eBoundLoc1, Acts::UnitConstants::mm}, + {Acts::eBoundPhi, 1.}, + {Acts::eBoundTheta, 1.}, + {Acts::eBoundQOverP, 1. / Acts::UnitConstants::GeV}, + {Acts::eBoundTime, Acts::UnitConstants::ns} + }}; + + CKFTracking::CKFTracking() { + } + + void CKFTracking::init(std::shared_ptr geo_svc, std::shared_ptr log) { + m_log = log; + m_acts_logger = eicrecon::getSpdlogLogger("CKF", m_log); + + m_geoSvc = geo_svc; + + m_BField = std::dynamic_pointer_cast(m_geoSvc->getFieldProvider()); + m_fieldctx = eicrecon::BField::BFieldVariant(m_BField); + + // eta bins, chi2 and #sourclinks per surface cutoffs + m_sourcelinkSelectorCfg = { + {Acts::GeometryIdentifier(), + {m_cfg.etaBins, m_cfg.chi2CutOff, + {m_cfg.numMeasurementsCutOff.begin(), m_cfg.numMeasurementsCutOff.end()} + } + }, + }; + m_trackFinderFunc = CKFTracking::makeCKFTrackingFunction(m_geoSvc->trackingGeometry(), m_BField, logger()); + } + + std::tuple< + std::vector, + std::vector + > + CKFTracking::process(const edm4eic::Measurement2DCollection& meas2Ds, + const edm4eic::TrackParametersCollection &init_trk_params) { + + + // create sourcelink and measurement containers + auto measurements = std::make_shared(); + + // need list here for stable addresses + std::list sourceLinkStorage; + ActsExamples::IndexSourceLinkContainer src_links; + src_links.reserve(meas2Ds.size()); + std::size_t hit_index = 0; + + + for (const auto& meas2D : meas2Ds) { + + // --follow example from ACTS to create source links + sourceLinkStorage.emplace_back(meas2D.getSurface(), hit_index); + ActsExamples::IndexSourceLink& sourceLink = sourceLinkStorage.back(); + // Add to output containers: + // index map and source link container are geometry-ordered. + // since the input is also geometry-ordered, new items can + // be added at the end. + src_links.insert(src_links.end(), sourceLink); + // --- + // Create ACTS measurements + Acts::Vector2 loc = Acts::Vector2::Zero(); + loc[Acts::eBoundLoc0] = meas2D.getLoc().a; + loc[Acts::eBoundLoc1] = meas2D.getLoc().b; + + + Acts::SquareMatrix2 cov = Acts::SquareMatrix2::Zero(); + cov(0, 0) = meas2D.getCovariance().xx; + cov(1, 1) = meas2D.getCovariance().yy; + cov(0, 1) = meas2D.getCovariance().xy; + cov(1, 0) = meas2D.getCovariance().xy; + +#if Acts_VERSION_MAJOR >= 36 + auto measurement = ActsExamples::makeFixedSizeMeasurement( + Acts::SourceLink{sourceLink}, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1); +#else + auto measurement = Acts::makeMeasurement( + Acts::SourceLink{sourceLink}, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1); +#endif + measurements->emplace_back(std::move(measurement)); + + hit_index++; + } + + ActsExamples::TrackParametersContainer acts_init_trk_params; + for (const auto& track_parameter: init_trk_params) { + + Acts::BoundVector params; + params(Acts::eBoundLoc0) = track_parameter.getLoc().a * Acts::UnitConstants::mm; // cylinder radius + params(Acts::eBoundLoc1) = track_parameter.getLoc().b * Acts::UnitConstants::mm; // cylinder length + params(Acts::eBoundPhi) = track_parameter.getPhi(); + params(Acts::eBoundTheta) = track_parameter.getTheta(); + params(Acts::eBoundQOverP) = track_parameter.getQOverP() / Acts::UnitConstants::GeV; + params(Acts::eBoundTime) = track_parameter.getTime() * Acts::UnitConstants::ns; + + double charge = std::copysign(1., track_parameter.getQOverP()); + + Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Zero(); + for (size_t i = 0; const auto& [a, x] : edm4eic_indexed_units) { + for (size_t j = 0; const auto& [b, y] : edm4eic_indexed_units) { + cov(a, b) = track_parameter.getCovariance()(i,j) * x * y; + ++j; + } + ++i; + } + + // Construct a perigee surface as the target surface + auto pSurface = Acts::Surface::makeShared(Acts::Vector3(0,0,0)); + + // Create parameters + acts_init_trk_params.emplace_back(pSurface, params, cov, Acts::ParticleHypothesis::pion()); + } + + //// Construct a perigee surface as the target surface + auto pSurface = Acts::Surface::makeShared(Acts::Vector3{0., 0., 0.}); + + ACTS_LOCAL_LOGGER(eicrecon::getSpdlogLogger("CKF", m_log, {"^No tracks found$"})); + + Acts::PropagatorPlainOptions pOptions; + pOptions.maxSteps = 10000; + + ActsExamples::PassThroughCalibrator pcalibrator; + ActsExamples::MeasurementCalibratorAdapter calibrator(pcalibrator, *measurements); + Acts::GainMatrixUpdater kfUpdater; +#if Acts_VERSION_MAJOR < 34 + Acts::GainMatrixSmoother kfSmoother; +#endif + Acts::MeasurementSelector measSel{m_sourcelinkSelectorCfg}; + + Acts::CombinatorialKalmanFilterExtensions + extensions; + extensions.calibrator.connect<&ActsExamples::MeasurementCalibratorAdapter::calibrate>( + &calibrator); + extensions.updater.connect< + &Acts::GainMatrixUpdater::operator()>( + &kfUpdater); +#if Acts_VERSION_MAJOR < 34 + extensions.smoother.connect< + &Acts::GainMatrixSmoother::operator()>( + &kfSmoother); +#endif + extensions.measurementSelector.connect< + &Acts::MeasurementSelector::select>( + &measSel); + + ActsExamples::IndexSourceLinkAccessor slAccessor; + slAccessor.container = &src_links; + Acts::SourceLinkAccessorDelegate + slAccessorDelegate; + slAccessorDelegate.connect<&ActsExamples::IndexSourceLinkAccessor::range>(&slAccessor); + + // Set the CombinatorialKalmanFilter options +#if Acts_VERSION_MAJOR < 34 + CKFTracking::TrackFinderOptions options( + m_geoctx, m_fieldctx, m_calibctx, slAccessorDelegate, + extensions, pOptions, &(*pSurface)); +#else + CKFTracking::TrackFinderOptions options( + m_geoctx, m_fieldctx, m_calibctx, slAccessorDelegate, + extensions, pOptions); +#endif + +#if Acts_VERSION_MAJOR >= 34 + Acts::Propagator, Acts::Navigator> extrapolator( + Acts::EigenStepper<>(m_BField), + Acts::Navigator({m_geoSvc->trackingGeometry()}, + logger().cloneWithSuffix("Navigator")), + logger().cloneWithSuffix("Propagator")); + + Acts::PropagatorOptions, + Acts::AbortList> + extrapolationOptions(m_geoctx, m_fieldctx); +#endif + + // Create track container + auto trackContainer = std::make_shared(); + auto trackStateContainer = std::make_shared(); + ActsExamples::TrackContainer acts_tracks(trackContainer, trackStateContainer); + + // Add seed number column + acts_tracks.addColumn("seed"); +#if Acts_VERSION_MAJOR >= 32 + Acts::ProxyAccessor seedNumber("seed"); +#else + Acts::TrackAccessor seedNumber("seed"); +#endif + + // Loop over seeds + for (std::size_t iseed = 0; iseed < acts_init_trk_params.size(); ++iseed) { + auto result = + (*m_trackFinderFunc)(acts_init_trk_params.at(iseed), options, acts_tracks); + + if (!result.ok()) { + m_log->debug("Track finding failed for seed {} with error {}", iseed, result.error()); + continue; + } + + // Set seed number for all found tracks + auto& tracksForSeed = result.value(); + for (auto& track : tracksForSeed) { + +#if Acts_VERSION_MAJOR >=34 + auto smoothingResult = Acts::smoothTrack(m_geoctx, track, logger()); + if (!smoothingResult.ok()) { + ACTS_ERROR("Smoothing for seed " + << iseed << " and track " << track.index() + << " failed with error " << smoothingResult.error()); + continue; + } + + auto extrapolationResult = Acts::extrapolateTrackToReferenceSurface( + track, *pSurface, extrapolator, extrapolationOptions, + Acts::TrackExtrapolationStrategy::firstOrLast, logger()); + if (!extrapolationResult.ok()) { + ACTS_ERROR("Extrapolation for seed " + << iseed << " and track " << track.index() + << " failed with error " << extrapolationResult.error()); + continue; + } +#endif + + seedNumber(track) = iseed; + } + } + + + // Move track states and track container to const containers + // NOTE Using the non-const containers leads to references to + // implicitly converted temporaries inside the Trajectories. + auto constTrackStateContainer = + std::make_shared( + std::move(*trackStateContainer)); + + auto constTrackContainer = + std::make_shared( + std::move(*trackContainer)); + + // FIXME JANA2 std::vector requires wrapping ConstTrackContainer, instead of: + //ConstTrackContainer constTracks(constTrackContainer, constTrackStateContainer); + std::vector constTracks_v; + constTracks_v.push_back( + new ActsExamples::ConstTrackContainer( + constTrackContainer, + constTrackStateContainer)); + auto& constTracks = *(constTracks_v.front()); + + // Seed number column accessor +#if Acts_VERSION_MAJOR >= 32 + const Acts::ConstProxyAccessor constSeedNumber("seed"); +#else + const Acts::ConstTrackAccessor constSeedNumber("seed"); +#endif + + // Prepare the output data with MultiTrajectory, per seed + std::vector acts_trajectories; + acts_trajectories.reserve(init_trk_params.size()); + + ActsExamples::Trajectories::IndexedParameters parameters; + std::vector tips; + + std::optional lastSeed; + for (const auto& track : constTracks) { + if (!lastSeed) { + lastSeed = constSeedNumber(track); + } + + if (constSeedNumber(track) != lastSeed.value()) { + // make copies and clear vectors + acts_trajectories.push_back(new ActsExamples::Trajectories( + constTracks.trackStateContainer(), + tips, parameters)); + + tips.clear(); + parameters.clear(); + } + + lastSeed = constSeedNumber(track); + + tips.push_back(track.tipIndex()); + parameters.emplace( + std::pair{track.tipIndex(), + ActsExamples::TrackParameters{track.referenceSurface().getSharedPtr(), + track.parameters(), track.covariance(), + track.particleHypothesis()}}); + } + + if (tips.empty()) { + m_log->info("Last trajectory is empty"); + } + + // last entry: move vectors + acts_trajectories.push_back(new ActsExamples::Trajectories( + constTracks.trackStateContainer(), + std::move(tips), std::move(parameters))); + + return std::make_tuple(std::move(acts_trajectories), std::move(constTracks_v)); + } + +} // namespace eicrecon diff --git a/source/tdis/tracking/CKFTracking.h b/source/tdis/tracking/CKFTracking.h new file mode 100644 index 0000000..df67aec --- /dev/null +++ b/source/tdis/tracking/CKFTracking.h @@ -0,0 +1,91 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "CKFTrackingConfig.h" +#include "DD4hepBField.h" +#include "algorithms/interfaces/WithPodConfig.h" + +class ActsGeometryProvider; + +namespace eicrecon { + +/** Fitting algorithm implementation . + * + * \ingroup tracking + */ + + class CKFTracking: public WithPodConfig { + public: + + + /// Find function that takes the above parameters + /// @note This is separated into a virtual interface to keep compilation units + /// small + class CKFTrackingFunction { + public: + virtual ~CKFTrackingFunction() = default; + + virtual TrackFinderResult operator()(const ActsExamples::TrackParameters&, + const TrackFinderOptions&, + ActsExamples::TrackContainer&) const = 0; + }; + + /// Create the track finder function implementation. + /// The magnetic field is intentionally given by-value since the variantresults + /// contains shared_ptr anyways. + static std::shared_ptr makeCKFTrackingFunction( + std::shared_ptr trackingGeometry, + std::shared_ptr magneticField, + const Acts::Logger& logger); + + CKFTracking(); + + void init(std::shared_ptr geo_svc, std::shared_ptr log); + + std::tuple< + std::vector, + std::vector + > + process(const edm4eic::Measurement2DCollection& meas2Ds, + const edm4eic::TrackParametersCollection &init_trk_params); + + private: + std::shared_ptr m_log; + std::shared_ptr m_acts_logger{nullptr}; + std::shared_ptr m_trackFinderFunc; + std::shared_ptr m_geoSvc; + + std::shared_ptr m_BField = nullptr; + Acts::GeometryContext m_geoctx; + Acts::CalibrationContext m_calibctx; + Acts::MagneticFieldContext m_fieldctx; + + Acts::MeasurementSelector::Config m_sourcelinkSelectorCfg; + + /// Private access to the logging instance + const Acts::Logger& logger() const { return *m_acts_logger; } + }; + +} // namespace eicrecon::Reco diff --git a/source/tdis/tracking/TruthTrackParameterFactory.h b/source/tdis/tracking/TruthTrackParameterFactory.h index 7cfa22f..1840251 100644 --- a/source/tdis/tracking/TruthTrackParameterFactory.h +++ b/source/tdis/tracking/TruthTrackParameterFactory.h @@ -2,96 +2,151 @@ #include #include +#include +#include +#include +#include "ActsGeometryService.h" #include "PadGeometryHelper.hpp" #include "podio_model/DigitizedMtpcMcHit.h" +#include "podio_model/DigitizedMtpcMcTrack.h" #include "podio_model/MutableTrackerHit.h" +#include "podio_model/TrackParametersCollection.h" #include "podio_model/TrackerHit.h" #include "podio_model/TrackerHitCollection.h" -namespace { - inline double get_resolution(const double pixel_size) { - constexpr const double sqrt_12 = 3.4641016151; - return pixel_size / sqrt_12; - } - inline double get_variance(const double pixel_size) { - const double res = get_resolution(pixel_size); - return res * res; - } -} // namespace namespace tdis::tracking { struct TruthTrackParameterFactory : public JOmniFactory { - PodioInput m_mc_hits_in {this, {"DigitizedMtpcMcHit"}}; - PodioOutput m_tracker_hits_out {this, "TrackerHit"}; + PodioInput m_mc_tracks_in {this, {"DigitizedMtpcMcTrack"}}; + PodioOutput m_trackers_out {this, "TruthTrackInitParameters"}; Service m_service_geometry{this}; + Service m_service_log{this}; + Parameter m_cfg_use_true_pos{this, "acts:use_true_position", true,"Use true hits xyz instead of digitized one"}; + Parameter m_cfg_momentum_smear{this, "acts:track_init:momentum_smear", 0.1,"GeV, Momentum smear for truth track initialization"}; + std::shared_ptr m_log; + void Configure() { m_service_geometry(); + m_log = m_service_log->logger("tracking:hit_reco"); } void ChangeRun(int32_t /*run_nr*/) { } + // Function to generate the next value in a normal distribution + double generateNormal(double mean, double stddev) { + // Create a random device and a generator + static std::random_device rd; + static std::mt19937 generator(rd()); + + // Create a normal distribution with the given mean and standard deviation + std::normal_distribution distribution(mean, stddev); + + // Generate and return the next value + return distribution(generator); + } + void Execute(int32_t /*run_nr*/, uint64_t /*evt_nr*/) { using namespace Acts::UnitLiterals; - auto rec_hits = std::make_unique(); + auto track_parameters = std::make_unique(); auto plane_positions = m_service_geometry->GetPlanePositions(); - for (auto mc_hit : *m_mc_hits_in()) { - const int plane = mc_hit.plane(); - const int ring = mc_hit.ring(); - const int pad = mc_hit.pad(); - const double z_to_gem = mc_hit.zToGem(); + // Loop over input particles + for (const auto& mc_track: *m_mc_tracks_in()) { - auto [pad_x,pad_y] = getPadCenter(ring, pad); - double plane_z = m_service_geometry().GetPlanePositions()[plane]; + // We need at least one hit + if(!mc_track.hits_size()) { + continue; + } + // require close to interaction vertex + auto v = mc_track.hits().at(0).truePosition(); // Use it as a vertex for now + double vx = v.x; + double vy = v.y; + double vz = v.z; - // Planes are set back to back like this: - // | || || || | - // This means that zToPlane is positive for even planes and negative for odd - // i.e. - // z = plane_z + zToPlane // for even - // z = plane_z - zToPlane // for odd - double calc_z = plane_z + (plane % 2 ? -z_to_gem : z_to_gem); + double magnitude = mc_track.momentum(); + double theta = mc_track.theta(); + double phi = mc_track.phi(); + const auto eta = -std::log(std::tan(theta/2)); + double px = magnitude * std::sin(theta) * std::cos(phi); + double py = magnitude * std::sin(theta) * std::sin(phi); + double pz = magnitude * std::cos(theta); - // Position - edm4hep::Vector3f position; - if(m_cfg_use_true_pos() && !std::isnan(mc_hit.truePosition().x)) { - position.x = mc_hit.truePosition().x; - position.y = mc_hit.truePosition().y; - position.z = mc_hit.truePosition().z; - } - else { - position.x = pad_x; - position.y = pad_y; - position.z = calc_z; + // require minimum momentum + const auto& p = mc_track.momentum(); + const auto pmag = std::hypot(px, py, pz); + + + + // modify initial momentum to avoid bleeding truth to results when fit fails + const auto pinit = pmag * generateNormal(1, m_cfg_momentum_smear()*Acts::UnitConstants::GeV); + + // define line surface for local position values + auto perigee = Acts::Surface::makeShared(Acts::Vector3(0,0,0)); + + // track particle back to transverse point-of-closest approach + // with respect to the defined line surface + auto linesurface_parameter = -(vx*px + vy*py)/(px*px + py*py); + + auto xpca = v.x + linesurface_parameter*px; + auto ypca = v.y + linesurface_parameter*py; + auto zpca = v.z + linesurface_parameter*pz; + + Acts::Vector3 global(xpca, ypca, zpca); + Acts::Vector3 direction(sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)); + + // convert from global to local coordinates using the defined line surface + auto local = perigee->globalToLocal(m_service_geometry->GetActsGeometryContext(), global, direction); + + if(!local.ok()) + { + m_log->error("skipping the track because globaltoLocal function failed"); + continue; } - // Covariance - double max_dimension = std::max(getPadApproxWidth(ring), getPadHight()); - double xy_variance = get_variance(max_dimension); - edm4eic::CovDiag3f cov{xy_variance, xy_variance, 1_cm}; // TODO correct Z variance - - uint32_t cell_id = 1000000*mc_hit.plane() + 1000*mc_hit.ring() + mc_hit.pad(); - - rec_hits->create( - cell_id, - position, - cov, - static_cast(mc_hit.time()), - static_cast(1_ns), // TODO correct time resolution - static_cast(mc_hit.adc()), - 0.0F); + Acts::Vector2 localpos = local.value(); + double charge = 1; // TODO ??? Proton? + + // Insert into edm4eic::TrackParameters, which uses numerical values in its specified units + auto track_parameter = track_parameters->create(); + track_parameter.type(-1); // type --> seed(-1) + track_parameter.loc({static_cast(localpos(0)), static_cast(localpos(1))}); // 2d location on surface [mm] + track_parameter.phi(phi); // phi [rad] + track_parameter.theta(theta); // theta [rad] + track_parameter.qOverP(charge / (pinit)); // Q/p [e/GeV] + track_parameter.time(mc_track.hits().at(0).time()); // time [ns] + edm4eic::Cov6f cov; + cov(0,0) = 1.0; // loc0 + cov(1,1) = 1.0; // loc1 + cov(2,2) = 0.05; // phi + cov(3,3) = 0.01; // theta + cov(4,4) = 0.1; // qOverP + cov(5,5) = 10e9; // time + track_parameter.covariance(cov); + + // // Debug output + // if (m_log->level() <= spdlog::level::debug) { + // m_log->debug("Invoke track finding seeded by truth particle with:"); + // m_log->debug(" p = {} GeV (smeared to {} GeV)", pmag / dd4hep::GeV, pinit / dd4hep::GeV); + // m_log->debug(" q = {}", charge); + // m_log->debug(" q/p = {} e/GeV (smeared to {} e/GeV)", charge / (pmag / dd4hep::GeV), charge / (pinit / dd4hep::GeV)); + // m_log->debug(" theta = {}", theta); + // m_log->debug(" phi = {}", phi); + // } } - m_tracker_hits_out() = std::move(rec_hits); + + + + m_trackers_out() = std::move(track_parameters); } }; } \ No newline at end of file