src/components/TrackParamVertexClusterInit.cpp

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