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