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