[EMCAL] implementation of number of local maxima variable

Author: fjonasALICE

Detectors/EMCAL/base/include/EMCALBase/ClusterFactory.h

@@ -401,6 +401,10 @@ class ClusterFactory
   /// in cell units
   void evalElipsAxis(gsl::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const;
 
+  ///
+  /// Calculate the number of local maxima in the cluster
+  void evalNExMax(gsl::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const;
+
   ///
   /// Time is set to the time of the digit with the maximum energy
   void evalTime(gsl::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const;

Detectors/EMCAL/base/include/EMCALBase/Geometry.h

@@ -429,6 +429,16 @@ class Geometry
   /// \return Position (0 - phi, 1 - eta) of the cell inside teh supermodule
   std::tuple<int, int> GetCellPhiEtaIndexInSModule(int supermoduleID, int moduleID, int phiInModule, int etaInModule) const;
 
+
+  /// \brief Get topological row and column of cell in SM (same as for clusteriser with artifical gaps)
+  /// \param supermoduleID super module number
+  /// \param moduleID module number
+  /// \param phiInModule index in phi direction in module
+  /// \param etaInModule index in phi direction in module
+  /// \return tuple with (row, column) of the cell, which is global numbering scheme
+  std::tuple<int,int> GetTopologicalRowColumn(int supermoduleID, int moduleID, int phiInModule, int etaInModule) const;
+
+
   /// \brief Adapt cell indices in supermodule to online indexing
   /// \param supermoduleID super module number of the channel/cell
   /// \param iphi row/phi cell index, modified for DCal

Detectors/EMCAL/base/src/ClusterFactory.cxx

@@ -120,6 +120,9 @@ o2::emcal::AnalysisCluster ClusterFactory<InputType>::buildCluster(int clusterIn
   evalElipsAxis(inputsIndices, clusterAnalysis);
   evalDispersion(inputsIndices, clusterAnalysis);
 
+  // evaluate number of local maxima
+  evalNExMax(inputsIndices, clusterAnalysis);
+
   evalCoreEnergy(inputsIndices, clusterAnalysis);
   evalTime(inputsIndices, clusterAnalysis);
 
@@ -489,6 +492,59 @@ void ClusterFactory<InputType>::evalCoreEnergy(gsl::span<const int> inputsIndice
   clusterAnalysis.setCoreEnergy(coreEnergy);
 }
 
+///
+/// Calculate the number of local maxima in the cluster
+//____________________________________________________________________________
+template <class InputType>
+void ClusterFactory<InputType>::evalNExMax(gsl::span<const int> inputsIndices, AnalysisCluster& clusterAnalysis) const
+{
+  // Pre-compute cell indices and energies for all cells in cluster to avoid multiple expensive geometry lookups
+  struct CellInfo {
+    int row;
+    int column;
+    double energy;
+  };
+  
+  std::vector<CellInfo> cellInfos;
+  cellInfos.reserve(inputsIndices.size());
+  
+  for (auto iInput : inputsIndices) {
+    auto [nSupMod, nModule, nIphi, nIeta] = mGeomPtr->GetCellIndex(mInputsContainer[iInput].getTower());
+    
+    // get a nice topological indexing that is done in exactly the same way as used by the clusterizer
+    // this way we can handle the shared cluster cases correctly
+    auto [row, column] = mGeomPtr->GetTopologicalRowColumn(nSupMod, nModule, nIphi, nIeta);
+    cellInfos.push_back({row, column, mInputsContainer[iInput].getEnergy()});
+  }
+  
+  // Now find local maxima using pre-computed data
+  int nExMax = 0;
+  for (size_t i = 0; i < cellInfos.size(); i++) {
+    // this cell is assumed to be local maximum unless we find a higher energy cell in the neighborhood
+    bool isExMax = true;
+    
+    // loop over all other cells in cluster
+    for (size_t j = 0; j < cellInfos.size(); j++) {
+      if (i == j) continue;
+      
+      // adjacent cell is any cell with adjacent phi or eta index
+      if (std::abs(cellInfos[i].row - cellInfos[j].row) <= 1 && 
+          std::abs(cellInfos[i].column - cellInfos[j].column) <= 1) {
+        
+        // if there is a cell with higher energy than the current cell, it is not a local maximum
+        if (cellInfos[j].energy > cellInfos[i].energy) {
+          isExMax = false;
+          break;
+        }
+      }
+    }
+    if (isExMax) {
+      nExMax++;
+    }
+  }
+  clusterAnalysis.setNExMax(nExMax);
+}
+
 ///
 /// Calculates the axis of the shower ellipsoid in eta and phi
 /// in cell units

Detectors/EMCAL/base/src/Geometry.cxx

@@ -1103,6 +1103,30 @@ std::tuple<int, int> Geometry::GetCellPhiEtaIndexInSModule(int supermoduleID, in
   return std::make_tuple(phiInSupermodule, etaInSupermodule);
 }
 
+std::tuple<int, int> Geometry::GetTopologicalRowColumn(int supermoduleID, int moduleID, int phiInModule, int etaInModule) const
+{
+  auto [iphi, ieta] = GetCellPhiEtaIndexInSModule(supermoduleID, moduleID, phiInModule, etaInModule);
+  int row = iphi;
+  int column = ieta;
+
+  // Add shifts wrt. supermodule and type of calorimeter
+  // NOTE:
+  // * Rows (phi) are arranged that one space is left empty between supermodules in phi
+  //   This is due to the physical gap that forbids clustering
+  // * For DCAL, there is an additional empty column between two supermodules in eta
+  //   Again, this is to account for the gap in DCAL
+
+  row += supermoduleID / 2 * (24 + 1);
+  // In DCAL, leave a gap between two SMs with same phi
+  if (!IsDCALSM(supermoduleID)) { // EMCAL
+    column += supermoduleID % 2 * 48;
+  } else {
+    column += supermoduleID % 2 * (48 + 1);
+  }
+
+  return std::make_tuple(row, column);
+}
+
 std::tuple<int, int> Geometry::ShiftOnlineToOfflineCellIndexes(Int_t supermoduleID, Int_t iphi, Int_t ieta) const
 {
   if (supermoduleID == 13 || supermoduleID == 15 || supermoduleID == 17) {