Commit 7074b891 authored by Chao Peng's avatar Chao Peng
Browse files

separate cluster reconstruction into a standalone component

parent a54eee1e
......@@ -36,6 +36,7 @@ gaudi_add_module(JugRecoPlugins
src/components/ParticlesFromTrackFit.cpp
src/components/CalorimeterIslandCluster.cpp
src/components/CrystalEndcapsReco.cpp
src/components/ClusterRecoCoG.cpp
LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO EICD::eicd DDRec Acts)
target_compile_options(JugRecoPlugins PRIVATE -Wno-suggest-override)
......
......@@ -2,7 +2,6 @@
* Island Clustering Algorithm for Calorimeter Blocks
* 1. group all the adjacent modules
* 2. split the groups between their local maxima with the energy deposit above <minClusterCenterEdep>
* 3. reconstruct the clustrers
*
* Author: Chao Peng (ANL), 09/27/2020
* References:
......@@ -39,7 +38,6 @@ class CalorimeterIslandCluster : public GaudiAlgorithm
public:
Gaudi::Property<double> m_groupRange{this, "groupRange", 1.8};
Gaudi::Property<double> m_minClusterCenterEdep{this, "minClusterCenterEdep", 50.0*MeV};
Gaudi::Property<double> m_logWeightThres{this, "logWeightThres", 4.2};
DataHandle<eic::CalorimeterHitCollection>
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::ClusterCollection>
......@@ -97,13 +95,6 @@ public:
<< "local maxima: " << maxima.hits_size() << endmsg;
}
// reconstruct hit position for the cluster
for (auto &cl : clusters) {
reconstruct(cl);
info() << cl.energy()/GeV << " GeV, (" << cl.position()[0]/mm << ", "
<< cl.position()[1]/mm << ", " << cl.position()[2]/mm << ")" << endmsg;
}
return StatusCode::SUCCESS;
}
......@@ -238,25 +229,6 @@ private:
}
return scoll;
}
void reconstruct(eic::Cluster cl) {
float totalE = 0.;
for (auto &hit : cl.hits()) {
totalE += hit.energy();
}
cl.energy(totalE);
// center of gravity with logarithmic weighting
float totalW = 0., x = 0., y = 0., z = 0.;
for (auto &hit : cl.hits()) {
float weight = m_logWeightThres + std::log(hit.energy()/totalE);
totalW += weight;
x += hit.position().x * weight;
y += hit.position().y * weight;
z += hit.position().z * weight;
}
cl.position() = std::array<float, 3>{x/totalW, y/totalW, z/totalW};
}
};
DECLARE_COMPONENT(CalorimeterIslandCluster)
......
/*
* Reconstruct the cluster with Center of Gravity method
* Logarithmic weighting is used for mimicing energy deposit in transverse direction
*
* Author: Chao Peng (ANL), 09/27/2020
*/
#include <algorithm>
#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/ClusterCollection.h"
using namespace Gaudi::Units;
namespace Jug::Reco {
class ClusterRecoCoG : public GaudiAlgorithm
{
public:
Gaudi::Property<double> m_logWeightThres{this, "logWeightThres", 4.2};
DataHandle<eic::ClusterCollection>
m_clusterCollection{"clusterCollection", Gaudi::DataHandle::Reader, this};
// ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
ClusterRecoCoG(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc)
{
declareProperty("clusterCollection", m_clusterCollection, "");
}
StatusCode initialize() override
{
if (GaudiAlgorithm::initialize().isFailure()) {
return StatusCode::FAILURE;
}
return StatusCode::SUCCESS;
}
StatusCode execute() override
{
// input collections
const auto &clusters = *m_clusterCollection.get();
// reconstruct hit position for the cluster
for (auto &cl : clusters) {
reconstruct(cl);
info() << cl.energy()/GeV << " GeV, (" << cl.position()[0]/mm << ", "
<< cl.position()[1]/mm << ", " << cl.position()[2]/mm << ")" << endmsg;
}
return StatusCode::SUCCESS;
}
private:
void reconstruct(eic::Cluster cl) {
float totalE = 0.;
for (auto &hit : cl.hits()) {
totalE += hit.energy();
}
cl.energy(totalE);
// center of gravity with logarithmic weighting
float totalW = 0., x = 0., y = 0., z = 0.;
for (auto &hit : cl.hits()) {
float weight = m_logWeightThres + std::log(hit.energy()/totalE);
totalW += weight;
x += hit.position().x * weight;
y += hit.position().y * weight;
z += hit.position().z * weight;
}
cl.position() = std::array<float, 3>{x/totalW, y/totalW, z/totalW};
}
};
DECLARE_COMPONENT(ClusterRecoCoG)
} // namespace Jug::Reco
......@@ -11,6 +11,7 @@ from Configurables import PodioInput
from Configurables import Jug__Digi__CrystalEndcapsDigi as CrystalEndcapsDigi
from Configurables import Jug__Reco__CrystalEndcapsReco as CrystalEndcapsReco
from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
podioinput = PodioInput("PodioReader", collections=["mcparticles","EcalHits"], OutputLevel=DEBUG)
emcaldigi = CrystalEndcapsDigi(inputHitCollection="EcalHits", outputHitCollection="RawDigiEcalHits")
......@@ -18,12 +19,13 @@ emcalreco = CrystalEndcapsReco(inputHitCollection="RawDigiEcalHits", outputHitCo
minModuleEdep=1.0*units.MeV)
emcalcluster = IslandCluster(inputHitCollection="RecoEcalHits", outputClusterCollection="EcalClusters",
minClusterCenterEdep=30*units.MeV, groupRange=2.0)
clusterreco = RecoCoG(clusterCollection="EcalClusters", logWeightThres=4.2)
out = PodioOutput("out", filename="reco_emcal_electrons_npsim.root")
out.outputCommands = ["keep EcalClusters"]
ApplicationMgr(
TopAlg = [podioinput, emcaldigi, emcalreco, emcalcluster, out],
TopAlg = [podioinput, emcaldigi, emcalreco, emcalcluster, clusterreco, out],
EvtSel = 'NONE',
EvtMax = 100,
ExtSvc = [podioevent],
......
Supports Markdown
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