Merge branch 'master' into 'master'

separate cluster reconstruction into a standalone component

See merge request !8

JugReco/CMakeLists.txt

@@ -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)

JugReco/src/components/CalorimeterIslandCluster.cpp

@@ -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)

JugReco/src/components/ClusterRecoCoG.cpp

+/*
+ *  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
+

JugReco/tests/options/example_crystalendcapsreco.py

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