make topoclustering work with new module-theta merged readout (#51)

* make topoclustering work with new module-theta merged readout * fix for theta readout * move ReadNoiseFromFileTool for theta-based readout to new tool in FCCee subdirectory * make inclusion of diagonal neighbours configurable * rename value of a property * check that tools are retrieved successfully
HEP-FCC · Oct 24, 2023 · bde81c8 · bde81c8
1 parent e1f62eb
commit bde81c8
Show file tree

Hide file tree

Showing 9 changed files with 1,292 additions and 18 deletions.
diff --git a/RecFCCeeCalorimeter/CMakeLists.txt b/RecFCCeeCalorimeter/CMakeLists.txt
@@ -12,6 +12,7 @@ gaudi_add_module(k4RecFCCeeCalorimeterPlugins
                       DD4hep::DDCore
                       EDM4HEP::edm4hep
                       FCCDetectors::DetSegmentation
+                      FCCDetectors::DetCommon
                       DD4hep::DDG4
                       ROOT::Core
                       ROOT::Hist

diff --git a/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.cpp b/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.cpp
@@ -117,8 +117,15 @@ StatusCode CaloTopoClusterFCCee::execute() {
             });
 
   std::map<uint, std::vector<std::pair<uint64_t, int>>> preClusterCollection;
-  CaloTopoClusterFCCee::buildingProtoCluster(m_neighbourSigma, m_lastNeighbourSigma, firstSeeds, m_allCells,
-                                        preClusterCollection);
+  StatusCode sc_buildProtoClusters =  CaloTopoClusterFCCee::buildingProtoCluster(m_neighbourSigma,
+										 m_lastNeighbourSigma,
+										 firstSeeds, m_allCells,
+										 preClusterCollection);
+  if (sc_buildProtoClusters.isFailure()) {
+    error() << "Unable to build the protoclusters!" << endmsg;
+    return StatusCode::FAILURE;
+  }
+
   // Build Clusters in edm
   debug() << "Building " << preClusterCollection.size() << " cluster." << endmsg;
   double checkTotEnergy = 0.;
@@ -132,10 +139,10 @@ StatusCode CaloTopoClusterFCCee::execute() {
     double energy = 0.;
     double deltaR = 0.;
     std::vector<double> posPhi (i.second.size());
-    std::vector<double> posEta (i.second.size());
+    std::vector<double> posTheta (i.second.size());
     std::vector<double> vecEnergy (i.second.size());
     double sumPhi = 0.;
-    double sumEta = 0.;
+    double sumTheta = 0.;
     std::map<int,int> system;
 
     for (auto pair : i.second) {
@@ -178,10 +185,10 @@ StatusCode CaloTopoClusterFCCee::execute() {
       posY += posCell.Y() * newCell.getEnergy();
       posZ += posCell.Z() * newCell.getEnergy();
       posPhi.push_back(posCell.Phi());
-      posEta.push_back(posCell.Eta());
+      posTheta.push_back(posCell.Theta());
       vecEnergy.push_back(newCell.getEnergy());
       sumPhi += posCell.Phi() * newCell.getEnergy();
-      sumEta += posCell.Eta() * newCell.getEnergy();
+      sumTheta += posCell.Theta() * newCell.getEnergy();
 
       cluster.addToHits(newCell);
       auto er = m_allCells.erase(cID);
@@ -193,10 +200,10 @@ StatusCode CaloTopoClusterFCCee::execute() {
     cluster.setPosition(edm4hep::Vector3f(posX / energy, posY / energy, posZ / energy));
     // store deltaR of cluster in time for the moment..
     sumPhi = sumPhi / energy;
-    sumEta = sumEta / energy;
+    sumTheta = sumTheta / energy;
     int counter = 0;
-    for (auto entryEta : posEta){
-      deltaR += sqrt(pow(entryEta-sumEta,2) + pow(posEta[counter]-sumPhi,2)) * vecEnergy[counter];
+    for (auto entryTheta : posTheta){
+      deltaR += sqrt(pow(entryTheta-sumTheta,2) + pow(posPhi[counter]-sumPhi,2)) * vecEnergy[counter];
       counter++;
     }
     cluster.addToShapeParameters(deltaR / energy);
@@ -208,7 +215,7 @@ StatusCode CaloTopoClusterFCCee::execute() {
       clusterWithMixedCells++;
 
     posPhi.clear();
-    posEta.clear();
+    posTheta.clear();
     vecEnergy.clear();
 
   }
@@ -225,12 +232,19 @@ void CaloTopoClusterFCCee::findingSeeds(const std::unordered_map<uint64_t, doubl
                                    std::vector<std::pair<uint64_t, double>>& aSeeds) {
   for (const auto& cell : aCells) {
     // retrieve the noise const and offset assigned to cell
-    //
     // first try to use the cache
     int system = m_decoder->get(cell.first, m_index_system);
     if (system == 4) { //ECal barrel
       int layer = m_decoder_ecal->get(cell.first, m_index_layer_ecal);
+
       double min_threshold = m_min_offset[layer] + m_min_noise[layer] * aNumSigma;
+
+      debug() << "m_min_offset[layer]   = " << m_min_offset[layer] << endmsg;
+      debug() << "m_min_noise[layer]   = " << m_min_noise[layer] << endmsg;
+      debug() << "aNumSigma   = " << aNumSigma << endmsg;
+      debug() << "min_threshold   = " << min_threshold << endmsg;
+      debug() << "abs(cell.second)   = " << abs(cell.second) << endmsg;
+
       if (abs(cell.second) < min_threshold) {
         // if we are below the minimum threshold for the full layer, no need to retrieve the exact value
         continue;
@@ -240,9 +254,10 @@ void CaloTopoClusterFCCee::findingSeeds(const std::unordered_map<uint64_t, doubl
     // we are above the minimum threshold of the layer. Let's see if we are above the threshold for this cell.
     double threshold = m_noiseTool->noiseOffset(cell.first) + ( m_noiseTool->noiseRMS(cell.first) * aNumSigma);
     if (msgLevel() <= MSG::VERBOSE){
-      verbose() << "noise offset    = " << m_noiseTool->noiseOffset(cell.first) << "GeV " << endmsg;
-      verbose() << "noise rms       = " << m_noiseTool->noiseRMS(cell.first) << "GeV " << endmsg;
-      verbose() << "seed threshold  = " << threshold << "GeV " << endmsg;
+      debug() << "noise offset    = " << m_noiseTool->noiseOffset(cell.first) << "GeV " << endmsg;
+      debug() << "noise rms       = " << m_noiseTool->noiseRMS(cell.first) << "GeV " << endmsg;
+      debug() << "seed threshold  = " << threshold << "GeV " << endmsg;
+      debug() << "======================================" << endmsg;
     }
     if (abs(cell.second) > threshold) {
       aSeeds.emplace_back(cell.first, cell.second);

diff --git a/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h b/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h
@@ -31,8 +31,7 @@ class Segmentation;
 }
 }
 
-/** @class CaloTopoClusterFCCeeAlgorithm Reconstruction/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h
- * CombinedCaloTopoCluster.h
+/** @class CaloTopoClusterFCCee k4RecCalorimeter/RecFCCeeCalorimeter/src/components/CaloTopoClusterFCCee.h
  *
  *  Algorithm building the topological clusters for the energy reconstruction, following ATLAS note
  *  ATL-LARG-PUB-2008-002.
@@ -126,15 +125,15 @@ class CaloTopoClusterFCCee : public GaudiAlgorithm {
   //ToolHandle<ICellPositionsTool> m_cellPositionsHFwdTool{"CellPositionsCaloDiscsTool", this};
 
   /// no segmentation used in HCal
-  Gaudi::Property<bool> m_noSegmentationHCalUsed{this, "noSegmentationHCal", true, "HCal Barrel readout without DD4hep eta-phi segmentation used."};
+  Gaudi::Property<bool> m_noSegmentationHCalUsed{this, "noSegmentationHCal", true, "HCal Barrel readout without DD4hep theta-module segmentation used."};
   /// Seed threshold in sigma
   Gaudi::Property<int> m_seedSigma{this, "seedSigma", 4, "number of sigma in noise threshold"};
   /// Neighbour threshold in sigma
   Gaudi::Property<int> m_neighbourSigma{this, "neighbourSigma", 2, "number of sigma in noise threshold"};
   /// Last neighbour threshold in sigma
   Gaudi::Property<int> m_lastNeighbourSigma{this, "lastNeighbourSigma", 0, "number of sigma in noise threshold"};
   /// Name of the electromagnetic calorimeter readout
-  Gaudi::Property<std::string> m_readoutName{this, "readoutName", "ECalBarrelPhiEta"};
+  Gaudi::Property<std::string> m_readoutName{this, "readoutName", "ECalBarrelModuleThetaMerged"};
 
   /// General decoder to encode the calorimeter sub-system to determine which positions tool to use
   dd4hep::DDSegmentation::BitFieldCoder* m_decoder = new dd4hep::DDSegmentation::BitFieldCoder("system:4");