Merge pull request #25633 from guitargeek/Egamma_5x5_shapes

Some EcalClusterTools maintainance and energy matrix in Egamma MVA Ntuplizers

RecoCaloTools/Navigation/interface/CaloRectangle.h

@@ -0,0 +1,112 @@
+#ifndef RecoCaloTools_Navigation_CaloRectangle_H
+#define RecoCaloTools_Navigation_CaloRectangle_H
+
+#include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
+#include "Geometry/CaloTopology/interface/CaloTopology.h"
+
+/*
+ * CaloRectangle is a class to create a rectangular range of DetIds around a
+ * central DetId. Meant to be used for range-based loops to calculate cluster
+ * shape variables.
+ */
+
+struct CaloRectangle {
+    const int iEtaOrIXMin;
+    const int iEtaOrIXMax;
+    const int iPhiOrIYMin;
+    const int iPhiOrIYMax;
+
+    template<class T>
+    auto operator()(T home, CaloTopology const& topology);
+
+};
+
+template<class T>
+T offsetBy(T start, CaloSubdetectorTopology const& topo, int dIEtaOrIX, int dIPhiOrIY)
+{
+    for(int i = 0; i < std::abs(dIEtaOrIX) && start != T(0); i++) {
+        start = dIEtaOrIX > 0 ? topo.goEast(start) : topo.goWest(start);
+    }
+
+    for(int i = 0; i < std::abs(dIPhiOrIY) && start != T(0); i++) {
+        start = dIPhiOrIY > 0 ? topo.goNorth(start) : topo.goSouth(start);
+    }
+    return start;
+}
+
+template<class T>
+class CaloRectangleRange {
+
+  public:
+
+    class Iterator {
+
+      public:
+
+        Iterator(T const& home, int iEtaOrIX, int iPhiOrIY, CaloRectangle const rectangle, CaloSubdetectorTopology const& topology)
+          : home_(home)
+          , rectangle_(rectangle)
+          , topology_(topology)
+          , iEtaOrIX_(iEtaOrIX)
+          , iPhiOrIY_(iPhiOrIY)
+        {}
+
+        Iterator& operator++() {
+            if(iPhiOrIY_ == rectangle_.iPhiOrIYMax) {
+                iPhiOrIY_ = rectangle_.iPhiOrIYMin;
+                iEtaOrIX_++;
+            } else ++iPhiOrIY_;
+            return *this;
+        }
+
+        int iEtaOrIX() const { return iEtaOrIX_; }
+        int iPhiOrIY() const { return iPhiOrIY_; }
+
+        bool operator==(Iterator const& other) const { return iEtaOrIX_ == other.iEtaOrIX() && iPhiOrIY_ == other.iPhiOrIY(); }
+        bool operator!=(Iterator const& other) const { return iEtaOrIX_ != other.iEtaOrIX() || iPhiOrIY_ != other.iPhiOrIY(); }
+
+        T operator*() const { return offsetBy(home_, topology_, iEtaOrIX_, iPhiOrIY_); }
+
+      private:
+
+        const T home_;
+
+        const CaloRectangle rectangle_;
+        CaloSubdetectorTopology const& topology_;
+
+        int iEtaOrIX_;
+        int iPhiOrIY_;
+    };
+
+  public:
+    CaloRectangleRange(CaloRectangle rectangle, T home, CaloTopology const& topology)
+      : home_(home)
+      , rectangle_(rectangle)
+      , topology_(*topology.getSubdetectorTopology(home))
+    {}
+
+   CaloRectangleRange(int size, T home, CaloTopology const& topology)
+      : home_(home)
+      , rectangle_{-size, size, -size, size}
+      , topology_(*topology.getSubdetectorTopology(home))
+    {}
+
+    auto begin() {
+        return Iterator(home_, rectangle_.iEtaOrIXMin, rectangle_.iPhiOrIYMin, rectangle_, topology_);
+    }
+    auto end() {
+        return Iterator(home_, rectangle_.iEtaOrIXMax + 1, rectangle_.iPhiOrIYMin, rectangle_, topology_);
+    }
+
+  private:
+    const T home_;
+    const CaloRectangle rectangle_;
+    CaloSubdetectorTopology const& topology_;
+};
+
+template<class T>
+auto CaloRectangle::operator()(T home, CaloTopology const& topology) {
+    return CaloRectangleRange<T>(*this, home, topology);
+}
+
+#endif

RecoEcal/EgammaClusterProducers/src/EcalDigiSelector.cc

@@ -134,14 +134,9 @@ void EcalDigiSelector::produce(edm::Event& evt, const edm::EventSetup& es)
            const CaloClusterPtr& bcref = clus1.seed();
            const BasicCluster *bc = bcref.get();
            //Get the maximum detid
-           std::pair<DetId, float> EDetty = 
-			   EcalClusterTools::getMaximum(*bc,rechits);
-           //get the 3x3 array centered on maximum detid.
-           std::vector<DetId> detvec = 
-			   EcalClusterTools::matrixDetId(topology,EDetty.first, -1, 1, -1, 1);
-           //Loop over the 3x3
-           for (int ik = 0;ik<int(detvec.size());++ik)
-             saveTheseDetIds.push_back(detvec[ik]);
+           DetId maxDetId = EcalClusterTools::getMaximum(*bc,rechits).first;
+           // Loop over the 3x3 array centered on maximum detid
+           for(DetId detId : CaloRectangleRange(1, maxDetId, *topology)) saveTheseDetIds.push_back(detId);
 
          }
          for (int detloop=0; detloop < int(saveTheseDetIds.size());++detloop){
@@ -196,13 +191,9 @@ void EcalDigiSelector::produce(edm::Event& evt, const edm::EventSetup& es)
            const CaloClusterPtr& bcref = clus1.seed();
            const BasicCluster *bc = bcref.get();
            //Get the maximum detid
-           std::pair<DetId, float> EDetty = EcalClusterTools::getMaximum(*bc,rechits);
-           //get the 3x3 array centered on maximum detid.
-           std::vector<DetId> detvec = 
-			   EcalClusterTools::matrixDetId(topology,EDetty.first, -1, 1, -1, 1);
-           //Loop over the 3x3
-           for (int ik = 0;ik<int(detvec.size());++ik)
-             saveTheseDetIds.insert(detvec[ik]);
+           DetId maxDetId = EcalClusterTools::getMaximum(*bc,rechits).first;
+           // Loop over the 3x3 array centered on maximum detid
+           for(DetId detId : CaloRectangleRange(1, maxDetId, *topology)) saveTheseDetIds.insert(detId);
          }
          int eecounter=0;
          for (EEDigiCollection::const_iterator blah = digis->begin();