acts/src/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 "JugBase/DataHandle.h"
0018 #include "JugBase/IGeoSvc.h"
0019 #include "JugTrack/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   using Clusters   = edm4eic::ClusterCollection;
0041   using VertexHits = edm4eic::TrackerHitCollection;
0042
0043   DataHandle<VertexHits> m_inputVertexHits{"inputVertexHits", Gaudi::DataHandle::Reader, this};
0044   DataHandle<Clusters> m_inputClusters{"inputClusters", Gaudi::DataHandle::Reader, this};
0045   DataHandle<TrackParametersContainer> m_outputInitialTrackParameters{"outputInitialTrackParameters",
0046                                                                       Gaudi::DataHandle::Writer, this};
0047   Gaudi::Property<double> m_maxHitRadius{this, "maxHitRadius", 40.0 * mm};
0048
0049 public:
0050   TrackParamVertexClusterInit(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
0051     declareProperty("inputVertexHits", m_inputVertexHits, "Vertex tracker hits");
0052     declareProperty("inputClusters", m_inputClusters, "Input clusters");
0053     declareProperty("outputInitialTrackParameters", m_outputInitialTrackParameters, "");
0054   }
0055
0056   StatusCode initialize() override {
0057     if (GaudiAlgorithm::initialize().isFailure()) {
0058       return StatusCode::FAILURE;
0059     }
0060     IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
0061     if (randSvc == nullptr) {
0062       return StatusCode::FAILURE;
0063     }
0064     return StatusCode::SUCCESS;
0065   }
0066
0067   StatusCode execute() override {
0068     // input collection
0069     const auto* const clusters = m_inputClusters.get();
0070     const auto* const vtx_hits = m_inputVertexHits.get();
0071     // Create output collections
0072     auto* init_trk_params = m_outputInitialTrackParameters.createAndPut();
0073
0074     double max_radius = m_maxHitRadius.value();
0075
0076     for (const auto& c : *clusters) {
0077
0078       using Acts::UnitConstants::GeV;
0079       using Acts::UnitConstants::MeV;
0080       using Acts::UnitConstants::mm;
0081       using Acts::UnitConstants::ns;
0082
0083       double p_cluster = c.getEnergy() * GeV;
0084       if (p_cluster / GeV < 0.1) {
0085         if (msgLevel(MSG::DEBUG)) {
0086           debug() << " skipping cluster with energy " << p_cluster / GeV << " GeV" << endmsg;
0087         }
0088         continue;
0089       }
0090
0091       auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3{0, 0, 0});
0092       for (const auto& t : *vtx_hits) {
0093
0094         double len = std::hypot(t.getPosition().x, t.getPosition().y, t.getPosition().z);
0095         if (len > max_radius) {
0096           continue;
0097         }
0098
0099         auto momentum = t.getPosition() * p_cluster / len;
0100
0101         Acts::BoundVector params;
0102         params(Acts::eBoundLoc0)   = 0.0 * mm;
0103         params(Acts::eBoundLoc1)   = 0.0 * mm;
0104         params(Acts::eBoundPhi)    = edm4hep::utils::angleAzimuthal(momentum);
0105         params(Acts::eBoundTheta)  = edm4hep::utils::anglePolar(momentum);
0106         params(Acts::eBoundQOverP) = 1 / p_cluster;
0107         params(Acts::eBoundTime)   = 0 * ns;
0108
0109         // debug() << "Invoke track finding seeded by truth particle with p = " << p / GeV << " GeV" << endmsg;
0110
0111         // add both charges to the track candidate...
0112         init_trk_params->push_back({pSurface, params, 1});
0113
0114         Acts::BoundVector params2;
0115         params2(Acts::eBoundLoc0)   = 0.0 * mm;
0116         params2(Acts::eBoundLoc1)   = 0.0 * mm;
0117         params2(Acts::eBoundPhi)    = edm4hep::utils::angleAzimuthal(momentum);
0118         params2(Acts::eBoundTheta)  = edm4hep::utils::anglePolar(momentum);
0119         params2(Acts::eBoundQOverP) = -1 / p_cluster;
0120         params2(Acts::eBoundTime)   = 0 * ns;
0121         init_trk_params->push_back({pSurface, params2, -1});
0122       }
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, 1,
0125       //                              std::make_optional(cov));
0126       // init_trk_params->emplace_back(Acts::Vector4(0 * mm, 0 * mm, 0 * mm, 0),
0127       //                              Acts::Vector3(c.x() * p / len, c.y() * p / len, c.z() * p / len), p, 0,
0128       //                              std::make_optional(cov));
0129       // debug() << "Invoke track finding seeded by truth particle with p = " << p_cluster/GeV  << " GeV" << endmsg;
0130     }
0131     return StatusCode::SUCCESS;
0132   }
0133 };
0134 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0135 DECLARE_COMPONENT(TrackParamVertexClusterInit)
0136
0137 } // namespace Jug::Reco