src/components/TrackProjector.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 wfan, Whitney Armstrong, Sylvester Joosten
0003
0004 #include <algorithm>
0005
0006 // Gaudi
0007 #include "GaudiAlg/GaudiAlgorithm.h"
0008 #include "GaudiKernel/ToolHandle.h"
0009 #include "GaudiAlg/Transformer.h"
0010 #include "GaudiAlg/GaudiTool.h"
0011 #include "GaudiKernel/RndmGenerators.h"
0012 #include "Gaudi/Property.h"
0013
0014 #include "DDRec/CellIDPositionConverter.h"
0015 #include "DDRec/SurfaceManager.h"
0016 #include "DDRec/Surface.h"
0017
0018 #include <k4FWCore/DataHandle.h>
0019 #include <k4Interface/IGeoSvc.h>
0020
0021 #include "Acts/EventData/MultiTrajectory.hpp"
0022 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
0023
0024 // Event Model related classes
0025 #include "edm4eic/EDM4eicVersion.h"
0026 #include "edm4eic/TrackerHitCollection.h"
0027 #include "edm4eic/TrackParametersCollection.h"
0028 #include "edm4eic/TrajectoryCollection.h"
0029 #include "edm4eic/TrackSegmentCollection.h"
0030 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0031 #include "ActsExamples/EventData/Track.hpp"
0032 #include "ActsExamples/EventData/Trajectories.hpp"
0033
0034 #include "Acts/Utilities/Helpers.hpp"
0035 #include "Acts/Geometry/GeometryIdentifier.hpp"
0036 #include "Acts/MagneticField/ConstantBField.hpp"
0037 #include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
0038 #include "Acts/Propagator/EigenStepper.hpp"
0039 #include "Acts/Surfaces/PerigeeSurface.hpp"
0040
0041 #include "edm4hep/utils/vector_utils.h"
0042
0043 #include <cmath>
0044
0045 namespace Jug::Reco {
0046
0047   /** Extrac the particles form fit trajectories.
0048    *
0049    * \ingroup tracking
0050    */
0051    class TrackProjector : public GaudiAlgorithm {
0052    private:
0053     DataHandle<ActsExamples::TrajectoriesContainer>        m_inputTrajectories{"inputTrajectories", Gaudi::DataHandle::Reader, this};
0054     DataHandle<edm4eic::TrackSegmentCollection> m_outputTrackSegments{"outputTrackSegments", Gaudi::DataHandle::Writer, this};
0055
0056     Gaudi::Property<unsigned int> m_firstInVolumeID{this, "firstInVolumeID", 0};
0057     Gaudi::Property<std::string> m_firstInVolumeName{this, "firstInVolumeName", ""};
0058     Gaudi::Property<float> m_firstSmallerThanZ{this, "firstSmallerThanZ", 0};
0059     Gaudi::Property<float> m_firstGreaterThanZ{this, "firstGreaterThanZ", 0};
0060     Gaudi::Property<float> m_firstGreaterThanR{this, "firstGreaterThanR", -1};
0061
0062     Acts::GeometryContext m_geoContext;
0063
0064     public:
0065     //  ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
0066     TrackProjector(const std::string& name, ISvcLocator* svcLoc)
0067         : GaudiAlgorithm(name, svcLoc) {
0068           declareProperty("inputTrajectories", m_inputTrajectories,"");
0069           declareProperty("outputTrackSegments", m_outputTrackSegments, "");
0070         }
0071
0072     StatusCode initialize() override {
0073       if (GaudiAlgorithm::initialize().isFailure())
0074         return StatusCode::FAILURE;
0075       return StatusCode::SUCCESS;
0076     }
0077
0078     StatusCode execute() override {
0079       // input collection
0080       const auto* const trajectories = m_inputTrajectories.get();
0081       // create output collections
0082       auto* track_segments = m_outputTrackSegments.createAndPut();
0083
0084       if (msgLevel(MSG::DEBUG)) {
0085         debug() << std::size(*trajectories) << " trajectories " << endmsg;
0086       }
0087
0088       // Loop over the trajectories
0089       for (const auto& traj : *trajectories) {
0090         // Get the entry index for the single trajectory
0091         // The trajectory entry indices and the multiTrajectory
0092         const auto& mj        = traj.multiTrajectory();
0093         const auto& trackTips = traj.tips();
0094         if (msgLevel(MSG::DEBUG)) {
0095           debug() << "# of elements in trackTips " <<trackTips.size() << endmsg;
0096         }
0097
0098         // Skip empty
0099         if (trackTips.empty()) {
0100           if (msgLevel(MSG::DEBUG)) {
0101             debug() << "Empty multiTrajectory." << endmsg;
0102           }
0103           continue;
0104         }
0105         auto& trackTip = trackTips.front();
0106
0107         // Collect the trajectory summary info
0108         auto trajState       = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
0109         int  m_nMeasurements = trajState.nMeasurements;
0110         int  m_nStates       = trajState.nStates;
0111         int  m_nCalibrated   = 0;
0112         if (msgLevel(MSG::DEBUG)) {
0113           debug() << "n measurement in trajectory " << m_nMeasurements << endmsg;
0114           debug() << "n state in trajectory " << m_nStates << endmsg;
0115         }
0116
0117         edm4eic::MutableTrackSegment track_segment;
0118
0119         // visit the track points
0120         mj.visitBackwards(trackTip, [&](auto&& trackstate) {
0121           // get volume info
0122           auto geoID = trackstate.referenceSurface().geometryId();
0123           auto volume = geoID.volume();
0124           auto layer = geoID.layer();
0125           if (trackstate.hasCalibrated()) {
0126             m_nCalibrated++;
0127           }
0128
0129           // get track state parameters and their covariances
0130           const auto& parameter = trackstate.predicted();
0131           const auto& covariance = trackstate.predictedCovariance();
0132
0133           // convert local to global
0134           auto global = trackstate.referenceSurface().localToGlobal(
0135             m_geoContext,
0136             {parameter[Acts::eBoundLoc0], parameter[Acts::eBoundLoc1]},
0137             {0, 0, 0}
0138           );
0139           // global position
0140           const decltype(edm4eic::TrackPoint::position) position {
0141             static_cast<float>(global.x()),
0142             static_cast<float>(global.y()),
0143             static_cast<float>(global.z())
0144           };
0145
0146           // local position
0147           const decltype(edm4eic::TrackParametersData::loc) loc {
0148             static_cast<float>(parameter[Acts::eBoundLoc0]),
0149             static_cast<float>(parameter[Acts::eBoundLoc1])
0150           };
0151           const edm4eic::Cov2f locError {
0152             static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)),
0153             static_cast<float>(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)),
0154             static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc1))
0155           };
0156           const decltype(edm4eic::TrackPoint::positionError) positionError{0, 0, 0};
0157           const decltype(edm4eic::TrackPoint::momentum) momentum = edm4hep::utils::sphericalToVector(
0158             static_cast<float>(1.0 / std::abs(parameter[Acts::eBoundQOverP])),
0159             static_cast<float>(parameter[Acts::eBoundTheta]),
0160             static_cast<float>(parameter[Acts::eBoundPhi])
0161           );
0162           const decltype(edm4eic::TrackPoint::momentumError) momentumError {
0163             static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundTheta)),
0164             static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundPhi)),
0165             static_cast<float>(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)),
0166             static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundPhi)),
0167             static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundQOverP)),
0168             static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundQOverP))
0169           };
0170           const float time{static_cast<float>(parameter(Acts::eBoundTime))};
0171           const float timeError{sqrt(static_cast<float>(covariance(Acts::eBoundTime, Acts::eBoundTime)))};
0172           const float theta(parameter[Acts::eBoundTheta]);
0173           const float phi(parameter[Acts::eBoundPhi]);
0174           const decltype(edm4eic::TrackPoint::directionError) directionError {
0175             static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundTheta)),
0176             static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundPhi)),
0177             static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundPhi))
0178           };
0179           const float pathLength = static_cast<float>(trackstate.pathLength());
0180           const float pathLengthError = 0;
0181
0182 #if EDM4EIC_VERSION_MAJOR >= 3
0183           uint64_t surface = 0; // trackstate.referenceSurface().geometryId().value(); FIXME - ASAN is not happy with this
0184           uint32_t system = 0;
0185 #endif
0186
0187           // Store track point
0188           track_segment.addToPoints({
0189 #if EDM4EIC_VERSION_MAJOR >= 3
0190             surface,
0191             system,
0192 #endif
0193             position,
0194             positionError,
0195             momentum,
0196             momentumError,
0197             time,
0198             timeError,
0199             theta,
0200             phi,
0201             directionError,
0202             pathLength,
0203             pathLengthError
0204           });
0205
0206           if (msgLevel(MSG::DEBUG)) {
0207             debug() << "******************************" << endmsg;
0208             debug() << "predicted variables: \n" << trackstate.predicted() << endmsg;
0209             debug() << "geoID = " << geoID << endmsg;
0210             debug() << "volume = " << volume << ", layer = " << layer << endmsg;
0211             debug() << "pathlength = " << pathLength << endmsg;
0212             debug() << "hasCalibrated = " << trackstate.hasCalibrated() << endmsg;
0213             debug() << "******************************" << endmsg;
0214           }
0215         });
0216
0217         if (msgLevel(MSG::DEBUG)) {
0218           debug() << "n calibrated state in trajectory " << m_nCalibrated << endmsg;
0219         }
0220
0221         // Add to output collection
0222         track_segments->push_back(track_segment);
0223       }
0224
0225       return StatusCode::SUCCESS;
0226     }
0227
0228   };
0229   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0230   DECLARE_COMPONENT(TrackProjector)
0231
0232 } // namespace Jug::Reco