From b18173f74b10b4dd976de42da653e0a8eaad63e6 Mon Sep 17 00:00:00 2001
From: Jihee Kim <jihee.kim@anl.gov>
Date: Sat, 3 Apr 2021 03:18:13 +0000
Subject: [PATCH] Simple Clustering
---
Examples/options/sampling_cluster2n1.py | 58 +++++++++
Examples/sampling_calo.py | 54 ++++++++
.../components/EcalTungstenSamplingDigi.cpp | 23 +++-
JugReco/CMakeLists.txt | 2 +
.../components/CalorimeterIslandCluster.cpp | 20 ++-
JugReco/src/components/ClusterRecoCoG.cpp | 12 +-
.../EcalTungstenSamplingCluster.cpp | 116 ++++++++++++++++++
.../components/EcalTungstenSamplingReco.cpp | 52 +++++---
.../src/components/SamplingECalHitsMerger.cpp | 110 +++++++++++++++++
9 files changed, 418 insertions(+), 29 deletions(-)
create mode 100644 Examples/options/sampling_cluster2n1.py
create mode 100644 Examples/sampling_calo.py
create mode 100644 JugReco/src/components/EcalTungstenSamplingCluster.cpp
create mode 100644 JugReco/src/components/SamplingECalHitsMerger.cpp
diff --git a/Examples/options/sampling_cluster2n1.py b/Examples/options/sampling_cluster2n1.py
new file mode 100644
index 00000000..125984f2
--- /dev/null
+++ b/Examples/options/sampling_cluster2n1.py
@@ -0,0 +1,58 @@
+from Gaudi.Configuration import *
+from GaudiKernel import SystemOfUnits as units
+
+from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+
+from Configurables import PodioInput
+from Configurables import Jug__Digi__EcalTungstenSamplingDigi as EcalTungstenSamplingDigi
+from Configurables import Jug__Reco__EcalTungstenSamplingReco as EcalTungstenSamplingReco
+from Configurables import Jug__Reco__SamplingECalHitsMerger as SamplingECalHitsMerger
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+
+geo_service = GeoSvc("GeoSvc", detectors=["../topside/test.xml"])
+podioevent = EICDataSvc("EventDataSvc", inputs=["barrel_electrons.root"], OutputLevel=DEBUG)
+
+podioinput = PodioInput("PodioReader", collections=["mcparticles", "EcalBarrelHits"], OutputLevel=DEBUG)
+emcaldigi = EcalTungstenSamplingDigi("ecal_digi",
+ inputHitCollection="EcalBarrelHits",
+ outputHitCollection="DigiEcalBarrelHits",
+ inputEnergyUnit=units.GeV,
+ inputTimeUnit=units.ns,
+ OutputLevel=DEBUG)
+emcalreco = EcalTungstenSamplingReco("ecal_reco",
+ inputHitCollection="DigiEcalBarrelHits",
+ outputHitCollection="RecoEcalBarrelHits",
+ OutputLevel=DEBUG)
+# readout id definition for barrel ecal
+# <id>system:8,barrel:3,module:4,layer:6,slice:5,x:32:-16,y:-16</id>
+# xy_merger sum layers/slices, masking (8+3+4, 8+3+4+5+6-1)
+xymerger = SamplingECalHitsMerger("ecal_xy_merger",
+ cellIDMaskRanges=[(15, 25)],
+ inputHitCollection="RecoEcalBarrelHits",
+ outputHitCollection="RecoEcalBarrelHitsXY")
+# xy_merger sum modules, masking (8+3+4+5+6, 8+3+4+5+6+32-1)
+zmerger = SamplingECalHitsMerger("ecal_z_merger",
+ cellIDMaskRanges=[(26, 57)],
+ inputHitCollection="RecoEcalBarrelHits",
+ outputHitCollection="RecoEcalBarrelHitsZ")
+emcalcluster = IslandCluster(inputHitCollection="RecoEcalBarrelHitsXY",
+ outputClusterCollection="EcalBarrelClusters",
+ minClusterCenterEdep=5.0*units.MeV,
+ splitCluster=False,
+ groupRange=5.0)
+clusterreco = RecoCoG(clusterCollection="EcalBarrelClusters", logWeightBase=6.2, OutputLevel=DEBUG)
+
+
+out = PodioOutput("out", filename="barrel_cluster.root")
+out.outputCommands = ["keep *"]
+
+ApplicationMgr(
+ TopAlg = [podioinput, emcaldigi, emcalreco, xymerger, zmerger, emcalcluster, clusterreco, out],
+ EvtSel = 'NONE',
+ EvtMax = 1000,
+ ExtSvc = [podioevent],
+ OutputLevel=DEBUG
+ )
+
diff --git a/Examples/sampling_calo.py b/Examples/sampling_calo.py
new file mode 100644
index 00000000..08d45fb6
--- /dev/null
+++ b/Examples/sampling_calo.py
@@ -0,0 +1,54 @@
+import numpy as np
+import pandas as pd
+import ROOT
+from ROOT import gROOT
+import argparse
+from matplotlib import pyplot as plt
+from matplotlib.patches import Ellipse
+import matplotlib.transforms as transforms
+
+parser = argparse.ArgumentParser(description='sampling calorimeter performance')
+parser.add_argument('file', type=str, help='path to root file')
+args = parser.parse_args()
+
+gROOT.Macro('rootlogon.C')
+# root dataframe
+rdf = ROOT.RDataFrame('events', args.file)
+
+rdf = rdf.Define('fraction', 'EcalBarrelClusters.energy/5000')\
+ .Define('r', 'sqrt(EcalBarrelClusters.position.x*EcalBarrelClusters.position.x + EcalBarrelClusters.position.y*EcalBarrelClusters.position.y)')\
+ .Define('z', 'EcalBarrelClusters.position.z')\
+ .Define('angle', 'acos(z/sqrt(r*r + z*z))/M_PI*180.')
+
+hist = rdf.Histo1D(ROOT.RDF.TH1DModel('energy', 'Sampling Fraction;Fraction;Counts', 280, 0.04, 0.6), 'fraction')
+c1 = ROOT.TCanvas('c1', '', 2560, 1440)
+hist.Fit('gaus')
+hist.Draw()
+c1.SaveAs('sampling_energy.png')
+
+c1 = ROOT.TCanvas('c1', '', 2560, 1440)
+hist = rdf.Histo1D(ROOT.RDF.TH1DModel('angle', ';Angle (deg);Counts', 180, 0, 180), 'angle')
+hist.Draw()
+c1.SaveAs('sampling_angle.png')
+
+c1 = ROOT.TCanvas('c1', '', 2560, 1440)
+hist = rdf.Histo2D(ROOT.RDF.TH2DModel('2d', ';Fraction;Angle (deg)', 40, 0.1, 0.5, 60, 60, 120), 'fraction', 'angle')
+hist.Draw('colz')
+c1.SaveAs('sampling_2d.png')
+
+f = ROOT.TFile(args.file)
+events = f.events
+
+for iev in range(events.GetEntries()):
+ events.GetEntry(iev)
+ xedep = []
+ for hit in events.RecoEcalBarrelHitsZ:
+ xedep.append((np.sqrt(hit.position.x**2 + hit.position.y**2), hit.energy))
+ df = pd.DataFrame(data=xedep, columns=['x', 'edep']).set_index('x', drop=True)
+ df.sort_index(inplace=True)
+ df.index -= 890
+ print(df.cumsum())
+
+fig, ax = plt.subplots()
+ax.plot(df.index, df.cumsum()['edep'].values)
+fig.savefig('test.png')
diff --git a/JugDigi/src/components/EcalTungstenSamplingDigi.cpp b/JugDigi/src/components/EcalTungstenSamplingDigi.cpp
index b77f6a0d..e36deb95 100644
--- a/JugDigi/src/components/EcalTungstenSamplingDigi.cpp
+++ b/JugDigi/src/components/EcalTungstenSamplingDigi.cpp
@@ -1,3 +1,4 @@
+// Digitize the simulation outputs from Ecal Tungsten Sampling Calorimeter
#include <algorithm>
#include <cmath>
@@ -15,6 +16,8 @@
#include "eicd/RawCalorimeterHitCollection.h"
#include "eicd/RawCalorimeterHitData.h"
+using namespace Gaudi::Units;
+
namespace Jug {
namespace Digi {
@@ -24,8 +27,15 @@ namespace Jug {
*/
class EcalTungstenSamplingDigi : public GaudiAlgorithm {
public:
- Gaudi::Property<double> m_energyResolution{this, "energyResolution", 0.11}; // 11%sqrt(E)
- Rndm::Numbers m_gaussDist;
+ Gaudi::Property<double> m_eRes{this, "energyResolution", 0.11}; // 11%sqrt(E)
+ Gaudi::Property<double> m_tRes{this, "timineResolution", 0.1*ns};
+ Gaudi::Property<double> m_eUnit{this, "inputEnergyUnit", GeV};
+ Gaudi::Property<double> m_tUnit{this, "inputTimeUnit", ns};
+ Gaudi::Property<int> m_capADC{this, "capacityADC", 8096};
+ Gaudi::Property<double> m_dyRangeADC{this, "DynamicRangeADC", 100*MeV};
+ Gaudi::Property<int> m_pedMeanADC{this, "pedestalMean", 400};
+ Gaudi::Property<double> m_pedSigmaADC{this, "pedestalSigma", 3.2};
+ Rndm::Numbers m_normDist;
DataHandle<dd4pod::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
this};
DataHandle<eic::RawCalorimeterHitCollection> m_outputHitCollection{"outputHitCollection",
@@ -43,7 +53,7 @@ namespace Jug {
if (GaudiAlgorithm::initialize().isFailure())
return StatusCode::FAILURE;
IRndmGenSvc* randSvc = svc<IRndmGenSvc>("RndmGenSvc", true);
- StatusCode sc = m_gaussDist.initialize(randSvc, Rndm::Gauss(0.0, m_energyResolution.value()));
+ StatusCode sc = m_normDist.initialize(randSvc, Rndm::Gauss(0.0, 1.0));
if (!sc.isSuccess()) {
return StatusCode::FAILURE;
}
@@ -59,11 +69,12 @@ namespace Jug {
auto rawhits = m_outputHitCollection.createAndPut();
eic::RawCalorimeterHitCollection* rawHitCollection = new eic::RawCalorimeterHitCollection();
for (const auto& ahit : *simhits) {
- double res = m_gaussDist() / sqrt(ahit.energyDeposit());
+ double res = m_normDist()*m_eRes / sqrt(ahit.energyDeposit()*m_eUnit/GeV);
+ double ped = m_pedMeanADC + m_normDist()*m_pedSigmaADC;
eic::RawCalorimeterHit rawhit(
(long long)ahit.cellID(),
- std::llround(ahit.energyDeposit() * (1. + res) * 1.0e6), // convert to keV integer
- (double)ahit.truth().time * 1.0e6);
+ std::llround(ped + ahit.energyDeposit()*(1. + res) * m_eUnit/m_dyRangeADC*m_capADC), // convert to ADC Value
+ (double)ahit.truth().time*m_tUnit/ns + m_normDist()*m_tRes/ns);
rawhits->push_back(rawhit);
}
return StatusCode::SUCCESS;
diff --git a/JugReco/CMakeLists.txt b/JugReco/CMakeLists.txt
index 184fdb12..675b92c6 100644
--- a/JugReco/CMakeLists.txt
+++ b/JugReco/CMakeLists.txt
@@ -46,11 +46,13 @@ gaudi_add_module(JugRecoPlugins
src/components/CalorimeterIslandCluster.cpp
src/components/CrystalEndcapsReco.cpp
src/components/EcalTungstenSamplingReco.cpp
+ src/components/EcalTungstenSamplingCluster.cpp
src/components/ClusterRecoCoG.cpp
src/components/PhotoMultiplierReco.cpp
src/components/PhotoRingClusters.cpp
src/components/FuzzyKClusters.cpp
src/components/EMCalReconstruction.cpp
+ src/components/SamplingECalHitsMerger.cpp
LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO EICD::eicd DDRec Acts)
target_compile_options(JugRecoPlugins PRIVATE -Wno-suggest-override)
diff --git a/JugReco/src/components/CalorimeterIslandCluster.cpp b/JugReco/src/components/CalorimeterIslandCluster.cpp
index e13ab6ea..d8a6e204 100644
--- a/JugReco/src/components/CalorimeterIslandCluster.cpp
+++ b/JugReco/src/components/CalorimeterIslandCluster.cpp
@@ -36,6 +36,7 @@ namespace Jug::Reco {
class CalorimeterIslandCluster : public GaudiAlgorithm {
public:
+ Gaudi::Property<bool> m_splitCluster{this, "splitCluster", true};
Gaudi::Property<double> m_groupRange{this, "groupRange", 1.8};
Gaudi::Property<double> m_minClusterCenterEdep{this, "minClusterCenterEdep", 50.0*MeV};
DataHandle<eic::CalorimeterHitCollection>
@@ -86,7 +87,7 @@ namespace Jug::Reco {
debug() << "we have " << groups.size() << " groups of hits" << endmsg;
for (auto &group : groups) {
- auto maxima = find_local_maxima(group);
+ auto maxima = find_maxima(group, !m_splitCluster);
split_group(group, maxima, clusters, split_hits);
debug() << "hits in a group: " << group.size() << ", "
<< "local maxima: " << maxima.hits_size() << endmsg;
@@ -120,9 +121,24 @@ private:
}
// find local maxima that above a certain threshold
- eic::Cluster find_local_maxima(const std::vector<eic::CalorimeterHit> &group) const
+ eic::Cluster find_maxima(const std::vector<eic::CalorimeterHit> &group, bool global = false) const
{
eic::Cluster maxima;
+ if (group.empty()) {
+ return maxima;
+ }
+
+ if (global) {
+ int mpos = 0;
+ for (size_t i = 0; i < group.size(); ++i) {
+ if (group[mpos].energy() < group[i].energy()) {
+ mpos = i;
+ }
+ }
+ maxima.addhits(group[mpos]);
+ return maxima;
+ }
+
for(auto &hit : group)
{
// not a qualified center
diff --git a/JugReco/src/components/ClusterRecoCoG.cpp b/JugReco/src/components/ClusterRecoCoG.cpp
index 786b4795..45c530d7 100644
--- a/JugReco/src/components/ClusterRecoCoG.cpp
+++ b/JugReco/src/components/ClusterRecoCoG.cpp
@@ -75,8 +75,8 @@ public:
auto hit = reconstruct(cl);
cl.energy(hit.energy());
cl.position(hit.position());
- // info() << cl.energy()/GeV << " GeV, (" << cl.position().x/mm << ", "
- // << cl.position().y/mm << ", " << cl.position().z/mm << ")" << endmsg;
+ debug() << cl.hits_size() << " hits: " << cl.energy()/GeV << " GeV, (" << cl.position().x/mm << ", "
+ << cl.position().y/mm << ", " << cl.position().z/mm << ")" << endmsg;
}
return StatusCode::SUCCESS;
@@ -95,6 +95,7 @@ private:
float totalE = 0., maxE = 0.;
auto centerID = cl.hits_begin()->cellID();
for (auto &hit : cl.hits()) {
+ // info() << "hit energy = " << hit.energy() << endmsg;
auto energy = hit.energy();
totalE += energy;
if (energy > maxE) {
@@ -109,11 +110,18 @@ private:
float tw = 0., x = 0., y = 0., z = 0.;
for (auto &hit : cl.hits()) {
// suppress low energy contributions
+ // info() << std::log(hit.energy()/totalE) << endmsg;
float w = std::max(0., m_logWeightBase + std::log(hit.energy()/totalE));
tw += w;
x += hit.local_x() * w;
y += hit.local_y() * w;
z += hit.local_z() * w;
+ /*
+ debug() << hit.cellID()
+ << "(" << hit.local_x() << ", " << hit.local_y() << ", " << hit.local_z() << "), "
+ << "(" << hit.x() << ", " << hit.y() << ", " << hit.z() << "), "
+ << endmsg;
+ */
}
res.local({x/tw, y/tw, z/tw + m_depthCorrection});
diff --git a/JugReco/src/components/EcalTungstenSamplingCluster.cpp b/JugReco/src/components/EcalTungstenSamplingCluster.cpp
new file mode 100644
index 00000000..c3810dbb
--- /dev/null
+++ b/JugReco/src/components/EcalTungstenSamplingCluster.cpp
@@ -0,0 +1,116 @@
+#include <algorithm>
+
+#include "Gaudi/Property.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiAlg/GaudiTool.h"
+#include "GaudiAlg/Transformer.h"
+#include "GaudiKernel/PhysicalConstants.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "DDRec/CellIDPositionConverter.h"
+#include "DDRec/Surface.h"
+#include "DDRec/SurfaceManager.h"
+
+// FCCSW
+#include "JugBase/DataHandle.h"
+#include "JugBase/IGeoSvc.h"
+
+// Event Model related classes
+#include "eicd/CalorimeterHitCollection.h"
+#include "eicd/ClusterCollection.h"
+
+using namespace Gaudi::Units;
+
+namespace Jug::Reco {
+
+ class EcalTungstenSamplingCluster : public GaudiAlgorithm {
+ public:
+ Gaudi::Property<double> m_minModuleEdep{this, "minModuleEdep", 0.5 * MeV};
+ DataHandle<eic::CalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
+ this};
+ DataHandle<eic::ClusterCollection> m_outputClusterCollection{"outputClusterCollection", Gaudi::DataHandle::Writer,
+ this};
+
+ /// Pointer to the geometry service
+ SmartIF<IGeoSvc> m_geoSvc;
+
+ // ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
+ EcalTungstenSamplingCluster(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc)
+ {
+ declareProperty("inputHitCollection", m_inputHitCollection, "");
+ declareProperty("outputClusterCollection", m_outputClusterCollection, "");
+ }
+
+ StatusCode initialize() override
+ {
+ if (GaudiAlgorithm::initialize().isFailure()) {
+ return StatusCode::FAILURE;
+ }
+ m_geoSvc = service("GeoSvc");
+ if (!m_geoSvc) {
+ error() << "Unable to locate Geometry Service. "
+ << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
+ return StatusCode::FAILURE;
+ }
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode execute() override
+ {
+ // input collections
+ const auto& hits = *m_inputHitCollection.get();
+ // Create output collections
+ auto& clusters = *m_outputClusterCollection.createAndPut();
+
+ // Loop over hits
+ std::vector<bool> visits(hits.size(), false);
+ double ref_pos_z = 0.0;
+ for (size_t i = 0; i < hits.size(); i++)
+ {
+ // Ignore hits that already visited
+ if (visits[i]) {
+ continue;
+ }
+ // Select above minimum energy deposit
+ if (hits[i].energy() > m_minModuleEdep/GeV) {
+ // Reference z position of hit
+ ref_pos_z = hits[i].position().z;
+ std::cout << "Before ref_pos_z: " << ref_pos_z << std::endl;
+ // Call function to add up energy deposit in the same layer based on z position
+ group_by_layer(i, ref_pos_z, hits, visits, clusters);
+ }
+ }
+ return StatusCode::SUCCESS;
+ }
+ private:
+ // Grouping hits by layers
+ void group_by_layer(int index, double ref_pos_z, const eic::CalorimeterHitCollection &hits, std::vector<bool> &visits, eic::ClusterCollection &clusters) const
+ {
+ visits[index] = true;
+ auto tot_edep = hits[index].energy();
+ double pos_x = hits[index].position().x;
+ double pos_y = hits[index].position().y;
+ double pos_z = hits[index].position().z;
+ double temp = ref_pos_z;
+ std::cout << "After ref_pos_z: " << temp << std::endl;
+ std::cout << "Reprint pos_z: " << pos_z << std::endl;
+ // Loop over hits to find hits on the same layer
+ for (size_t i = 0; i < hits.size(); i++) {
+ if(visits[i]) {
+ continue;
+ }
+ // Add up energy deposit based on the same z position and above energy threshold
+ if((double)hits[i].position().z == ref_pos_z && hits[i].energy() > m_minModuleEdep/GeV) {
+ tot_edep += hits[i].energy();
+ visits[i] = true;
+ }
+ }
+ // Save info as a cluster
+ // TODO: position x and y determination
+ clusters.push_back(eic::Cluster{tot_edep,{pos_x,pos_y,pos_z}, {0,0,0,0,0,0}, 0, 0});
+ return;
+ }
+ };
+ DECLARE_COMPONENT(EcalTungstenSamplingCluster)
+} // namespace Jug::Reco
diff --git a/JugReco/src/components/EcalTungstenSamplingReco.cpp b/JugReco/src/components/EcalTungstenSamplingReco.cpp
index 23e20b2a..5bedd82f 100644
--- a/JugReco/src/components/EcalTungstenSamplingReco.cpp
+++ b/JugReco/src/components/EcalTungstenSamplingReco.cpp
@@ -1,4 +1,8 @@
+// Reconstruct digitized outputs fof Ecal Tungsten Sampling Calorimeter
+// It is exactly the reverse step of JugDigi/src/components/EcalTungstenSamplingDigi.cpp
+
#include <algorithm>
+#include <bitset>
#include "Gaudi/Property.h"
#include "GaudiAlg/GaudiAlgorithm.h"
@@ -25,8 +29,11 @@ using namespace Gaudi::Units;
namespace Jug::Reco {
class EcalTungstenSamplingReco : public GaudiAlgorithm {
public:
- Gaudi::Property<double> m_samplingFraction{this, "samplingFraction", 0.25};
- Gaudi::Property<double> m_minModuleEdep{this, "minModuleEdep", 0.5 * MeV};
+ Gaudi::Property<int> m_capADC{this, "capacityADC", 8096};
+ Gaudi::Property<double> m_dyRangeADC{this, "DynamicRangeADC", 100*MeV};
+ Gaudi::Property<int> m_pedMeanADC{this, "pedestalMean", 400};
+ Gaudi::Property<double> m_pedSigmaADC{this, "pedestalSigma", 3.2};
+ Gaudi::Property<double> m_thresholdADC{this, "thresholdFactor", 3.0};
DataHandle<eic::RawCalorimeterHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader,
this};
DataHandle<eic::CalorimeterHitCollection> m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer,
@@ -64,24 +71,31 @@ namespace Jug::Reco {
// energy time reconstruction
for (auto& rh : rawhits) {
- float energy = rh.amplitude() / 1.0e6; // convert keV -> GeV
- if (energy >= (m_minModuleEdep / GeV)) {
- float time = rh.timeStamp() / 1.0e6; // ns
- auto id = rh.cellID();
- // global positions
- auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
- // local positions
- auto pos = m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position();
- // cell dimension
- auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(id);
- hits.push_back(eic::CalorimeterHit{id,
- energy/m_samplingFraction,
- time,
- {gpos.x() / dd4hep::mm, gpos.y() / dd4hep::mm, gpos.z() / dd4hep::mm},
- {pos.x() / dd4hep::mm, pos.y() / dd4hep::mm, pos.z() / dd4hep::mm},
- {dim[0] / dd4hep::mm, dim[1] / dd4hep::mm, 0.0},
- 0});
+ // did not pass the threshold
+ if ((rh.amplitude() - m_pedMeanADC) < m_thresholdADC*m_pedSigmaADC) {
+ continue;
}
+ float energy = (rh.amplitude() - m_pedMeanADC) / (float) m_capADC * m_dyRangeADC; // convert ADC -> energy
+ float time = rh.timeStamp(); // ns
+ auto id = rh.cellID();
+ // global positions
+ auto gpos = m_geoSvc->cellIDPositionConverter()->position(id);
+ // local positions
+ auto volman = m_geoSvc->detector()->volumeManager();
+ auto alignment = volman.lookupDetector(id).nominal();
+ auto pos = alignment.worldToLocal(dd4hep::Position(gpos.x(), gpos.y(), gpos.z()));
+ // auto pos = m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().position();
+ // cell dimension
+ auto dim = m_geoSvc->cellIDPositionConverter()->cellDimensions(id);
+ // info() << std::bitset<64>(id) << "\n"
+ // << m_geoSvc->cellIDPositionConverter()->findContext(id)->volumePlacement().volIDs().str() << endmsg;
+ hits.push_back(eic::CalorimeterHit{id,
+ energy,
+ time,
+ {gpos.x() / dd4hep::mm, gpos.y() / dd4hep::mm, gpos.z() / dd4hep::mm},
+ {pos.x() / dd4hep::mm, pos.y() / dd4hep::mm, pos.z() / dd4hep::mm},
+ {dim[0] / dd4hep::mm, dim[1] / dd4hep::mm, 0.0},
+ 0});
}
return StatusCode::SUCCESS;
diff --git a/JugReco/src/components/SamplingECalHitsMerger.cpp b/JugReco/src/components/SamplingECalHitsMerger.cpp
new file mode 100644
index 00000000..c7926085
--- /dev/null
+++ b/JugReco/src/components/SamplingECalHitsMerger.cpp
@@ -0,0 +1,110 @@
+/*
+ * A clustering algorithm to reduce 3D clustering to 2D (x, y) + 1D (depth)
+ * 2D clustering is formed by summing all layers in the same module
+ * 1D clustering is formed by summing all modules in the same layer
+ *
+ * Author: Chao Peng (ANL), 03/31/2021
+ */
+#include <bitset>
+#include <algorithm>
+#include <unordered_map>
+
+#include "Gaudi/Property.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiAlg/Transformer.h"
+#include "GaudiAlg/GaudiTool.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include "GaudiKernel/PhysicalConstants.h"
+
+#include "DDRec/CellIDPositionConverter.h"
+#include "DDRec/SurfaceManager.h"
+#include "DDRec/Surface.h"
+
+// FCCSW
+#include "JugBase/DataHandle.h"
+#include "JugBase/IGeoSvc.h"
+
+// Event Model related classes
+#include "eicd/CalorimeterHitCollection.h"
+
+using namespace Gaudi::Units;
+
+namespace Jug::Reco {
+
+class SamplingECalHitsMerger : public GaudiAlgorithm {
+public:
+ Gaudi::Property<std::vector<std::pair<int, int>>>
+ u_cellIDMaskRanges{this, "cellIDMaskRanges", {{0, 31}}};
+ DataHandle<eic::CalorimeterHitCollection>
+ m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<eic::CalorimeterHitCollection>
+ m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
+ int64_t id_mask;
+
+ // ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
+ SamplingECalHitsMerger(const std::string& name, ISvcLocator* svcLoc)
+ : GaudiAlgorithm(name, svcLoc)
+ {
+ declareProperty("inputHitCollection", m_inputHitCollection, "");
+ declareProperty("outputHitCollection", m_outputHitCollection, "");
+ }
+
+ StatusCode initialize() override
+ {
+ if (GaudiAlgorithm::initialize().isFailure()) {
+ return StatusCode::FAILURE;
+ }
+
+ // build masks from input
+ id_mask = 0;
+ for (auto &p : u_cellIDMaskRanges) {
+ debug() << "masking bit " << p.first << " - " << p.second << endmsg;
+ for (int64_t k = p.first; k <= p.second; ++k) {
+ id_mask |= (int64_t(1) << k);
+ }
+ }
+ id_mask = ~id_mask;
+ debug() << "cellID mask = " << std::bitset<64>(id_mask) << endmsg;
+
+ return StatusCode::SUCCESS;
+ }
+
+ StatusCode execute() override
+ {
+ // input collections
+ const auto &hits = *m_inputHitCollection.get();
+ // Create output collections
+ auto &mhits = *m_outputHitCollection.createAndPut();
+
+ // sum energies that has the same id
+ std::unordered_map<long long, size_t> merge_map;
+ for (auto &h : hits) {
+ auto id = (h.cellID() & id_mask);
+ // debug() << h.cellID() << " - " << std::bitset<64>(h.cellID()) << endmsg;
+ auto it = merge_map.find(id);
+ if (it == merge_map.end()) {
+ merge_map[id] = mhits.size();
+ mhits.push_back(h.clone());
+ debug() << mhits[mhits.size() - 1].cellID() << " - " << std::bitset<64>(id) << endmsg;
+ } else {
+ mhits[it->second].energy(mhits[it->second].energy() + h.energy());
+ }
+ }
+
+ for (auto &h : mhits) {
+ debug() << h.cellID() << ": " << h.energy() << endmsg;
+ }
+
+ debug() << "Size before = " << hits.size() << ", after = " << mhits.size() << endmsg;
+
+ return StatusCode::SUCCESS;
+ }
+
+
+}; // class SamplingECalHitsMerger
+
+DECLARE_COMPONENT(SamplingECalHitsMerger)
+
+} // namespace Jug::Reco
+
--
GitLab