separate reconstruction from clustering

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)

JugReco/src/components/CalorimeterIslandCluster.cpp

 /*
  *  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;
         }
 

JugReco/src/components/CrystalEndcapsReco.cpp

+#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
+
+

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],