src/components/TrackParamClusterInit.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Sylvester Joosten
0003
0004 #include <cmath>
0005 // Gaudi
0006 #include "Gaudi/Property.h"
0007 #include "GaudiAlg/GaudiAlgorithm.h"
0008 #include "GaudiAlg/GaudiTool.h"
0009 #include "GaudiAlg/Transformer.h"
0010 #include "GaudiKernel/RndmGenerators.h"
0011 #include "GaudiKernel/ToolHandle.h"
0012
0013 #include "Acts/Definitions/Common.hpp"
0014 #include "Acts/Definitions/Units.hpp"
0015 #include <k4FWCore/DataHandle.h>
0016 #include <k4Interface/IGeoSvc.h>
0017 #include "ActsExamples/EventData/Track.hpp"
0018
0019 #include "edm4eic/ClusterCollection.h"
0020 #include "edm4eic/TrackerHitCollection.h"
0021 #include "edm4hep/utils/vector_utils.h"
0022
0023 #include "Acts/Surfaces/PerigeeSurface.hpp"
0024
0025 ///// (Reconstructed) track parameters e.g. close to the vertex.
0026 // using TrackParameters = Acts::CurvilinearTrackParameters;
0027
0028 ///// Container of reconstructed track states for multiple tracks.
0029 // using TrackParametersContainer = std::vector<TrackParameters>;
0030
0031 ///// MultiTrajectory definition
0032 // using Trajectory = Acts::MultiTrajectory<SourceLink>;
0033
0034 ///// Container for the truth fitting/finding track(s)
0035 // using TrajectoryContainer = std::vector<SimMultiTrajectory>;
0036
0037 namespace Jug::Reco {
0038
0039 /** Initial Track parameters from MC truth.
0040  *
0041  *  TrackParmetersContainer
0042  *
0043  *  \ingroup tracking
0044  */
0045 class TrackParamClusterInit : public GaudiAlgorithm {
0046 private:
0047   DataHandle<edm4eic::ClusterCollection> m_inputClusters{"inputClusters", Gaudi::DataHandle::Reader, this};
0048   DataHandle<ActsExamples::TrackParametersContainer> m_outputInitialTrackParameters{"outputInitialTrackParameters",
0049                                                                       Gaudi::DataHandle::Writer, this};
0050
0051 public:
0052   TrackParamClusterInit(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
0053     declareProperty("inputClusters", m_inputClusters, "Input clusters");
0054     declareProperty("outputInitialTrackParameters", m_outputInitialTrackParameters, "");
0055   }
0056
0057   StatusCode initialize() override {
0058     if (GaudiAlgorithm::initialize().isFailure()) {
0059       return StatusCode::FAILURE;
0060     }
0061     IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
0062     if (randSvc == nullptr) {
0063       return StatusCode::FAILURE;
0064     }
0065     return StatusCode::SUCCESS;
0066   }
0067
0068   StatusCode execute() override {
0069     // input collection
0070     const auto* const clusters = m_inputClusters.get();
0071     // Create output collections
0072     auto* init_trk_params = m_outputInitialTrackParameters.createAndPut();
0073
0074     for (const auto& c : *clusters) {
0075
0076       using Acts::UnitConstants::GeV;
0077       using Acts::UnitConstants::MeV;
0078       using Acts::UnitConstants::mm;
0079       using Acts::UnitConstants::ns;
0080
0081       double p = c.getEnergy() * GeV;
0082       if (p < 0.1 * GeV) {
0083         continue;
0084       }
0085       double len    = edm4hep::utils::magnitude(c.getPosition());
0086       auto momentum = c.getPosition() * p / len;
0087
0088       Acts::BoundVector params;
0089       params(Acts::eBoundLoc0)   = 0.0 * mm;
0090       params(Acts::eBoundLoc1)   = 0.0 * mm;
0091       params(Acts::eBoundPhi)    = edm4hep::utils::angleAzimuthal(momentum);
0092       params(Acts::eBoundTheta)  = edm4hep::utils::anglePolar(momentum);
0093       params(Acts::eBoundQOverP) = 1 / p;
0094       params(Acts::eBoundTime)   = 0 * ns;
0095
0096       auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0, 0, 0});
0097
0098       if (msgLevel(MSG::DEBUG)) {
0099         debug() << "Invoke track finding seeded by truth particle with p = " << p / GeV << " GeV" << endmsg;
0100       }
0101
0102       // add both charges to the track candidate...
0103       init_trk_params->push_back({pSurface, params, {}, Acts::ParticleHypothesis::pion()});
0104
0105       Acts::BoundVector params2;
0106       params2(Acts::eBoundLoc0)   = 0.0 * mm;
0107       params2(Acts::eBoundLoc1)   = 0.0 * mm;
0108       params2(Acts::eBoundPhi)    = edm4hep::utils::angleAzimuthal(momentum);
0109       params2(Acts::eBoundTheta)  = edm4hep::utils::anglePolar(momentum);
0110       params2(Acts::eBoundQOverP) = -1 / p;
0111       params2(Acts::eBoundTime)   = 0 * ns;
0112       init_trk_params->push_back({pSurface, params2, {}, Acts::ParticleHypothesis::pion()});
0113
0114       // acts v1.2.0:
0115       // init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
0116       //                              Acts::Vector3(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, -1,
0117       //                              std::make_optional(cov));
0118       // debug() << init_trk_params->back() << endmsg;
0119       // init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
0120       //                              Acts::Vector3(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 1,
0121       //                              std::make_optional(cov));
0122       ////debug() << init_trk_params->back() << endmsg;
0123       // init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
0124       //                              Acts::Vector3(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 0,
0125       //                              std::make_optional(cov));
0126     }
0127     return StatusCode::SUCCESS;
0128   }
0129 };
0130 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0131 DECLARE_COMPONENT(TrackParamClusterInit)
0132
0133 } // namespace Jug::Reco