From 2b5af88c3da0fb6531f7462eb2e4c4c3e01caeb6 Mon Sep 17 00:00:00 2001
From: Sylvester Joosten <sjoosten@anl.gov>
Date: Sun, 22 Aug 2021 01:08:00 +0000
Subject: [PATCH] Added simple class to associate truth PID with tracking info
---
JugReco/CMakeLists.txt | 1 +
.../src/components/ParticleWithTruthPID.cpp | 92 +++++++++++++++++++
2 files changed, 93 insertions(+)
create mode 100644 JugReco/src/components/ParticleWithTruthPID.cpp
diff --git a/JugReco/CMakeLists.txt b/JugReco/CMakeLists.txt
index b90538d8..4b0b0c7b 100644
--- a/JugReco/CMakeLists.txt
+++ b/JugReco/CMakeLists.txt
@@ -22,6 +22,7 @@ gaudi_add_module(JugRecoPlugins
src/components/EcalTungstenSamplingCluster.cpp
src/components/EnergyPositionClusterMerger.cpp
src/components/ClusterRecoCoG.cpp
+ src/components/ParticleWithTruthPID.cpp
src/components/PhotoMultiplierReco.cpp
src/components/PhotoRingClusters.cpp
src/components/FuzzyKClusters.cpp
diff --git a/JugReco/src/components/ParticleWithTruthPID.cpp b/JugReco/src/components/ParticleWithTruthPID.cpp
new file mode 100644
index 00000000..dccd4453
--- /dev/null
+++ b/JugReco/src/components/ParticleWithTruthPID.cpp
@@ -0,0 +1,92 @@
+#include "Gaudi/Algorithm.h"
+#include "GaudiAlg/GaudiTool.h"
+#include "GaudiAlg/Producer.h"
+#include "GaudiAlg/Transformer.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include <algorithm>
+#include <cmath>
+
+#include "JugBase/DataHandle.h"
+#include "JugBase/UniqueID.h"
+
+// Event Model related classes
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "eicd/ReconstructedParticleCollection.h"
+#include "eicd/TrackParametersCollection.h"
+#include "eicd/VectorPolar.h"
+
+namespace Jug::Rec {
+
+class ParticleWithTruthPID : public GaudiAlgorithm, AlgorithmIDMixin<> {
+public:
+ DataHandle<dd4pod::Geant4ParticleCollection> m_inputTruthCollection{"inputMC", Gaudi::DataHandle::Reader, this};
+ DataHandle<eic::TrackParametersCollection> m_inputTrackCollection{"inputTracks", Gaudi::DataHandle::Reader, this};
+ DataHandle<eic::ReconstructedParticleCollection> m_outputParticleCollection{"outputParticles",
+ Gaudi::DataHandle::Writer, this};
+ ParticleWithTruthPID(const std::string& name, ISvcLocator* svcLoc)
+ : GaudiAlgorithm(name, svcLoc), AlgorithmIDMixin(name, info()) {
+ declareProperty("inputMC", m_inputTruthCollection, "mcparticles");
+ declareProperty("inputTracks", m_inputTrackCollection, "outputTrackParameters");
+ declareProperty("outputParticles", m_outputParticleCollection, "ReconstructedParticles");
+ }
+ StatusCode initialize() override {
+ if (GaudiAlgorithm::initialize().isFailure())
+ return StatusCode::FAILURE;
+
+ return StatusCode::SUCCESS;
+ }
+ StatusCode execute() override {
+ // input collection
+ const auto& mc = *(m_inputTruthCollection.get());
+ const auto& tracks = *(m_inputTrackCollection.get());
+ auto& part = *(m_outputParticleCollection.createAndPut());
+
+ std::vector<bool> consumed(tracks.size(), false);
+ int ID = 0;
+ for (const auto& p : mc) {
+ if (p.genStatus() != 1) {
+ continue;
+ }
+ // check if we find a good match
+ int best_match = -1;
+ double best_delta = std::numeric_limits<double>::max();
+ for (int it = 0; it < tracks.size(); ++it) {
+ if (consumed[it]) {
+ continue;
+ }
+ const auto& trk = tracks[it];
+ eic::VectorPolar mom{1.0 / std::abs(trk.qOverP()), trk.direction().theta, trk.direction().phi};
+ double delta = std::hypot(p.ps().x - mom.x(), p.ps().y - mom.y(), p.ps().z - mom.z());
+ if (delta < best_delta) {
+ best_match = it;
+ best_delta = delta;
+ }
+ }
+ if (best_match > 0) {
+ consumed[best_match] = true;
+ const auto& trk = tracks[best_match];
+ eic::VectorPolar mom{1.0 / std::abs(trk.qOverP()), trk.direction().theta, trk.direction().phi};
+ eic::ReconstructedParticle rec_part{ID++, // index
+ {mom.x(), mom.y(), mom.z()}, // momentum
+ {0., 0., 0.}, // vertex
+ 0., // time
+ p.pdgID(), // PID
+ static_cast<int16_t>(0), // Status
+ static_cast<int16_t>(std::copysign(1., trk.qOverP())), // charge
+ algorithmID(), // Algorithm type
+ 1., // particle weight
+ mom.mag(), // 3-momentum magnitude [GeV]
+ std::hypot(mom.mag(), p.mass()), // energy [GeV]
+ static_cast<float>(p.mass())}; // mass [GeV]
+ part.push_back(rec_part);
+ }
+ }
+
+ return StatusCode::SUCCESS;
+ }
+};
+
+DECLARE_COMPONENT(ParticleWithTruthPID)
+
+} // namespace Jug::Rec
+
--
GitLab