File indexing completed on 2024-09-28 07:02:59
0001
0002
0003
0004 #include <Acts/Definitions/TrackParametrization.hpp>
0005 #include <Acts/EventData/MultiTrajectoryHelpers.hpp>
0006 #include <Acts/Geometry/GeometryIdentifier.hpp>
0007 #include <Acts/Utilities/UnitVectors.hpp>
0008 #include <ActsExamples/EventData/Trajectories.hpp>
0009 #include <edm4eic/Cov2f.h>
0010 #include <edm4eic/Cov3f.h>
0011 #include <edm4eic/TrackParametersCollection.h>
0012 #include <edm4eic/TrackPoint.h>
0013 #include <edm4eic/TrackSegmentCollection.h>
0014 #include <edm4hep/Vector3f.h>
0015 #include <edm4hep/utils/vector_utils.h>
0016 #include <fmt/core.h>
0017 #include <fmt/ostream.h>
0018 #include <spdlog/logger.h>
0019 #include <stdint.h>
0020 #include <cmath>
0021 #include <iterator>
0022 #include <utility>
0023
0024 #include "TrackProjector.h"
0025 #include "extensions/spdlog/SpdlogFormatters.h" // IWYU pragma: keep
0026
0027 #if FMT_VERSION >= 90000
0028 template<> struct fmt::formatter<Acts::GeometryIdentifier> : fmt::ostream_formatter {};
0029 #endif
0030
0031 namespace eicrecon {
0032
0033 void
0034 TrackProjector::init(std::shared_ptr<const ActsGeometryProvider> geo_svc, std::shared_ptr<spdlog::logger> logger) {
0035 m_log = logger;
0036 m_geo_provider = geo_svc;
0037 }
0038
0039
0040 std::unique_ptr<edm4eic::TrackSegmentCollection> TrackProjector::execute(std::vector<const ActsExamples::Trajectories *> trajectories) {
0041
0042
0043 auto track_segments = std::make_unique<edm4eic::TrackSegmentCollection>();
0044 m_log->debug("Track projector event process. Num of input trajectories: {}", std::size(trajectories));
0045
0046
0047 for (const auto &traj: trajectories) {
0048
0049
0050 const auto &mj = traj->multiTrajectory();
0051 const auto &trackTips = traj->tips();
0052 m_log->debug("------ Trajectory ------");
0053 m_log->debug(" Num of elements in trackTips {}", trackTips.size());
0054
0055
0056 if (trackTips.empty()) {
0057 m_log->debug(" Empty multiTrajectory.");
0058 continue;
0059 }
0060 const auto &trackTip = trackTips.front();
0061
0062
0063 auto trajState = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
0064 int m_nMeasurements = trajState.nMeasurements;
0065 int m_nStates = trajState.nStates;
0066 int m_nCalibrated = 0;
0067 m_log->debug(" Num measurement in trajectory {}", m_nMeasurements);
0068 m_log->debug(" Num state in trajectory {}", m_nStates);
0069
0070 auto track_segment = track_segments->create();
0071
0072
0073 mj.visitBackwards(trackTip, [&](auto &&trackstate) {
0074
0075 auto geoID = trackstate.referenceSurface().geometryId();
0076 auto volume = geoID.volume();
0077 auto layer = geoID.layer();
0078
0079 if (trackstate.hasCalibrated()) {
0080 m_nCalibrated++;
0081 }
0082
0083
0084 const auto &boundParams = trackstate.predicted();
0085 const auto &boundCov = trackstate.predictedCovariance();
0086
0087
0088 auto global = trackstate.referenceSurface().localToGlobal(
0089 m_geo_provider->getActsGeometryContext(),
0090 {boundParams[Acts::eBoundLoc0], boundParams[Acts::eBoundLoc1]},
0091 Acts::makeDirectionFromPhiTheta(
0092 boundParams[Acts::eBoundPhi],
0093 boundParams[Acts::eBoundTheta]
0094 )
0095 );
0096
0097 #if Acts_VERSION_MAJOR >= 34
0098 auto freeParams = Acts::transformBoundToFreeParameters(
0099 trackstate.referenceSurface(),
0100 m_geo_provider->getActsGeometryContext(),
0101 boundParams);
0102 auto jacobian = trackstate.referenceSurface().boundToFreeJacobian(
0103 m_geo_provider->getActsGeometryContext(),
0104 freeParams.template segment<3>(Acts::eFreePos0),
0105 freeParams.template segment<3>(Acts::eFreeDir0)
0106 );
0107 #else
0108 auto jacobian = trackstate.referenceSurface().boundToFreeJacobian(
0109 m_geo_provider->getActsGeometryContext(),
0110 boundParams
0111 );
0112 #endif
0113 auto freeCov = jacobian * boundCov * jacobian.transpose();
0114
0115
0116 const decltype(edm4eic::TrackPoint::position) position{
0117 static_cast<float>(global.x()),
0118 static_cast<float>(global.y()),
0119 static_cast<float>(global.z())
0120 };
0121
0122
0123 const decltype(edm4eic::TrackParametersData::loc) loc{
0124 static_cast<float>(boundParams[Acts::eBoundLoc0]),
0125 static_cast<float>(boundParams[Acts::eBoundLoc1])
0126 };
0127 const edm4eic::Cov2f locError{
0128 static_cast<float>(boundCov(Acts::eBoundLoc0, Acts::eBoundLoc0)),
0129 static_cast<float>(boundCov(Acts::eBoundLoc1, Acts::eBoundLoc1)),
0130 static_cast<float>(boundCov(Acts::eBoundLoc0, Acts::eBoundLoc1))
0131 };
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
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 );
0147 const decltype(edm4eic::TrackPoint::momentumError) momentumError{
0148 static_cast<float>(boundCov(Acts::eBoundTheta, Acts::eBoundTheta)),
0149 static_cast<float>(boundCov(Acts::eBoundPhi, Acts::eBoundPhi)),
0150 static_cast<float>(boundCov(Acts::eBoundQOverP, Acts::eBoundQOverP)),
0151 static_cast<float>(boundCov(Acts::eBoundTheta, Acts::eBoundPhi)),
0152 static_cast<float>(boundCov(Acts::eBoundTheta, Acts::eBoundQOverP)),
0153 static_cast<float>(boundCov(Acts::eBoundPhi, Acts::eBoundQOverP))
0154 };
0155 const float time{static_cast<float>(boundParams(Acts::eBoundTime))};
0156 const float timeError{sqrt(static_cast<float>(boundCov(Acts::eBoundTime, Acts::eBoundTime)))};
0157 const float theta(boundParams[Acts::eBoundTheta]);
0158 const float phi(boundParams[Acts::eBoundPhi]);
0159 const decltype(edm4eic::TrackPoint::directionError) directionError{
0160 static_cast<float>(boundCov(Acts::eBoundTheta, Acts::eBoundTheta)),
0161 static_cast<float>(boundCov(Acts::eBoundPhi, Acts::eBoundPhi)),
0162 static_cast<float>(boundCov(Acts::eBoundTheta, Acts::eBoundPhi))
0163 };
0164 const float 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
0171 track_segment.addToPoints({
0172 surface,
0173 system,
0174 position,
0175 positionError,
0176 momentum,
0177 momentumError,
0178 time,
0179 timeError,
0180 theta,
0181 phi,
0182 directionError,
0183 pathLength,
0184 pathLengthError
0185 });
0186
0187
0188 m_log->debug(" ******************************");
0189 m_log->debug(" position: {}", position);
0190 m_log->debug(" positionError: {}", positionError);
0191 m_log->debug(" momentum: {}", momentum);
0192 m_log->debug(" momentumError: {}", momentumError);
0193 m_log->debug(" time: {}", time);
0194 m_log->debug(" timeError: {}", timeError);
0195 m_log->debug(" theta: {}", theta);
0196 m_log->debug(" phi: {}", phi);
0197 m_log->debug(" directionError: {}", directionError);
0198 m_log->debug(" pathLength: {}", pathLength);
0199 m_log->debug(" pathLengthError: {}", pathLengthError);
0200 m_log->debug(" geoID = {}", geoID);
0201 m_log->debug(" volume = {}, layer = {}", volume, layer);
0202 m_log->debug(" pathlength = {}", pathLength);
0203 m_log->debug(" hasCalibrated = {}", trackstate.hasCalibrated());
0204 m_log->debug(" ******************************");
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214 });
0215
0216 m_log->debug(" Num calibrated state in trajectory {}", m_nCalibrated);
0217 m_log->debug("------ end of trajectory process ------");
0218 }
0219
0220 m_log->debug("END OF Track projector event process");
0221 return std::move(track_segments);
0222 }
0223
0224
0225 }