File indexing completed on 2025-07-10 08:57:28
0001
0002
0003
0004 #include <cmath>
0005
0006 #include "Gaudi/Algorithm.h"
0007 #include "GaudiKernel/ToolHandle.h"
0008 #include "GaudiKernel/RndmGenerators.h"
0009 #include "Gaudi/Property.h"
0010
0011 #include <k4FWCore/DataHandle.h>
0012 #include <k4Interface/IGeoSvc.h>
0013 #include "JugBase/IParticleSvc.h"
0014 #include "ActsExamples/EventData/Track.hpp"
0015
0016
0017
0018
0019
0020 #include "edm4eic/EDM4eicVersion.h"
0021 #include "edm4hep/MCParticleCollection.h"
0022 #include "edm4eic/TrackParametersCollection.h"
0023 #include "Math/Vector3D.h"
0024
0025
0026 namespace Jug::Reco {
0027
0028
0029
0030
0031
0032
0033
0034
0035 class TruthTrackSeeding : public Gaudi::Algorithm {
0036 private:
0037 mutable DataHandle<edm4hep::MCParticleCollection> m_inputMCParticles{"inputMCParticles", Gaudi::DataHandle::Reader,
0038 this};
0039 mutable DataHandle<edm4eic::TrackParametersCollection> m_outputTrackParameters{"outputTrackParameters",
0040 Gaudi::DataHandle::Writer, this};
0041 SmartIF<IParticleSvc> m_pidSvc;
0042
0043 public:
0044 TruthTrackSeeding(const std::string& name, ISvcLocator* svcLoc) : Gaudi::Algorithm(name, svcLoc) {
0045 declareProperty("inputMCParticles", m_inputMCParticles, "mcparticle truth data from npsim");
0046 declareProperty("outputTrackParameters", m_outputTrackParameters, "Output initial track parameters");
0047 }
0048
0049 StatusCode initialize() override {
0050 if (Gaudi::Algorithm::initialize().isFailure()) {
0051 return StatusCode::FAILURE;
0052 }
0053 IRndmGenSvc* randSvc = Gaudi::svcLocator()->service<IRndmGenSvc>("RndmGenSvc", true);
0054 if (randSvc == nullptr) {
0055 return StatusCode::FAILURE;
0056 }
0057 m_pidSvc = service("ParticleSvc");
0058 if (!m_pidSvc) {
0059 error() << "Unable to locate Particle Service. "
0060 << "Make sure you have ParticleSvc in the configuration."
0061 << endmsg;
0062 return StatusCode::FAILURE;
0063 }
0064 return StatusCode::SUCCESS;
0065 }
0066
0067 StatusCode execute(const EventContext&) const override {
0068
0069 const auto* const mcparts = m_inputMCParticles.get();
0070
0071 auto* init_trk_params = m_outputTrackParameters.createAndPut();
0072
0073 for(const auto& part : *mcparts) {
0074
0075
0076 if(part.getGeneratorStatus() != 1 ) {
0077 continue;
0078 }
0079
0080 const auto& pvec = part.getMomentum();
0081 const auto p = std::hypot(pvec.x, pvec.y, pvec.z);
0082 const auto phi = std::atan2(pvec.x, pvec.y);
0083 const auto theta = std::atan2(std::hypot(pvec.x, pvec.y), pvec.z);
0084
0085
0086
0087
0088 const auto charge = static_cast<float>(m_pidSvc->particle(part.getPDG()).charge);
0089 if (abs(charge) < std::numeric_limits<double>::epsilon()) {
0090 continue;
0091 }
0092
0093 const auto q_over_p = charge / p;
0094
0095 auto params = init_trk_params->create();
0096 params.setType(-1);
0097 params.setLoc({0.0F, 0.0F});
0098 params.setTheta(theta);
0099 params.setPhi(phi);
0100 params.setQOverP(q_over_p * .05F);
0101 params.setTime(part.getTime());
0102 #if EDM4EIC_VERSION_MAJOR >= 5
0103 edm4eic::Cov6f cov;
0104 cov(0,0) = 0.1;
0105 cov(0,1) = 0.1;
0106 cov(1,1) = 0.1;
0107 cov(2,2) = 0.1;
0108 cov(3,3) = 0.1;
0109 cov(4,4) = 0.1;
0110 cov(5,5) = 0.1;
0111 params.setCovariance(cov);
0112 #else
0113 params.setLocError({0.1, 0.1, 0.1});
0114 params.setMomentumError({0.1, 0.1, 0.1});
0115 params.setTimeError(0.1);
0116 params.setCharge(charge);
0117 #endif
0118
0119
0120
0121
0122
0123 if (msgLevel(MSG::DEBUG)) {
0124 debug() << "Invoke track finding seeded by truth particle with p = " << p << " GeV" << endmsg;
0125 debug() << " charge = " << charge << endmsg;
0126 debug() << " q/p = " << charge / p << endmsg;
0127 }
0128 }
0129 return StatusCode::SUCCESS;
0130 }
0131 };
0132
0133 DECLARE_COMPONENT(TruthTrackSeeding)
0134
0135 }