algorithms/tracking/TrackProjector.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2025 wfan, Whitney Armstrong, Sylvester Joosten, Dmitry Kalinkin
0003
0004 #include <Acts/Definitions/TrackParametrization.hpp>
0005 #include <Acts/EventData/MultiTrajectoryHelpers.hpp>
0006 #if Acts_VERSION_MAJOR >= 34
0007 #include <Acts/EventData/TransformationHelpers.hpp>
0008 #endif
0009 #include <Acts/Geometry/GeometryIdentifier.hpp>
0010 #include <Acts/Utilities/UnitVectors.hpp>
0011 #include <ActsExamples/EventData/Trajectories.hpp>
0012 #include <algorithms/service.h>
0013 #include <edm4eic/Cov2f.h>
0014 #include <edm4eic/Cov3f.h>
0015 #include <edm4eic/TrackParametersCollection.h>
0016 #include <edm4eic/TrackPoint.h>
0017 #include <edm4eic/TrackSegmentCollection.h>
0018 #include <edm4hep/Vector3f.h>
0019 #include <edm4hep/utils/vector_utils.h>
0020 #include <fmt/core.h>
0021 #include <fmt/ostream.h>
0022 #include <cmath>
0023 #include <cstddef>
0024 #include <cstdint>
0025 #include <gsl/pointers>
0026 #include <iterator>
0027
0028 #include "TrackProjector.h"
0029 #include "algorithms/interfaces/ActsSvc.h"
0030 #include "extensions/spdlog/SpdlogFormatters.h" // IWYU pragma: keep
0031
0032 #if FMT_VERSION >= 90000
0033 template <> struct fmt::formatter<Acts::GeometryIdentifier> : fmt::ostream_formatter {};
0034 #endif // FMT_VERSION >= 90000
0035
0036 namespace eicrecon {
0037
0038 void TrackProjector::init() {
0039   auto& serviceSvc = algorithms::ServiceSvc::instance();
0040   m_geo_provider   = serviceSvc.service<algorithms::ActsSvc>("ActsSvc")->acts_geometry_provider();
0041 }
0042
0043 void TrackProjector::process(const Input& input, const Output& output) const {
0044   const auto [acts_trajectories, tracks] = input;
0045   auto [track_segments]                  = output;
0046
0047   debug("Track projector event process. Num of input trajectories: {}",
0048         std::size(acts_trajectories));
0049
0050   // Loop over the trajectories
0051   for (std::size_t i = 0; const auto& traj : acts_trajectories) {
0052     // Get the entry index for the single trajectory
0053     // The trajectory entry indices and the multiTrajectory
0054     const auto& mj        = traj->multiTrajectory();
0055     const auto& trackTips = traj->tips();
0056     debug("------ Trajectory ------");
0057     debug("  Num of elements in trackTips {}", trackTips.size());
0058
0059     // Skip empty
0060     if (trackTips.empty()) {
0061       debug("  Empty multiTrajectory.");
0062       continue;
0063     }
0064     const auto& trackTip = trackTips.front();
0065
0066     // Collect the trajectory summary info
0067     auto trajState      = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
0068     int m_nMeasurements = trajState.nMeasurements;
0069     int m_nStates       = trajState.nStates;
0070     int m_nCalibrated   = 0;
0071     debug("  Num measurement in trajectory {}", m_nMeasurements);
0072     debug("  Num state in trajectory {}", m_nStates);
0073
0074     auto track_segment = track_segments->create();
0075
0076     // corresponding track
0077     if (tracks->size() == acts_trajectories.size()) {
0078       trace("track segment connected to track {}", i);
0079       track_segment.setTrack((*tracks)[i]);
0080       ++i;
0081     }
0082
0083     // visit the track points
0084     mj.visitBackwards(trackTip, [&](auto&& trackstate) {
0085       // get volume info
0086       auto geoID  = trackstate.referenceSurface().geometryId();
0087       auto volume = geoID.volume();
0088       auto layer  = geoID.layer();
0089
0090       if (trackstate.hasCalibrated()) {
0091         m_nCalibrated++;
0092       }
0093
0094       // get track state bound parameters and their boundCovs
0095       const auto& boundParams = trackstate.predicted();
0096       const auto& boundCov    = trackstate.predictedCovariance();
0097
0098       // convert local to global
0099       auto global = trackstate.referenceSurface().localToGlobal(
0100           m_geo_provider->getActsGeometryContext(),
0101           {boundParams[Acts::eBoundLoc0], boundParams[Acts::eBoundLoc1]},
0102           Acts::makeDirectionFromPhiTheta(boundParams[Acts::eBoundPhi],
0103                                           boundParams[Acts::eBoundTheta]));
0104
0105 #if Acts_VERSION_MAJOR >= 34
0106       auto freeParams = Acts::transformBoundToFreeParameters(
0107           trackstate.referenceSurface(), m_geo_provider->getActsGeometryContext(), boundParams);
0108       auto jacobian = trackstate.referenceSurface().boundToFreeJacobian(
0109           m_geo_provider->getActsGeometryContext(), freeParams.template segment<3>(Acts::eFreePos0),
0110           freeParams.template segment<3>(Acts::eFreeDir0));
0111 #else
0112                 auto jacobian = trackstate.referenceSurface().boundToFreeJacobian(
0113                         m_geo_provider->getActsGeometryContext(),
0114                         boundParams
0115                 );
0116 #endif
0117       auto freeCov = jacobian * boundCov * jacobian.transpose();
0118
0119       // global position
0120       const decltype(edm4eic::TrackPoint::position) position{static_cast<float>(global.x()),
0121                                                              static_cast<float>(global.y()),
0122                                                              static_cast<float>(global.z())};
0123
0124       // local position
0125       const decltype(edm4eic::TrackParametersData::loc) loc{
0126           static_cast<float>(boundParams[Acts::eBoundLoc0]),
0127           static_cast<float>(boundParams[Acts::eBoundLoc1])};
0128       const edm4eic::Cov2f locError{
0129           static_cast<float>(boundCov(Acts::eBoundLoc0, Acts::eBoundLoc0)),
0130           static_cast<float>(boundCov(Acts::eBoundLoc1, Acts::eBoundLoc1)),
0131           static_cast<float>(boundCov(Acts::eBoundLoc0, Acts::eBoundLoc1))};
0132       const decltype(edm4eic::TrackPoint::positionError) positionError{
0133           static_cast<float>(freeCov(Acts::eFreePos0, Acts::eFreePos0)),
0134           static_cast<float>(freeCov(Acts::eFreePos1, Acts::eFreePos1)),
0135           static_cast<float>(freeCov(Acts::eFreePos2, Acts::eFreePos2)),
0136           static_cast<float>(freeCov(Acts::eFreePos0, Acts::eFreePos1)),
0137           static_cast<float>(freeCov(Acts::eFreePos0, Acts::eFreePos2)),
0138           static_cast<float>(freeCov(Acts::eFreePos1, Acts::eFreePos2)),
0139       };
0140
0141       // momentum
0142       const decltype(edm4eic::TrackPoint::momentum) momentum = edm4hep::utils::sphericalToVector(
0143           static_cast<float>(1.0 / std::abs(boundParams[Acts::eBoundQOverP])),
0144           static_cast<float>(boundParams[Acts::eBoundTheta]),
0145           static_cast<float>(boundParams[Acts::eBoundPhi]));
0146       const decltype(edm4eic::TrackPoint::momentumError) momentumError{
0147           static_cast<float>(boundCov(Acts::eBoundTheta, Acts::eBoundTheta)),
0148           static_cast<float>(boundCov(Acts::eBoundPhi, Acts::eBoundPhi)),
0149           static_cast<float>(boundCov(Acts::eBoundQOverP, Acts::eBoundQOverP)),
0150           static_cast<float>(boundCov(Acts::eBoundTheta, Acts::eBoundPhi)),
0151           static_cast<float>(boundCov(Acts::eBoundTheta, Acts::eBoundQOverP)),
0152           static_cast<float>(boundCov(Acts::eBoundPhi, Acts::eBoundQOverP))};
0153       const float time{static_cast<float>(boundParams(Acts::eBoundTime))};
0154       const float timeError{static_cast<float>(sqrt(boundCov(Acts::eBoundTime, Acts::eBoundTime)))};
0155       const float theta(boundParams[Acts::eBoundTheta]);
0156       const float phi(boundParams[Acts::eBoundPhi]);
0157       const decltype(edm4eic::TrackPoint::directionError) directionError{
0158           static_cast<float>(boundCov(Acts::eBoundTheta, Acts::eBoundTheta)),
0159           static_cast<float>(boundCov(Acts::eBoundPhi, Acts::eBoundPhi)),
0160           static_cast<float>(boundCov(Acts::eBoundTheta, Acts::eBoundPhi))};
0161       const auto pathLength       = static_cast<float>(trackstate.pathLength());
0162       const float pathLengthError = 0;
0163
0164       uint64_t surface = trackstate.referenceSurface().geometryId().value();
0165       uint32_t system  = 0;
0166
0167       // Store track point
0168       track_segment.addToPoints({surface, system, position, positionError, momentum, momentumError,
0169                                  time, timeError, theta, phi, directionError, pathLength,
0170                                  pathLengthError});
0171
0172       debug("  ******************************");
0173       debug("    position: {}", position);
0174       debug("    positionError: {}", positionError);
0175       debug("    momentum: {}", momentum);
0176       debug("    momentumError: {}", momentumError);
0177       debug("    time: {}", time);
0178       debug("    timeError: {}", timeError);
0179       debug("    theta: {}", theta);
0180       debug("    phi: {}", phi);
0181       debug("    directionError: {}", directionError);
0182       debug("    pathLength: {}", pathLength);
0183       debug("    pathLengthError: {}", pathLengthError);
0184       debug("    geoID = {}", geoID);
0185       debug("    volume = {}, layer = {}", volume, layer);
0186       debug("    pathlength = {}", pathLength);
0187       debug("    hasCalibrated = {}", trackstate.hasCalibrated());
0188       debug("  ******************************");
0189
0190       // Local position on the reference surface.
0191       //debug("boundParams[eBoundLoc0] = {}", boundParams[Acts::eBoundLoc0]);
0192       //debug("boundParams[eBoundLoc1] = {}", boundParams[Acts::eBoundLoc1]);
0193       //debug("boundParams[eBoundPhi] = {}", boundParams[Acts::eBoundPhi]);
0194       //debug("boundParams[eBoundTheta] = {}", boundParams[Acts::eBoundTheta]);
0195       //debug("boundParams[eBoundQOverP] = {}", boundParams[Acts::eBoundQOverP]);
0196       //debug("boundParams[eBoundTime] = {}", boundParams[Acts::eBoundTime]);
0197       //debug("predicted variables: {}", trackstate.predicted());
0198     });
0199
0200     debug("  Num calibrated state in trajectory {}", m_nCalibrated);
0201     debug("------ end of trajectory process ------");
0202   }
0203
0204   debug("END OF Track projector event process");
0205 }
0206
0207 } // namespace eicrecon