From a54eee1e2d0b620ccec874e6f99a349d1fbf85d2 Mon Sep 17 00:00:00 2001
From: Chao <cpeng@anl.gov>
Date: Sun, 27 Sep 2020 17:44:03 -0500
Subject: [PATCH] separate reconstruction from clustering
---
JugReco/CMakeLists.txt | 1 +
.../components/CalorimeterIslandCluster.cpp | 38 ++------
JugReco/src/components/CrystalEndcapsReco.cpp | 86 +++++++++++++++++++
.../options/example_crystalendcapsreco.py | 11 ++-
4 files changed, 103 insertions(+), 33 deletions(-)
create mode 100644 JugReco/src/components/CrystalEndcapsReco.cpp
diff --git a/JugReco/CMakeLists.txt b/JugReco/CMakeLists.txt
index 9bed9225..29500669 100644
--- a/JugReco/CMakeLists.txt
+++ b/JugReco/CMakeLists.txt
@@ -35,6 +35,7 @@ gaudi_add_module(JugRecoPlugins
src/components/TrackParamTruthInit.cpp
src/components/ParticlesFromTrackFit.cpp
src/components/CalorimeterIslandCluster.cpp
+ src/components/CrystalEndcapsReco.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 778d0c2e..8d716a1e 100644
--- a/JugReco/src/components/CalorimeterIslandCluster.cpp
+++ b/JugReco/src/components/CalorimeterIslandCluster.cpp
@@ -1,6 +1,6 @@
/*
* Island Clustering Algorithm for Calorimeter Blocks
- * 1. group all the adjacent modules with the energy deposit above <minModuleEdep>
+ * 1. group all the adjacent modules
* 2. split the groups between their local maxima with the energy deposit above <minClusterCenterEdep>
* 3. reconstruct the clustrers
*
@@ -29,7 +29,6 @@
// Event Model related classes
#include "eicd/CalorimeterHitCollection.h"
-#include "eicd/RawCalorimeterHitCollection.h"
#include "eicd/ClusterCollection.h"
using namespace Gaudi::Units;
@@ -38,10 +37,10 @@ namespace Jug::Reco {
class CalorimeterIslandCluster : public GaudiAlgorithm
{
public:
- Gaudi::Property<double> m_minModuleEdep{this, "minModuleEdep", 0.5*MeV};
+ 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::RawCalorimeterHitCollection>
+ DataHandle<eic::CalorimeterHitCollection>
m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
DataHandle<eic::ClusterCollection>
m_outputClusterCollection{"outputClusterCollection", Gaudi::DataHandle::Writer, this};
@@ -73,26 +72,10 @@ public:
StatusCode execute() override
{
// input collections
- const auto &rawhits = *m_inputHitCollection.get();
+ const auto &hits = *m_inputHitCollection.get();
// Create output collections
auto &clusters = *m_outputClusterCollection.createAndPut();
- info() << "we have " << rawhits.size() << " raw hits" << endmsg;
-
- // energy time reconstruction
- eic::CalorimeterHitCollection hits;
- for (auto &rh : rawhits) {
- float energy = rh.amplitude()/100.*MeV;
- if (energy >= m_minModuleEdep) {
- float time = rh.timeStamp();
- auto pos = m_geoSvc->cellIDPositionConverter()->position(rh.cellID0());
- hits.push_back(eic::CalorimeterHit{
- rh.cellID0(), rh.cellID1(), energy, time, {pos.X(), pos.Y(), pos.Z()}, 0
- });
- }
- }
- info() << "we have " << hits.size() << " hits" << endmsg;
-
// group neighboring hits
std::vector<bool> visits(hits.size(), false);
eic::ClusterCollection groups;
@@ -136,12 +119,12 @@ private:
// info() << std::abs(pos1.x - pos2.x) << ", " << (dim1[0] + dim2[0])/2. << ", "
// << std::abs(pos1.y - pos2.y) << ", " << (dim1[1] + dim2[1])/2. << endmsg;
- return (std::abs(pos1.x - pos2.x) <= (dim1[0] + dim2[0])/1.3) &&
- (std::abs(pos1.y - pos2.y) <= (dim1[1] + dim2[1])/1.3);
+ return (std::abs(pos1.x - pos2.x) <= (dim1[0] + dim2[0])/2.*m_groupRange) &&
+ (std::abs(pos1.y - pos2.y) <= (dim1[1] + dim2[1])/2.*m_groupRange);
}
// grouping function with Depth-First Search
- void dfs_group(eic::Cluster group, int idx, eic::CalorimeterHitCollection &hits, std::vector<bool> &visits)
+ void dfs_group(eic::Cluster group, int idx, const eic::CalorimeterHitCollection &hits, std::vector<bool> &visits)
{
auto hit = hits[idx];
group.addhits(hit);
@@ -203,12 +186,7 @@ private:
if (maxima.hits_size() == 0) {
return scoll;
} else if (maxima.hits_size() == 1) {
- auto cl = clusters.create();
- for (auto hit : group.hits()) {
- auto shit = hit.clone();
- cl.addhits(shit);
- scoll.push_back(shit);
- }
+ clusters.push_back(group.clone());
return scoll;
}
diff --git a/JugReco/src/components/CrystalEndcapsReco.cpp b/JugReco/src/components/CrystalEndcapsReco.cpp
new file mode 100644
index 00000000..8e152e85
--- /dev/null
+++ b/JugReco/src/components/CrystalEndcapsReco.cpp
@@ -0,0 +1,86 @@
+#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/CalorimeterHitCollection.h"
+#include "eicd/RawCalorimeterHitCollection.h"
+
+using namespace Gaudi::Units;
+
+namespace Jug::Reco {
+class CrystalEndcapsReco : public GaudiAlgorithm
+{
+public:
+ Gaudi::Property<double> m_minModuleEdep{this, "minModuleEdep", 0.5*MeV};
+ DataHandle<eic::RawCalorimeterHitCollection>
+ m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<eic::CalorimeterHitCollection>
+ m_outputHitCollection{"outputHitCollection", Gaudi::DataHandle::Writer, this};
+ /// Pointer to the geometry service
+ SmartIF<IGeoSvc> m_geoSvc;
+
+ // ill-formed: using GaudiAlgorithm::GaudiAlgorithm;
+ CrystalEndcapsReco(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;
+ }
+ 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 &rawhits = *m_inputHitCollection.get();
+ // Create output collections
+ auto &hits = *m_outputHitCollection.createAndPut();
+
+ // energy time reconstruction
+ for (auto &rh : rawhits) {
+ float energy = rh.amplitude()/100.*MeV;
+ if (energy >= m_minModuleEdep) {
+ float time = rh.timeStamp();
+ auto pos = m_geoSvc->cellIDPositionConverter()->position(rh.cellID0());
+ hits.push_back(eic::CalorimeterHit{
+ rh.cellID0(), rh.cellID1(), energy, time, {pos.X(), pos.Y(), pos.Z()}, 0
+ });
+ }
+ }
+
+ return StatusCode::SUCCESS;
+ }
+};
+
+DECLARE_COMPONENT(CrystalEndcapsReco)
+
+} // namespace Jug::Reco
+
+
diff --git a/JugReco/tests/options/example_crystalendcapsreco.py b/JugReco/tests/options/example_crystalendcapsreco.py
index e091016a..2ab6ff0e 100644
--- a/JugReco/tests/options/example_crystalendcapsreco.py
+++ b/JugReco/tests/options/example_crystalendcapsreco.py
@@ -2,23 +2,28 @@ from Gaudi.Configuration import *
from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+from GaudiKernel import SystemOfUnits as units
geo_service = GeoSvc("GeoSvc", detectors=["../NPDet/src/GenericDetectors/calorimeters/compact/Crystal_example.xml"])
podioevent = EICDataSvc("EventDataSvc", inputs=["output_emcal_electrons_npsim.root"], OutputLevel=DEBUG)
from Configurables import PodioInput
-from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
from Configurables import Jug__Digi__CrystalEndcapsDigi as CrystalEndcapsDigi
+from Configurables import Jug__Reco__CrystalEndcapsReco as CrystalEndcapsReco
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
podioinput = PodioInput("PodioReader", collections=["mcparticles","EcalHits"], OutputLevel=DEBUG)
emcaldigi = CrystalEndcapsDigi(inputHitCollection="EcalHits", outputHitCollection="RawDigiEcalHits")
-emcalcluster = IslandCluster(inputHitCollection="RawDigiEcalHits", outputClusterCollection="EcalClusters")
+emcalreco = CrystalEndcapsReco(inputHitCollection="RawDigiEcalHits", outputHitCollection="RecoEcalHits",
+ minModuleEdep=1.0*units.MeV)
+emcalcluster = IslandCluster(inputHitCollection="RecoEcalHits", outputClusterCollection="EcalClusters",
+ minClusterCenterEdep=30*units.MeV, groupRange=2.0)
out = PodioOutput("out", filename="reco_emcal_electrons_npsim.root")
out.outputCommands = ["keep EcalClusters"]
ApplicationMgr(
- TopAlg = [podioinput, emcaldigi, emcalcluster, out],
+ TopAlg = [podioinput, emcaldigi, emcalreco, emcalcluster, out],
EvtSel = 'NONE',
EvtMax = 100,
ExtSvc = [podioevent],
--
GitLab