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