acts/src/ParticlesFromTrackFit.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Sylvester Joosten, Wouter Deconinck
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/ReconstructedParticleCollection.h"
0026 #include "edm4eic/TrackerHitCollection.h"
0027 #include "edm4eic/TrackParametersCollection.h"
0028 #include "JugTrack/IndexSourceLink.hpp"
0029 #include "JugTrack/Track.hpp"
0030 #include "JugTrack/Trajectories.hpp"
0031
0032 #include "Acts/Utilities/Helpers.hpp"
0033
0034 #include "edm4hep/utils/vector_utils.h"
0035
0036 #include <cmath>
0037
0038 namespace Jug::Reco {
0039
0040   /** Extract the particles form fit trajectories.
0041    *
0042    * \ingroup tracking
0043    */
0044    class ParticlesFromTrackFit : public GaudiAlgorithm {
0045    private:
0046     DataHandle<TrajectoriesContainer>     m_inputTrajectories{"inputTrajectories", Gaudi::DataHandle::Reader, this};
0047     DataHandle<edm4eic::ReconstructedParticleCollection> m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer, this};
0048     DataHandle<edm4eic::TrackParametersCollection> m_outputTrackParameters{"outputTrackParameters", Gaudi::DataHandle::Writer, this};
0049
0050    public:
0051     //  ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
0052     ParticlesFromTrackFit(const std::string& name, ISvcLocator* svcLoc)
0053         : GaudiAlgorithm(name, svcLoc) {
0054           declareProperty("inputTrajectories", m_inputTrajectories,"");
0055           declareProperty("outputParticles", m_outputParticles, "");
0056           declareProperty("outputTrackParameters", m_outputTrackParameters, "Acts Track Parameters");
0057         }
0058
0059     StatusCode initialize() override {
0060       if (GaudiAlgorithm::initialize().isFailure()) {
0061         return StatusCode::FAILURE;
0062       }
0063       return StatusCode::SUCCESS;
0064     }
0065
0066     StatusCode execute() override {
0067       // input collection
0068       const auto* const trajectories = m_inputTrajectories.get();
0069       // create output collections
0070       auto* rec_parts = m_outputParticles.createAndPut();
0071       auto* track_pars = m_outputTrackParameters.createAndPut();
0072
0073       if (msgLevel(MSG::DEBUG)) {
0074         debug() << std::size(*trajectories) << " trajectories " << endmsg;
0075       }
0076
0077       // Loop over the trajectories
0078         for (const auto& traj : *trajectories) {
0079
0080           // Get the entry index for the single trajectory
0081           // The trajectory entry indices and the multiTrajectory
0082           const auto& mj        = traj.multiTrajectory();
0083           const auto& trackTips = traj.tips();
0084           if (trackTips.empty()) {
0085             if (msgLevel(MSG::DEBUG)) {
0086               debug() << "Empty multiTrajectory." << endmsg;
0087             }
0088             continue;
0089           }
0090
0091           const auto& trackTip = trackTips.front();
0092
0093           // Collect the trajectory summary info
0094           auto trajState       = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
0095           //int  m_nMeasurements = trajState.nMeasurements;
0096           //int  m_nStates       = trajState.nStates;
0097
0098           // Get the fitted track parameter
0099           //
0100           if (traj.hasTrackParameters(trackTip)) {
0101             const auto& boundParam = traj.trackParameters(trackTip);
0102             const auto& parameter  = boundParam.parameters();
0103             const auto& covariance = *boundParam.covariance();
0104             if (msgLevel(MSG::DEBUG)) {
0105               debug() << "loc 0 = " << parameter[Acts::eBoundLoc0] << endmsg;
0106               debug() << "loc 1 = " << parameter[Acts::eBoundLoc1] << endmsg;
0107               debug() << "phi   = " << parameter[Acts::eBoundPhi] << endmsg;
0108               debug() << "theta = " << parameter[Acts::eBoundTheta] << endmsg;
0109               debug() << "q/p   = " << parameter[Acts::eBoundQOverP] << endmsg;
0110               debug() << "p     = " << 1.0 / parameter[Acts::eBoundQOverP] << endmsg;
0111
0112               debug() << "err phi = " << sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)) << endmsg;
0113               debug() << "err th  = " << sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)) << endmsg;
0114               debug() << "err q/p = " << sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)) << endmsg;
0115
0116               debug() << " chi2 = " << trajState.chi2Sum << endmsg;
0117             }
0118
0119             const decltype(edm4eic::TrackParametersData::loc) loc {
0120               static_cast<float>(parameter[Acts::eBoundLoc0]),
0121               static_cast<float>(parameter[Acts::eBoundLoc1])
0122             };
0123             const decltype(edm4eic::TrackParametersData::momentumError) momentumError {
0124               static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundTheta)),
0125               static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundPhi)),
0126               static_cast<float>(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)),
0127               static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundPhi)),
0128               static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundQOverP)),
0129               static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundQOverP))};
0130             const decltype(edm4eic::TrackParametersData::locError) locError {
0131               static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)),
0132               static_cast<float>(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)),
0133               static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc1))};
0134             const float timeError{sqrt(static_cast<float>(covariance(Acts::eBoundTime, Acts::eBoundTime)))};
0135
0136             edm4eic::TrackParameters pars{
0137               0, // type: track head --> 0
0138               loc,
0139               locError,
0140               static_cast<float>(parameter[Acts::eBoundTheta]),
0141               static_cast<float>(parameter[Acts::eBoundPhi]),
0142               static_cast<float>(parameter[Acts::eBoundQOverP]),
0143               momentumError,
0144               static_cast<float>(parameter[Acts::eBoundTime]),
0145               timeError,
0146               static_cast<float>(boundParam.charge())};
0147             track_pars->push_back(pars);
0148           }
0149
0150           auto tsize = trackTips.size();
0151           if (msgLevel(MSG::DEBUG)) {
0152             debug() << "# fitted parameters : " << tsize << endmsg;
0153           }
0154           if (tsize == 0) {
0155             continue;
0156           }
0157
0158           mj.visitBackwards(tsize - 1, [&](auto&& trackstate) {
0159             // debug() << trackstate.hasPredicted() << endmsg;
0160             // debug() << trackstate.predicted() << endmsg;
0161             auto params = trackstate.predicted(); //<< endmsg;
0162
0163             double p0 = (1.0 / params[Acts::eBoundQOverP]) / Acts::UnitConstants::GeV;
0164             if (msgLevel(MSG::DEBUG)) {
0165               debug() << "track predicted p = " << p0 << " GeV" << endmsg;
0166             }
0167             if (std::abs(p0) > 500) {
0168               if (msgLevel(MSG::DEBUG)) {
0169                 debug() << "skipping" << endmsg;
0170               }
0171               return;
0172             }
0173
0174             auto rec_part = rec_parts->create();
0175             rec_part.setMomentum(
0176               edm4hep::utils::sphericalToVector(
0177                 1.0 / std::abs(params[Acts::eBoundQOverP]),
0178                 params[Acts::eBoundTheta],
0179                 params[Acts::eBoundPhi])
0180             );
0181             rec_part.setCharge(static_cast<int16_t>(std::copysign(1., params[Acts::eBoundQOverP])));
0182           });
0183       }
0184
0185       return StatusCode::SUCCESS;
0186     }
0187
0188   };
0189   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0190   DECLARE_COMPONENT(ParticlesFromTrackFit)
0191
0192 } // namespace Jug::Reco