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