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