Project CMSSW displayed by LXR CMSSW_14_2_X_2024-10-23-2300/​RecoEcal/​EgammaCoreTools/​interface/​EcalClusterTools.h	Source navigation Diff markup Identifier search General search
	Version: CMSSW_14_2_X_2024-10-23-2300 CMSSW_14_2_X_2024-10-22-2300 CMSSW_14_2_X_2024-10-21-2300 CMSSW_14_2_X_2024-10-20-2300 CMSSW_14_2_X_2024-10-18-2300 CMSSW_14_2_X_2024-10-17-2300 CMSSW_14_2_X_2024-10-16-2300 Revert   Change

File indexing completed on 2024-04-06 12:24:43

 #ifndef RecoEcal_EgammaCoreTools_EcalClusterTools_h
 #define RecoEcal_EgammaCoreTools_EcalClusterTools_h
 
 /** \class EcalClusterTools
  *  
  * various cluster tools (e.g. cluster shapes)
  *
  * \author Federico Ferri
  *
  * editing author: M.B. Anderson
  * 
  *
  */
 
 /** Note about "noZS" and the new templated version of EcalClusterTools
  *
  * "noZS" means "no zero-suppression" and means that the rechit quantity is 
  * computed without fractions and always including the full matrix of crystals
  * about the seed.
  *
  * To access the noZS variables use the "noZS::EcalClusterTools"
  * The interface at this point is the same as the standard tools, except the
  * behavior is altered to provide the output as described above.
  *
  * This was achieved by templating the EcalClusterTools.
  * Do not use EcalClusterToolsT<> directly.
  *
  */
 
 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
 //#include "DataFormats/Math/interface/Point3D.h"
 #include "DataFormats/Math/interface/Vector3D.h"
 //includes for ShowerShape function to work
 #include <vector>
 #include <cmath>
 #include <TMath.h>
 #include <TMatrixT.h>
 #include <TMatrixD.h>
 #include <TVectorT.h>
 #include <TVectorD.h>
 
 #include "DataFormats/DetId/interface/DetId.h"
 #include "DataFormats/EcalDetId/interface/EBDetId.h"
 #include "DataFormats/EcalDetId/interface/EEDetId.h"
 #include "RecoCaloTools/Navigation/interface/CaloRectangle.h"
 
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 
 #include "CLHEP/Geometry/Transform3D.h"
 
 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
 #include "Geometry/CaloTopology/interface/CaloTopology.h"
 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
 #include "RecoEcal/EgammaCoreTools/interface/EgammaLocalCovParamDefaults.h"
 
 class DetId;
 class CaloTopology;
 class CaloGeometry;
 
 struct Cluster2ndMoments {
   // major and minor cluster moments wrt principale axes:
   float sMaj;
   float sMin;
   // angle between sMaj and phi:
   float alpha;
   Cluster2ndMoments() : sMaj(0.), sMin(0.), alpha(0.) {}
 };
 
 template <bool noZS>
 class EcalClusterToolsT {
 public:
   // various energies in the matrix nxn surrounding the maximum energy crystal of the input cluster
   //we use an eta/phi coordinate system rather than phi/eta
   //note e3x2 does not have a definate eta/phi geometry, it takes the maximum 3x2 block containing the
   //seed regardless of whether that 3 in eta or phi
   static float e1x3(const reco::BasicCluster &cluster,
                     const EcalRecHitCollection *recHits,
                     const CaloTopology *topology);
 
   static float e3x1(const reco::BasicCluster &cluster,
                     const EcalRecHitCollection *recHits,
                     const CaloTopology *topology);
 
   static float e1x5(const reco::BasicCluster &cluster,
                     const EcalRecHitCollection *recHits,
                     const CaloTopology *topology);
 
   static float e5x1(const reco::BasicCluster &cluster,
                     const EcalRecHitCollection *recHits,
                     const CaloTopology *topology);
 
   static float e2x2(const reco::BasicCluster &cluster,
                     const EcalRecHitCollection *recHits,
                     const CaloTopology *topology);
 
   static float e3x2(const reco::BasicCluster &cluster,
                     const EcalRecHitCollection *recHits,
                     const CaloTopology *topology);
 
   static float e3x3(const reco::BasicCluster &cluster,
                     const EcalRecHitCollection *recHits,
                     const CaloTopology *topology);
 
   static float e4x4(const reco::BasicCluster &cluster,
                     const EcalRecHitCollection *recHits,
                     const CaloTopology *topology);
 
   static float e5x5(const reco::BasicCluster &cluster,
                     const EcalRecHitCollection *recHits,
                     const CaloTopology *topology);
   static int n5x5(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology);
 
   // energy in the 2x5 strip right of the max crystal (does not contain max crystal)
   // 2 crystals wide in eta, 5 wide in phi.
   static float e2x5Right(const reco::BasicCluster &cluster,
                          const EcalRecHitCollection *recHits,
                          const CaloTopology *topology);
   // energy in the 2x5 strip left of the max crystal (does not contain max crystal)
 
   static float e2x5Left(const reco::BasicCluster &cluster,
                         const EcalRecHitCollection *recHits,
                         const CaloTopology *topology);
   // energy in the 5x2 strip above the max crystal (does not contain max crystal)
   // 5 crystals wide in eta, 2 wide in phi.
 
   static float e2x5Top(const reco::BasicCluster &cluster,
                        const EcalRecHitCollection *recHits,
                        const CaloTopology *topology);
   // energy in the 5x2 strip below the max crystal (does not contain max crystal)
 
   static float e2x5Bottom(const reco::BasicCluster &cluster,
                           const EcalRecHitCollection *recHits,
                           const CaloTopology *topology);
   // energy in a 2x5 strip containing the seed (max) crystal.
   // 2 crystals wide in eta, 5 wide in phi.
   // it is the maximum of either (1x5left + 1x5center) or (1x5right + 1x5center)
   static float e2x5Max(const reco::BasicCluster &cluster,
                        const EcalRecHitCollection *recHits,
                        const CaloTopology *topology);
 
   // energies in the crystal left, right, top, bottom w.r.t. to the most energetic crystal
   static float eLeft(const reco::BasicCluster &cluster,
                      const EcalRecHitCollection *recHits,
                      const CaloTopology *topology);
 
   static float eRight(const reco::BasicCluster &cluster,
                       const EcalRecHitCollection *recHits,
                       const CaloTopology *topology);
 
   static float eTop(const reco::BasicCluster &cluster,
                     const EcalRecHitCollection *recHits,
                     const CaloTopology *topology);
 
   static float eBottom(const reco::BasicCluster &cluster,
                        const EcalRecHitCollection *recHits,
                        const CaloTopology *topology);
   // the energy of the most energetic crystal in the cluster
 
   static float eMax(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
 
   // the energy of the second most energetic crystal in the cluster
   static float e2nd(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
 
   // get the DetId and the energy of the maximum energy crystal of the input cluster
   static std::pair<DetId, float> getMaximum(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
 
   static std::vector<float> energyBasketFractionEta(const reco::BasicCluster &cluster,
                                                     const EcalRecHitCollection *recHits);
 
   static std::vector<float> energyBasketFractionPhi(const reco::BasicCluster &cluster,
                                                     const EcalRecHitCollection *recHits);
 
   // return a vector v with v[0] = etaLat, v[1] = phiLat, v[2] = lat
   static std::vector<float> lat(const reco::BasicCluster &cluster,
                                 const EcalRecHitCollection *recHits,
                                 const CaloGeometry *geometry,
                                 bool logW = true,
                                 float w0 = 4.7);
 
   // return an array v with v[0] = covEtaEta, v[1] = covEtaPhi, v[2] = covPhiPhi
   static std::array<float, 3> covariances(const reco::BasicCluster &cluster,
                                           const EcalRecHitCollection *recHits,
                                           const CaloTopology *topology,
                                           const CaloGeometry *geometry,
                                           float w0 = 4.7);
 
   // return an array v with v[0] = covIEtaIEta, v[1] = covIEtaIPhi, v[2] = covIPhiIPhi
   //this function calculates differences in eta/phi in units of crystals not global eta/phi
   //this is gives better performance in the crack regions of the calorimeter but gives otherwise identical results to covariances function
   //   except that it doesnt need an eta based correction funtion in the endcap
   //it is multipled by an approprate crystal size to ensure it gives similar values to covariances(...)
   //
   //Warning: covIEtaIEta has been studied by egamma, but so far covIPhiIPhi hasnt been studied extensively so there could be a bug in
   //         the covIPhiIEta or covIPhiIPhi calculations. I dont think there is but as it hasnt been heavily used, there might be one
   static std::array<float, 3> localCovariances(const reco::BasicCluster &cluster,
                                                const EcalRecHitCollection *recHits,
                                                const CaloTopology *topology,
                                                float w0 = EgammaLocalCovParamDefaults::kRelEnCut,
                                                const EcalPFRecHitThresholds *thresholds = nullptr,
                                                float multEB = 0.0,
                                                float multEE = 0.0);
 
   static std::array<float, 3> scLocalCovariances(const reco::SuperCluster &cluster,
                                                  const EcalRecHitCollection *recHits,
                                                  const CaloTopology *topology,
                                                  float w0 = 4.7);
 
   // cluster second moments with respect to principal axes:
   static Cluster2ndMoments cluster2ndMoments(const reco::BasicCluster &basicCluster,
                                              const EcalRecHitCollection &recHits,
                                              double phiCorrectionFactor = 0.8,
                                              double w0 = 4.7,
                                              bool useLogWeights = true);
 
   static Cluster2ndMoments cluster2ndMoments(const reco::SuperCluster &superCluster,
                                              const EcalRecHitCollection &recHits,
                                              double phiCorrectionFactor = 0.8,
                                              double w0 = 4.7,
                                              bool useLogWeights = true);
   static Cluster2ndMoments cluster2ndMoments(const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs,
                                              double phiCorrectionFactor = 0.8,
                                              double w0 = 4.7,
                                              bool useLogWeights = true);
 
   static double zernike20(const reco::BasicCluster &cluster,
                           const EcalRecHitCollection *recHits,
                           const CaloGeometry *geometry,
                           double R0 = 6.6,
                           bool logW = true,
                           float w0 = 4.7);
   static double zernike42(const reco::BasicCluster &cluster,
                           const EcalRecHitCollection *recHits,
                           const CaloGeometry *geometry,
                           double R0 = 6.6,
                           bool logW = true,
                           float w0 = 4.7);
 
   // get the detId's of a matrix centered in the maximum energy crystal = (0,0)
   // the size is specified by ixMin, ixMax, iyMin, iyMax in unit of crystals
   static std::vector<DetId> matrixDetId(const CaloTopology *topology, DetId id, CaloRectangle rectangle);
   static std::vector<DetId> matrixDetId(const CaloTopology *topology, DetId id, int size) {
     return matrixDetId(topology, id, CaloRectangle{-size, size, -size, size});
   }
 
   // get the energy deposited in a matrix centered in the maximum energy crystal = (0,0)
   // the size is specified by ixMin, ixMax, iyMin, iyMax in unit of crystals
   static float matrixEnergy(const reco::BasicCluster &cluster,
                             const EcalRecHitCollection *recHits,
                             const CaloTopology *topology,
                             DetId id,
                             CaloRectangle rectangle);
   static int matrixSize(const reco::BasicCluster &cluster,
                         const EcalRecHitCollection *recHits,
                         const CaloTopology *topology,
                         DetId id,
                         CaloRectangle rectangle);
 
   static float getFraction(const std::vector<std::pair<DetId, float>> &v_id, DetId id);
   // get the DetId and the energy of the maximum energy crystal in a vector of DetId
   static std::pair<DetId, float> getMaximum(const std::vector<std::pair<DetId, float>> &v_id,
                                             const EcalRecHitCollection *recHits);
 
   // get the energy of a DetId, return 0 if the DetId is not in the collection
   static float recHitEnergy(DetId id, const EcalRecHitCollection *recHits);
 
   //Shower shape variables return vector <Roundness, Angle> of a photon
   static std::vector<float> roundnessBarrelSuperClusters(const reco::SuperCluster &superCluster,
                                                          const EcalRecHitCollection &recHits,
                                                          int weightedPositionMethod = 0,
                                                          float energyThreshold = 0.0);
   static std::vector<float> roundnessBarrelSuperClustersUserExtended(const reco::SuperCluster &superCluster,
                                                                      const EcalRecHitCollection &recHits,
                                                                      int ieta_delta = 0,
                                                                      int iphi_delta = 0,
                                                                      float energyRHThresh = 0.00000,
                                                                      int weightedPositionMethod = 0);
   static std::vector<float> roundnessSelectedBarrelRecHits(
       const std::vector<std::pair<const EcalRecHit *, float>> &rhVector, int weightedPositionMethod = 0);
 
   //works out the number of staturated crystals in 5x5
   static int nrSaturatedCrysIn5x5(const DetId &id, const EcalRecHitCollection *recHits, const CaloTopology *topology);
 
 private:
   struct EcalClusterEnergyDeposition {
     double deposited_energy;
     double r;
     double phi;
   };
 
   static std::vector<EcalClusterEnergyDeposition> getEnergyDepTopology(const reco::BasicCluster &cluster,
                                                                        const EcalRecHitCollection *recHits,
                                                                        const CaloGeometry *geometry,
                                                                        bool logW,
                                                                        float w0);
 
   static math::XYZVector meanClusterPosition(const reco::BasicCluster &cluster,
                                              const EcalRecHitCollection *recHits,
                                              const CaloTopology *topology,
                                              const CaloGeometry *geometry);
 
   //return energy weighted mean distance from the seed crystal in number of crystals
   //<iEta,iPhi>, iPhi is not defined for endcap and is returned as zero
   static std::pair<float, float> mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster,
                                                                  const EcalRecHitCollection *recHits,
                                                                  const CaloTopology *topology);
 
   static std::pair<float, float> mean5x5PositionInXY(const reco::BasicCluster &cluster,
                                                      const EcalRecHitCollection *recHits,
                                                      const CaloTopology *topology);
 
   static double f00(double r) { return 1; }
   static double f11(double r) { return r; }
   static double f20(double r) { return 2.0 * r * r - 1.0; }
   static double f22(double r) { return r * r; }
   static double f31(double r) { return 3.0 * r * r * r - 2.0 * r; }
   static double f33(double r) { return r * r * r; }
   static double f40(double r) { return 6.0 * r * r * r * r - 6.0 * r * r + 1.0; }
   static double f42(double r) { return 4.0 * r * r * r * r - 3.0 * r * r; }
   static double f44(double r) { return r * r * r * r; }
   static double f51(double r) { return 10.0 * pow(r, 5) - 12.0 * pow(r, 3) + 3.0 * r; }
   static double f53(double r) { return 5.0 * pow(r, 5) - 4.0 * pow(r, 3); }
   static double f55(double r) { return pow(r, 5); }
 
   static double absZernikeMoment(const reco::BasicCluster &cluster,
                                  const EcalRecHitCollection *recHits,
                                  const CaloGeometry *geometry,
                                  int n,
                                  int m,
                                  double R0,
                                  bool logW,
                                  float w0);
   static double fast_AbsZernikeMoment(const reco::BasicCluster &cluster,
                                       const EcalRecHitCollection *recHits,
                                       const CaloGeometry *geometry,
                                       int n,
                                       int m,
                                       double R0,
                                       bool logW,
                                       float w0);
   static double calc_AbsZernikeMoment(const reco::BasicCluster &cluster,
                                       const EcalRecHitCollection *recHits,
                                       const CaloGeometry *geometry,
                                       int n,
                                       int m,
                                       double R0,
                                       bool logW,
                                       float w0);
 
   static double factorial(int n) {
     double res = 1.;
     for (int i = 2; i <= n; ++i)
       res *= i;
     return res;
   }
 
   //useful functions for the localCovariances function
   static float getIEta(const DetId &id);
   static float getIPhi(const DetId &id);
   static float getNormedIX(const DetId &id);
   static float getNormedIY(const DetId &id);
   static float getDPhiEndcap(const DetId &crysId, float meanX, float meanY);
   static float getNrCrysDiffInEta(const DetId &crysId, const DetId &orginId);
   static float getNrCrysDiffInPhi(const DetId &crysId, const DetId &orginId);
 
   //useful functions for showerRoundnessBarrel function
   static int deltaIEta(int seed_ieta, int rh_ieta);
   static int deltaIPhi(int seed_iphi, int rh_iphi);
   static std::vector<int> getSeedPosition(const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs);
   static float getSumEnergy(const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs);
   static float computeWeight(float eRH, float energyTotal, int weightedPositionMethod);
 };
 
 // implementation
 template <bool noZS>
 float EcalClusterToolsT<noZS>::getFraction(const std::vector<std::pair<DetId, float>> &v_id, DetId id) {
   if (noZS)
     return 1.0;
   float frac = 0.0;
   for (size_t i = 0; i < v_id.size(); ++i) {
     if (v_id[i].first.rawId() == id.rawId()) {
       frac = v_id[i].second;
       break;
     }
   }
   return frac;
 }
 
 template <bool noZS>
 std::pair<DetId, float> EcalClusterToolsT<noZS>::getMaximum(const std::vector<std::pair<DetId, float>> &v_id,
                                                             const EcalRecHitCollection *recHits) {
   float max = 0;
   DetId id(0);
   for (size_t i = 0; i < v_id.size(); ++i) {
     float energy = recHitEnergy(v_id[i].first, recHits) * (noZS ? 1.0 : v_id[i].second);
     if (energy > max) {
       max = energy;
       id = v_id[i].first;
     }
   }
   return std::pair<DetId, float>(id, max);
 }
 
 template <bool noZS>
 std::pair<DetId, float> EcalClusterToolsT<noZS>::getMaximum(const reco::BasicCluster &cluster,
                                                             const EcalRecHitCollection *recHits) {
   return getMaximum(cluster.hitsAndFractions(), recHits);
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::recHitEnergy(DetId id, const EcalRecHitCollection *recHits) {
   if (id == DetId(0)) {
     return 0;
   } else {
     EcalRecHitCollection::const_iterator it = recHits->find(id);
     if (it != recHits->end()) {
       if (noZS && (it->checkFlag(EcalRecHit::kTowerRecovered) || it->checkFlag(EcalRecHit::kWeird) ||
                    (it->detid().subdetId() == EcalBarrel && it->checkFlag(EcalRecHit::kDiWeird))))
         return 0.0;
       else
         return it->energy();
     } else {
       //throw cms::Exception("EcalRecHitNotFound") << "The recHit corresponding to the DetId" << id.rawId() << " not found in the EcalRecHitCollection";
       // the recHit is not in the collection (hopefully zero suppressed)
       return 0;
     }
   }
   return 0;
 }
 
 // Returns the energy in a rectangle of crystals
 // specified in eta by ixMin and ixMax
 //       and in phi by iyMin and iyMax
 //
 // Reference picture (X=seed crystal)
 //    iy ___________
 //     2 |_|_|_|_|_|
 //     1 |_|_|_|_|_|
 //     0 |_|_|X|_|_|
 //    -1 |_|_|_|_|_|
 //    -2 |_|_|_|_|_|
 //      -2 -1 0 1 2 ix
 template <bool noZS>
 float EcalClusterToolsT<noZS>::matrixEnergy(const reco::BasicCluster &cluster,
                                             const EcalRecHitCollection *recHits,
                                             const CaloTopology *topology,
                                             DetId id,
                                             CaloRectangle rectangle) {
   float energy = 0;
   auto const &vId = cluster.hitsAndFractions();
 
   for (auto const &detId : rectangle(id, *topology)) {
     energy += recHitEnergy(detId, recHits) * getFraction(vId, detId);
   }
 
   return energy;
 }
 
 template <bool noZS>
 int EcalClusterToolsT<noZS>::matrixSize(const reco::BasicCluster &cluster,
                                         const EcalRecHitCollection *recHits,
                                         const CaloTopology *topology,
                                         DetId id,
                                         CaloRectangle rectangle) {
   int result = 0;
 
   for (auto const &detId : rectangle(id, *topology)) {
     const float energy = recHitEnergy(detId, recHits);
     const float frac = getFraction(cluster.hitsAndFractions(), detId);
     if (energy * frac > 0)
       result++;
   }
   return result;
 }
 
 template <bool noZS>
 std::vector<DetId> EcalClusterToolsT<noZS>::matrixDetId(const CaloTopology *topology,
                                                         DetId id,
                                                         CaloRectangle rectangle) {
   std::vector<DetId> v;
   for (auto const &detId : rectangle(id, *topology)) {
     if (detId != DetId(0))
       v.push_back(detId);
   }
   return v;
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e2x2(const reco::BasicCluster &cluster,
                                     const EcalRecHitCollection *recHits,
                                     const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   std::list<float> energies;
   float max_E = matrixEnergy(cluster, recHits, topology, id, {-1, 0, -1, 0});
   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-1, 0, 0, 1}));
   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {0, 1, 0, 1}));
   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {0, 1, -1, 0}));
   return max_E;
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e3x2(const reco::BasicCluster &cluster,
                                     const EcalRecHitCollection *recHits,
                                     const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   float max_E = matrixEnergy(cluster, recHits, topology, id, {-1, 1, -1, 0});
   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {0, 1, -1, 1}));
   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-1, 1, 0, 1}));
   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-1, 0, -1, 1}));
   return max_E;
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e3x3(const reco::BasicCluster &cluster,
                                     const EcalRecHitCollection *recHits,
                                     const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {-1, 1, -1, 1});
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e4x4(const reco::BasicCluster &cluster,
                                     const EcalRecHitCollection *recHits,
                                     const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   float max_E = matrixEnergy(cluster, recHits, topology, id, {-1, 2, -2, 1});
   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-2, 1, -2, 1}));
   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-2, 1, -1, 2}));
   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-1, 2, -1, 2}));
   return max_E;
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e5x5(const reco::BasicCluster &cluster,
                                     const EcalRecHitCollection *recHits,
                                     const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {-2, 2, -2, 2});
 }
 
 template <bool noZS>
 int EcalClusterToolsT<noZS>::n5x5(const reco::BasicCluster &cluster,
                                   const EcalRecHitCollection *recHits,
                                   const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixSize(cluster, recHits, topology, id, {-2, 2, -2, 2});
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::eMax(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits) {
   return getMaximum(cluster.hitsAndFractions(), recHits).second;
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e2nd(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits) {
   std::vector<float> energies;
   const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
   energies.reserve(v_id.size());
   if (v_id.size() < 2)
     return 0;
   for (size_t i = 0; i < v_id.size(); ++i) {
     energies.push_back(recHitEnergy(v_id[i].first, recHits) * (noZS ? 1.0 : v_id[i].second));
   }
   std::partial_sort(energies.begin(), energies.begin() + 2, energies.end(), std::greater<float>());
   return energies[1];
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e2x5Right(const reco::BasicCluster &cluster,
                                          const EcalRecHitCollection *recHits,
                                          const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {1, 2, -2, 2});
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e2x5Left(const reco::BasicCluster &cluster,
                                         const EcalRecHitCollection *recHits,
                                         const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {-2, -1, -2, 2});
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e2x5Top(const reco::BasicCluster &cluster,
                                        const EcalRecHitCollection *recHits,
                                        const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {-2, 2, 1, 2});
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e2x5Bottom(const reco::BasicCluster &cluster,
                                           const EcalRecHitCollection *recHits,
                                           const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {-2, 2, -2, -1});
 }
 
 // Energy in 2x5 strip containing the max crystal.
 // Adapted from code by Sam Harper
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e2x5Max(const reco::BasicCluster &cluster,
                                        const EcalRecHitCollection *recHits,
                                        const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
 
   // 1x5 strip left of seed
   float left = matrixEnergy(cluster, recHits, topology, id, {-1, -1, -2, 2});
   // 1x5 strip right of seed
   float right = matrixEnergy(cluster, recHits, topology, id, {1, 1, -2, 2});
   // 1x5 strip containing seed
   float centre = matrixEnergy(cluster, recHits, topology, id, {0, 0, -2, 2});
 
   // Return the maximum of (left+center) or (right+center) strip
   return left > right ? left + centre : right + centre;
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e1x5(const reco::BasicCluster &cluster,
                                     const EcalRecHitCollection *recHits,
                                     const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {0, 0, -2, 2});
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e5x1(const reco::BasicCluster &cluster,
                                     const EcalRecHitCollection *recHits,
                                     const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {-2, 2, 0, 0});
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e1x3(const reco::BasicCluster &cluster,
                                     const EcalRecHitCollection *recHits,
                                     const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {0, 0, -1, 1});
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::e3x1(const reco::BasicCluster &cluster,
                                     const EcalRecHitCollection *recHits,
                                     const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {-1, 1, 0, 0});
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::eLeft(const reco::BasicCluster &cluster,
                                      const EcalRecHitCollection *recHits,
                                      const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {-1, -1, 0, 0});
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::eRight(const reco::BasicCluster &cluster,
                                       const EcalRecHitCollection *recHits,
                                       const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {1, 1, 0, 0});
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::eTop(const reco::BasicCluster &cluster,
                                     const EcalRecHitCollection *recHits,
                                     const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {0, 0, 1, 1});
 }
 
 template <bool noZS>
 float EcalClusterToolsT<noZS>::eBottom(const reco::BasicCluster &cluster,
                                        const EcalRecHitCollection *recHits,
                                        const CaloTopology *topology) {
   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
   return matrixEnergy(cluster, recHits, topology, id, {0, 0, -1, -1});
 }
 
 template <bool noZS>
 std::vector<float> EcalClusterToolsT<noZS>::energyBasketFractionEta(const reco::BasicCluster &cluster,
                                                                     const EcalRecHitCollection *recHits) {
   std::vector<float> basketFraction(2 * EBDetId::kModulesPerSM);
   float clusterEnergy = cluster.energy();
   const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
   if (v_id[0].first.subdetId() != EcalBarrel) {
     edm::LogWarning("EcalClusterToolsT<noZS>::energyBasketFractionEta")
         << "Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel "
            "basic-clusters. Returning empty vector.";
     return basketFraction;
   }
   for (size_t i = 0; i < v_id.size(); ++i) {
     basketFraction[EBDetId(v_id[i].first).im() - 1 + EBDetId(v_id[i].first).positiveZ() * EBDetId::kModulesPerSM] +=
         recHitEnergy(v_id[i].first, recHits) * v_id[i].second / clusterEnergy;
   }
   std::sort(basketFraction.rbegin(), basketFraction.rend());
   return basketFraction;
 }
 
 template <bool noZS>
 std::vector<float> EcalClusterToolsT<noZS>::energyBasketFractionPhi(const reco::BasicCluster &cluster,
                                                                     const EcalRecHitCollection *recHits) {
   std::vector<float> basketFraction(2 * (EBDetId::MAX_IPHI / EBDetId::kCrystalsInPhi));
   float clusterEnergy = cluster.energy();
   const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
   if (v_id[0].first.subdetId() != EcalBarrel) {
     edm::LogWarning("EcalClusterToolsT<noZS>::energyBasketFractionPhi")
         << "Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel "
            "basic-clusters. Returning empty vector.";
     return basketFraction;
   }
   for (size_t i = 0; i < v_id.size(); ++i) {
     basketFraction[(EBDetId(v_id[i].first).iphi() - 1) / EBDetId::kCrystalsInPhi +
                    EBDetId(v_id[i].first).positiveZ() * EBDetId::kTowersInPhi] +=
         recHitEnergy(v_id[i].first, recHits) * (noZS ? 1.0 : v_id[i].second) / clusterEnergy;
   }
   std::sort(basketFraction.rbegin(), basketFraction.rend());
   return basketFraction;
 }
 
 template <bool noZS>
 std::vector<typename EcalClusterToolsT<noZS>::EcalClusterEnergyDeposition>
 EcalClusterToolsT<noZS>::getEnergyDepTopology(const reco::BasicCluster &cluster,
                                               const EcalRecHitCollection *recHits,
                                               const CaloGeometry *geometry,
                                               bool logW,
                                               float w0) {
   std::vector<typename EcalClusterToolsT<noZS>::EcalClusterEnergyDeposition> energyDistribution;
   // init a map of the energy deposition centered on the
   // cluster centroid. This is for momenta calculation only.
   CLHEP::Hep3Vector clVect(cluster.position().x(), cluster.position().y(), cluster.position().z());
   CLHEP::Hep3Vector clDir(clVect);
   clDir *= 1.0 / clDir.mag();
   // in the transverse plane, axis perpendicular to clusterDir
   CLHEP::Hep3Vector theta_axis(clDir.y(), -clDir.x(), 0.0);
   theta_axis *= 1.0 / theta_axis.mag();
   CLHEP::Hep3Vector phi_axis = theta_axis.cross(clDir);
 
   const std::vector<std::pair<DetId, float>> &clusterDetIds = cluster.hitsAndFractions();
 
   EcalClusterEnergyDeposition clEdep;
   EcalRecHit testEcalRecHit;
   std::vector<std::pair<DetId, float>>::const_iterator posCurrent;
   // loop over crystals
   for (posCurrent = clusterDetIds.begin(); posCurrent != clusterDetIds.end(); ++posCurrent) {
     EcalRecHitCollection::const_iterator itt = recHits->find((*posCurrent).first);
     testEcalRecHit = *itt;
 
     if (((*posCurrent).first != DetId(0)) && (recHits->find((*posCurrent).first) != recHits->end())) {
       clEdep.deposited_energy = testEcalRecHit.energy() * (noZS ? 1.0 : (*posCurrent).second);
       // if logarithmic weight is requested, apply cut on minimum energy of the recHit
       if (logW) {
         //double w0 = parameterMap_.find("W0")->second;
 
         double weight = std::max(0.0, w0 + log(std::abs(clEdep.deposited_energy) / cluster.energy()));
         if (weight == 0) {
           LogDebug("ClusterShapeAlgo") << "Crystal has insufficient energy: E = " << clEdep.deposited_energy
                                        << " GeV; skipping... ";
           continue;
         } else
           LogDebug("ClusterShapeAlgo") << "===> got crystal. Energy = " << clEdep.deposited_energy << " GeV. ";
       }
       DetId id_ = (*posCurrent).first;
       const CaloSubdetectorGeometry *geo = geometry->getSubdetectorGeometry(id_);
       auto this_cell = geo->getGeometry(id_);
       const GlobalPoint &cellPos = this_cell->getPosition();
       CLHEP::Hep3Vector gblPos(cellPos.x(), cellPos.y(), cellPos.z());  //surface position?
       // Evaluate the distance from the cluster centroid
       CLHEP::Hep3Vector diff = gblPos - clVect;
       // Important: for the moment calculation, only the "lateral distance" is important
       // "lateral distance" r_i = distance of the digi position from the axis Origin-Cluster Center
       // ---> subtract the projection on clDir
       CLHEP::Hep3Vector DigiVect = diff - diff.dot(clDir) * clDir;
       clEdep.r = DigiVect.mag();
       LogDebug("ClusterShapeAlgo") << "E = " << clEdep.deposited_energy << "\tdiff = " << diff.mag()
                                    << "\tr = " << clEdep.r;
       clEdep.phi = DigiVect.angle(theta_axis);
       if (DigiVect.dot(phi_axis) < 0)
         clEdep.phi = 2 * M_PI - clEdep.phi;
       energyDistribution.push_back(clEdep);
     }
   }
   return energyDistribution;
 }
 
 template <bool noZS>
 std::vector<float> EcalClusterToolsT<noZS>::lat(const reco::BasicCluster &cluster,
                                                 const EcalRecHitCollection *recHits,
                                                 const CaloGeometry *geometry,
                                                 bool logW,
                                                 float w0) {
   std::vector<EcalClusterToolsT::EcalClusterEnergyDeposition> energyDistribution =
       getEnergyDepTopology(cluster, recHits, geometry, logW, w0);
 
   std::vector<float> lat;
   double r, redmoment = 0;
   double phiRedmoment = 0;
   double etaRedmoment = 0;
   int n, n1, n2, tmp;
   int clusterSize = energyDistribution.size();
   float etaLat_, phiLat_, lat_;
   if (clusterSize < 3) {
     etaLat_ = 0.0;
     lat_ = 0.0;
     lat.push_back(0.);
     lat.push_back(0.);
     lat.push_back(0.);
     return lat;
   }
 
   n1 = 0;
   n2 = 1;
   if (energyDistribution[1].deposited_energy > energyDistribution[0].deposited_energy) {
     tmp = n2;
     n2 = n1;
     n1 = tmp;
   }
   for (int i = 2; i < clusterSize; i++) {
     n = i;
     if (energyDistribution[i].deposited_energy > energyDistribution[n1].deposited_energy) {
       tmp = n2;
       n2 = n1;
       n1 = i;
       n = tmp;
     } else {
       if (energyDistribution[i].deposited_energy > energyDistribution[n2].deposited_energy) {
         tmp = n2;
         n2 = i;
         n = tmp;
       }
     }
 
     r = energyDistribution[n].r;
     redmoment += r * r * energyDistribution[n].deposited_energy;
     double rphi = r * cos(energyDistribution[n].phi);
     phiRedmoment += rphi * rphi * energyDistribution[n].deposited_energy;
     double reta = r * sin(energyDistribution[n].phi);
     etaRedmoment += reta * reta * energyDistribution[n].deposited_energy;
   }
   double e1 = energyDistribution[n1].deposited_energy;
   double e2 = energyDistribution[n2].deposited_energy;
 
   lat_ = redmoment / (redmoment + 2.19 * 2.19 * (e1 + e2));
   phiLat_ = phiRedmoment / (phiRedmoment + 2.19 * 2.19 * (e1 + e2));
   etaLat_ = etaRedmoment / (etaRedmoment + 2.19 * 2.19 * (e1 + e2));
 
   lat.push_back(etaLat_);
   lat.push_back(phiLat_);
   lat.push_back(lat_);
   return lat;
 }
 
 template <bool noZS>
 math::XYZVector EcalClusterToolsT<noZS>::meanClusterPosition(const reco::BasicCluster &cluster,
                                                              const EcalRecHitCollection *recHits,
                                                              const CaloTopology *topology,
                                                              const CaloGeometry *geometry) {
   // find mean energy position of a 5x5 cluster around the maximum
   math::XYZVector meanPosition(0.0, 0.0, 0.0);
   const std::vector<std::pair<DetId, float>> &hsAndFs = cluster.hitsAndFractions();
   std::vector<DetId> v_id = matrixDetId(topology, getMaximum(cluster, recHits).first, 2);
   for (const std::pair<DetId, float> &hitAndFrac : hsAndFs) {
     for (std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it) {
       if (hitAndFrac.first != *it && !noZS)
         continue;
       const CaloSubdetectorGeometry *geo = geometry->getSubdetectorGeometry(*it);
       GlobalPoint positionGP = geo->getGeometry(*it)->getPosition();
       math::XYZVector position(positionGP.x(), positionGP.y(), positionGP.z());
       meanPosition = meanPosition + recHitEnergy(*it, recHits) * position * hitAndFrac.second;
     }
     if (noZS)
       break;
   }
   return meanPosition / e5x5(cluster, recHits, topology);
 }
 
 //returns mean energy weighted eta/phi in crystals from the seed
 //iPhi is not defined for endcap and is returned as zero
 //return <eta,phi>
 //we have an issue in working out what to do for negative energies
 //I (Sam Harper) think it makes sense to ignore crystals with E<0 in the calculation as they are ignored
 //in the sigmaIEtaIEta calculation (well they arent fully ignored, they do still contribute to the e5x5 sum
 //in the sigmaIEtaIEta calculation but not here)
 template <bool noZS>
 std::pair<float, float> EcalClusterToolsT<noZS>::mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster,
                                                                                  const EcalRecHitCollection *recHits,
                                                                                  const CaloTopology *topology) {
   DetId seedId = getMaximum(cluster, recHits).first;
   float meanDEta = 0.;
   float meanDPhi = 0.;
   float energySum = 0.;
 
   const std::vector<std::pair<DetId, float>> &hsAndFs = cluster.hitsAndFractions();
   std::vector<DetId> v_id = matrixDetId(topology, seedId, 2);
   for (const std::pair<DetId, float> &hAndF : hsAndFs) {
     for (std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it) {
       if (hAndF.first != *it && !noZS)
         continue;
       float energy = recHitEnergy(*it, recHits) * hAndF.second;
       if (energy < 0.)
         continue;  //skipping negative energy crystals
       meanDEta += energy * getNrCrysDiffInEta(*it, seedId);
       meanDPhi += energy * getNrCrysDiffInPhi(*it, seedId);
       energySum += energy;
     }
     if (noZS)
       break;
   }
   meanDEta /= energySum;
   meanDPhi /= energySum;
   return std::pair<float, float>(meanDEta, meanDPhi);
 }
 
 //returns mean energy weighted x/y in normalised crystal coordinates
 //only valid for endcap, returns 0,0 for barrel
 //we have an issue in working out what to do for negative energies
 //I (Sam Harper) think it makes sense to ignore crystals with E<0 in the calculation as they are ignored
 //in the sigmaIEtaIEta calculation (well they arent fully ignored, they do still contribute to the e5x5 sum
 //in the sigmaIEtaIEta calculation but not here)
 template <bool noZS>
 std::pair<float, float> EcalClusterToolsT<noZS>::mean5x5PositionInXY(const reco::BasicCluster &cluster,
                                                                      const EcalRecHitCollection *recHits,
                                                                      const CaloTopology *topology) {
   DetId seedId = getMaximum(cluster, recHits).first;
 
   std::pair<float, float> meanXY(0., 0.);
   if (seedId.subdetId() == EcalBarrel)
     return meanXY;
 
   float energySum = 0.;
 
   const std::vector<std::pair<DetId, float>> &hsAndFs = cluster.hitsAndFractions();
   std::vector<DetId> v_id = matrixDetId(topology, seedId, 2);
   for (const std::pair<DetId, float> &hAndF : hsAndFs) {
     for (std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it) {
       if (hAndF.first != *it && !noZS)
         continue;
       float energy = recHitEnergy(*it, recHits) * hAndF.second;
       if (energy < 0.)
         continue;  //skipping negative energy crystals
       meanXY.first += energy * getNormedIX(*it);
       meanXY.second += energy * getNormedIY(*it);
       energySum += energy;
     }
     if (noZS)
       break;
   }
   meanXY.first /= energySum;
   meanXY.second /= energySum;
   return meanXY;
 }
 
 template <bool noZS>
 std::array<float, 3> EcalClusterToolsT<noZS>::covariances(const reco::BasicCluster &cluster,
                                                           const EcalRecHitCollection *recHits,
                                                           const CaloTopology *topology,
                                                           const CaloGeometry *geometry,
                                                           float w0) {
   float e_5x5 = e5x5(cluster, recHits, topology);
   float covEtaEta, covEtaPhi, covPhiPhi;
   if (e_5x5 >= 0.) {
     //double w0_ = parameterMap_.find("W0")->second;
     const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
     math::XYZVector meanPosition = meanClusterPosition(cluster, recHits, topology, geometry);
 
     // now we can calculate the covariances
     double numeratorEtaEta = 0;
     double numeratorEtaPhi = 0;
     double numeratorPhiPhi = 0;
     double denominator = 0;
 
     DetId id = getMaximum(v_id, recHits).first;
     CaloRectangle rectangle{-2, 2, -2, 2};
     for (auto const &detId : rectangle(id, *topology)) {
       float frac = getFraction(v_id, detId);
       float energy = recHitEnergy(detId, recHits) * frac;
 
       if (energy <= 0)
         continue;
 
       const CaloSubdetectorGeometry *geo = geometry->getSubdetectorGeometry(detId);
       GlobalPoint position = geo->getGeometry(detId)->getPosition();
 
       double dPhi = position.phi() - meanPosition.phi();
       if (dPhi > +Geom::pi()) {
         dPhi = Geom::twoPi() - dPhi;
       }
       if (dPhi < -Geom::pi()) {
         dPhi = Geom::twoPi() + dPhi;
       }
 
       double dEta = position.eta() - meanPosition.eta();
       double w = 0.;
       w = std::max(0.0f, w0 + std::log(energy / e_5x5));
 
       denominator += w;
       numeratorEtaEta += w * dEta * dEta;
       numeratorEtaPhi += w * dEta * dPhi;
       numeratorPhiPhi += w * dPhi * dPhi;
     }
 
     if (denominator != 0.0) {
       covEtaEta = numeratorEtaEta / denominator;
       covEtaPhi = numeratorEtaPhi / denominator;
       covPhiPhi = numeratorPhiPhi / denominator;
     } else {
       covEtaEta = 999.9;
       covEtaPhi = 999.9;
       covPhiPhi = 999.9;
     }
 
   } else {
     // Warn the user if there was no energy in the cells and return zeroes.
     //       std::cout << "\ClusterShapeAlgo::Calculate_Covariances:  no energy in supplied cells.\n";
     covEtaEta = 0;
     covEtaPhi = 0;
     covPhiPhi = 0;
   }
   std::array<float, 3> v{{covEtaEta, covEtaPhi, covPhiPhi}};
   return v;
 }
 
 //for covIEtaIEta,covIEtaIPhi and covIPhiIPhi are defined but only covIEtaIEta has been actively studied
 //instead of using absolute eta/phi it counts crystals normalised so that it gives identical results to normal covariances except near the cracks where of course its better
 //it also does not require any eta correction function in the endcap
 //it is multipled by an approprate crystal size to ensure it gives similar values to covariances(...)
 template <bool noZS>
 std::array<float, 3> EcalClusterToolsT<noZS>::localCovariances(const reco::BasicCluster &cluster,
                                                                const EcalRecHitCollection *recHits,
                                                                const CaloTopology *topology,
                                                                float w0,
                                                                const EcalPFRecHitThresholds *thresholds,
                                                                float multEB,
                                                                float multEE) {
   float e_5x5 = e5x5(cluster, recHits, topology);
   float covEtaEta, covEtaPhi, covPhiPhi;
 
   if (e_5x5 >= 0.) {
     //double w0_ = parameterMap_.find("W0")->second;
     const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
     std::pair<float, float> mean5x5PosInNrCrysFromSeed = mean5x5PositionInLocalCrysCoord(cluster, recHits, topology);
     std::pair<float, float> mean5x5XYPos = mean5x5PositionInXY(cluster, recHits, topology);
 
     // now we can calculate the covariances
     double numeratorEtaEta = 0;
     double numeratorEtaPhi = 0;
     double numeratorPhiPhi = 0;
     double denominator = 0;
 
     //these allow us to scale the localCov by the crystal size
     //so that the localCovs have the same average value as the normal covs
     const double barrelCrysSize = 0.01745;  //approximate size of crystal in eta,phi in barrel
     const double endcapCrysSize = 0.0447;   //the approximate crystal size sigmaEtaEta was corrected to in the endcap
 
     DetId seedId = getMaximum(v_id, recHits).first;
 
     bool isBarrel = seedId.subdetId() == EcalBarrel;
     const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
     float mult = isBarrel ? multEB : multEE;  // we will multiply PF RecHit threshold by mult.
                                               // mult = 1 should work reasonably well for noise cleaning.
     // dedicated studies showed mult=1.25 works best for Run3 in endcap, where noise is high.
     // If no noise cleaning is intended then put mult=0.
 
     CaloRectangle rectangle{-2, 2, -2, 2};
     for (auto const &detId : rectangle(seedId, *topology)) {
       float frac = getFraction(v_id, detId);
       float energy = recHitEnergy(detId, recHits) * frac;
 
       if ((thresholds == nullptr) && (mult != 0.0)) {
         throw cms::Exception("EmptyPFRechHitThresColl")
             << "In EcalClusterTools::localCovariances, if EcalPFRecHitThresholds==nulptr, then multEB and multEE "
                "should be 0 as well.";
       }
       float rhThres = 0.0;
       if (thresholds != nullptr) {
         rhThres = (*thresholds)[detId];  // access PFRechit thresholds for noise cleaning
       }
 
       if (energy <= (rhThres * mult))
         continue;
 
       float dEta = getNrCrysDiffInEta(detId, seedId) - mean5x5PosInNrCrysFromSeed.first;
       float dPhi = 0;
 
       if (isBarrel)
         dPhi = getNrCrysDiffInPhi(detId, seedId) - mean5x5PosInNrCrysFromSeed.second;
       else
         dPhi = getDPhiEndcap(detId, mean5x5XYPos.first, mean5x5XYPos.second);
 
       double w = std::max(0.0f, w0 + std::log(energy / e_5x5));
 
       denominator += w;
       numeratorEtaEta += w * dEta * dEta;
       numeratorEtaPhi += w * dEta * dPhi;
       numeratorPhiPhi += w * dPhi * dPhi;
     }
 
     //multiplying by crysSize to make the values compariable to normal covariances
     if (denominator != 0.0) {
       covEtaEta = crysSize * crysSize * numeratorEtaEta / denominator;
       covEtaPhi = crysSize * crysSize * numeratorEtaPhi / denominator;
       covPhiPhi = crysSize * crysSize * numeratorPhiPhi / denominator;
     } else {
       covEtaEta = 999.9;
       covEtaPhi = 999.9;
       covPhiPhi = 999.9;
     }
 
   } else {
     // Warn the user if there was no energy in the cells and return zeroes.
     //       std::cout << "\ClusterShapeAlgo::Calculate_Covariances:  no energy in supplied cells.\n";
     covEtaEta = 0;
     covEtaPhi = 0;
     covPhiPhi = 0;
   }
   std::array<float, 3> v{{covEtaEta, covEtaPhi, covPhiPhi}};
   return v;
 }
 
 template <bool noZS>
 double EcalClusterToolsT<noZS>::zernike20(const reco::BasicCluster &cluster,
                                           const EcalRecHitCollection *recHits,
                                           const CaloGeometry *geometry,
                                           double R0,
                                           bool logW,
                                           float w0) {
   return absZernikeMoment(cluster, recHits, geometry, 2, 0, R0, logW, w0);
 }
 
 template <bool noZS>
 double EcalClusterToolsT<noZS>::zernike42(const reco::BasicCluster &cluster,
                                           const EcalRecHitCollection *recHits,
                                           const CaloGeometry *geometry,
                                           double R0,
                                           bool logW,
                                           float w0) {
   return absZernikeMoment(cluster, recHits, geometry, 4, 2, R0, logW, w0);
 }
 
 template <bool noZS>
 double EcalClusterToolsT<noZS>::absZernikeMoment(const reco::BasicCluster &cluster,
                                                  const EcalRecHitCollection *recHits,
                                                  const CaloGeometry *geometry,
                                                  int n,
                                                  int m,
                                                  double R0,
                                                  bool logW,
                                                  float w0) {
   // 1. Check if n,m are correctly
   if ((m > n) || ((n - m) % 2 != 0) || (n < 0) || (m < 0))
     return -1;
 
   // 2. Check if n,R0 are within validity Range :
   // n>20 or R0<2.19cm  just makes no sense !
   if ((n > 20) || (R0 <= 2.19))
     return -1;
   if (n <= 5)
     return fast_AbsZernikeMoment(cluster, recHits, geometry, n, m, R0, logW, w0);
   else
     return calc_AbsZernikeMoment(cluster, recHits, geometry, n, m, R0, logW, w0);
 }
 
 template <bool noZS>
 double EcalClusterToolsT<noZS>::fast_AbsZernikeMoment(const reco::BasicCluster &cluster,
                                                       const EcalRecHitCollection *recHits,
                                                       const CaloGeometry *geometry,
                                                       int n,
                                                       int m,
                                                       double R0,
                                                       bool logW,
                                                       float w0) {
   double r, ph, e, Re = 0, Im = 0;
   double TotalEnergy = cluster.energy();
   int index = (n / 2) * (n / 2) + (n / 2) + m;
   std::vector<EcalClusterEnergyDeposition> energyDistribution =
       getEnergyDepTopology(cluster, recHits, geometry, logW, w0);
   int clusterSize = energyDistribution.size();
   if (clusterSize < 3)
     return 0.0;
 
   for (int i = 0; i < clusterSize; i++) {
     r = energyDistribution[i].r / R0;
     if (r < 1) {
       std::vector<double> pol;
       pol.push_back(f00(r));
       pol.push_back(f11(r));
       pol.push_back(f20(r));
       pol.push_back(f22(r));
       pol.push_back(f31(r));
       pol.push_back(f33(r));
       pol.push_back(f40(r));
       pol.push_back(f42(r));
       pol.push_back(f44(r));
       pol.push_back(f51(r));
       pol.push_back(f53(r));
       pol.push_back(f55(r));
       ph = (energyDistribution[i]).phi;
       e = energyDistribution[i].deposited_energy;
       Re = Re + e / TotalEnergy * pol[index] * cos((double)m * ph);
       Im = Im - e / TotalEnergy * pol[index] * sin((double)m * ph);
     }
   }
   return sqrt(Re * Re + Im * Im);
 }
 
 template <bool noZS>
 double EcalClusterToolsT<noZS>::calc_AbsZernikeMoment(const reco::BasicCluster &cluster,
                                                       const EcalRecHitCollection *recHits,
                                                       const CaloGeometry *geometry,
                                                       int n,
                                                       int m,
                                                       double R0,
                                                       bool logW,
                                                       float w0) {
   double r, ph, e, Re = 0, Im = 0, f_nm;
   double TotalEnergy = cluster.energy();
   std::vector<EcalClusterEnergyDeposition> energyDistribution =
       getEnergyDepTopology(cluster, recHits, geometry, logW, w0);
   int clusterSize = energyDistribution.size();
   if (clusterSize < 3)
     return 0.0;
 
   for (int i = 0; i < clusterSize; ++i) {
     r = energyDistribution[i].r / R0;
     if (r < 1) {
       ph = energyDistribution[i].phi;
       e = energyDistribution[i].deposited_energy;
       f_nm = 0;
       for (int s = 0; s <= (n - m) / 2; s++) {
         if (s % 2 == 0) {
           f_nm = f_nm + factorial(n - s) * pow(r, (double)(n - 2 * s)) /
                             (factorial(s) * factorial((n + m) / 2 - s) * factorial((n - m) / 2 - s));
         } else {
           f_nm = f_nm - factorial(n - s) * pow(r, (double)(n - 2 * s)) /
                             (factorial(s) * factorial((n + m) / 2 - s) * factorial((n - m) / 2 - s));
         }
       }
       Re = Re + e / TotalEnergy * f_nm * cos((double)m * ph);
       Im = Im - e / TotalEnergy * f_nm * sin((double)m * ph);
     }
   }
   return sqrt(Re * Re + Im * Im);
 }
 
 //returns the crystal 'eta' from the det id
 //it is defined as the number of crystals from the centre in the eta direction
 //for the barrel with its eta/phi geometry it is always integer
 //for the endcap it is fractional due to the x/y geometry
 template <bool noZS>
 float EcalClusterToolsT<noZS>::getIEta(const DetId &id) {
   if (id.det() == DetId::Ecal) {
     if (id.subdetId() == EcalBarrel) {
       EBDetId ebId(id);
       return ebId.ieta();
     } else if (id.subdetId() == EcalEndcap) {
       float iXNorm = getNormedIX(id);
       float iYNorm = getNormedIY(id);
 
       return std::sqrt(iXNorm * iXNorm + iYNorm * iYNorm);
     }
   }
   return 0.;
 }
 
 //returns the crystal 'phi' from the det id
 //it is defined as the number of crystals from the centre in the phi direction
 //for the barrel with its eta/phi geometry it is always integer
 //for the endcap it is not defined
 template <bool noZS>
 float EcalClusterToolsT<noZS>::getIPhi(const DetId &id) {
   if (id.det() == DetId::Ecal) {
     if (id.subdetId() == EcalBarrel) {
       EBDetId ebId(id);
       return ebId.iphi();
     }
   }
   return 0.;
 }
 
 //want to map 1=-50,50=-1,51=1 and 100 to 50 so sub off one if zero or neg
 template <bool noZS>
 float EcalClusterToolsT<noZS>::getNormedIX(const DetId &id) {
   if (id.det() == DetId::Ecal && id.subdetId() == EcalEndcap) {
     EEDetId eeId(id);
     int iXNorm = eeId.ix() - 50;
     if (iXNorm <= 0)
       iXNorm--;
     return iXNorm;
   }
   return 0;
 }
 
 //want to map 1=-50,50=-1,51=1 and 100 to 50 so sub off one if zero or neg
 template <bool noZS>
 float EcalClusterToolsT<noZS>::getNormedIY(const DetId &id) {
   if (id.det() == DetId::Ecal && id.subdetId() == EcalEndcap) {
     EEDetId eeId(id);
     int iYNorm = eeId.iy() - 50;
     if (iYNorm <= 0)
       iYNorm--;
     return iYNorm;
   }
   return 0;
 }
 
 //nr crystals crysId is away from orgin id in eta
 template <bool noZS>
 float EcalClusterToolsT<noZS>::getNrCrysDiffInEta(const DetId &crysId, const DetId &orginId) {
   float crysIEta = getIEta(crysId);
   float orginIEta = getIEta(orginId);
   bool isBarrel = orginId.subdetId() == EcalBarrel;
 
   float nrCrysDiff = crysIEta - orginIEta;
 
   //no iEta=0 in barrel, so if go from positive to negative
   //need to reduce abs(detEta) by 1
   if (isBarrel) {
     if (crysIEta * orginIEta < 0) {  // -1 to 1 transition
       if (crysIEta > 0)
         nrCrysDiff--;
       else
         nrCrysDiff++;
     }
   }
   return nrCrysDiff;
 }
 
 //nr crystals crysId is away from orgin id in phi
 template <bool noZS>
 float EcalClusterToolsT<noZS>::getNrCrysDiffInPhi(const DetId &crysId, const DetId &orginId) {
   float crysIPhi = getIPhi(crysId);
   float orginIPhi = getIPhi(orginId);
   bool isBarrel = orginId.subdetId() == EcalBarrel;
 
   float nrCrysDiff = crysIPhi - orginIPhi;
 
   if (isBarrel) {  //if barrel, need to map into 0-180
     if (nrCrysDiff > +180) {
       nrCrysDiff = nrCrysDiff - 360;
     }
     if (nrCrysDiff < -180) {
       nrCrysDiff = nrCrysDiff + 360;
     }
   }
   return nrCrysDiff;
 }
 
 //nr crystals crysId is away from mean phi in 5x5 in phi
 template <bool noZS>
 float EcalClusterToolsT<noZS>::getDPhiEndcap(const DetId &crysId, float meanX, float meanY) {
   float iXNorm = getNormedIX(crysId);
   float iYNorm = getNormedIY(crysId);
 
   float hitLocalR2 = (iXNorm - meanX) * (iXNorm - meanX) + (iYNorm - meanY) * (iYNorm - meanY);
   float hitR2 = iXNorm * iXNorm + iYNorm * iYNorm;
   float meanR2 = meanX * meanX + meanY * meanY;
   float hitR = sqrt(hitR2);
   float meanR = sqrt(meanR2);
 
   float tmp = (hitR2 + meanR2 - hitLocalR2) / (2 * hitR * meanR);
   if (tmp < -1)
     tmp = -1;
   if (tmp > 1)
     tmp = 1;
   float phi = acos(tmp);
   float dPhi = hitR * phi;
 
   return dPhi;
 }
 
 template <bool noZS>
 std::array<float, 3> EcalClusterToolsT<noZS>::scLocalCovariances(const reco::SuperCluster &cluster,
                                                                  const EcalRecHitCollection *recHits,
                                                                  const CaloTopology *topology,
                                                                  float w0) {
   const reco::BasicCluster bcluster = *(cluster.seed());
 
   float e_5x5 = e5x5(bcluster, recHits, topology);
   float covEtaEta, covEtaPhi, covPhiPhi;
 
   if (e_5x5 >= 0.) {
     const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
     std::pair<float, float> mean5x5PosInNrCrysFromSeed = mean5x5PositionInLocalCrysCoord(bcluster, recHits, topology);
     std::pair<float, float> mean5x5XYPos = mean5x5PositionInXY(cluster, recHits, topology);
     // now we can calculate the covariances
     double numeratorEtaEta = 0;
     double numeratorEtaPhi = 0;
     double numeratorPhiPhi = 0;
     double denominator = 0;
 
     const double barrelCrysSize = 0.01745;  //approximate size of crystal in eta,phi in barrel
     const double endcapCrysSize = 0.0447;   //the approximate crystal size sigmaEtaEta was corrected to in the endcap
 
     DetId seedId = getMaximum(v_id, recHits).first;
     bool isBarrel = seedId.subdetId() == EcalBarrel;
 
     const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
 
     for (size_t i = 0; i < v_id.size(); ++i) {
       float frac = getFraction(v_id, v_id[i].first);
       float energy = recHitEnergy(v_id[i].first, recHits) * frac;
 
       if (energy <= 0)
         continue;
 
       float dEta = getNrCrysDiffInEta(v_id[i].first, seedId) - mean5x5PosInNrCrysFromSeed.first;
       float dPhi = 0;
       if (isBarrel)
         dPhi = getNrCrysDiffInPhi(v_id[i].first, seedId) - mean5x5PosInNrCrysFromSeed.second;
       else
         dPhi = getDPhiEndcap(v_id[i].first, mean5x5XYPos.first, mean5x5XYPos.second);
 
       double w = 0.;
       w = std::max(0.0f, w0 + std::log(energy / e_5x5));
 
       denominator += w;
       numeratorEtaEta += w * dEta * dEta;
       numeratorEtaPhi += w * dEta * dPhi;
       numeratorPhiPhi += w * dPhi * dPhi;
     }
 
     //multiplying by crysSize to make the values compariable to normal covariances
     if (denominator != 0.0) {
       covEtaEta = crysSize * crysSize * numeratorEtaEta / denominator;
       covEtaPhi = crysSize * crysSize * numeratorEtaPhi / denominator;
       covPhiPhi = crysSize * crysSize * numeratorPhiPhi / denominator;
     } else {
       covEtaEta = 999.9;
       covEtaPhi = 999.9;
       covPhiPhi = 999.9;
     }
 
   } else {
     // Warn the user if there was no energy in the cells and return zeroes.
     // std::cout << "\ClusterShapeAlgo::Calculate_Covariances:  no energy in supplied cells.\n";
     covEtaEta = 0;
     covEtaPhi = 0;
     covPhiPhi = 0;
   }
 
   std::array<float, 3> v{{covEtaEta, covEtaPhi, covPhiPhi}};
   return v;
 }
 
 // compute cluster second moments with respect to principal axes (eigenvectors of sEtaEta, sPhiPhi, sEtaPhi matrix)
 // store also angle alpha between major axis and phi.
 // takes into account shower elongation in phi direction due to magnetic field effect:
 // default value of 0.8 ensures sMaj = sMin for unconverted photons
 // (if phiCorrectionFactor=1 sMaj > sMin and alpha=0 also for unconverted photons)
 template <bool noZS>
 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments(const reco::BasicCluster &basicCluster,
                                                              const EcalRecHitCollection &recHits,
                                                              double phiCorrectionFactor,
                                                              double w0,
                                                              bool useLogWeights) {
   if (recHits.empty())
     return Cluster2ndMoments();
 
   std::vector<std::pair<const EcalRecHit *, float>> RH_ptrs_fracs;
 
   const std::vector<std::pair<DetId, float>> &myHitsPair = basicCluster.hitsAndFractions();
 
   for (unsigned int i = 0; i < myHitsPair.size(); i++) {
     //get pointer to recHit object
     EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
     if (myRH != recHits.end()) {
       RH_ptrs_fracs.push_back(std::make_pair(&(*myRH), myHitsPair[i].second));
     }
   }
 
   return EcalClusterToolsT<noZS>::cluster2ndMoments(RH_ptrs_fracs, phiCorrectionFactor, w0, useLogWeights);
 }
 
 template <bool noZS>
 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments(const reco::SuperCluster &superCluster,
                                                              const EcalRecHitCollection &recHits,
                                                              double phiCorrectionFactor,
                                                              double w0,
                                                              bool useLogWeights) {
   return EcalClusterToolsT<noZS>::cluster2ndMoments(
       *(superCluster.seed()), recHits, phiCorrectionFactor, w0, useLogWeights);
 }
 
 template <bool noZS>
 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments(
     const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs,
     double phiCorrectionFactor,
     double w0,
     bool useLogWeights) {
   if (RH_ptrs_fracs.empty())
     return Cluster2ndMoments();
 
   double mid_eta(0), mid_phi(0), mid_x(0), mid_y(0);
 
   double Etot = EcalClusterToolsT<noZS>::getSumEnergy(RH_ptrs_fracs);
 
   double max_phi = -10.;
   double min_phi = 100.;
 
   std::vector<double> etaDetId;
   std::vector<double> phiDetId;
   std::vector<double> xDetId;
   std::vector<double> yDetId;
   std::vector<double> wiDetId;
 
   unsigned int nCry = 0;
   double denominator = 0.;
   bool isBarrel(true);
 
   // loop over rechits and compute weights:
   for (std::vector<std::pair<const EcalRecHit *, float>>::const_iterator rhf_ptr = RH_ptrs_fracs.begin();
        rhf_ptr != RH_ptrs_fracs.end();
        rhf_ptr++) {
     const EcalRecHit *rh_ptr = rhf_ptr->first;
 
     //get iEta, iPhi
     double temp_eta(0), temp_phi(0), temp_x(0), temp_y(0);
     isBarrel = rh_ptr->detid().subdetId() == EcalBarrel;
 
     if (isBarrel) {
       temp_eta = (getIEta(rh_ptr->detid()) > 0. ? getIEta(rh_ptr->detid()) + 84.5 : getIEta(rh_ptr->detid()) + 85.5);
       temp_phi = getIPhi(rh_ptr->detid()) - 0.5;
     } else {
       temp_eta = getIEta(rh_ptr->detid());
       temp_x = getNormedIX(rh_ptr->detid());
       temp_y = getNormedIY(rh_ptr->detid());
     }
 
     double temp_ene = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
 
     double temp_wi = ((useLogWeights) ? std::max(0.0, w0 + std::log(std::abs(temp_ene) / Etot)) : temp_ene);
 
     if (temp_phi > max_phi)
       max_phi = temp_phi;
     if (temp_phi < min_phi)
       min_phi = temp_phi;
     etaDetId.push_back(temp_eta);
     phiDetId.push_back(temp_phi);
     xDetId.push_back(temp_x);
     yDetId.push_back(temp_y);
     wiDetId.push_back(temp_wi);
     denominator += temp_wi;
     nCry++;
   }
 
   if (isBarrel) {
     // correct phi wrap-around:
     if (max_phi == 359.5 && min_phi == 0.5) {
       for (unsigned int i = 0; i < nCry; i++) {
         if (phiDetId[i] - 179. > 0.)
           phiDetId[i] -= 360.;
         mid_phi += phiDetId[i] * wiDetId[i];
         mid_eta += etaDetId[i] * wiDetId[i];
       }
     } else {
       for (unsigned int i = 0; i < nCry; i++) {
         mid_phi += phiDetId[i] * wiDetId[i];
         mid_eta += etaDetId[i] * wiDetId[i];
       }
     }
   } else {
     for (unsigned int i = 0; i < nCry; i++) {
       mid_eta += etaDetId[i] * wiDetId[i];
       mid_x += xDetId[i] * wiDetId[i];
       mid_y += yDetId[i] * wiDetId[i];
     }
   }
 
   mid_eta /= denominator;
   mid_phi /= denominator;
   mid_x /= denominator;
   mid_y /= denominator;
 
   // See = sigma eta eta
   // Spp = (B field corrected) sigma phi phi
   // See = (B field corrected) sigma eta phi
   double See = 0.;
   double Spp = 0.;
   double Sep = 0.;
   double deta(0), dphi(0);
   // compute (phi-corrected) covariance matrix:
   for (unsigned int i = 0; i < nCry; i++) {
     if (isBarrel) {
       deta = etaDetId[i] - mid_eta;
       dphi = phiDetId[i] - mid_phi;
     } else {
       deta = etaDetId[i] - mid_eta;
       float hitLocalR2 = (xDetId[i] - mid_x) * (xDetId[i] - mid_x) + (yDetId[i] - mid_y) * (yDetId[i] - mid_y);
       float hitR2 = xDetId[i] * xDetId[i] + yDetId[i] * yDetId[i];
       float meanR2 = mid_x * mid_x + mid_y * mid_y;
       float hitR = sqrt(hitR2);
       float meanR = sqrt(meanR2);
       float phi = acos((hitR2 + meanR2 - hitLocalR2) / (2 * hitR * meanR));
       dphi = hitR * phi;
     }
     See += (wiDetId[i] * deta * deta) / denominator;
     Spp += phiCorrectionFactor * (wiDetId[i] * dphi * dphi) / denominator;
     Sep += sqrt(phiCorrectionFactor) * (wiDetId[i] * deta * dphi) / denominator;
   }
 
   Cluster2ndMoments returnMoments;
 
   // compute matrix eigenvalues:
   returnMoments.sMaj = ((See + Spp) + sqrt((See - Spp) * (See - Spp) + 4. * Sep * Sep)) / 2.;
   returnMoments.sMin = ((See + Spp) - sqrt((See - Spp) * (See - Spp) + 4. * Sep * Sep)) / 2.;
 
   returnMoments.alpha = atan((See - Spp + sqrt((Spp - See) * (Spp - See) + 4. * Sep * Sep)) / (2. * Sep));
 
   return returnMoments;
 }
 
 //compute shower shapes: roundness and angle in a vector. Roundness is 0th element, Angle is 1st element.
 //description: uses classical mechanics inertia tensor.
 //             roundness is smaller_eValue/larger_eValue
 //             angle is the angle from the iEta axis to the smallest eVector (a.k.a. the shower's elongated axis)
 // this function uses only recHits belonging to a SC above energyThreshold (default 0)
 // you can select linear weighting = energy_recHit/total_energy         (weightedPositionMethod=0) default
 // or log weighting = max( 0.0, 4.2 + log(energy_recHit/total_energy) ) (weightedPositionMethod=1)
 template <bool noZS>
 std::vector<float> EcalClusterToolsT<noZS>::roundnessBarrelSuperClusters(
     const reco::SuperCluster &superCluster,
     const EcalRecHitCollection &recHits,
     int weightedPositionMethod,
     float energyThreshold) {  //int positionWeightingMethod=0){
   std::vector<std::pair<const EcalRecHit *, float>> RH_ptrs_fracs;
   const std::vector<std::pair<DetId, float>> &myHitsPair = superCluster.hitsAndFractions();
   for (unsigned int i = 0; i < myHitsPair.size(); ++i) {
     //get pointer to recHit object
     EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
     if (myRH != recHits.end() && myRH->energy() * (noZS ? 1.0 : myHitsPair[i].second) > energyThreshold) {
       //require rec hit to have positive energy
       RH_ptrs_fracs.push_back(std::make_pair(&(*myRH), myHitsPair[i].second));
     }
   }
   std::vector<float> temp =
       EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(RH_ptrs_fracs, weightedPositionMethod);
   return temp;
 }
 
 // this function uses all recHits within specified window ( with default values ieta_delta=2, iphi_delta=5) around SC's highest recHit.
 // recHits must pass an energy threshold "energyRHThresh" (default 0)
 // you can select linear weighting = energy_recHit/total_energy         (weightedPositionMethod=0)
 // or log weighting = max( 0.0, 4.2 + log(energy_recHit/total_energy) ) (weightedPositionMethod=1)
 template <bool noZS>
 std::vector<float> EcalClusterToolsT<noZS>::roundnessBarrelSuperClustersUserExtended(
     const reco::SuperCluster &superCluster,
     const EcalRecHitCollection &recHits,
     int ieta_delta,
     int iphi_delta,
     float energyRHThresh,
     int weightedPositionMethod) {
   std::vector<std::pair<const EcalRecHit *, float>> RH_ptrs_fracs;
   const std::vector<std::pair<DetId, float>> &myHitsPair = superCluster.hitsAndFractions();
   for (unsigned int i = 0; i < myHitsPair.size(); ++i) {
     //get pointer to recHit object
     EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
     if (myRH != recHits.end() && myRH->energy() * (noZS ? 1.0 : myHitsPair[i].second) > energyRHThresh)
       RH_ptrs_fracs.push_back(std::make_pair(&(*myRH), myHitsPair[i].second));
   }
 
   std::vector<int> seedPosition = EcalClusterToolsT<noZS>::getSeedPosition(RH_ptrs_fracs);
 
   for (EcalRecHitCollection::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++) {
     EBDetId EBdetIdi(rh->detid());
     float the_fraction = 0;
     //if(rh != recHits.end())
     bool inEtaWindow = (abs(deltaIEta(seedPosition[0], EBdetIdi.ieta())) <= ieta_delta);
     bool inPhiWindow = (abs(deltaIPhi(seedPosition[1], EBdetIdi.iphi())) <= iphi_delta);
     bool passEThresh = (rh->energy() > energyRHThresh);
     bool alreadyCounted = false;
 
     // figure out if the rechit considered now is already inside the SC
     bool is_SCrh_inside_recHits = false;
     for (unsigned int i = 0; i < myHitsPair.size(); i++) {
       EcalRecHitCollection::const_iterator SCrh = recHits.find(myHitsPair[i].first);
       if (SCrh != recHits.end()) {
         the_fraction = myHitsPair[i].second;
         is_SCrh_inside_recHits = true;
         if (rh->detid() == SCrh->detid())
           alreadyCounted = true;
       }
     }  //for loop over SC's recHits
 
     if (is_SCrh_inside_recHits && !alreadyCounted && passEThresh && inEtaWindow && inPhiWindow) {
       RH_ptrs_fracs.push_back(std::make_pair(&(*rh), the_fraction));
     }
 
   }  //for loop over rh
   return EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(RH_ptrs_fracs, weightedPositionMethod);
 }
 
 // this function computes the roundness and angle variables for vector of pointers to recHits you pass it
 // you can select linear weighting = energy_recHit/total_energy         (weightedPositionMethod=0)
 // or log weighting = max( 0.0, 4.2 + log(energy_recHit/total_energy) ) (weightedPositionMethod=1)
 template <bool noZS>
 std::vector<float> EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(
     const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs,
     int weightedPositionMethod) {  //int weightedPositionMethod = 0){
   //positionWeightingMethod = 0 linear weighting, 1 log energy weighting
 
   std::vector<float> shapes;  // this is the returning vector
 
   //make sure photon has more than one crystal; else roundness and angle suck
   if (RH_ptrs_fracs.size() < 2) {
     shapes.push_back(-3);
     shapes.push_back(-3);
     return shapes;
   }
 
   //Find highest E RH (Seed) and save info, compute sum total energy used
   std::vector<int> seedPosition = EcalClusterToolsT<noZS>::getSeedPosition(RH_ptrs_fracs);  // *recHits);
   int tempInt = seedPosition[0];
   if (tempInt < 0)
     tempInt++;
   float energyTotal = EcalClusterToolsT<noZS>::getSumEnergy(RH_ptrs_fracs);
 
   //1st loop over rechits: compute new weighted center position in coordinates centered on seed
   float centerIEta = 0.;
   float centerIPhi = 0.;
   float denominator = 0.;
 
   for (std::vector<std::pair<const EcalRecHit *, float>>::const_iterator rhf_ptr = RH_ptrs_fracs.begin();
        rhf_ptr != RH_ptrs_fracs.end();
        rhf_ptr++) {
     const EcalRecHit *rh_ptr = rhf_ptr->first;
     //get iEta, iPhi
     EBDetId EBdetIdi(rh_ptr->detid());
     if (fabs(energyTotal) < 0.0001) {
       // don't /0, bad!
       shapes.push_back(-2);
       shapes.push_back(-2);
       return shapes;
     }
     float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
     float weight = 0;
     if (std::abs(weightedPositionMethod) < 0.0001) {  //linear
       weight = rh_energy / energyTotal;
     } else {  //logrithmic
       weight = std::max(0.0, 4.2 + log(rh_energy / energyTotal));
     }
     denominator += weight;
     centerIEta += weight * deltaIEta(seedPosition[0], EBdetIdi.ieta());
     centerIPhi += weight * deltaIPhi(seedPosition[1], EBdetIdi.iphi());
   }
   if (fabs(denominator) < 0.0001) {
     // don't /0, bad!
     shapes.push_back(-2);
     shapes.push_back(-2);
     return shapes;
   }
   centerIEta = centerIEta / denominator;
   centerIPhi = centerIPhi / denominator;
 
   //2nd loop over rechits: compute inertia tensor
   TMatrixDSym inertia(2);  //initialize 2d inertia tensor
   double inertia00 = 0.;
   double inertia01 = 0.;  // = inertia10 b/c matrix should be symmetric
   double inertia11 = 0.;
   for (std::vector<std::pair<const EcalRecHit *, float>>::const_iterator rhf_ptr = RH_ptrs_fracs.begin();
        rhf_ptr != RH_ptrs_fracs.end();
        rhf_ptr++) {
     const EcalRecHit *rh_ptr = rhf_ptr->first;
     //get iEta, iPhi
     EBDetId EBdetIdi(rh_ptr->detid());
 
     if (fabs(energyTotal) < 0.0001) {
       // don't /0, bad!
       shapes.push_back(-2);
       shapes.push_back(-2);
       return shapes;
     }
     float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
     float weight = 0;
     if (std::abs(weightedPositionMethod) < 0.0001) {  //linear
       weight = rh_energy / energyTotal;
     } else {  //logrithmic
       weight = std::max(0.0, 4.2 + log(rh_energy / energyTotal));
     }
 
     float ieta_rh_to_center = deltaIEta(seedPosition[0], EBdetIdi.ieta()) - centerIEta;
     float iphi_rh_to_center = deltaIPhi(seedPosition[1], EBdetIdi.iphi()) - centerIPhi;
 
     inertia00 += weight * iphi_rh_to_center * iphi_rh_to_center;
     inertia01 -= weight * iphi_rh_to_center * ieta_rh_to_center;
     inertia11 += weight * ieta_rh_to_center * ieta_rh_to_center;
   }
 
   inertia[0][0] = inertia00;
   inertia[0][1] = inertia01;  // use same number here
   inertia[1][0] = inertia01;  // and here to insure symmetry
   inertia[1][1] = inertia11;
 
   //step 1 find principal axes of inertia
   TMatrixD eVectors(2, 2);
   TVectorD eValues(2);
   //std::cout<<"EcalClusterToolsT<noZS>::showerRoundness- about to compute eVectors"<<std::endl;
   eVectors = inertia.EigenVectors(eValues);  //ordered highest eV to lowest eV (I checked!)
   //and eVectors are in columns of matrix! I checked!
   //and they are normalized to 1
 
   //step 2 select eta component of smaller eVal's eVector
   TVectorD smallerAxis(2);          //easiest to spin SC on this axis (smallest eVal)
   smallerAxis[0] = eVectors[0][1];  //row,col  //eta component
   smallerAxis[1] = eVectors[1][1];  //phi component
 
   //step 3 compute interesting quatities
   Double_t temp = fabs(smallerAxis[0]);  // closer to 1 ->beamhalo, closer to 0 something else
   if (fabs(eValues[0]) < 0.0001) {
     // don't /0, bad!
     shapes.push_back(-2);
     shapes.push_back(-2);
     return shapes;
   }
 
   float Roundness = eValues[1] / eValues[0];
   float Angle = acos(temp);
 
   if (-0.00001 < Roundness && Roundness < 0)
     Roundness = 0.;
   if (-0.00001 < Angle && Angle < 0)
     Angle = 0.;
 
   shapes.push_back(Roundness);
   shapes.push_back(Angle);
   return shapes;
 }
 
 template <bool noZS>
 int EcalClusterToolsT<noZS>::nrSaturatedCrysIn5x5(const DetId &id,
                                                   const EcalRecHitCollection *recHits,
                                                   const CaloTopology *topology) {
   int nrSat = 0;
   CaloRectangle rectangle{-2, 2, -2, 2};
   for (auto const &detId : rectangle(id, *topology)) {
     auto recHitIt = recHits->find(detId);
     if (recHitIt != recHits->end() && recHitIt->checkFlag(EcalRecHit::kSaturated)) {
       nrSat++;
     }
   }
   return nrSat;
 }
 
 //private functions useful for roundnessBarrelSuperClusters etc.
 //compute delta iphi between a seed and a particular recHit
 //iphi [1,360]
 //safe gaurds are put in to ensure the difference is between [-180,180]
 template <bool noZS>
 int EcalClusterToolsT<noZS>::deltaIPhi(int seed_iphi, int rh_iphi) {
   int rel_iphi = rh_iphi - seed_iphi;
   // take care of cyclic variable iphi [1,360]
   if (rel_iphi > 180)
     rel_iphi = rel_iphi - 360;
   if (rel_iphi < -180)
     rel_iphi = rel_iphi + 360;
   return rel_iphi;
 }
 
 //compute delta ieta between a seed and a particular recHit
 //ieta [-85,-1] and [1,85]
 //safe gaurds are put in to shift the negative ieta +1 to make an ieta=0 so differences are computed correctly
 template <bool noZS>
 int EcalClusterToolsT<noZS>::deltaIEta(int seed_ieta, int rh_ieta) {
   // get rid of the fact that there is no ieta=0
   if (seed_ieta < 0)
     seed_ieta++;
   if (rh_ieta < 0)
     rh_ieta++;
   int rel_ieta = rh_ieta - seed_ieta;
   return rel_ieta;
 }
 
 //return (ieta,iphi) of highest energy recHit of the recHits passed to this function
 template <bool noZS>
 std::vector<int> EcalClusterToolsT<noZS>::getSeedPosition(
     const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs) {
   std::vector<int> seedPosition;
   float eSeedRH = 0;
   int iEtaSeedRH = 0;
   int iPhiSeedRH = 0;
 
   for (auto const &rhf : RH_ptrs_fracs) {
     const EcalRecHit *rh_ptr = rhf.first;
     //get iEta, iPhi
     EBDetId EBdetIdi(rh_ptr->detid());
     float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf.second);
 
     if (eSeedRH < rh_energy) {
       eSeedRH = rh_energy;
       iEtaSeedRH = EBdetIdi.ieta();
       iPhiSeedRH = EBdetIdi.iphi();
     }
 
   }  // for loop
 
   seedPosition.push_back(iEtaSeedRH);
   seedPosition.push_back(iPhiSeedRH);
   return seedPosition;
 }
 
 // return the total energy of the recHits passed to this function
 template <bool noZS>
 float EcalClusterToolsT<noZS>::getSumEnergy(const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs) {
   float sumE = 0.;
   for (const auto &hAndF : RH_ptrs_fracs) {
     sumE += hAndF.first->energy() * (noZS ? 1.0 : hAndF.second);
   }
   return sumE;
 }
 
 typedef EcalClusterToolsT<false> EcalClusterTools;
 
 namespace noZS {
   typedef EcalClusterToolsT<true> EcalClusterTools;
 }
 
 #endif

This page was automatically generated by the 2.2.1 LXR engine. The LXR team