From 5082cb206ed7e46e9a1ada62ded9689c5608712d Mon Sep 17 00:00:00 2001
From: Sylvester Joosten <sjoosten@anl.gov>
Date: Sat, 21 Aug 2021 22:18:39 +0000
Subject: [PATCH] Added algorithm to merge cluster info (e.g. imaging and scfi
info)
---
JugBase/JugBase/UniqueID.h | 3 +-
JugReco/CMakeLists.txt | 1 +
JugReco/src/components/CalorimeterHitReco.cpp | 4 +-
.../CalorimeterHitsEtaPhiProjector.cpp | 4 +-
.../src/components/CalorimeterHitsMerger.cpp | 4 +-
JugReco/src/components/ClusterRecoCoG.cpp | 4 +-
.../EnergyPositionClusterMerger.cpp | 130 ++++++++++++++++++
JugReco/src/components/ImagingClusterReco.cpp | 4 +-
JugReco/src/components/ImagingPixelMerger.cpp | 4 +-
JugReco/src/components/ImagingPixelReco.cpp | 4 +-
.../src/components/PhotoMultiplierReco.cpp | 4 +-
JugReco/src/components/SimpleClustering.cpp | 4 +-
.../components/TrackerHitReconstruction.cpp | 4 +-
13 files changed, 153 insertions(+), 21 deletions(-)
create mode 100644 JugReco/src/components/EnergyPositionClusterMerger.cpp
diff --git a/JugBase/JugBase/UniqueID.h b/JugBase/JugBase/UniqueID.h
index 3c9989d2..977249ff 100644
--- a/JugBase/JugBase/UniqueID.h
+++ b/JugBase/JugBase/UniqueID.h
@@ -6,6 +6,7 @@
//
// Both in function and mixin form (latter is useful for algorithms)
+#include <cstdint>
#include <functional>
#include <limits>
#include <string>
@@ -20,7 +21,7 @@ template <class Integer> Integer uniqueID(const std::string& s) {
return static_cast<Integer>(fullID & max);
}
-template <class Integer> class AlgorithmIDMixin {
+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";
diff --git a/JugReco/CMakeLists.txt b/JugReco/CMakeLists.txt
index 68acc233..b90538d8 100644
--- a/JugReco/CMakeLists.txt
+++ b/JugReco/CMakeLists.txt
@@ -20,6 +20,7 @@ gaudi_add_module(JugRecoPlugins
src/components/CrystalEndcapsReco.cpp
src/components/EcalTungstenSamplingReco.cpp
src/components/EcalTungstenSamplingCluster.cpp
+ src/components/EnergyPositionClusterMerger.cpp
src/components/ClusterRecoCoG.cpp
src/components/PhotoMultiplierReco.cpp
src/components/PhotoRingClusters.cpp
diff --git a/JugReco/src/components/CalorimeterHitReco.cpp b/JugReco/src/components/CalorimeterHitReco.cpp
index 31ebdc89..00da050a 100644
--- a/JugReco/src/components/CalorimeterHitReco.cpp
+++ b/JugReco/src/components/CalorimeterHitReco.cpp
@@ -37,7 +37,7 @@ namespace Jug::Reco {
* Reconstruct digitized outputs, paired with Jug::Digi::CalorimeterHitDigi
* \ingroup reco
*/
- class CalorimeterHitReco : public GaudiAlgorithm, AlgorithmIDMixin<int32_t> {
+ class CalorimeterHitReco : public GaudiAlgorithm, AlgorithmIDMixin<> {
public:
// length unit from dd4hep, should be fixed
@@ -81,7 +81,7 @@ namespace Jug::Reco {
CalorimeterHitReco(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
- , AlgorithmIDMixin(name, info())
+ , AlgorithmIDMixin<>(name, info())
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
diff --git a/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp b/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp
index cbfabc37..953e209b 100644
--- a/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp
+++ b/JugReco/src/components/CalorimeterHitsEtaPhiProjector.cpp
@@ -57,7 +57,7 @@ namespace Jug::Reco {
*
* \ingroup reco
*/
- class CalorimeterHitsEtaPhiProjector : public GaudiAlgorithm, AlgorithmIDMixin<int32_t> {
+ class CalorimeterHitsEtaPhiProjector : public GaudiAlgorithm, AlgorithmIDMixin<> {
public:
Gaudi::Property<std::vector<double>> u_gridSizes{this, "gridSizes", {0.001, 0.001*rad}};
DataHandle<eic::CalorimeterHitCollection> m_inputHitCollection{
@@ -69,7 +69,7 @@ namespace Jug::Reco {
CalorimeterHitsEtaPhiProjector(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
- , AlgorithmIDMixin(name, info())
+ , AlgorithmIDMixin<>(name, info())
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
diff --git a/JugReco/src/components/CalorimeterHitsMerger.cpp b/JugReco/src/components/CalorimeterHitsMerger.cpp
index 156c7206..872a7d9f 100644
--- a/JugReco/src/components/CalorimeterHitsMerger.cpp
+++ b/JugReco/src/components/CalorimeterHitsMerger.cpp
@@ -43,7 +43,7 @@ namespace Jug::Reco {
*
* \ingroup reco
*/
- class CalorimeterHitsMerger : public GaudiAlgorithm, AlgorithmIDMixin<int32_t> {
+ class CalorimeterHitsMerger : public GaudiAlgorithm, AlgorithmIDMixin<> {
public:
Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
Gaudi::Property<std::string> m_readout{this, "readoutClass", ""};
@@ -61,7 +61,7 @@ namespace Jug::Reco {
CalorimeterHitsMerger(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
- , AlgorithmIDMixin(name, info())
+ , AlgorithmIDMixin<>(name, info())
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
diff --git a/JugReco/src/components/ClusterRecoCoG.cpp b/JugReco/src/components/ClusterRecoCoG.cpp
index 3cb873cb..0b9d58d6 100644
--- a/JugReco/src/components/ClusterRecoCoG.cpp
+++ b/JugReco/src/components/ClusterRecoCoG.cpp
@@ -63,7 +63,7 @@ namespace Jug::Reco {
*
* \ingroup reco
*/
- class ClusterRecoCoG : public GaudiAlgorithm, AlgorithmIDMixin<int32_t> {
+ class ClusterRecoCoG : public GaudiAlgorithm, AlgorithmIDMixin<> {
public:
Gaudi::Property<double> m_sampFrac{this, "samplingFraction", 1.0};
Gaudi::Property<double> m_logWeightBase{this, "logWeightBase", 3.6};
@@ -85,7 +85,7 @@ namespace Jug::Reco {
ClusterRecoCoG(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
- , AlgorithmIDMixin(name, info())
+ , AlgorithmIDMixin<>(name, info())
{
declareProperty("inputHitCollection", m_inputHits, "");
declareProperty("inputProtoClusterCollection", m_inputProto, "");
diff --git a/JugReco/src/components/EnergyPositionClusterMerger.cpp b/JugReco/src/components/EnergyPositionClusterMerger.cpp
new file mode 100644
index 00000000..5f7af366
--- /dev/null
+++ b/JugReco/src/components/EnergyPositionClusterMerger.cpp
@@ -0,0 +1,130 @@
+#include <limits>
+// Gaudi
+#include "Gaudi/Property.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiAlg/GaudiTool.h"
+#include "GaudiAlg/Transformer.h"
+#include "GaudiKernel/PhysicalConstants.h"
+
+#include "JugBase/DataHandle.h"
+#include "JugBase/UniqueID.h"
+
+// Event Model related classes
+#include "eicd/ClusterCollection.h"
+#include "eicd/MergedClusterRelationsCollection.h"
+
+using namespace Gaudi::Units;
+
+namespace Jug::Reco {
+
+/** Simple algorithm to merge the energy measurement from cluster1 with the position
+ * measurement of cluster2 (in case matching clusters are found). If not, it will
+ * propagate the raw cluster from cluster1 or cluster2
+ *
+ * Matching occurs based on the cluster phi, z and E variables, with tolerances
+ * defined in the options file. A negative tolerance effectively disables
+ * a check. The energy tolerance is defined as a relative number (e.g. .1)
+ *
+ * In case of ambiguity the closest cluster is merged.
+ *
+ * \ingroup reco
+ */
+class EnergyPositionClusterMerger : public GaudiAlgorithm, public AlgorithmIDMixin<> {
+public:
+ // Input
+ DataHandle<eic::ClusterCollection> m_energyClusters{"energyClusters", Gaudi::DataHandle::Reader, this};
+ DataHandle<eic::ClusterCollection> m_positionClusters{"positionClusters", Gaudi::DataHandle::Reader, this};
+ // Output
+ DataHandle<eic::ClusterCollection> m_outputClusters{"outputClusters", Gaudi::DataHandle::Writer, this};
+ DataHandle<eic::MergedClusterRelationsCollection> m_relations{"outputCluster", Gaudi::DataHandle::Writer,
+ this}; // namespace Jug::Reco
+ // Negative values mean the tolerance check is disabled
+ Gaudi::Property<double> m_zToleranceUnits{this, "zTolerance", -1 * cm};
+ Gaudi::Property<double> m_phiToleranceUnits{this, "phiTolerance", 20 * degree};
+ Gaudi::Property<double> m_energyRelTolerance{this, "energyRelTolerance", 0.3};
+ // Unitless (GeV/mm/ns/rad) versions of these tolerances
+ double m_zTolerance;
+ double m_phiTolerance;
+
+public:
+ EnergyPositionClusterMerger(const std::string& name, ISvcLocator* svcLoc)
+ : GaudiAlgorithm(name, svcLoc), AlgorithmIDMixin(name, info()) {
+ declareProperty("energyClusters", m_energyClusters, "Cluster collection with good energy precision");
+ declareProperty("positionClusters", m_positionClusters, "Cluster collection with good position precision");
+ }
+
+ StatusCode initialize() override {
+ m_zTolerance = m_zToleranceUnits / mm;
+ m_phiTolerance = m_phiTolerance / rad;
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode execute() override {
+ // input
+ const auto& e_clus = *(m_energyClusters.get());
+ const auto& pos_clus = *(m_positionClusters.get());
+ // output
+ auto& merged = *(m_outputClusters.createAndPut());
+ auto& relations = *(m_relations.createAndPut());
+
+ std::vector<bool> consumed(e_clus.size(), false);
+
+ int idx = 0; // merged cluster ID
+ // use position clusters as starting point
+ for (const auto& pc : pos_clus) {
+ // check if we find a good match
+ int best_match = -1;
+ double best_delta = std::numeric_limits<double>::max();
+ for (size_t ie = 0; ie < e_clus.size(); ++ie) {
+ if (consumed[ie]) {
+ continue;
+ }
+ const auto& ec = e_clus[ie];
+ // 1. stop if not within tolerance
+ 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)) {
+ continue;
+ }
+ // --> if we get here, we have a match within tolerance. Now treat the case
+ // where we have multiple matches. In this case take the one with the closest
+ // energies.
+ // 2. best match?
+ const double delta = fabs(pc.energy() - ec.energy());
+ if (delta < best_delta) {
+ best_delta = delta;
+ best_match = ie;
+ }
+ }
+ // Create a merged cluster if we find a good match
+ if (best_match >= 0) {
+ const auto& ec = e_clus[best_match];
+ auto new_clus = merged.create();
+ new_clus.ID(idx++);
+ new_clus.source(algorithmID());
+ new_clus.energy(ec.energy());
+ new_clus.energyError(ec.energyError());
+ new_clus.time(pc.time());
+ new_clus.nhits(pc.nhits() + ec.nhits());
+ new_clus.position(pc.position());
+ new_clus.positionError(pc.positionError());
+ new_clus.radius(pc.radius());
+ new_clus.skewness(pc.skewness());
+ auto rel = relations.create();
+ rel.clusterID(new_clus.ID());
+ rel.size(2);
+ rel.parent()[0] = {pc.source(), pc.ID()};
+ rel.parent()[1] = {ec.source(), ec.ID()};
+ // label our energy cluster as consumed
+ consumed[best_match] = true;
+ }
+ }
+ // That's all!
+
+ return StatusCode::SUCCESS;
+ }
+};
+DECLARE_COMPONENT(EnergyPositionClusterMerger)
+
+} // namespace Jug::Reco
diff --git a/JugReco/src/components/ImagingClusterReco.cpp b/JugReco/src/components/ImagingClusterReco.cpp
index 29c7c255..66611c8b 100644
--- a/JugReco/src/components/ImagingClusterReco.cpp
+++ b/JugReco/src/components/ImagingClusterReco.cpp
@@ -46,7 +46,7 @@ namespace Jug::Reco {
*
* \ingroup reco
*/
-class ImagingClusterReco : public GaudiAlgorithm, AlgorithmIDMixin<int32_t> {
+class ImagingClusterReco : public GaudiAlgorithm, AlgorithmIDMixin<> {
public:
Gaudi::Property<double> m_sampFrac{this, "samplingFraction", 1.0};
Gaudi::Property<int> m_trackStopLayer{this, "trackStopLayer", 9};
@@ -63,7 +63,7 @@ public:
ImagingClusterReco(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
- , AlgorithmIDMixin(name, info()) {
+ , AlgorithmIDMixin<>(name, info()) {
declareProperty("inputProtoClusterCollection", m_inputProtoClusterCollection, "");
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputLayerCollection", m_outputLayerCollection, "");
diff --git a/JugReco/src/components/ImagingPixelMerger.cpp b/JugReco/src/components/ImagingPixelMerger.cpp
index da8cea79..102c2237 100644
--- a/JugReco/src/components/ImagingPixelMerger.cpp
+++ b/JugReco/src/components/ImagingPixelMerger.cpp
@@ -56,7 +56,7 @@ namespace Jug::Reco {
*
* \ingroup reco
*/
- class ImagingPixelMerger : public GaudiAlgorithm, AlgorithmIDMixin<int32_t> {
+ class ImagingPixelMerger : public GaudiAlgorithm, AlgorithmIDMixin<> {
public:
Gaudi::Property<int> m_nHits{this, "numberOfHits", 20};
Gaudi::Property<int> m_nLayers{this, "numberOfLayers", 20};
@@ -69,7 +69,7 @@ namespace Jug::Reco {
ImagingPixelMerger(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
- , AlgorithmIDMixin(name, info())
+ , AlgorithmIDMixin<>(name, info())
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
diff --git a/JugReco/src/components/ImagingPixelReco.cpp b/JugReco/src/components/ImagingPixelReco.cpp
index e3945f39..74ebe43b 100644
--- a/JugReco/src/components/ImagingPixelReco.cpp
+++ b/JugReco/src/components/ImagingPixelReco.cpp
@@ -36,7 +36,7 @@ namespace Jug::Reco {
*
* \ingroup reco
*/
- class ImagingPixelReco : public GaudiAlgorithm, AlgorithmIDMixin<int32_t> {
+ class ImagingPixelReco : public GaudiAlgorithm, AlgorithmIDMixin<> {
public:
// geometry service
Gaudi::Property<std::string> m_geoSvcName{this, "geoServiceName", "GeoSvc"};
@@ -71,7 +71,7 @@ namespace Jug::Reco {
ImagingPixelReco(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
- , AlgorithmIDMixin(name, info())
+ , AlgorithmIDMixin<>(name, info())
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
diff --git a/JugReco/src/components/PhotoMultiplierReco.cpp b/JugReco/src/components/PhotoMultiplierReco.cpp
index eb4b44db..3183b4f1 100644
--- a/JugReco/src/components/PhotoMultiplierReco.cpp
+++ b/JugReco/src/components/PhotoMultiplierReco.cpp
@@ -40,7 +40,7 @@ namespace Jug::Reco {
*
* \ingroup reco
*/
- class PhotoMultiplierReco : public GaudiAlgorithm, AlgorithmIDMixin<int32_t> {
+ class PhotoMultiplierReco : public GaudiAlgorithm, AlgorithmIDMixin<> {
public:
DataHandle<eic::RawPMTHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::PMTHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
@@ -54,7 +54,7 @@ namespace Jug::Reco {
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
PhotoMultiplierReco(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
- , AlgorithmIDMixin(name, info())
+ , AlgorithmIDMixin<>(name, info())
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
diff --git a/JugReco/src/components/SimpleClustering.cpp b/JugReco/src/components/SimpleClustering.cpp
index f758dd07..a86635a2 100644
--- a/JugReco/src/components/SimpleClustering.cpp
+++ b/JugReco/src/components/SimpleClustering.cpp
@@ -30,7 +30,7 @@ namespace Jug::Reco {
*
* \ingroup reco
*/
- class SimpleClustering : public GaudiAlgorithm, AlgorithmIDMixin<int32_t> {
+ class SimpleClustering : public GaudiAlgorithm, AlgorithmIDMixin<> {
public:
using RecHits = eic::CalorimeterHitCollection;
using ProtoClusters = eic::ProtoClusterCollection;
@@ -48,7 +48,7 @@ namespace Jug::Reco {
SimpleClustering(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
- , AlgorithmIDMixin(name, info()) {
+ , AlgorithmIDMixin<>(name, info()) {
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputProtoClusterCollection", m_outputClusters, "Output proto clusters");
declareProperty("outputClusterCollection", m_outputClusters, "Output clusters");
diff --git a/JugReco/src/components/TrackerHitReconstruction.cpp b/JugReco/src/components/TrackerHitReconstruction.cpp
index 739c4ad2..d497b491 100644
--- a/JugReco/src/components/TrackerHitReconstruction.cpp
+++ b/JugReco/src/components/TrackerHitReconstruction.cpp
@@ -31,7 +31,7 @@ namespace Jug {
*
* \ingroup reco
*/
- class TrackerHitReconstruction : public GaudiAlgorithm, AlgorithmIDMixin<int32_t> {
+ class TrackerHitReconstruction : public GaudiAlgorithm, AlgorithmIDMixin<> {
public:
Gaudi::Property<double> m_timeResolution{this, "timeResolution", 10};
Rndm::Numbers m_gaussDist;
@@ -46,7 +46,7 @@ namespace Jug {
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
TrackerHitReconstruction(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
- , AlgorithmIDMixin(name, info())
+ , AlgorithmIDMixin<>(name, info())
{
declareProperty("inputHitCollection", m_inputHitCollection, "");
declareProperty("outputHitCollection", m_outputHitCollection, "");
--
GitLab