src/components/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 <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/ReconstructedParticleCollection.h"
0027 #include "edm4eic/TrackerHitCollection.h"
0028 #include "edm4eic/TrackParametersCollection.h"
0029 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0030 #include "ActsExamples/EventData/Track.hpp"
0031 #include "ActsExamples/EventData/Trajectories.hpp"
0032
0033 #include "Acts/Utilities/Helpers.hpp"
0034
0035 #include "edm4hep/utils/vector_utils.h"
0036
0037 #include <cmath>
0038
0039 namespace Jug::Reco {
0040
0041   #if EDM4EIC_VERSION_MAJOR >= 5
0042     // This array relates the Acts and EDM4eic covariance matrices, including
0043     // the unit conversion to get from Acts units into EDM4eic units.
0044     //
0045     // Note: std::map is not constexpr, so we use a constexpr std::array
0046     // std::array initialization need double braces since arrays are aggregates
0047     // ref: https://en.cppreference.com/w/cpp/language/aggregate_initialization
0048     static constexpr std::array<std::pair<Acts::BoundIndices, double>, 6> edm4eic_indexed_units{{
0049       {Acts::eBoundLoc0, Acts::UnitConstants::mm},
0050       {Acts::eBoundLoc1, Acts::UnitConstants::mm},
0051       {Acts::eBoundTheta, 1.},
0052       {Acts::eBoundPhi, 1.},
0053       {Acts::eBoundQOverP, 1. / Acts::UnitConstants::GeV},
0054       {Acts::eBoundTime, Acts::UnitConstants::ns}
0055     }};
0056   #endif
0057
0058   /** Extract the particles form fit trajectories.
0059    *
0060    * \ingroup tracking
0061    */
0062    class ParticlesFromTrackFit : public GaudiAlgorithm {
0063    private:
0064     DataHandle<ActsExamples::TrajectoriesContainer>     m_inputTrajectories{"inputTrajectories", Gaudi::DataHandle::Reader, this};
0065     DataHandle<edm4eic::ReconstructedParticleCollection> m_outputParticles{"outputParticles", Gaudi::DataHandle::Writer, this};
0066     DataHandle<edm4eic::TrackParametersCollection> m_outputTrackParameters{"outputTrackParameters", Gaudi::DataHandle::Writer, this};
0067
0068    public:
0069     //  ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
0070     ParticlesFromTrackFit(const std::string& name, ISvcLocator* svcLoc)
0071         : GaudiAlgorithm(name, svcLoc) {
0072           declareProperty("inputTrajectories", m_inputTrajectories,"");
0073           declareProperty("outputParticles", m_outputParticles, "");
0074           declareProperty("outputTrackParameters", m_outputTrackParameters, "Acts Track Parameters");
0075         }
0076
0077     StatusCode initialize() override {
0078       if (GaudiAlgorithm::initialize().isFailure()) {
0079         return StatusCode::FAILURE;
0080       }
0081       return StatusCode::SUCCESS;
0082     }
0083
0084     StatusCode execute() override {
0085       // input collection
0086       const auto* const trajectories = m_inputTrajectories.get();
0087       // create output collections
0088       auto* rec_parts = m_outputParticles.createAndPut();
0089       auto* track_pars = m_outputTrackParameters.createAndPut();
0090
0091       if (msgLevel(MSG::DEBUG)) {
0092         debug() << std::size(*trajectories) << " trajectories " << endmsg;
0093       }
0094
0095       // Loop over the trajectories
0096         for (const auto& traj : *trajectories) {
0097
0098           // Get the entry index for the single trajectory
0099           // The trajectory entry indices and the multiTrajectory
0100           const auto& mj        = traj.multiTrajectory();
0101           const auto& trackTips = traj.tips();
0102           if (trackTips.empty()) {
0103             if (msgLevel(MSG::DEBUG)) {
0104               debug() << "Empty multiTrajectory." << endmsg;
0105             }
0106             continue;
0107           }
0108
0109           const auto& trackTip = trackTips.front();
0110
0111           // Collect the trajectory summary info
0112           auto trajState       = Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
0113           //int  m_nMeasurements = trajState.nMeasurements;
0114           //int  m_nStates       = trajState.nStates;
0115
0116           // Get the fitted track parameter
0117           //
0118           if (traj.hasTrackParameters(trackTip)) {
0119             const auto& boundParam = traj.trackParameters(trackTip);
0120             const auto& parameter  = boundParam.parameters();
0121             const auto& covariance = *boundParam.covariance();
0122             if (msgLevel(MSG::DEBUG)) {
0123               debug() << "loc 0 = " << parameter[Acts::eBoundLoc0] << endmsg;
0124               debug() << "loc 1 = " << parameter[Acts::eBoundLoc1] << endmsg;
0125               debug() << "phi   = " << parameter[Acts::eBoundPhi] << endmsg;
0126               debug() << "theta = " << parameter[Acts::eBoundTheta] << endmsg;
0127               debug() << "q/p   = " << parameter[Acts::eBoundQOverP] << endmsg;
0128               debug() << "p     = " << 1.0 / parameter[Acts::eBoundQOverP] << endmsg;
0129
0130               debug() << "err phi = " << sqrt(covariance(Acts::eBoundPhi, Acts::eBoundPhi)) << endmsg;
0131               debug() << "err th  = " << sqrt(covariance(Acts::eBoundTheta, Acts::eBoundTheta)) << endmsg;
0132               debug() << "err q/p = " << sqrt(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)) << endmsg;
0133
0134               debug() << " chi2 = " << trajState.chi2Sum << endmsg;
0135             }
0136
0137             auto pars = track_pars->create();
0138             pars.setType(0); // type: track head --> 0
0139             pars.setLoc({
0140               static_cast<float>(parameter[Acts::eBoundLoc0]),
0141               static_cast<float>(parameter[Acts::eBoundLoc1])
0142             });
0143             pars.setTheta(static_cast<float>(parameter[Acts::eBoundTheta]));
0144             pars.setPhi(static_cast<float>(parameter[Acts::eBoundPhi]));
0145             pars.setQOverP(static_cast<float>(parameter[Acts::eBoundQOverP]));
0146             pars.setTime(static_cast<float>(parameter[Acts::eBoundTime]));
0147             #if EDM4EIC_VERSION_MAJOR >= 5
0148               edm4eic::Cov6f cov;
0149               for (size_t i = 0; const auto& [a, x] : edm4eic_indexed_units) {
0150                 for (size_t j = 0; const auto& [b, y] : edm4eic_indexed_units) {
0151                   // FIXME why not pars.getCovariance()(i,j) = covariance(a,b) / x / y;
0152                   cov(i,j) = covariance(a,b) / x / y;
0153                 }
0154               }
0155               pars.setCovariance(cov);
0156             #else
0157               pars.setCharge(static_cast<float>(boundParam.charge()));
0158               pars.setLocError({
0159                     static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)),
0160                     static_cast<float>(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)),
0161                     static_cast<float>(covariance(Acts::eBoundLoc0, Acts::eBoundLoc1))
0162                 });
0163               pars.setMomentumError({
0164                     static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundTheta)),
0165                     static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundPhi)),
0166                     static_cast<float>(covariance(Acts::eBoundQOverP, Acts::eBoundQOverP)),
0167                     static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundPhi)),
0168                     static_cast<float>(covariance(Acts::eBoundTheta, Acts::eBoundQOverP)),
0169                     static_cast<float>(covariance(Acts::eBoundPhi, Acts::eBoundQOverP))
0170                 });
0171               pars.setTimeError(sqrt(static_cast<float>(covariance(Acts::eBoundTime, Acts::eBoundTime))));
0172             #endif
0173           }
0174
0175           auto tsize = trackTips.size();
0176           if (msgLevel(MSG::DEBUG)) {
0177             debug() << "# fitted parameters : " << tsize << endmsg;
0178           }
0179           if (tsize == 0) {
0180             continue;
0181           }
0182
0183           mj.visitBackwards(tsize - 1, [&](auto&& trackstate) {
0184             // debug() << trackstate.hasPredicted() << endmsg;
0185             // debug() << trackstate.predicted() << endmsg;
0186             auto params = trackstate.predicted(); //<< endmsg;
0187
0188             double p0 = (1.0 / params[Acts::eBoundQOverP]) / Acts::UnitConstants::GeV;
0189             if (msgLevel(MSG::DEBUG)) {
0190               debug() << "track predicted p = " << p0 << " GeV" << endmsg;
0191             }
0192             if (std::abs(p0) > 500) {
0193               if (msgLevel(MSG::DEBUG)) {
0194                 debug() << "skipping" << endmsg;
0195               }
0196               return;
0197             }
0198
0199             auto rec_part = rec_parts->create();
0200             rec_part.setMomentum(
0201               edm4hep::utils::sphericalToVector(
0202                 1.0 / std::abs(params[Acts::eBoundQOverP]),
0203                 params[Acts::eBoundTheta],
0204                 params[Acts::eBoundPhi])
0205             );
0206             rec_part.setCharge(static_cast<int16_t>(std::copysign(1., params[Acts::eBoundQOverP])));
0207           });
0208       }
0209
0210       return StatusCode::SUCCESS;
0211     }
0212
0213   };
0214   // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0215   DECLARE_COMPONENT(ParticlesFromTrackFit)
0216
0217 } // namespace Jug::Reco