Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2024-09-27 07:03:01 

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 wfan, Whitney Armstrong, Sylvester Joosten
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 // FMT_VERSION >= 90000
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         // create output collections
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         // Loop over the trajectories
0047         for (const auto &traj: trajectories) {
0048             // Get the entry index for the single trajectory
0049             // The trajectory entry indices and the multiTrajectory
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             // Skip empty
0056             if (trackTips.empty()) {
0057                 m_log->debug("  Empty multiTrajectory.");
0058                 continue;
0059             }
0060             const auto &trackTip = trackTips.front();
0061 
0062             // Collect the trajectory summary info
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             // visit the track points
0073             mj.visitBackwards(trackTip, [&](auto &&trackstate) {
0074                 // get volume info
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                 // get track state bound parameters and their boundCovs
0084                 const auto &boundParams = trackstate.predicted();
0085                 const auto &boundCov = trackstate.predictedCovariance();
0086 
0087                 // convert local to global
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                 // global position
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                 // local position
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                 // 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                 );
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                 // Store track point
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                 // Local position on the reference surface.
0207                 //m_log->debug("boundParams[eBoundLoc0] = {}", boundParams[Acts::eBoundLoc0]);
0208                 //m_log->debug("boundParams[eBoundLoc1] = {}", boundParams[Acts::eBoundLoc1]);
0209                 //m_log->debug("boundParams[eBoundPhi] = {}", boundParams[Acts::eBoundPhi]);
0210                 //m_log->debug("boundParams[eBoundTheta] = {}", boundParams[Acts::eBoundTheta]);
0211                 //m_log->debug("boundParams[eBoundQOverP] = {}", boundParams[Acts::eBoundQOverP]);
0212                 //m_log->debug("boundParams[eBoundTime] = {}", boundParams[Acts::eBoundTime]);
0213                 //m_log->debug("predicted variables: {}", trackstate.predicted());
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 } // namespace eicrecon::Reco
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!