Skip to content

GitLab

Commit b18173f7 authored by Jihee Kim's avatar Jihee Kim
Browse files

Simple Clustering

parent 93d6e22a
Pipeline #7921 passed with stages
in 8 minutes and 43 seconds
Hide whitespace changes
Inline Side-by-side
Showing with 418 additions and 29 deletions
Examples/options/sampling_cluster2n1.py 0 → 100644
View file @ b18173f7
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
)
Examples/sampling_calo.py 0 → 100644
View file @ b18173f7
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')
JugDigi/src/components/EcalTungstenSamplingDigi.cpp
View file @ b18173f7
// 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;
......
JugReco/CMakeLists.txt
View file @ b18173f7
......@@ -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)
......
JugReco/src/components/CalorimeterIslandCluster.cpp
View file @ b18173f7
......@@ -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
......
JugReco/src/components/ClusterRecoCoG.cpp
View file @ b18173f7
......@@ -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});
......
JugReco/src/components/EcalTungstenSamplingCluster.cpp 0 → 100644
View file @ b18173f7
#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
JugReco/src/components/EcalTungstenSamplingReco.cpp
View file @ b18173f7
// 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;
......
JugReco/src/components/SamplingECalHitsMerger.cpp 0 → 100644
View file @ b18173f7
/*
* 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
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment