From f486383374d8d96b5956a145536697dbf2abae3d Mon Sep 17 00:00:00 2001
From: Sylvester Joosten <sjoosten@anl.gov>
Date: Sun, 22 Aug 2021 20:30:16 +0000
Subject: [PATCH] Added debug info and fixed some minor inconsistencies
---
JugBase/JugBase/UniqueID.h | 2 +-
JugReco/CMakeLists.txt | 2 +-
.../src/components/CalorimeterHitsMerger.cpp | 4 ++-
.../components/CalorimeterIslandCluster.cpp | 21 +++++++-----
JugReco/src/components/ClusterRecoCoG.cpp | 14 +++++---
.../EnergyPositionClusterMerger.cpp | 22 +++++++++++-
JugReco/src/components/ImagingClusterReco.cpp | 13 ++++---
JugReco/src/components/ImagingTopoCluster.cpp | 23 ++++++-------
...TruthPID.cpp => ParticlesWithTruthPID.cpp} | 34 +++++++++++++------
JugReco/src/components/SimpleClustering.cpp | 15 +++++---
.../src/components/TrackingHitsCollector2.cpp | 8 +++--
11 files changed, 108 insertions(+), 50 deletions(-)
rename JugReco/src/components/{ParticleWithTruthPID.cpp => ParticlesWithTruthPID.cpp} (74%)
diff --git a/JugBase/JugBase/UniqueID.h b/JugBase/JugBase/UniqueID.h
index 977249ff..e21e4c9e 100644
--- a/JugBase/JugBase/UniqueID.h
+++ b/JugBase/JugBase/UniqueID.h
@@ -24,7 +24,7 @@ template <class Integer> Integer uniqueID(const std::string& s) {
template <class Integer = int32_t> class AlgorithmIDMixin {
public:
AlgorithmIDMixin(const std::string& name, MsgStream& out) : m_id{uniqueID<Integer>(name)} {
- out << "Unique ID associated with '" << name << "': " << m_id << "\n";
+ out << "Unique ID associated with '" << name << "': " << m_id << endmsg;
}
AlgorithmIDMixin() = delete;
AlgorithmIDMixin(const AlgorithmIDMixin&) = delete;
diff --git a/JugReco/CMakeLists.txt b/JugReco/CMakeLists.txt
index 4b0b0c7b..068ea53e 100644
--- a/JugReco/CMakeLists.txt
+++ b/JugReco/CMakeLists.txt
@@ -22,7 +22,7 @@ gaudi_add_module(JugRecoPlugins
src/components/EcalTungstenSamplingCluster.cpp
src/components/EnergyPositionClusterMerger.cpp
src/components/ClusterRecoCoG.cpp
- src/components/ParticleWithTruthPID.cpp
+ src/components/ParticlesWithTruthPID.cpp
src/components/PhotoMultiplierReco.cpp
src/components/PhotoRingClusters.cpp
src/components/FuzzyKClusters.cpp
diff --git a/JugReco/src/components/CalorimeterHitsMerger.cpp b/JugReco/src/components/CalorimeterHitsMerger.cpp
index 872a7d9f..174ce91c 100644
--- a/JugReco/src/components/CalorimeterHitsMerger.cpp
+++ b/JugReco/src/components/CalorimeterHitsMerger.cpp
@@ -166,7 +166,9 @@ namespace Jug::Reco {
href.dimension()});
}
- debug() << "Size before = " << inputs.size() << ", after = " << outputs.size() << endmsg;
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << "Size before = " << inputs.size() << ", after = " << outputs.size() << endmsg;
+ }
return StatusCode::SUCCESS;
}
diff --git a/JugReco/src/components/CalorimeterIslandCluster.cpp b/JugReco/src/components/CalorimeterIslandCluster.cpp
index 93e6a569..59fe127f 100644
--- a/JugReco/src/components/CalorimeterIslandCluster.cpp
+++ b/JugReco/src/components/CalorimeterIslandCluster.cpp
@@ -182,11 +182,14 @@ namespace Jug::Reco {
std::vector<bool> visits(hits.size(), false);
for (size_t i = 0; i < hits.size(); ++i) {
- debug() << fmt::format("hit {:d}: energy = {:.4f} MeV, local = ({:.4f}, {:.4f}) mm, "
- "global=({:.4f}, {:.4f}, {:.4f}) mm, layer = {:d}, sector = {:d}.",
- i, hits[i].energy()*1000., hits[i].local().x, hits[i].local().y,
- hits[i].position().x, hits[i].position().y, hits[i].position().z,
- hits[i].layer(), hits[i].sector()) << endmsg;
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << fmt::format("hit {:d}: energy = {:.4f} MeV, local = ({:.4f}, {:.4f}) mm, "
+ "global=({:.4f}, {:.4f}, {:.4f}) mm, layer = {:d}, sector = {:d}.",
+ i, hits[i].energy() * 1000., hits[i].local().x, hits[i].local().y,
+ hits[i].position().x, hits[i].position().y, hits[i].position().z, hits[i].layer(),
+ hits[i].sector())
+ << endmsg;
+ }
// already in a group
if (visits[i]) {
continue;
@@ -203,9 +206,11 @@ namespace Jug::Reco {
}
auto maxima = find_maxima(group, !m_splitCluster.value());
split_group(group, maxima, clusterID, proto);
- debug() << "hits in a group: " << group.size() << ", "
- << "local maxima: " << maxima.size() << endmsg;
- debug() << "total number of clusters so far: " << clusterID << ", " << endmsg;
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << "hits in a group: " << group.size() << ", "
+ << "local maxima: " << maxima.size() << endmsg;
+ debug() << "total number of clusters so far: " << clusterID << ", " << endmsg;
+ }
}
return StatusCode::SUCCESS;
diff --git a/JugReco/src/components/ClusterRecoCoG.cpp b/JugReco/src/components/ClusterRecoCoG.cpp
index 0b9d58d6..2edc71c9 100644
--- a/JugReco/src/components/ClusterRecoCoG.cpp
+++ b/JugReco/src/components/ClusterRecoCoG.cpp
@@ -159,8 +159,10 @@ namespace Jug::Reco {
for (const auto& [idx, hit_info] : cluster_map) {
auto cl = reconstruct(hit_info, idx);
- debug() << cl.nhits() << " hits: " << cl.energy() / GeV << " GeV, (" << cl.position().x / mm
- << ", " << cl.position().y / mm << ", " << cl.position().z / mm << ")" << endmsg;
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << cl.nhits() << " hits: " << cl.energy() / GeV << " GeV, (" << cl.position().x / mm << ", "
+ << cl.position().y / mm << ", " << cl.position().z / mm << ")" << endmsg;
+ }
clusters.push_back(cl);
info.push_back({cl.ID(), cl.position(), cl.position().eta()});
}
@@ -186,7 +188,9 @@ namespace Jug::Reco {
cl.nhits(hit_info.size());
// no hits
- debug() << "hit size = " << hit_info.size() << endmsg;
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << "hit size = " << hit_info.size() << endmsg;
+ }
if (hit_info.empty()) {
return cl;
}
@@ -196,7 +200,9 @@ namespace Jug::Reco {
float maxE = 0.;
auto time = hit_info[0].second.time();
for (const auto& [proto, hit] : hit_info) {
- debug() << "hit energy = " << hit.energy() << " hit weight: " << proto.weight() << endmsg;
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << "hit energy = " << hit.energy() << " hit weight: " << proto.weight() << endmsg;
+ }
auto energy = hit.energy() * proto.weight();
totalE += energy;
if (energy > maxE) {
diff --git a/JugReco/src/components/EnergyPositionClusterMerger.cpp b/JugReco/src/components/EnergyPositionClusterMerger.cpp
index 5f7af366..39aa263d 100644
--- a/JugReco/src/components/EnergyPositionClusterMerger.cpp
+++ b/JugReco/src/components/EnergyPositionClusterMerger.cpp
@@ -1,4 +1,7 @@
#include <limits>
+#include <numbers>
+
+#include <fmt/format.h>
// Gaudi
#include "Gaudi/Property.h"
#include "GaudiAlg/GaudiAlgorithm.h"
@@ -81,10 +84,15 @@ public:
}
const auto& ec = e_clus[ie];
// 1. stop if not within tolerance
+ // (make sure to handle rollover of phi properly)
+ double dphi = pc.position().phi() - ec.position().phi();
+ if (std::abs(dphi) > M_PI) {
+ dphi = std::abs(dphi) - M_PI;
+ }
if ((m_energyRelTolerance > 0 &&
(ec.energy() <= 0 || fabs((pc.energy() - ec.energy()) / ec.energy()) > m_energyRelTolerance)) ||
(m_zTolerance > 0 && fabs(pc.position().z - ec.position().z) > m_zTolerance) ||
- (m_phiTolerance > 0 && fabs(pc.position().phi() - ec.position().phi()) > m_phiTolerance)) {
+ (m_phiTolerance > 0 && dphi > m_phiTolerance)) {
continue;
}
// --> if we get here, we have a match within tolerance. Now treat the case
@@ -118,6 +126,18 @@ public:
rel.parent()[1] = {ec.source(), ec.ID()};
// label our energy cluster as consumed
consumed[best_match] = true;
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << fmt::format("Matched position cluster {} with energy cluster {}\n", pc.ID(), ec.ID()) << endmsg;
+ debug() << fmt::format(" - Position cluster: (E: {}, phi: {}, z: {})", pc.energy(), pc.position().phi(),
+ pc.position().z)
+ << endmsg;
+ debug() << fmt::format(" - Energy cluster: (E: {}, phi: {}, z: {})", ec.energy(), ec.position().phi(),
+ ec.position().z)
+ << endmsg;
+ debug() << fmt::format(" ---> Merged cluster: (E: {}, phi: {}, z: {})", new_clus.energy(),
+ new_clus.position().phi(), new_clus.position().z)
+ << endmsg;
+ }
}
}
// That's all!
diff --git a/JugReco/src/components/ImagingClusterReco.cpp b/JugReco/src/components/ImagingClusterReco.cpp
index 66611c8b..85abd8d0 100644
--- a/JugReco/src/components/ImagingClusterReco.cpp
+++ b/JugReco/src/components/ImagingClusterReco.cpp
@@ -125,11 +125,14 @@ public:
// debug output
int idx = 0;
- for (const auto& cl : clusters) {
- debug() << fmt::format("Cluster {:d}: Edep = {:.3f} MeV, Dir = ({:.3f}, {:.3f}) deg", cl.ID(), cl.energy() * 1000.,
- info[idx].direction().theta / M_PI * 180., info[idx].direction().phi / M_PI * 180.)
- << endmsg;
- idx += 1;
+ if (msgLevel(MSG::DEBUG)) {
+ for (const auto& cl : clusters) {
+ debug() << fmt::format("Cluster {:d}: Edep = {:.3f} MeV, Dir = ({:.3f}, {:.3f}) deg", cl.ID(),
+ cl.energy() * 1000., info[idx].direction().theta / M_PI * 180.,
+ info[idx].direction().phi / M_PI * 180.)
+ << endmsg;
+ idx += 1;
+ }
}
return StatusCode::SUCCESS;
diff --git a/JugReco/src/components/ImagingTopoCluster.cpp b/JugReco/src/components/ImagingTopoCluster.cpp
index 4c70715b..2943a8e4 100644
--- a/JugReco/src/components/ImagingTopoCluster.cpp
+++ b/JugReco/src/components/ImagingTopoCluster.cpp
@@ -132,13 +132,12 @@ namespace Jug::Reco {
std::vector<bool> visits(hits.size(), false);
std::vector<std::vector<eic::ConstCalorimeterHit>> groups;
for (size_t i = 0; i < hits.size(); ++i) {
- /*
- debug() << fmt::format("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {})",
- i + 1,
- hits[i].local().x, hits[i].local().y, hits[i].z(),
- hits[i].position().x, hits[i].position().y, hits[i].position().z)
- << endmsg;
- */
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << fmt::format("hit {:d}: local position = ({}, {}, {}), global position = ({}, {}, {})", i + 1,
+ hits[i].local().x, hits[i].local().y, hits[i].position().z, hits[i].position().x,
+ hits[i].position().y, hits[i].position().z)
+ << endmsg;
+ }
// already in a group, or not energetic enough to form a cluster
if (visits[i] || hits[i].energy() < minClusterCenterEdep) {
continue;
@@ -147,12 +146,12 @@ namespace Jug::Reco {
// create a new group, and group all the neighboring hits
dfs_group(groups.back(), i, hits, visits);
}
- debug() << "found " << groups.size() << " potential clusters (groups of hits)" << endmsg;
- /*
- for (size_t i = 0; i < groups.size(); ++i) {
- debug() << fmt::format("group {}: {} hits", i, groups[i].size()) << endmsg;
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << "found " << groups.size() << " potential clusters (groups of hits)" << endmsg;
+ for (size_t i = 0; i < groups.size(); ++i) {
+ debug() << fmt::format("group {}: {} hits", i, groups[i].size()) << endmsg;
+ }
}
- */
// form clusters
size_t clusterID = 0;
diff --git a/JugReco/src/components/ParticleWithTruthPID.cpp b/JugReco/src/components/ParticlesWithTruthPID.cpp
similarity index 74%
rename from JugReco/src/components/ParticleWithTruthPID.cpp
rename to JugReco/src/components/ParticlesWithTruthPID.cpp
index dccd4453..8a6d0f19 100644
--- a/JugReco/src/components/ParticleWithTruthPID.cpp
+++ b/JugReco/src/components/ParticlesWithTruthPID.cpp
@@ -1,10 +1,13 @@
+#include <algorithm>
+#include <cmath>
+
+#include <fmt/format.h>
+
#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"
@@ -17,13 +20,15 @@
namespace Jug::Rec {
-class ParticleWithTruthPID : public GaudiAlgorithm, AlgorithmIDMixin<> {
+class ParticlesWithTruthPID : 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<dd4pod::Geant4ParticleCollection> m_inputTruthCollection{"inputMCParticles", Gaudi::DataHandle::Reader,
+ this};
+ DataHandle<eic::TrackParametersCollection> m_inputTrackCollection{"inputTrackParameters", Gaudi::DataHandle::Reader,
+ this};
DataHandle<eic::ReconstructedParticleCollection> m_outputParticleCollection{"outputParticles",
Gaudi::DataHandle::Writer, this};
- ParticleWithTruthPID(const std::string& name, ISvcLocator* svcLoc)
+ ParticlesWithTruthPID(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc), AlgorithmIDMixin(name, info()) {
declareProperty("inputMC", m_inputTruthCollection, "mcparticles");
declareProperty("inputTracks", m_inputTrackCollection, "outputTrackParameters");
@@ -50,7 +55,7 @@ public:
// 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) {
+ for (size_t it = 0; it < tracks.size(); ++it) {
if (consumed[it]) {
continue;
}
@@ -75,10 +80,17 @@ public:
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]
+ mom.mag(), // 3-momentum magnitude [GeV]
+ static_cast<float>(std::hypot(mom.mag(), p.mass())), // energy [GeV]
+ static_cast<float>(p.mass())}; // mass [GeV]
part.push_back(rec_part);
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << fmt::format("Matched track {} with MC particle {}\n", trk.ID(), p.ID()) << endmsg;
+ debug() << fmt::format(" - Track: (mom: {}, theta: {}, phi: {})", mom.mag(), mom.theta, mom.phi) << endmsg;
+ debug() << fmt::format(" - MC particle: (mom: {}, theta: {}, phi: {}, type: {}", p.ps().mag(),
+ p.ps().theta(), p.ps().phi(), p.pdgID())
+ << endmsg;
+ }
}
}
@@ -86,7 +98,7 @@ public:
}
};
-DECLARE_COMPONENT(ParticleWithTruthPID)
+DECLARE_COMPONENT(ParticlesWithTruthPID)
} // namespace Jug::Rec
diff --git a/JugReco/src/components/SimpleClustering.cpp b/JugReco/src/components/SimpleClustering.cpp
index a86635a2..97844dca 100644
--- a/JugReco/src/components/SimpleClustering.cpp
+++ b/JugReco/src/components/SimpleClustering.cpp
@@ -93,7 +93,9 @@ namespace Jug::Reco {
the_hits.push_back(h);
}
- debug() << " max_dist = " << max_dist << endmsg;
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << " max_dist = " << max_dist << endmsg;
+ }
while (have_ref && ref_hit.energy() > min_energy) {
@@ -107,11 +109,14 @@ namespace Jug::Reco {
remaining_hits.push_back(h); }
}
- debug() << " cluster size " << cluster_hits.size() << endmsg;
double total_energy =
std::accumulate(std::begin(cluster_hits), std::end(cluster_hits), 0.0,
[](double t, const eic::ConstCalorimeterHit& h1) { return (t + h1.energy()); });
- debug() << " total_energy = " << total_energy << endmsg;
+
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << " total_energy = " << total_energy << endmsg;
+ debug() << " cluster size " << cluster_hits.size() << endmsg;
+ }
eic::Cluster cl;
cl.ID(clusters.size());
@@ -137,7 +142,9 @@ namespace Jug::Reco {
remaining_hits.clear();
}
}
- debug() << " clusters found: " << clusters.size() << endmsg;
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << " clusters found: " << clusters.size() << endmsg;
+ }
return StatusCode::SUCCESS;
}
diff --git a/JugReco/src/components/TrackingHitsCollector2.cpp b/JugReco/src/components/TrackingHitsCollector2.cpp
index c187b339..a218baf8 100644
--- a/JugReco/src/components/TrackingHitsCollector2.cpp
+++ b/JugReco/src/components/TrackingHitsCollector2.cpp
@@ -51,10 +51,14 @@ namespace Jug::Reco {
StatusCode execute() override
{
auto outputHits = m_trackingHits.createAndPut();
- debug() << "execute collector" << endmsg;
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << "execute collector" << endmsg;
+ }
for(const auto& hits: m_hitCollections) {
const eic::TrackerHitCollection* hitCol = hits->get();
- debug() << "col n hits: " << hitCol->size() << endmsg;
+ if (msgLevel(MSG::DEBUG)) {
+ debug() << "col n hits: " << hitCol->size() << endmsg;
+ }
for (const auto& ahit : *hitCol) {
outputHits->push_back(ahit.clone());
}
--
GitLab