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