diff --git a/src/algorithms/tracking/ActsToTracks.cc b/src/algorithms/tracking/ActsToTracks.cc new file mode 100644 index 0000000000..0bfeee2bd4 --- /dev/null +++ b/src/algorithms/tracking/ActsToTracks.cc @@ -0,0 +1,195 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Whitney Armstrong, Wouter Deconinck, Dmitry Romanov, Shujie Li, Dmitry Kalinkin + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "ActsToTracks.h" + +namespace eicrecon { + +// 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} +}}; + +void ActsToTracks::init() { +} + +void ActsToTracks::process(const Input& input, const Output& output) const { + const auto [meas2Ds, acts_trajectories] = input; + auto [trajectories, track_parameters, tracks] = output; + + // Loop over trajectories + for (const auto traj : acts_trajectories) { + // The trajectory entry indices and the multiTrajectory + const auto& trackTips = traj->tips(); + const auto& mj = traj->multiTrajectory(); + if (trackTips.empty()) { + warning("Empty multiTrajectory."); + continue; + } + + // Loop over all trajectories in a multiTrajectory + // FIXME: we only retain the first trackTips entry + for (auto trackTip : decltype(trackTips){trackTips.front()}) { + // Collect the trajectory summary info + auto trajectoryState = + Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip); + + // Check if the reco track has fitted track parameters + if (not traj->hasTrackParameters(trackTip)) { + warning( + "No fitted track parameters for trajectory with entry index = {}", + trackTip); + continue; + } + + // Create trajectory + auto trajectory = trajectories->create(); + trajectory.setNMeasurements(trajectoryState.nMeasurements); + trajectory.setNStates(trajectoryState.nStates); + trajectory.setNOutliers(trajectoryState.nOutliers); + trajectory.setNHoles(trajectoryState.nHoles); + trajectory.setNSharedHits(trajectoryState.nSharedHits); + + debug("trajectory state, measurement, outlier, hole: {} {} {} {}", + trajectoryState.nStates, + trajectoryState.nMeasurements, + trajectoryState.nOutliers, + trajectoryState.nHoles); + + for (const auto& measurementChi2 : trajectoryState.measurementChi2) { + trajectory.addToMeasurementChi2(measurementChi2); + } + + for (const auto& outlierChi2 : trajectoryState.outlierChi2) { + trajectory.addToOutlierChi2(outlierChi2); + } + + // Get the fitted track parameter + const auto& boundParam = traj->trackParameters(trackTip); + const auto& parameter = boundParam.parameters(); + const auto& covariance = *boundParam.covariance(); + + auto pars = track_parameters->create(); + pars.setType(0); // type: track head --> 0 + pars.setLoc({ + static_cast(parameter[Acts::eBoundLoc0]), + static_cast(parameter[Acts::eBoundLoc1]) + }); + pars.setTheta(static_cast(parameter[Acts::eBoundTheta])); + pars.setPhi(static_cast(parameter[Acts::eBoundPhi])); + pars.setQOverP(static_cast(parameter[Acts::eBoundQOverP])); + pars.setTime(static_cast(parameter[Acts::eBoundTime])); + edm4eic::Cov6f cov; + for (size_t i = 0; const auto& [a, x] : edm4eic_indexed_units) { + for (size_t j = 0; const auto& [b, y] : edm4eic_indexed_units) { + // FIXME why not pars.getCovariance()(i,j) = covariance(a,b) / x / y; + cov(i,j) = covariance(a,b) / x / y; + } + } + pars.setCovariance(cov); + + trajectory.addToTrackParameters(pars); + + // Fill tracks + auto track = tracks->create(); + track.setType( // Flag that defines the type of track + pars.getType() + ); + track.setPosition( // Track 3-position at the vertex + edm4hep::Vector3f() + ); + track.setMomentum( // Track 3-momentum at the vertex [GeV] + edm4hep::Vector3f() + ); + track.setPositionMomentumCovariance( // Covariance matrix in basis [x,y,z,px,py,pz] + edm4eic::Cov6f() + ); + track.setTime( // Track time at the vertex [ns] + static_cast(parameter[Acts::eBoundTime]) + ); + track.setTimeError( // Error on the track vertex time + sqrt(static_cast(covariance(Acts::eBoundTime, Acts::eBoundTime))) + ); + track.setCharge( // Particle charge + std::copysign(1., parameter[Acts::eBoundQOverP]) + ); + track.setChi2(trajectoryState.chi2Sum); // Total chi2 + track.setNdf(trajectoryState.NDF); // Number of degrees of freedom + track.setPdg( // PDG particle ID hypothesis + boundParam.particleHypothesis().absolutePdg() + ); + track.setTrajectory(trajectory); // Trajectory of this track + + // save measurement2d to good measurements or outliers according to srclink index + // fix me: ideally, this should be integrated into multitrajectoryhelper + // fix me: should say "OutlierMeasurements" instead of "OutlierHits" etc + mj.visitBackwards(trackTip, [&](const auto& state) { + + auto geoID = state.referenceSurface().geometryId().value(); + auto typeFlags = state.typeFlags(); + + // find the associated hit (2D measurement) with state sourcelink index + // fix me: calibrated or not? + if (state.hasUncalibratedSourceLink()) { + + std::size_t srclink_index = state.getUncalibratedSourceLink().template get().index(); + + // no hit on this state/surface, skip + if (typeFlags.test(Acts::TrackStateFlag::HoleFlag)) { + debug("No hit found on geo id={}", geoID); + + } else { + auto meas2D = (*meas2Ds) [srclink_index]; + if (typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) { + track.addToMeasurements(meas2D); + trajectory.addToMeasurements_deprecated(meas2D); + debug("Measurement on geo id={}, index={}, loc={},{}", + geoID, srclink_index, meas2D.getLoc().a, meas2D.getLoc().b); + + } + else if (typeFlags.test(Acts::TrackStateFlag::OutlierFlag)) { + trajectory.addToOutliers_deprecated(meas2D); + debug("Outlier on geo id={}, index={}, loc={},{}", + geoID, srclink_index, meas2D.getLoc().a, meas2D.getLoc().b); + + } + } + } + + }); + + } + } +} + +} diff --git a/src/algorithms/tracking/ActsToTracks.h b/src/algorithms/tracking/ActsToTracks.h new file mode 100644 index 0000000000..95bf88ec9f --- /dev/null +++ b/src/algorithms/tracking/ActsToTracks.h @@ -0,0 +1,32 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Dmitry Kalinkin + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace eicrecon { + +using ActsToTracksAlgorithm = + algorithms::Algorithm< + algorithms::Input>, + algorithms::Output + >; + +class ActsToTracks : public ActsToTracksAlgorithm { +public: + ActsToTracks(std::string_view name) : ActsToTracksAlgorithm{name, {"inputActsTrajectories", "inputActsConstTrackContainer"}, {"outputTrajectoryCollection", "outputTrackParametersCollection", "outputTrackCollection"}, "Converts ACTS trajectories to EDM4eic"} {}; + + void init() final; + void process(const Input&, const Output&) const final; +}; + +} diff --git a/src/algorithms/tracking/CKFTracking.cc b/src/algorithms/tracking/CKFTracking.cc index 10d8402900..de29d15f98 100644 --- a/src/algorithms/tracking/CKFTracking.cc +++ b/src/algorithms/tracking/CKFTracking.cc @@ -9,12 +9,10 @@ #include #include #include -#include #include #include #include #include -#include #include #include #include @@ -32,12 +30,9 @@ #include #include #include -#include #include -#include #include #include -#include #include #include #include @@ -94,9 +89,6 @@ namespace eicrecon { } std::tuple< - std::unique_ptr, - std::unique_ptr, - std::unique_ptr, std::vector, std::vector > @@ -172,10 +164,6 @@ namespace eicrecon { acts_init_trk_params.emplace_back(pSurface, params, cov, Acts::ParticleHypothesis::pion()); } - auto trajectories = std::make_unique(); - auto track_parameters = std::make_unique(); - auto tracks = std::make_unique(); - //// Construct a perigee surface as the target surface auto pSurface = Acts::Surface::makeShared(Acts::Vector3{0., 0., 0.}); @@ -308,152 +296,7 @@ namespace eicrecon { constTracks.trackStateContainer(), std::move(tips), std::move(parameters))); - - // Loop over trajectories - for (const auto* traj : acts_trajectories) { - // The trajectory entry indices and the multiTrajectory - const auto& trackTips = traj->tips(); - const auto& mj = traj->multiTrajectory(); - if (trackTips.empty()) { - m_log->warn("Empty multiTrajectory."); - continue; - } - - // Loop over all trajectories in a multiTrajectory - // FIXME: we only retain the first trackTips entry - for (auto trackTip : decltype(trackTips){trackTips.front()}) { - // Collect the trajectory summary info - auto trajectoryState = - Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip); - - // Check if the reco track has fitted track parameters - if (not traj->hasTrackParameters(trackTip)) { - m_log->warn( - "No fitted track parameters for trajectory with entry index = {}", - trackTip); - continue; - } - - // Create trajectory - auto trajectory = trajectories->create(); - trajectory.setNMeasurements(trajectoryState.nMeasurements); - trajectory.setNStates(trajectoryState.nStates); - trajectory.setNOutliers(trajectoryState.nOutliers); - trajectory.setNHoles(trajectoryState.nHoles); - trajectory.setNSharedHits(trajectoryState.nSharedHits); - - m_log->debug("trajectory state, measurement, outlier, hole: {} {} {} {}", - trajectoryState.nStates, - trajectoryState.nMeasurements, - trajectoryState.nOutliers, - trajectoryState.nHoles); - - for (const auto& measurementChi2 : trajectoryState.measurementChi2) { - trajectory.addToMeasurementChi2(measurementChi2); - } - - for (const auto& outlierChi2 : trajectoryState.outlierChi2) { - trajectory.addToOutlierChi2(outlierChi2); - } - - // Get the fitted track parameter - const auto& boundParam = traj->trackParameters(trackTip); - const auto& parameter = boundParam.parameters(); - const auto& covariance = *boundParam.covariance(); - - auto pars = track_parameters->create(); - pars.setType(0); // type: track head --> 0 - pars.setLoc({ - static_cast(parameter[Acts::eBoundLoc0]), - static_cast(parameter[Acts::eBoundLoc1]) - }); - pars.setTheta(static_cast(parameter[Acts::eBoundTheta])); - pars.setPhi(static_cast(parameter[Acts::eBoundPhi])); - pars.setQOverP(static_cast(parameter[Acts::eBoundQOverP])); - pars.setTime(static_cast(parameter[Acts::eBoundTime])); - edm4eic::Cov6f cov; - for (size_t i = 0; const auto& [a, x] : edm4eic_indexed_units) { - for (size_t j = 0; const auto& [b, y] : edm4eic_indexed_units) { - // FIXME why not pars.getCovariance()(i,j) = covariance(a,b) / x / y; - cov(i,j) = covariance(a,b) / x / y; - } - } - pars.setCovariance(cov); - - trajectory.addToTrackParameters(pars); - - // Fill tracks - auto track = tracks->create(); - track.setType( // Flag that defines the type of track - pars.getType() - ); - track.setPosition( // Track 3-position at the vertex - edm4hep::Vector3f() - ); - track.setMomentum( // Track 3-momentum at the vertex [GeV] - edm4hep::Vector3f() - ); - track.setPositionMomentumCovariance( // Covariance matrix in basis [x,y,z,px,py,pz] - edm4eic::Cov6f() - ); - track.setTime( // Track time at the vertex [ns] - static_cast(parameter[Acts::eBoundTime]) - ); - track.setTimeError( // Error on the track vertex time - sqrt(static_cast(covariance(Acts::eBoundTime, Acts::eBoundTime))) - ); - track.setCharge( // Particle charge - std::copysign(1., parameter[Acts::eBoundQOverP]) - ); - track.setChi2(trajectoryState.chi2Sum); // Total chi2 - track.setNdf(trajectoryState.NDF); // Number of degrees of freedom - track.setPdg( // PDG particle ID hypothesis - boundParam.particleHypothesis().absolutePdg() - ); - track.setTrajectory(trajectory); // Trajectory of this track - - // save measurement2d to good measurements or outliers according to srclink index - // fix me: ideally, this should be integrated into multitrajectoryhelper - // fix me: should say "OutlierMeasurements" instead of "OutlierHits" etc - mj.visitBackwards(trackTip, [&](const auto& state) { - - auto geoID = state.referenceSurface().geometryId().value(); - auto typeFlags = state.typeFlags(); - - // find the associated hit (2D measurement) with state sourcelink index - // fix me: calibrated or not? - if (state.hasUncalibratedSourceLink()) { - - std::size_t srclink_index = state.getUncalibratedSourceLink().template get().index(); - - // no hit on this state/surface, skip - if (typeFlags.test(Acts::TrackStateFlag::HoleFlag)) { - m_log->debug("No hit found on geo id={}", geoID); - - } else { - auto meas2D = meas2Ds[srclink_index]; - if (typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) { - track.addToMeasurements(meas2D); - trajectory.addToMeasurements_deprecated(meas2D); - m_log->debug("Measurement on geo id={}, index={}, loc={},{}", - geoID, srclink_index, meas2D.getLoc().a, meas2D.getLoc().b); - - } - else if (typeFlags.test(Acts::TrackStateFlag::OutlierFlag)) { - trajectory.addToOutliers_deprecated(meas2D); - m_log->debug("Outlier on geo id={}, index={}, loc={},{}", - geoID, srclink_index, meas2D.getLoc().a, meas2D.getLoc().b); - - } - } - } - - }); - - } - } - - return std::make_tuple(std::move(trajectories), std::move(track_parameters), std::move(tracks), std::move(acts_trajectories), std::move(constTracks_v)); + return std::make_tuple(std::move(acts_trajectories), std::move(constTracks_v)); } } // namespace eicrecon diff --git a/src/algorithms/tracking/CKFTracking.h b/src/algorithms/tracking/CKFTracking.h index 3b9b847eb9..6373488b79 100644 --- a/src/algorithms/tracking/CKFTracking.h +++ b/src/algorithms/tracking/CKFTracking.h @@ -17,9 +17,7 @@ #include #include #include -#include #include -#include #include #include #include @@ -73,9 +71,6 @@ namespace eicrecon { void init(std::shared_ptr geo_svc, std::shared_ptr log); std::tuple< - std::unique_ptr, - std::unique_ptr, - std::unique_ptr, std::vector, std::vector > diff --git a/src/global/tracking/ActsToTracks_factory.h b/src/global/tracking/ActsToTracks_factory.h new file mode 100644 index 0000000000..8fdc30a3ce --- /dev/null +++ b/src/global/tracking/ActsToTracks_factory.h @@ -0,0 +1,46 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright 2024, Dmitry Kalinkin + +#pragma once + +#include + +#include "algorithms/tracking/ActsToTracks.h" +#include "extensions/jana/JOmniFactory.h" + +namespace eicrecon { + +class ActsToTracks_factory + : public JOmniFactory { +public: + using AlgoT = eicrecon::ActsToTracks; + +private: + std::unique_ptr m_algo; + + PodioInput m_measurements_input {this}; + Input m_acts_trajectories_input {this}; + PodioOutput m_trajectories_output {this}; + PodioOutput m_parameters_output {this}; + PodioOutput m_tracks_output {this}; + +public: + void Configure() { + m_algo = std::make_unique(this->GetPrefix()); + m_algo->level((algorithms::LogLevel)logger()->level()); + m_algo->init(); + }; + + void ChangeRun(int64_t run_number) {}; + + void Process(int64_t run_number, uint64_t event_number) { + std::vector> acts_trajectories_input; + for (auto acts_traj : m_acts_trajectories_input()) { + acts_trajectories_input.push_back(acts_traj); + } + m_algo->process({m_measurements_input(), acts_trajectories_input}, + {m_trajectories_output().get(), m_parameters_output().get(), m_tracks_output().get()}); + } +}; + +} // namespace eicrecon diff --git a/src/global/tracking/CKFTracking_factory.h b/src/global/tracking/CKFTracking_factory.h index 803e4e63e4..91f2e7c2fd 100644 --- a/src/global/tracking/CKFTracking_factory.h +++ b/src/global/tracking/CKFTracking_factory.h @@ -30,9 +30,6 @@ class CKFTracking_factory : PodioInput m_parameters_input {this}; PodioInput m_measurements_input {this}; - PodioOutput m_trajectories_output {this}; - PodioOutput m_parameters_output {this}; - PodioOutput m_tracks_output {this}; Output m_acts_trajectories_output {this}; Output m_acts_tracks_output {this}; @@ -54,9 +51,6 @@ class CKFTracking_factory : void Process(int64_t run_number, uint64_t event_number) { std::tie( - m_trajectories_output(), - m_parameters_output(), - m_tracks_output(), m_acts_trajectories_output(), m_acts_tracks_output() ) = m_algo->process( diff --git a/src/global/tracking/tracking.cc b/src/global/tracking/tracking.cc index 7f3070afd4..fa37a317d8 100644 --- a/src/global/tracking/tracking.cc +++ b/src/global/tracking/tracking.cc @@ -12,17 +12,18 @@ #include #include -#include "CKFTrackingConfig.h" +#include "ActsToTracks.h" +#include "ActsToTracks_factory.h" #include "CKFTracking_factory.h" #include "IterativeVertexFinder_factory.h" -#include "TracksToParticlesConfig.h" -#include "TracksToParticles_factory.h" #include "TrackParamTruthInit_factory.h" #include "TrackProjector_factory.h" #include "TrackPropagationConfig.h" #include "TrackPropagation_factory.h" #include "TrackSeeding_factory.h" #include "TrackerMeasurementFromHits_factory.h" +#include "TracksToParticlesConfig.h" +#include "TracksToParticles_factory.h" #include "extensions/jana/JOmniFactoryGeneratorT.h" #include "factories/meta/CollectionCollector_factory.h" #include "services/geometry/dd4hep/DD4hep_service.h" @@ -87,12 +88,23 @@ void InitPlugin(JApplication *app) { "InitTrackParams", "CentralTrackerMeasurements" }, + { + "CentralCKFActsTrajectories", + "CentralCKFActsTracks", + }, + app + )); + + app->Add(new JOmniFactoryGeneratorT( + "CentralCKFTracks", + { + "CentralTrackerMeasurements", + "CentralCKFActsTrajectories", + }, { "CentralCKFTrajectories", "CentralCKFTrackParameters", "CentralCKFTracks", - "CentralCKFActsTrajectories", - "CentralCKFActsTracks", }, app )); @@ -111,12 +123,23 @@ void InitPlugin(JApplication *app) { "CentralTrackSeedingResults", "CentralTrackerMeasurements" }, + { + "CentralCKFSeededActsTrajectories", + "CentralCKFSeededActsTracks", + }, + app + )); + + app->Add(new JOmniFactoryGeneratorT( + "CentralCKFSeededTracks", + { + "CentralTrackerMeasurements", + "CentralCKFSeededActsTrajectories", + }, { "CentralCKFSeededTrajectories", "CentralCKFSeededTrackParameters", "CentralCKFSeededTracks", - "CentralCKFSeededActsTrajectories", - "CentralCKFSeededActsTracks", }, app ));