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