Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-12 22:39:43

0001 #ifndef RecoEcal_EgammaCoreTools_EcalClusterTools_h
0002 #define RecoEcal_EgammaCoreTools_EcalClusterTools_h
0003 
0004 /** \class EcalClusterTools
0005  *  
0006  * various cluster tools (e.g. cluster shapes)
0007  *
0008  * \author Federico Ferri
0009  *
0010  * editing author: M.B. Anderson
0011  * 
0012  *
0013  */
0014 
0015 /** Note about "noZS" and the new templated version of EcalClusterTools
0016  *
0017  * "noZS" means "no zero-suppression" and means that the rechit quantity is 
0018  * computed without fractions and always including the full matrix of crystals
0019  * about the seed.
0020  *
0021  * To access the noZS variables use the "noZS::EcalClusterTools"
0022  * The interface at this point is the same as the standard tools, except the
0023  * behavior is altered to provide the output as described above.
0024  *
0025  * This was achieved by templating the EcalClusterTools.
0026  * Do not use EcalClusterToolsT<> directly.
0027  *
0028  */
0029 
0030 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0031 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0032 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0033 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0034 //#include "DataFormats/Math/interface/Point3D.h"
0035 #include "DataFormats/Math/interface/Vector3D.h"
0036 //includes for ShowerShape function to work
0037 #include <vector>
0038 #include <cmath>
0039 #include <TMath.h>
0040 #include <TMatrixT.h>
0041 #include <TMatrixD.h>
0042 #include <TVectorT.h>
0043 #include <TVectorD.h>
0044 
0045 #include "DataFormats/DetId/interface/DetId.h"
0046 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0047 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0048 #include "RecoCaloTools/Navigation/interface/CaloRectangle.h"
0049 
0050 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0051 
0052 #include "CLHEP/Geometry/Transform3D.h"
0053 
0054 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0055 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0056 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0057 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0058 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0059 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0060 #include "RecoEcal/EgammaCoreTools/interface/EgammaLocalCovParamDefaults.h"
0061 
0062 class DetId;
0063 class CaloTopology;
0064 class CaloGeometry;
0065 
0066 struct Cluster2ndMoments {
0067   // major and minor cluster moments wrt principale axes:
0068   float sMaj;
0069   float sMin;
0070   // angle between sMaj and phi:
0071   float alpha;
0072   Cluster2ndMoments() : sMaj(0.), sMin(0.), alpha(0.) {}
0073 };
0074 
0075 template <bool noZS>
0076 class EcalClusterToolsT {
0077 public:
0078   // various energies in the matrix nxn surrounding the maximum energy crystal of the input cluster
0079   //we use an eta/phi coordinate system rather than phi/eta
0080   //note e3x2 does not have a definate eta/phi geometry, it takes the maximum 3x2 block containing the
0081   //seed regardless of whether that 3 in eta or phi
0082   static float e1x3(const reco::BasicCluster &cluster,
0083                     const EcalRecHitCollection *recHits,
0084                     const CaloTopology *topology);
0085 
0086   static float e3x1(const reco::BasicCluster &cluster,
0087                     const EcalRecHitCollection *recHits,
0088                     const CaloTopology *topology);
0089 
0090   static float e1x5(const reco::BasicCluster &cluster,
0091                     const EcalRecHitCollection *recHits,
0092                     const CaloTopology *topology);
0093 
0094   static float e5x1(const reco::BasicCluster &cluster,
0095                     const EcalRecHitCollection *recHits,
0096                     const CaloTopology *topology);
0097 
0098   static float e2x2(const reco::BasicCluster &cluster,
0099                     const EcalRecHitCollection *recHits,
0100                     const CaloTopology *topology);
0101 
0102   static float e3x2(const reco::BasicCluster &cluster,
0103                     const EcalRecHitCollection *recHits,
0104                     const CaloTopology *topology);
0105 
0106   static float e3x3(const reco::BasicCluster &cluster,
0107                     const EcalRecHitCollection *recHits,
0108                     const CaloTopology *topology);
0109 
0110   static float e4x4(const reco::BasicCluster &cluster,
0111                     const EcalRecHitCollection *recHits,
0112                     const CaloTopology *topology);
0113 
0114   static float e5x5(const reco::BasicCluster &cluster,
0115                     const EcalRecHitCollection *recHits,
0116                     const CaloTopology *topology);
0117   static int n5x5(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology);
0118 
0119   // energy in the 2x5 strip right of the max crystal (does not contain max crystal)
0120   // 2 crystals wide in eta, 5 wide in phi.
0121   static float e2x5Right(const reco::BasicCluster &cluster,
0122                          const EcalRecHitCollection *recHits,
0123                          const CaloTopology *topology);
0124   // energy in the 2x5 strip left of the max crystal (does not contain max crystal)
0125 
0126   static float e2x5Left(const reco::BasicCluster &cluster,
0127                         const EcalRecHitCollection *recHits,
0128                         const CaloTopology *topology);
0129   // energy in the 5x2 strip above the max crystal (does not contain max crystal)
0130   // 5 crystals wide in eta, 2 wide in phi.
0131 
0132   static float e2x5Top(const reco::BasicCluster &cluster,
0133                        const EcalRecHitCollection *recHits,
0134                        const CaloTopology *topology);
0135   // energy in the 5x2 strip below the max crystal (does not contain max crystal)
0136 
0137   static float e2x5Bottom(const reco::BasicCluster &cluster,
0138                           const EcalRecHitCollection *recHits,
0139                           const CaloTopology *topology);
0140   // energy in a 2x5 strip containing the seed (max) crystal.
0141   // 2 crystals wide in eta, 5 wide in phi.
0142   // it is the maximum of either (1x5left + 1x5center) or (1x5right + 1x5center)
0143   static float e2x5Max(const reco::BasicCluster &cluster,
0144                        const EcalRecHitCollection *recHits,
0145                        const CaloTopology *topology);
0146 
0147   // energies in the crystal left, right, top, bottom w.r.t. to the most energetic crystal
0148   static float eLeft(const reco::BasicCluster &cluster,
0149                      const EcalRecHitCollection *recHits,
0150                      const CaloTopology *topology);
0151 
0152   static float eRight(const reco::BasicCluster &cluster,
0153                       const EcalRecHitCollection *recHits,
0154                       const CaloTopology *topology);
0155 
0156   static float eTop(const reco::BasicCluster &cluster,
0157                     const EcalRecHitCollection *recHits,
0158                     const CaloTopology *topology);
0159 
0160   static float eBottom(const reco::BasicCluster &cluster,
0161                        const EcalRecHitCollection *recHits,
0162                        const CaloTopology *topology);
0163   // the energy of the most energetic crystal in the cluster
0164 
0165   static float eMax(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
0166 
0167   // the energy of the second most energetic crystal in the cluster
0168   static float e2nd(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
0169 
0170   // get the DetId and the energy of the maximum energy crystal of the input cluster
0171   static std::pair<DetId, float> getMaximum(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
0172 
0173   static std::vector<float> energyBasketFractionEta(const reco::BasicCluster &cluster,
0174                                                     const EcalRecHitCollection *recHits);
0175 
0176   static std::vector<float> energyBasketFractionPhi(const reco::BasicCluster &cluster,
0177                                                     const EcalRecHitCollection *recHits);
0178 
0179   // return a vector v with v[0] = etaLat, v[1] = phiLat, v[2] = lat
0180   static std::vector<float> lat(const reco::BasicCluster &cluster,
0181                                 const EcalRecHitCollection *recHits,
0182                                 const CaloGeometry *geometry,
0183                                 bool logW = true,
0184                                 float w0 = 4.7);
0185 
0186   // return an array v with v[0] = covEtaEta, v[1] = covEtaPhi, v[2] = covPhiPhi
0187   static std::array<float, 3> covariances(const reco::BasicCluster &cluster,
0188                                           const EcalRecHitCollection *recHits,
0189                                           const CaloTopology *topology,
0190                                           const CaloGeometry *geometry,
0191                                           float w0 = 4.7);
0192 
0193   // return an array v with v[0] = covIEtaIEta, v[1] = covIEtaIPhi, v[2] = covIPhiIPhi
0194   //this function calculates differences in eta/phi in units of crystals not global eta/phi
0195   //this is gives better performance in the crack regions of the calorimeter but gives otherwise identical results to covariances function
0196   //   except that it doesnt need an eta based correction funtion in the endcap
0197   //it is multipled by an approprate crystal size to ensure it gives similar values to covariances(...)
0198   //
0199   //Warning: covIEtaIEta has been studied by egamma, but so far covIPhiIPhi hasnt been studied extensively so there could be a bug in
0200   //         the covIPhiIEta or covIPhiIPhi calculations. I dont think there is but as it hasnt been heavily used, there might be one
0201   static std::array<float, 3> localCovariances(const reco::BasicCluster &cluster,
0202                                                const EcalRecHitCollection *recHits,
0203                                                const CaloTopology *topology,
0204                                                float w0 = EgammaLocalCovParamDefaults::kRelEnCut,
0205                                                const EcalPFRecHitThresholds *thresholds = nullptr,
0206                                                float multEB = 0.0,
0207                                                float multEE = 0.0);
0208 
0209   static std::array<float, 3> scLocalCovariances(const reco::SuperCluster &cluster,
0210                                                  const EcalRecHitCollection *recHits,
0211                                                  const CaloTopology *topology,
0212                                                  float w0 = 4.7);
0213 
0214   // cluster second moments with respect to principal axes:
0215   static Cluster2ndMoments cluster2ndMoments(const reco::BasicCluster &basicCluster,
0216                                              const EcalRecHitCollection &recHits,
0217                                              double phiCorrectionFactor = 0.8,
0218                                              double w0 = 4.7,
0219                                              bool useLogWeights = true);
0220 
0221   static Cluster2ndMoments cluster2ndMoments(const reco::SuperCluster &superCluster,
0222                                              const EcalRecHitCollection &recHits,
0223                                              double phiCorrectionFactor = 0.8,
0224                                              double w0 = 4.7,
0225                                              bool useLogWeights = true);
0226   static Cluster2ndMoments cluster2ndMoments(const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs,
0227                                              double phiCorrectionFactor = 0.8,
0228                                              double w0 = 4.7,
0229                                              bool useLogWeights = true);
0230 
0231   static double zernike20(const reco::BasicCluster &cluster,
0232                           const EcalRecHitCollection *recHits,
0233                           const CaloGeometry *geometry,
0234                           double R0 = 6.6,
0235                           bool logW = true,
0236                           float w0 = 4.7);
0237   static double zernike42(const reco::BasicCluster &cluster,
0238                           const EcalRecHitCollection *recHits,
0239                           const CaloGeometry *geometry,
0240                           double R0 = 6.6,
0241                           bool logW = true,
0242                           float w0 = 4.7);
0243 
0244   // get the detId's of a matrix centered in the maximum energy crystal = (0,0)
0245   // the size is specified by ixMin, ixMax, iyMin, iyMax in unit of crystals
0246   static std::vector<DetId> matrixDetId(const CaloTopology *topology, DetId id, CaloRectangle rectangle);
0247   static std::vector<DetId> matrixDetId(const CaloTopology *topology, DetId id, int size) {
0248     return matrixDetId(topology, id, CaloRectangle{-size, size, -size, size});
0249   }
0250 
0251   // get the energy deposited in a matrix centered in the maximum energy crystal = (0,0)
0252   // the size is specified by ixMin, ixMax, iyMin, iyMax in unit of crystals
0253   static float matrixEnergy(const reco::BasicCluster &cluster,
0254                             const EcalRecHitCollection *recHits,
0255                             const CaloTopology *topology,
0256                             DetId id,
0257                             CaloRectangle rectangle);
0258   static int matrixSize(const reco::BasicCluster &cluster,
0259                         const EcalRecHitCollection *recHits,
0260                         const CaloTopology *topology,
0261                         DetId id,
0262                         CaloRectangle rectangle);
0263 
0264   static float getFraction(const std::vector<std::pair<DetId, float>> &v_id, DetId id);
0265   // get the DetId and the energy of the maximum energy crystal in a vector of DetId
0266   static std::pair<DetId, float> getMaximum(const std::vector<std::pair<DetId, float>> &v_id,
0267                                             const EcalRecHitCollection *recHits);
0268 
0269   // get the energy of a DetId, return 0 if the DetId is not in the collection
0270   static float recHitEnergy(DetId id, const EcalRecHitCollection *recHits);
0271 
0272   //Shower shape variables return vector <Roundness, Angle> of a photon
0273   static std::vector<float> roundnessBarrelSuperClusters(const reco::SuperCluster &superCluster,
0274                                                          const EcalRecHitCollection &recHits,
0275                                                          int weightedPositionMethod = 0,
0276                                                          float energyThreshold = 0.0);
0277   static std::vector<float> roundnessBarrelSuperClustersUserExtended(const reco::SuperCluster &superCluster,
0278                                                                      const EcalRecHitCollection &recHits,
0279                                                                      int ieta_delta = 0,
0280                                                                      int iphi_delta = 0,
0281                                                                      float energyRHThresh = 0.00000,
0282                                                                      int weightedPositionMethod = 0);
0283   static std::vector<float> roundnessSelectedBarrelRecHits(
0284       const std::vector<std::pair<const EcalRecHit *, float>> &rhVector, int weightedPositionMethod = 0);
0285 
0286   //works out the number of staturated crystals in 5x5
0287   static int nrSaturatedCrysIn5x5(const DetId &id, const EcalRecHitCollection *recHits, const CaloTopology *topology);
0288 
0289 private:
0290   struct EcalClusterEnergyDeposition {
0291     double deposited_energy;
0292     double r;
0293     double phi;
0294   };
0295 
0296   static std::vector<EcalClusterEnergyDeposition> getEnergyDepTopology(const reco::BasicCluster &cluster,
0297                                                                        const EcalRecHitCollection *recHits,
0298                                                                        const CaloGeometry *geometry,
0299                                                                        bool logW,
0300                                                                        float w0);
0301 
0302   static math::XYZVector meanClusterPosition(const reco::BasicCluster &cluster,
0303                                              const EcalRecHitCollection *recHits,
0304                                              const CaloTopology *topology,
0305                                              const CaloGeometry *geometry);
0306 
0307   //return energy weighted mean distance from the seed crystal in number of crystals
0308   //<iEta,iPhi>, iPhi is not defined for endcap and is returned as zero
0309   static std::pair<float, float> mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster,
0310                                                                  const EcalRecHitCollection *recHits,
0311                                                                  const CaloTopology *topology);
0312 
0313   static std::pair<float, float> mean5x5PositionInXY(const reco::BasicCluster &cluster,
0314                                                      const EcalRecHitCollection *recHits,
0315                                                      const CaloTopology *topology);
0316 
0317   static double f00(double r) { return 1; }
0318   static double f11(double r) { return r; }
0319   static double f20(double r) { return 2.0 * r * r - 1.0; }
0320   static double f22(double r) { return r * r; }
0321   static double f31(double r) { return 3.0 * r * r * r - 2.0 * r; }
0322   static double f33(double r) { return r * r * r; }
0323   static double f40(double r) { return 6.0 * r * r * r * r - 6.0 * r * r + 1.0; }
0324   static double f42(double r) { return 4.0 * r * r * r * r - 3.0 * r * r; }
0325   static double f44(double r) { return r * r * r * r; }
0326   static double f51(double r) { return 10.0 * pow(r, 5) - 12.0 * pow(r, 3) + 3.0 * r; }
0327   static double f53(double r) { return 5.0 * pow(r, 5) - 4.0 * pow(r, 3); }
0328   static double f55(double r) { return pow(r, 5); }
0329 
0330   static double absZernikeMoment(const reco::BasicCluster &cluster,
0331                                  const EcalRecHitCollection *recHits,
0332                                  const CaloGeometry *geometry,
0333                                  int n,
0334                                  int m,
0335                                  double R0,
0336                                  bool logW,
0337                                  float w0);
0338   static double fast_AbsZernikeMoment(const reco::BasicCluster &cluster,
0339                                       const EcalRecHitCollection *recHits,
0340                                       const CaloGeometry *geometry,
0341                                       int n,
0342                                       int m,
0343                                       double R0,
0344                                       bool logW,
0345                                       float w0);
0346   static double calc_AbsZernikeMoment(const reco::BasicCluster &cluster,
0347                                       const EcalRecHitCollection *recHits,
0348                                       const CaloGeometry *geometry,
0349                                       int n,
0350                                       int m,
0351                                       double R0,
0352                                       bool logW,
0353                                       float w0);
0354 
0355   static double factorial(int n) {
0356     double res = 1.;
0357     for (int i = 2; i <= n; ++i)
0358       res *= i;
0359     return res;
0360   }
0361 
0362   //useful functions for the localCovariances function
0363   static float getIEta(const DetId &id);
0364   static float getIPhi(const DetId &id);
0365   static float getNormedIX(const DetId &id);
0366   static float getNormedIY(const DetId &id);
0367   static float getDPhiEndcap(const DetId &crysId, float meanX, float meanY);
0368   static float getNrCrysDiffInEta(const DetId &crysId, const DetId &orginId);
0369   static float getNrCrysDiffInPhi(const DetId &crysId, const DetId &orginId);
0370 
0371   //useful functions for showerRoundnessBarrel function
0372   static int deltaIEta(int seed_ieta, int rh_ieta);
0373   static int deltaIPhi(int seed_iphi, int rh_iphi);
0374   static std::vector<int> getSeedPosition(const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs);
0375   static float getSumEnergy(const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs);
0376   static float computeWeight(float eRH, float energyTotal, int weightedPositionMethod);
0377 };
0378 
0379 // implementation
0380 template <bool noZS>
0381 float EcalClusterToolsT<noZS>::getFraction(const std::vector<std::pair<DetId, float>> &v_id, DetId id) {
0382   if (noZS)
0383     return 1.0;
0384   float frac = 0.0;
0385   for (size_t i = 0; i < v_id.size(); ++i) {
0386     if (v_id[i].first.rawId() == id.rawId()) {
0387       frac = v_id[i].second;
0388       break;
0389     }
0390   }
0391   return frac;
0392 }
0393 
0394 template <bool noZS>
0395 std::pair<DetId, float> EcalClusterToolsT<noZS>::getMaximum(const std::vector<std::pair<DetId, float>> &v_id,
0396                                                             const EcalRecHitCollection *recHits) {
0397   float max = 0;
0398   DetId id(0);
0399   for (size_t i = 0; i < v_id.size(); ++i) {
0400     float energy = recHitEnergy(v_id[i].first, recHits) * (noZS ? 1.0 : v_id[i].second);
0401     if (energy > max) {
0402       max = energy;
0403       id = v_id[i].first;
0404     }
0405   }
0406   return std::pair<DetId, float>(id, max);
0407 }
0408 
0409 template <bool noZS>
0410 std::pair<DetId, float> EcalClusterToolsT<noZS>::getMaximum(const reco::BasicCluster &cluster,
0411                                                             const EcalRecHitCollection *recHits) {
0412   return getMaximum(cluster.hitsAndFractions(), recHits);
0413 }
0414 
0415 template <bool noZS>
0416 float EcalClusterToolsT<noZS>::recHitEnergy(DetId id, const EcalRecHitCollection *recHits) {
0417   if (id == DetId(0)) {
0418     return 0;
0419   } else {
0420     EcalRecHitCollection::const_iterator it = recHits->find(id);
0421     if (it != recHits->end()) {
0422       if (noZS && (it->checkFlag(EcalRecHit::kTowerRecovered) || it->checkFlag(EcalRecHit::kWeird) ||
0423                    (it->detid().subdetId() == EcalBarrel && it->checkFlag(EcalRecHit::kDiWeird))))
0424         return 0.0;
0425       else
0426         return it->energy();
0427     } else {
0428       //throw cms::Exception("EcalRecHitNotFound") << "The recHit corresponding to the DetId" << id.rawId() << " not found in the EcalRecHitCollection";
0429       // the recHit is not in the collection (hopefully zero suppressed)
0430       return 0;
0431     }
0432   }
0433   return 0;
0434 }
0435 
0436 // Returns the energy in a rectangle of crystals
0437 // specified in eta by ixMin and ixMax
0438 //       and in phi by iyMin and iyMax
0439 //
0440 // Reference picture (X=seed crystal)
0441 //    iy ___________
0442 //     2 |_|_|_|_|_|
0443 //     1 |_|_|_|_|_|
0444 //     0 |_|_|X|_|_|
0445 //    -1 |_|_|_|_|_|
0446 //    -2 |_|_|_|_|_|
0447 //      -2 -1 0 1 2 ix
0448 template <bool noZS>
0449 float EcalClusterToolsT<noZS>::matrixEnergy(const reco::BasicCluster &cluster,
0450                                             const EcalRecHitCollection *recHits,
0451                                             const CaloTopology *topology,
0452                                             DetId id,
0453                                             CaloRectangle rectangle) {
0454   float energy = 0;
0455   auto const &vId = cluster.hitsAndFractions();
0456 
0457   for (auto const &detId : rectangle(id, *topology)) {
0458     energy += recHitEnergy(detId, recHits) * getFraction(vId, detId);
0459   }
0460 
0461   return energy;
0462 }
0463 
0464 template <bool noZS>
0465 int EcalClusterToolsT<noZS>::matrixSize(const reco::BasicCluster &cluster,
0466                                         const EcalRecHitCollection *recHits,
0467                                         const CaloTopology *topology,
0468                                         DetId id,
0469                                         CaloRectangle rectangle) {
0470   int result = 0;
0471 
0472   for (auto const &detId : rectangle(id, *topology)) {
0473     const float energy = recHitEnergy(detId, recHits);
0474     const float frac = getFraction(cluster.hitsAndFractions(), detId);
0475     if (energy * frac > 0)
0476       result++;
0477   }
0478   return result;
0479 }
0480 
0481 template <bool noZS>
0482 std::vector<DetId> EcalClusterToolsT<noZS>::matrixDetId(const CaloTopology *topology,
0483                                                         DetId id,
0484                                                         CaloRectangle rectangle) {
0485   std::vector<DetId> v;
0486   for (auto const &detId : rectangle(id, *topology)) {
0487     if (detId != DetId(0))
0488       v.push_back(detId);
0489   }
0490   return v;
0491 }
0492 
0493 template <bool noZS>
0494 float EcalClusterToolsT<noZS>::e2x2(const reco::BasicCluster &cluster,
0495                                     const EcalRecHitCollection *recHits,
0496                                     const CaloTopology *topology) {
0497   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0498   std::list<float> energies;
0499   float max_E = matrixEnergy(cluster, recHits, topology, id, {-1, 0, -1, 0});
0500   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-1, 0, 0, 1}));
0501   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {0, 1, 0, 1}));
0502   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {0, 1, -1, 0}));
0503   return max_E;
0504 }
0505 
0506 template <bool noZS>
0507 float EcalClusterToolsT<noZS>::e3x2(const reco::BasicCluster &cluster,
0508                                     const EcalRecHitCollection *recHits,
0509                                     const CaloTopology *topology) {
0510   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0511   float max_E = matrixEnergy(cluster, recHits, topology, id, {-1, 1, -1, 0});
0512   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {0, 1, -1, 1}));
0513   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-1, 1, 0, 1}));
0514   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-1, 0, -1, 1}));
0515   return max_E;
0516 }
0517 
0518 template <bool noZS>
0519 float EcalClusterToolsT<noZS>::e3x3(const reco::BasicCluster &cluster,
0520                                     const EcalRecHitCollection *recHits,
0521                                     const CaloTopology *topology) {
0522   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0523   return matrixEnergy(cluster, recHits, topology, id, {-1, 1, -1, 1});
0524 }
0525 
0526 template <bool noZS>
0527 float EcalClusterToolsT<noZS>::e4x4(const reco::BasicCluster &cluster,
0528                                     const EcalRecHitCollection *recHits,
0529                                     const CaloTopology *topology) {
0530   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0531   float max_E = matrixEnergy(cluster, recHits, topology, id, {-1, 2, -2, 1});
0532   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-2, 1, -2, 1}));
0533   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-2, 1, -1, 2}));
0534   max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-1, 2, -1, 2}));
0535   return max_E;
0536 }
0537 
0538 template <bool noZS>
0539 float EcalClusterToolsT<noZS>::e5x5(const reco::BasicCluster &cluster,
0540                                     const EcalRecHitCollection *recHits,
0541                                     const CaloTopology *topology) {
0542   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0543   return matrixEnergy(cluster, recHits, topology, id, {-2, 2, -2, 2});
0544 }
0545 
0546 template <bool noZS>
0547 int EcalClusterToolsT<noZS>::n5x5(const reco::BasicCluster &cluster,
0548                                   const EcalRecHitCollection *recHits,
0549                                   const CaloTopology *topology) {
0550   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0551   return matrixSize(cluster, recHits, topology, id, {-2, 2, -2, 2});
0552 }
0553 
0554 template <bool noZS>
0555 float EcalClusterToolsT<noZS>::eMax(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits) {
0556   return getMaximum(cluster.hitsAndFractions(), recHits).second;
0557 }
0558 
0559 template <bool noZS>
0560 float EcalClusterToolsT<noZS>::e2nd(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits) {
0561   std::vector<float> energies;
0562   const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
0563   energies.reserve(v_id.size());
0564   if (v_id.size() < 2)
0565     return 0;
0566   for (size_t i = 0; i < v_id.size(); ++i) {
0567     energies.push_back(recHitEnergy(v_id[i].first, recHits) * (noZS ? 1.0 : v_id[i].second));
0568   }
0569   std::partial_sort(energies.begin(), energies.begin() + 2, energies.end(), std::greater<float>());
0570   return energies[1];
0571 }
0572 
0573 template <bool noZS>
0574 float EcalClusterToolsT<noZS>::e2x5Right(const reco::BasicCluster &cluster,
0575                                          const EcalRecHitCollection *recHits,
0576                                          const CaloTopology *topology) {
0577   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0578   return matrixEnergy(cluster, recHits, topology, id, {1, 2, -2, 2});
0579 }
0580 
0581 template <bool noZS>
0582 float EcalClusterToolsT<noZS>::e2x5Left(const reco::BasicCluster &cluster,
0583                                         const EcalRecHitCollection *recHits,
0584                                         const CaloTopology *topology) {
0585   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0586   return matrixEnergy(cluster, recHits, topology, id, {-2, -1, -2, 2});
0587 }
0588 
0589 template <bool noZS>
0590 float EcalClusterToolsT<noZS>::e2x5Top(const reco::BasicCluster &cluster,
0591                                        const EcalRecHitCollection *recHits,
0592                                        const CaloTopology *topology) {
0593   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0594   return matrixEnergy(cluster, recHits, topology, id, {-2, 2, 1, 2});
0595 }
0596 
0597 template <bool noZS>
0598 float EcalClusterToolsT<noZS>::e2x5Bottom(const reco::BasicCluster &cluster,
0599                                           const EcalRecHitCollection *recHits,
0600                                           const CaloTopology *topology) {
0601   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0602   return matrixEnergy(cluster, recHits, topology, id, {-2, 2, -2, -1});
0603 }
0604 
0605 // Energy in 2x5 strip containing the max crystal.
0606 // Adapted from code by Sam Harper
0607 template <bool noZS>
0608 float EcalClusterToolsT<noZS>::e2x5Max(const reco::BasicCluster &cluster,
0609                                        const EcalRecHitCollection *recHits,
0610                                        const CaloTopology *topology) {
0611   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0612 
0613   // 1x5 strip left of seed
0614   float left = matrixEnergy(cluster, recHits, topology, id, {-1, -1, -2, 2});
0615   // 1x5 strip right of seed
0616   float right = matrixEnergy(cluster, recHits, topology, id, {1, 1, -2, 2});
0617   // 1x5 strip containing seed
0618   float centre = matrixEnergy(cluster, recHits, topology, id, {0, 0, -2, 2});
0619 
0620   // Return the maximum of (left+center) or (right+center) strip
0621   return left > right ? left + centre : right + centre;
0622 }
0623 
0624 template <bool noZS>
0625 float EcalClusterToolsT<noZS>::e1x5(const reco::BasicCluster &cluster,
0626                                     const EcalRecHitCollection *recHits,
0627                                     const CaloTopology *topology) {
0628   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0629   return matrixEnergy(cluster, recHits, topology, id, {0, 0, -2, 2});
0630 }
0631 
0632 template <bool noZS>
0633 float EcalClusterToolsT<noZS>::e5x1(const reco::BasicCluster &cluster,
0634                                     const EcalRecHitCollection *recHits,
0635                                     const CaloTopology *topology) {
0636   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0637   return matrixEnergy(cluster, recHits, topology, id, {-2, 2, 0, 0});
0638 }
0639 
0640 template <bool noZS>
0641 float EcalClusterToolsT<noZS>::e1x3(const reco::BasicCluster &cluster,
0642                                     const EcalRecHitCollection *recHits,
0643                                     const CaloTopology *topology) {
0644   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0645   return matrixEnergy(cluster, recHits, topology, id, {0, 0, -1, 1});
0646 }
0647 
0648 template <bool noZS>
0649 float EcalClusterToolsT<noZS>::e3x1(const reco::BasicCluster &cluster,
0650                                     const EcalRecHitCollection *recHits,
0651                                     const CaloTopology *topology) {
0652   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0653   return matrixEnergy(cluster, recHits, topology, id, {-1, 1, 0, 0});
0654 }
0655 
0656 template <bool noZS>
0657 float EcalClusterToolsT<noZS>::eLeft(const reco::BasicCluster &cluster,
0658                                      const EcalRecHitCollection *recHits,
0659                                      const CaloTopology *topology) {
0660   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0661   return matrixEnergy(cluster, recHits, topology, id, {-1, -1, 0, 0});
0662 }
0663 
0664 template <bool noZS>
0665 float EcalClusterToolsT<noZS>::eRight(const reco::BasicCluster &cluster,
0666                                       const EcalRecHitCollection *recHits,
0667                                       const CaloTopology *topology) {
0668   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0669   return matrixEnergy(cluster, recHits, topology, id, {1, 1, 0, 0});
0670 }
0671 
0672 template <bool noZS>
0673 float EcalClusterToolsT<noZS>::eTop(const reco::BasicCluster &cluster,
0674                                     const EcalRecHitCollection *recHits,
0675                                     const CaloTopology *topology) {
0676   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0677   return matrixEnergy(cluster, recHits, topology, id, {0, 0, 1, 1});
0678 }
0679 
0680 template <bool noZS>
0681 float EcalClusterToolsT<noZS>::eBottom(const reco::BasicCluster &cluster,
0682                                        const EcalRecHitCollection *recHits,
0683                                        const CaloTopology *topology) {
0684   DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0685   return matrixEnergy(cluster, recHits, topology, id, {0, 0, -1, -1});
0686 }
0687 
0688 template <bool noZS>
0689 std::vector<float> EcalClusterToolsT<noZS>::energyBasketFractionEta(const reco::BasicCluster &cluster,
0690                                                                     const EcalRecHitCollection *recHits) {
0691   std::vector<float> basketFraction(2 * EBDetId::kModulesPerSM);
0692   float clusterEnergy = cluster.energy();
0693   const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
0694   if (v_id[0].first.subdetId() != EcalBarrel) {
0695     edm::LogWarning("EcalClusterToolsT<noZS>::energyBasketFractionEta")
0696         << "Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel "
0697            "basic-clusters. Returning empty vector.";
0698     return basketFraction;
0699   }
0700   for (size_t i = 0; i < v_id.size(); ++i) {
0701     basketFraction[EBDetId(v_id[i].first).im() - 1 + EBDetId(v_id[i].first).positiveZ() * EBDetId::kModulesPerSM] +=
0702         recHitEnergy(v_id[i].first, recHits) * v_id[i].second / clusterEnergy;
0703   }
0704   std::sort(basketFraction.rbegin(), basketFraction.rend());
0705   return basketFraction;
0706 }
0707 
0708 template <bool noZS>
0709 std::vector<float> EcalClusterToolsT<noZS>::energyBasketFractionPhi(const reco::BasicCluster &cluster,
0710                                                                     const EcalRecHitCollection *recHits) {
0711   std::vector<float> basketFraction(2 * (EBDetId::MAX_IPHI / EBDetId::kCrystalsInPhi));
0712   float clusterEnergy = cluster.energy();
0713   const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
0714   if (v_id[0].first.subdetId() != EcalBarrel) {
0715     edm::LogWarning("EcalClusterToolsT<noZS>::energyBasketFractionPhi")
0716         << "Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel "
0717            "basic-clusters. Returning empty vector.";
0718     return basketFraction;
0719   }
0720   for (size_t i = 0; i < v_id.size(); ++i) {
0721     basketFraction[(EBDetId(v_id[i].first).iphi() - 1) / EBDetId::kCrystalsInPhi +
0722                    EBDetId(v_id[i].first).positiveZ() * EBDetId::kTowersInPhi] +=
0723         recHitEnergy(v_id[i].first, recHits) * (noZS ? 1.0 : v_id[i].second) / clusterEnergy;
0724   }
0725   std::sort(basketFraction.rbegin(), basketFraction.rend());
0726   return basketFraction;
0727 }
0728 
0729 template <bool noZS>
0730 std::vector<typename EcalClusterToolsT<noZS>::EcalClusterEnergyDeposition>
0731 EcalClusterToolsT<noZS>::getEnergyDepTopology(const reco::BasicCluster &cluster,
0732                                               const EcalRecHitCollection *recHits,
0733                                               const CaloGeometry *geometry,
0734                                               bool logW,
0735                                               float w0) {
0736   std::vector<typename EcalClusterToolsT<noZS>::EcalClusterEnergyDeposition> energyDistribution;
0737   // init a map of the energy deposition centered on the
0738   // cluster centroid. This is for momenta calculation only.
0739   CLHEP::Hep3Vector clVect(cluster.position().x(), cluster.position().y(), cluster.position().z());
0740   CLHEP::Hep3Vector clDir(clVect);
0741   clDir *= 1.0 / clDir.mag();
0742   // in the transverse plane, axis perpendicular to clusterDir
0743   CLHEP::Hep3Vector theta_axis(clDir.y(), -clDir.x(), 0.0);
0744   theta_axis *= 1.0 / theta_axis.mag();
0745   CLHEP::Hep3Vector phi_axis = theta_axis.cross(clDir);
0746 
0747   const std::vector<std::pair<DetId, float>> &clusterDetIds = cluster.hitsAndFractions();
0748 
0749   EcalClusterEnergyDeposition clEdep;
0750   EcalRecHit testEcalRecHit;
0751   std::vector<std::pair<DetId, float>>::const_iterator posCurrent;
0752   // loop over crystals
0753   for (posCurrent = clusterDetIds.begin(); posCurrent != clusterDetIds.end(); ++posCurrent) {
0754     EcalRecHitCollection::const_iterator itt = recHits->find((*posCurrent).first);
0755     testEcalRecHit = *itt;
0756 
0757     if (((*posCurrent).first != DetId(0)) && (recHits->find((*posCurrent).first) != recHits->end())) {
0758       clEdep.deposited_energy = testEcalRecHit.energy() * (noZS ? 1.0 : (*posCurrent).second);
0759       // if logarithmic weight is requested, apply cut on minimum energy of the recHit
0760       if (logW) {
0761         //double w0 = parameterMap_.find("W0")->second;
0762 
0763         double weight = std::max(0.0, w0 + log(std::abs(clEdep.deposited_energy) / cluster.energy()));
0764         if (weight == 0) {
0765           LogDebug("ClusterShapeAlgo") << "Crystal has insufficient energy: E = " << clEdep.deposited_energy
0766                                        << " GeV; skipping... ";
0767           continue;
0768         } else
0769           LogDebug("ClusterShapeAlgo") << "===> got crystal. Energy = " << clEdep.deposited_energy << " GeV. ";
0770       }
0771       DetId id_ = (*posCurrent).first;
0772       const CaloSubdetectorGeometry *geo = geometry->getSubdetectorGeometry(id_);
0773       auto this_cell = geo->getGeometry(id_);
0774       const GlobalPoint &cellPos = this_cell->getPosition();
0775       CLHEP::Hep3Vector gblPos(cellPos.x(), cellPos.y(), cellPos.z());  //surface position?
0776       // Evaluate the distance from the cluster centroid
0777       CLHEP::Hep3Vector diff = gblPos - clVect;
0778       // Important: for the moment calculation, only the "lateral distance" is important
0779       // "lateral distance" r_i = distance of the digi position from the axis Origin-Cluster Center
0780       // ---> subtract the projection on clDir
0781       CLHEP::Hep3Vector DigiVect = diff - diff.dot(clDir) * clDir;
0782       clEdep.r = DigiVect.mag();
0783       LogDebug("ClusterShapeAlgo") << "E = " << clEdep.deposited_energy << "\tdiff = " << diff.mag()
0784                                    << "\tr = " << clEdep.r;
0785       clEdep.phi = DigiVect.angle(theta_axis);
0786       if (DigiVect.dot(phi_axis) < 0)
0787         clEdep.phi = 2 * M_PI - clEdep.phi;
0788       energyDistribution.push_back(clEdep);
0789     }
0790   }
0791   return energyDistribution;
0792 }
0793 
0794 template <bool noZS>
0795 std::vector<float> EcalClusterToolsT<noZS>::lat(const reco::BasicCluster &cluster,
0796                                                 const EcalRecHitCollection *recHits,
0797                                                 const CaloGeometry *geometry,
0798                                                 bool logW,
0799                                                 float w0) {
0800   std::vector<EcalClusterToolsT::EcalClusterEnergyDeposition> energyDistribution =
0801       getEnergyDepTopology(cluster, recHits, geometry, logW, w0);
0802 
0803   std::vector<float> lat;
0804   double r, redmoment = 0;
0805   double phiRedmoment = 0;
0806   double etaRedmoment = 0;
0807   int n, n1, n2, tmp;
0808   int clusterSize = energyDistribution.size();
0809   float etaLat_, phiLat_, lat_;
0810   if (clusterSize < 3) {
0811     etaLat_ = 0.0;
0812     lat_ = 0.0;
0813     lat.push_back(0.);
0814     lat.push_back(0.);
0815     lat.push_back(0.);
0816     return lat;
0817   }
0818 
0819   n1 = 0;
0820   n2 = 1;
0821   if (energyDistribution[1].deposited_energy > energyDistribution[0].deposited_energy) {
0822     tmp = n2;
0823     n2 = n1;
0824     n1 = tmp;
0825   }
0826   for (int i = 2; i < clusterSize; i++) {
0827     n = i;
0828     if (energyDistribution[i].deposited_energy > energyDistribution[n1].deposited_energy) {
0829       tmp = n2;
0830       n2 = n1;
0831       n1 = i;
0832       n = tmp;
0833     } else {
0834       if (energyDistribution[i].deposited_energy > energyDistribution[n2].deposited_energy) {
0835         tmp = n2;
0836         n2 = i;
0837         n = tmp;
0838       }
0839     }
0840 
0841     r = energyDistribution[n].r;
0842     redmoment += r * r * energyDistribution[n].deposited_energy;
0843     double rphi = r * cos(energyDistribution[n].phi);
0844     phiRedmoment += rphi * rphi * energyDistribution[n].deposited_energy;
0845     double reta = r * sin(energyDistribution[n].phi);
0846     etaRedmoment += reta * reta * energyDistribution[n].deposited_energy;
0847   }
0848   double e1 = energyDistribution[n1].deposited_energy;
0849   double e2 = energyDistribution[n2].deposited_energy;
0850 
0851   lat_ = redmoment / (redmoment + 2.19 * 2.19 * (e1 + e2));
0852   phiLat_ = phiRedmoment / (phiRedmoment + 2.19 * 2.19 * (e1 + e2));
0853   etaLat_ = etaRedmoment / (etaRedmoment + 2.19 * 2.19 * (e1 + e2));
0854 
0855   lat.push_back(etaLat_);
0856   lat.push_back(phiLat_);
0857   lat.push_back(lat_);
0858   return lat;
0859 }
0860 
0861 template <bool noZS>
0862 math::XYZVector EcalClusterToolsT<noZS>::meanClusterPosition(const reco::BasicCluster &cluster,
0863                                                              const EcalRecHitCollection *recHits,
0864                                                              const CaloTopology *topology,
0865                                                              const CaloGeometry *geometry) {
0866   // find mean energy position of a 5x5 cluster around the maximum
0867   math::XYZVector meanPosition(0.0, 0.0, 0.0);
0868   const std::vector<std::pair<DetId, float>> &hsAndFs = cluster.hitsAndFractions();
0869   std::vector<DetId> v_id = matrixDetId(topology, getMaximum(cluster, recHits).first, 2);
0870   for (const std::pair<DetId, float> &hitAndFrac : hsAndFs) {
0871     for (std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it) {
0872       if (hitAndFrac.first != *it && !noZS)
0873         continue;
0874       const CaloSubdetectorGeometry *geo = geometry->getSubdetectorGeometry(*it);
0875       GlobalPoint positionGP = geo->getGeometry(*it)->getPosition();
0876       math::XYZVector position(positionGP.x(), positionGP.y(), positionGP.z());
0877       meanPosition = meanPosition + recHitEnergy(*it, recHits) * position * hitAndFrac.second;
0878     }
0879     if (noZS)
0880       break;
0881   }
0882   return meanPosition / e5x5(cluster, recHits, topology);
0883 }
0884 
0885 //returns mean energy weighted eta/phi in crystals from the seed
0886 //iPhi is not defined for endcap and is returned as zero
0887 //return <eta,phi>
0888 //we have an issue in working out what to do for negative energies
0889 //I (Sam Harper) think it makes sense to ignore crystals with E<0 in the calculation as they are ignored
0890 //in the sigmaIEtaIEta calculation (well they arent fully ignored, they do still contribute to the e5x5 sum
0891 //in the sigmaIEtaIEta calculation but not here)
0892 template <bool noZS>
0893 std::pair<float, float> EcalClusterToolsT<noZS>::mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster,
0894                                                                                  const EcalRecHitCollection *recHits,
0895                                                                                  const CaloTopology *topology) {
0896   DetId seedId = getMaximum(cluster, recHits).first;
0897   float meanDEta = 0.;
0898   float meanDPhi = 0.;
0899   float energySum = 0.;
0900 
0901   const std::vector<std::pair<DetId, float>> &hsAndFs = cluster.hitsAndFractions();
0902   std::vector<DetId> v_id = matrixDetId(topology, seedId, 2);
0903   for (const std::pair<DetId, float> &hAndF : hsAndFs) {
0904     for (std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it) {
0905       if (hAndF.first != *it && !noZS)
0906         continue;
0907       float energy = recHitEnergy(*it, recHits) * hAndF.second;
0908       if (energy < 0.)
0909         continue;  //skipping negative energy crystals
0910       meanDEta += energy * getNrCrysDiffInEta(*it, seedId);
0911       meanDPhi += energy * getNrCrysDiffInPhi(*it, seedId);
0912       energySum += energy;
0913     }
0914     if (noZS)
0915       break;
0916   }
0917   meanDEta /= energySum;
0918   meanDPhi /= energySum;
0919   return std::pair<float, float>(meanDEta, meanDPhi);
0920 }
0921 
0922 //returns mean energy weighted x/y in normalised crystal coordinates
0923 //only valid for endcap, returns 0,0 for barrel
0924 //we have an issue in working out what to do for negative energies
0925 //I (Sam Harper) think it makes sense to ignore crystals with E<0 in the calculation as they are ignored
0926 //in the sigmaIEtaIEta calculation (well they arent fully ignored, they do still contribute to the e5x5 sum
0927 //in the sigmaIEtaIEta calculation but not here)
0928 template <bool noZS>
0929 std::pair<float, float> EcalClusterToolsT<noZS>::mean5x5PositionInXY(const reco::BasicCluster &cluster,
0930                                                                      const EcalRecHitCollection *recHits,
0931                                                                      const CaloTopology *topology) {
0932   DetId seedId = getMaximum(cluster, recHits).first;
0933 
0934   std::pair<float, float> meanXY(0., 0.);
0935   if (seedId.subdetId() == EcalBarrel)
0936     return meanXY;
0937 
0938   float energySum = 0.;
0939 
0940   const std::vector<std::pair<DetId, float>> &hsAndFs = cluster.hitsAndFractions();
0941   std::vector<DetId> v_id = matrixDetId(topology, seedId, 2);
0942   for (const std::pair<DetId, float> &hAndF : hsAndFs) {
0943     for (std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it) {
0944       if (hAndF.first != *it && !noZS)
0945         continue;
0946       float energy = recHitEnergy(*it, recHits) * hAndF.second;
0947       if (energy < 0.)
0948         continue;  //skipping negative energy crystals
0949       meanXY.first += energy * getNormedIX(*it);
0950       meanXY.second += energy * getNormedIY(*it);
0951       energySum += energy;
0952     }
0953     if (noZS)
0954       break;
0955   }
0956   meanXY.first /= energySum;
0957   meanXY.second /= energySum;
0958   return meanXY;
0959 }
0960 
0961 template <bool noZS>
0962 std::array<float, 3> EcalClusterToolsT<noZS>::covariances(const reco::BasicCluster &cluster,
0963                                                           const EcalRecHitCollection *recHits,
0964                                                           const CaloTopology *topology,
0965                                                           const CaloGeometry *geometry,
0966                                                           float w0) {
0967   float e_5x5 = e5x5(cluster, recHits, topology);
0968   float covEtaEta, covEtaPhi, covPhiPhi;
0969   if (e_5x5 >= 0.) {
0970     //double w0_ = parameterMap_.find("W0")->second;
0971     const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
0972     math::XYZVector meanPosition = meanClusterPosition(cluster, recHits, topology, geometry);
0973 
0974     // now we can calculate the covariances
0975     double numeratorEtaEta = 0;
0976     double numeratorEtaPhi = 0;
0977     double numeratorPhiPhi = 0;
0978     double denominator = 0;
0979 
0980     DetId id = getMaximum(v_id, recHits).first;
0981     CaloRectangle rectangle{-2, 2, -2, 2};
0982     for (auto const &detId : rectangle(id, *topology)) {
0983       float frac = getFraction(v_id, detId);
0984       float energy = recHitEnergy(detId, recHits) * frac;
0985 
0986       if (energy <= 0)
0987         continue;
0988 
0989       const CaloSubdetectorGeometry *geo = geometry->getSubdetectorGeometry(detId);
0990       GlobalPoint position = geo->getGeometry(detId)->getPosition();
0991 
0992       double dPhi = position.phi() - meanPosition.phi();
0993       if (dPhi > +Geom::pi()) {
0994         dPhi = Geom::twoPi() - dPhi;
0995       }
0996       if (dPhi < -Geom::pi()) {
0997         dPhi = Geom::twoPi() + dPhi;
0998       }
0999 
1000       double dEta = position.eta() - meanPosition.eta();
1001       double w = 0.;
1002       w = std::max(0.0f, w0 + std::log(energy / e_5x5));
1003 
1004       denominator += w;
1005       numeratorEtaEta += w * dEta * dEta;
1006       numeratorEtaPhi += w * dEta * dPhi;
1007       numeratorPhiPhi += w * dPhi * dPhi;
1008     }
1009 
1010     if (denominator != 0.0) {
1011       covEtaEta = numeratorEtaEta / denominator;
1012       covEtaPhi = numeratorEtaPhi / denominator;
1013       covPhiPhi = numeratorPhiPhi / denominator;
1014     } else {
1015       covEtaEta = 999.9;
1016       covEtaPhi = 999.9;
1017       covPhiPhi = 999.9;
1018     }
1019 
1020   } else {
1021     // Warn the user if there was no energy in the cells and return zeroes.
1022     //       std::cout << "\ClusterShapeAlgo::Calculate_Covariances:  no energy in supplied cells.\n";
1023     covEtaEta = 0;
1024     covEtaPhi = 0;
1025     covPhiPhi = 0;
1026   }
1027   std::array<float, 3> v{{covEtaEta, covEtaPhi, covPhiPhi}};
1028   return v;
1029 }
1030 
1031 //for covIEtaIEta,covIEtaIPhi and covIPhiIPhi are defined but only covIEtaIEta has been actively studied
1032 //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
1033 //it also does not require any eta correction function in the endcap
1034 //it is multipled by an approprate crystal size to ensure it gives similar values to covariances(...)
1035 template <bool noZS>
1036 std::array<float, 3> EcalClusterToolsT<noZS>::localCovariances(const reco::BasicCluster &cluster,
1037                                                                const EcalRecHitCollection *recHits,
1038                                                                const CaloTopology *topology,
1039                                                                float w0,
1040                                                                const EcalPFRecHitThresholds *thresholds,
1041                                                                float multEB,
1042                                                                float multEE) {
1043   float e_5x5 = e5x5(cluster, recHits, topology);
1044   float covEtaEta, covEtaPhi, covPhiPhi;
1045 
1046   if (e_5x5 >= 0.) {
1047     //double w0_ = parameterMap_.find("W0")->second;
1048     const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
1049     std::pair<float, float> mean5x5PosInNrCrysFromSeed = mean5x5PositionInLocalCrysCoord(cluster, recHits, topology);
1050     std::pair<float, float> mean5x5XYPos = mean5x5PositionInXY(cluster, recHits, topology);
1051 
1052     // now we can calculate the covariances
1053     double numeratorEtaEta = 0;
1054     double numeratorEtaPhi = 0;
1055     double numeratorPhiPhi = 0;
1056     double denominator = 0;
1057 
1058     //these allow us to scale the localCov by the crystal size
1059     //so that the localCovs have the same average value as the normal covs
1060     const double barrelCrysSize = 0.01745;  //approximate size of crystal in eta,phi in barrel
1061     const double endcapCrysSize = 0.0447;   //the approximate crystal size sigmaEtaEta was corrected to in the endcap
1062 
1063     DetId seedId = getMaximum(v_id, recHits).first;
1064 
1065     bool isBarrel = seedId.subdetId() == EcalBarrel;
1066     const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1067     float mult = isBarrel ? multEB : multEE;  // we will multiply PF RecHit threshold by mult.
1068                                               // mult = 1 should work reasonably well for noise cleaning.
1069     // dedicated studies showed mult=1.25 works best for Run3 in endcap, where noise is high.
1070     // If no noise cleaning is intended then put mult=0.
1071 
1072     CaloRectangle rectangle{-2, 2, -2, 2};
1073     for (auto const &detId : rectangle(seedId, *topology)) {
1074       float frac = getFraction(v_id, detId);
1075       float energy = recHitEnergy(detId, recHits) * frac;
1076 
1077       if ((thresholds == nullptr) && (mult != 0.0)) {
1078         throw cms::Exception("EmptyPFRechHitThresColl")
1079             << "In EcalClusterTools::localCovariances, if EcalPFRecHitThresholds==nulptr, then multEB and multEE "
1080                "should be 0 as well.";
1081       }
1082       float rhThres = 0.0;
1083       if (thresholds != nullptr) {
1084         rhThres = (*thresholds)[detId];  // access PFRechit thresholds for noise cleaning
1085       }
1086 
1087       if (energy <= (rhThres * mult))
1088         continue;
1089 
1090       float dEta = getNrCrysDiffInEta(detId, seedId) - mean5x5PosInNrCrysFromSeed.first;
1091       float dPhi = 0;
1092 
1093       if (isBarrel)
1094         dPhi = getNrCrysDiffInPhi(detId, seedId) - mean5x5PosInNrCrysFromSeed.second;
1095       else
1096         dPhi = getDPhiEndcap(detId, mean5x5XYPos.first, mean5x5XYPos.second);
1097 
1098       double w = std::max(0.0f, w0 + std::log(energy / e_5x5));
1099 
1100       denominator += w;
1101       numeratorEtaEta += w * dEta * dEta;
1102       numeratorEtaPhi += w * dEta * dPhi;
1103       numeratorPhiPhi += w * dPhi * dPhi;
1104     }
1105 
1106     //multiplying by crysSize to make the values compariable to normal covariances
1107     if (denominator != 0.0) {
1108       covEtaEta = crysSize * crysSize * numeratorEtaEta / denominator;
1109       covEtaPhi = crysSize * crysSize * numeratorEtaPhi / denominator;
1110       covPhiPhi = crysSize * crysSize * numeratorPhiPhi / denominator;
1111     } else {
1112       covEtaEta = 999.9;
1113       covEtaPhi = 999.9;
1114       covPhiPhi = 999.9;
1115     }
1116 
1117   } else {
1118     // Warn the user if there was no energy in the cells and return zeroes.
1119     //       std::cout << "\ClusterShapeAlgo::Calculate_Covariances:  no energy in supplied cells.\n";
1120     covEtaEta = 0;
1121     covEtaPhi = 0;
1122     covPhiPhi = 0;
1123   }
1124   std::array<float, 3> v{{covEtaEta, covEtaPhi, covPhiPhi}};
1125   return v;
1126 }
1127 
1128 template <bool noZS>
1129 double EcalClusterToolsT<noZS>::zernike20(const reco::BasicCluster &cluster,
1130                                           const EcalRecHitCollection *recHits,
1131                                           const CaloGeometry *geometry,
1132                                           double R0,
1133                                           bool logW,
1134                                           float w0) {
1135   return absZernikeMoment(cluster, recHits, geometry, 2, 0, R0, logW, w0);
1136 }
1137 
1138 template <bool noZS>
1139 double EcalClusterToolsT<noZS>::zernike42(const reco::BasicCluster &cluster,
1140                                           const EcalRecHitCollection *recHits,
1141                                           const CaloGeometry *geometry,
1142                                           double R0,
1143                                           bool logW,
1144                                           float w0) {
1145   return absZernikeMoment(cluster, recHits, geometry, 4, 2, R0, logW, w0);
1146 }
1147 
1148 template <bool noZS>
1149 double EcalClusterToolsT<noZS>::absZernikeMoment(const reco::BasicCluster &cluster,
1150                                                  const EcalRecHitCollection *recHits,
1151                                                  const CaloGeometry *geometry,
1152                                                  int n,
1153                                                  int m,
1154                                                  double R0,
1155                                                  bool logW,
1156                                                  float w0) {
1157   // 1. Check if n,m are correctly
1158   if ((m > n) || ((n - m) % 2 != 0) || (n < 0) || (m < 0))
1159     return -1;
1160 
1161   // 2. Check if n,R0 are within validity Range :
1162   // n>20 or R0<2.19cm  just makes no sense !
1163   if ((n > 20) || (R0 <= 2.19))
1164     return -1;
1165   if (n <= 5)
1166     return fast_AbsZernikeMoment(cluster, recHits, geometry, n, m, R0, logW, w0);
1167   else
1168     return calc_AbsZernikeMoment(cluster, recHits, geometry, n, m, R0, logW, w0);
1169 }
1170 
1171 template <bool noZS>
1172 double EcalClusterToolsT<noZS>::fast_AbsZernikeMoment(const reco::BasicCluster &cluster,
1173                                                       const EcalRecHitCollection *recHits,
1174                                                       const CaloGeometry *geometry,
1175                                                       int n,
1176                                                       int m,
1177                                                       double R0,
1178                                                       bool logW,
1179                                                       float w0) {
1180   double r, ph, e, Re = 0, Im = 0;
1181   double TotalEnergy = cluster.energy();
1182   int index = (n / 2) * (n / 2) + (n / 2) + m;
1183   std::vector<EcalClusterEnergyDeposition> energyDistribution =
1184       getEnergyDepTopology(cluster, recHits, geometry, logW, w0);
1185   int clusterSize = energyDistribution.size();
1186   if (clusterSize < 3)
1187     return 0.0;
1188 
1189   for (int i = 0; i < clusterSize; i++) {
1190     r = energyDistribution[i].r / R0;
1191     if (r < 1) {
1192       std::vector<double> pol;
1193       pol.push_back(f00(r));
1194       pol.push_back(f11(r));
1195       pol.push_back(f20(r));
1196       pol.push_back(f22(r));
1197       pol.push_back(f31(r));
1198       pol.push_back(f33(r));
1199       pol.push_back(f40(r));
1200       pol.push_back(f42(r));
1201       pol.push_back(f44(r));
1202       pol.push_back(f51(r));
1203       pol.push_back(f53(r));
1204       pol.push_back(f55(r));
1205       ph = (energyDistribution[i]).phi;
1206       e = energyDistribution[i].deposited_energy;
1207       Re = Re + e / TotalEnergy * pol[index] * cos((double)m * ph);
1208       Im = Im - e / TotalEnergy * pol[index] * sin((double)m * ph);
1209     }
1210   }
1211   return sqrt(Re * Re + Im * Im);
1212 }
1213 
1214 template <bool noZS>
1215 double EcalClusterToolsT<noZS>::calc_AbsZernikeMoment(const reco::BasicCluster &cluster,
1216                                                       const EcalRecHitCollection *recHits,
1217                                                       const CaloGeometry *geometry,
1218                                                       int n,
1219                                                       int m,
1220                                                       double R0,
1221                                                       bool logW,
1222                                                       float w0) {
1223   double r, ph, e, Re = 0, Im = 0, f_nm;
1224   double TotalEnergy = cluster.energy();
1225   std::vector<EcalClusterEnergyDeposition> energyDistribution =
1226       getEnergyDepTopology(cluster, recHits, geometry, logW, w0);
1227   int clusterSize = energyDistribution.size();
1228   if (clusterSize < 3)
1229     return 0.0;
1230 
1231   for (int i = 0; i < clusterSize; ++i) {
1232     r = energyDistribution[i].r / R0;
1233     if (r < 1) {
1234       ph = energyDistribution[i].phi;
1235       e = energyDistribution[i].deposited_energy;
1236       f_nm = 0;
1237       for (int s = 0; s <= (n - m) / 2; s++) {
1238         if (s % 2 == 0) {
1239           f_nm = f_nm + factorial(n - s) * pow(r, (double)(n - 2 * s)) /
1240                             (factorial(s) * factorial((n + m) / 2 - s) * factorial((n - m) / 2 - s));
1241         } else {
1242           f_nm = f_nm - factorial(n - s) * pow(r, (double)(n - 2 * s)) /
1243                             (factorial(s) * factorial((n + m) / 2 - s) * factorial((n - m) / 2 - s));
1244         }
1245       }
1246       Re = Re + e / TotalEnergy * f_nm * cos((double)m * ph);
1247       Im = Im - e / TotalEnergy * f_nm * sin((double)m * ph);
1248     }
1249   }
1250   return sqrt(Re * Re + Im * Im);
1251 }
1252 
1253 //returns the crystal 'eta' from the det id
1254 //it is defined as the number of crystals from the centre in the eta direction
1255 //for the barrel with its eta/phi geometry it is always integer
1256 //for the endcap it is fractional due to the x/y geometry
1257 template <bool noZS>
1258 float EcalClusterToolsT<noZS>::getIEta(const DetId &id) {
1259   if (id.det() == DetId::Ecal) {
1260     if (id.subdetId() == EcalBarrel) {
1261       EBDetId ebId(id);
1262       return ebId.ieta();
1263     } else if (id.subdetId() == EcalEndcap) {
1264       float iXNorm = getNormedIX(id);
1265       float iYNorm = getNormedIY(id);
1266 
1267       return std::sqrt(iXNorm * iXNorm + iYNorm * iYNorm);
1268     }
1269   }
1270   return 0.;
1271 }
1272 
1273 //returns the crystal 'phi' from the det id
1274 //it is defined as the number of crystals from the centre in the phi direction
1275 //for the barrel with its eta/phi geometry it is always integer
1276 //for the endcap it is not defined
1277 template <bool noZS>
1278 float EcalClusterToolsT<noZS>::getIPhi(const DetId &id) {
1279   if (id.det() == DetId::Ecal) {
1280     if (id.subdetId() == EcalBarrel) {
1281       EBDetId ebId(id);
1282       return ebId.iphi();
1283     }
1284   }
1285   return 0.;
1286 }
1287 
1288 //want to map 1=-50,50=-1,51=1 and 100 to 50 so sub off one if zero or neg
1289 template <bool noZS>
1290 float EcalClusterToolsT<noZS>::getNormedIX(const DetId &id) {
1291   if (id.det() == DetId::Ecal && id.subdetId() == EcalEndcap) {
1292     EEDetId eeId(id);
1293     int iXNorm = eeId.ix() - 50;
1294     if (iXNorm <= 0)
1295       iXNorm--;
1296     return iXNorm;
1297   }
1298   return 0;
1299 }
1300 
1301 //want to map 1=-50,50=-1,51=1 and 100 to 50 so sub off one if zero or neg
1302 template <bool noZS>
1303 float EcalClusterToolsT<noZS>::getNormedIY(const DetId &id) {
1304   if (id.det() == DetId::Ecal && id.subdetId() == EcalEndcap) {
1305     EEDetId eeId(id);
1306     int iYNorm = eeId.iy() - 50;
1307     if (iYNorm <= 0)
1308       iYNorm--;
1309     return iYNorm;
1310   }
1311   return 0;
1312 }
1313 
1314 //nr crystals crysId is away from orgin id in eta
1315 template <bool noZS>
1316 float EcalClusterToolsT<noZS>::getNrCrysDiffInEta(const DetId &crysId, const DetId &orginId) {
1317   float crysIEta = getIEta(crysId);
1318   float orginIEta = getIEta(orginId);
1319   bool isBarrel = orginId.subdetId() == EcalBarrel;
1320 
1321   float nrCrysDiff = crysIEta - orginIEta;
1322 
1323   //no iEta=0 in barrel, so if go from positive to negative
1324   //need to reduce abs(detEta) by 1
1325   if (isBarrel) {
1326     if (crysIEta * orginIEta < 0) {  // -1 to 1 transition
1327       if (crysIEta > 0)
1328         nrCrysDiff--;
1329       else
1330         nrCrysDiff++;
1331     }
1332   }
1333   return nrCrysDiff;
1334 }
1335 
1336 //nr crystals crysId is away from orgin id in phi
1337 template <bool noZS>
1338 float EcalClusterToolsT<noZS>::getNrCrysDiffInPhi(const DetId &crysId, const DetId &orginId) {
1339   float crysIPhi = getIPhi(crysId);
1340   float orginIPhi = getIPhi(orginId);
1341   bool isBarrel = orginId.subdetId() == EcalBarrel;
1342 
1343   float nrCrysDiff = crysIPhi - orginIPhi;
1344 
1345   if (isBarrel) {  //if barrel, need to map into 0-180
1346     if (nrCrysDiff > +180) {
1347       nrCrysDiff = nrCrysDiff - 360;
1348     }
1349     if (nrCrysDiff < -180) {
1350       nrCrysDiff = nrCrysDiff + 360;
1351     }
1352   }
1353   return nrCrysDiff;
1354 }
1355 
1356 //nr crystals crysId is away from mean phi in 5x5 in phi
1357 template <bool noZS>
1358 float EcalClusterToolsT<noZS>::getDPhiEndcap(const DetId &crysId, float meanX, float meanY) {
1359   float iXNorm = getNormedIX(crysId);
1360   float iYNorm = getNormedIY(crysId);
1361 
1362   float hitLocalR2 = (iXNorm - meanX) * (iXNorm - meanX) + (iYNorm - meanY) * (iYNorm - meanY);
1363   float hitR2 = iXNorm * iXNorm + iYNorm * iYNorm;
1364   float meanR2 = meanX * meanX + meanY * meanY;
1365   float hitR = sqrt(hitR2);
1366   float meanR = sqrt(meanR2);
1367 
1368   float tmp = (hitR2 + meanR2 - hitLocalR2) / (2 * hitR * meanR);
1369   if (tmp < -1)
1370     tmp = -1;
1371   if (tmp > 1)
1372     tmp = 1;
1373   float phi = acos(tmp);
1374   float dPhi = hitR * phi;
1375 
1376   return dPhi;
1377 }
1378 
1379 template <bool noZS>
1380 std::array<float, 3> EcalClusterToolsT<noZS>::scLocalCovariances(const reco::SuperCluster &cluster,
1381                                                                  const EcalRecHitCollection *recHits,
1382                                                                  const CaloTopology *topology,
1383                                                                  float w0) {
1384   const reco::BasicCluster bcluster = *(cluster.seed());
1385 
1386   float e_5x5 = e5x5(bcluster, recHits, topology);
1387   float covEtaEta, covEtaPhi, covPhiPhi;
1388 
1389   if (e_5x5 >= 0.) {
1390     const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
1391     std::pair<float, float> mean5x5PosInNrCrysFromSeed = mean5x5PositionInLocalCrysCoord(bcluster, recHits, topology);
1392     std::pair<float, float> mean5x5XYPos = mean5x5PositionInXY(cluster, recHits, topology);
1393     // now we can calculate the covariances
1394     double numeratorEtaEta = 0;
1395     double numeratorEtaPhi = 0;
1396     double numeratorPhiPhi = 0;
1397     double denominator = 0;
1398 
1399     const double barrelCrysSize = 0.01745;  //approximate size of crystal in eta,phi in barrel
1400     const double endcapCrysSize = 0.0447;   //the approximate crystal size sigmaEtaEta was corrected to in the endcap
1401 
1402     DetId seedId = getMaximum(v_id, recHits).first;
1403     bool isBarrel = seedId.subdetId() == EcalBarrel;
1404 
1405     const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1406 
1407     for (size_t i = 0; i < v_id.size(); ++i) {
1408       float frac = getFraction(v_id, v_id[i].first);
1409       float energy = recHitEnergy(v_id[i].first, recHits) * frac;
1410 
1411       if (energy <= 0)
1412         continue;
1413 
1414       float dEta = getNrCrysDiffInEta(v_id[i].first, seedId) - mean5x5PosInNrCrysFromSeed.first;
1415       float dPhi = 0;
1416       if (isBarrel)
1417         dPhi = getNrCrysDiffInPhi(v_id[i].first, seedId) - mean5x5PosInNrCrysFromSeed.second;
1418       else
1419         dPhi = getDPhiEndcap(v_id[i].first, mean5x5XYPos.first, mean5x5XYPos.second);
1420 
1421       double w = 0.;
1422       w = std::max(0.0f, w0 + std::log(energy / e_5x5));
1423 
1424       denominator += w;
1425       numeratorEtaEta += w * dEta * dEta;
1426       numeratorEtaPhi += w * dEta * dPhi;
1427       numeratorPhiPhi += w * dPhi * dPhi;
1428     }
1429 
1430     //multiplying by crysSize to make the values compariable to normal covariances
1431     if (denominator != 0.0) {
1432       covEtaEta = crysSize * crysSize * numeratorEtaEta / denominator;
1433       covEtaPhi = crysSize * crysSize * numeratorEtaPhi / denominator;
1434       covPhiPhi = crysSize * crysSize * numeratorPhiPhi / denominator;
1435     } else {
1436       covEtaEta = 999.9;
1437       covEtaPhi = 999.9;
1438       covPhiPhi = 999.9;
1439     }
1440 
1441   } else {
1442     // Warn the user if there was no energy in the cells and return zeroes.
1443     // std::cout << "\ClusterShapeAlgo::Calculate_Covariances:  no energy in supplied cells.\n";
1444     covEtaEta = 0;
1445     covEtaPhi = 0;
1446     covPhiPhi = 0;
1447   }
1448 
1449   std::array<float, 3> v{{covEtaEta, covEtaPhi, covPhiPhi}};
1450   return v;
1451 }
1452 
1453 // compute cluster second moments with respect to principal axes (eigenvectors of sEtaEta, sPhiPhi, sEtaPhi matrix)
1454 // store also angle alpha between major axis and phi.
1455 // takes into account shower elongation in phi direction due to magnetic field effect:
1456 // default value of 0.8 ensures sMaj = sMin for unconverted photons
1457 // (if phiCorrectionFactor=1 sMaj > sMin and alpha=0 also for unconverted photons)
1458 template <bool noZS>
1459 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments(const reco::BasicCluster &basicCluster,
1460                                                              const EcalRecHitCollection &recHits,
1461                                                              double phiCorrectionFactor,
1462                                                              double w0,
1463                                                              bool useLogWeights) {
1464   if (recHits.empty())
1465     return Cluster2ndMoments();
1466 
1467   std::vector<std::pair<const EcalRecHit *, float>> RH_ptrs_fracs;
1468 
1469   const std::vector<std::pair<DetId, float>> &myHitsPair = basicCluster.hitsAndFractions();
1470 
1471   for (unsigned int i = 0; i < myHitsPair.size(); i++) {
1472     //get pointer to recHit object
1473     EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
1474     if (myRH != recHits.end()) {
1475       RH_ptrs_fracs.push_back(std::make_pair(&(*myRH), myHitsPair[i].second));
1476     }
1477   }
1478 
1479   return EcalClusterToolsT<noZS>::cluster2ndMoments(RH_ptrs_fracs, phiCorrectionFactor, w0, useLogWeights);
1480 }
1481 
1482 template <bool noZS>
1483 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments(const reco::SuperCluster &superCluster,
1484                                                              const EcalRecHitCollection &recHits,
1485                                                              double phiCorrectionFactor,
1486                                                              double w0,
1487                                                              bool useLogWeights) {
1488   return EcalClusterToolsT<noZS>::cluster2ndMoments(
1489       *(superCluster.seed()), recHits, phiCorrectionFactor, w0, useLogWeights);
1490 }
1491 
1492 template <bool noZS>
1493 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments(
1494     const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs,
1495     double phiCorrectionFactor,
1496     double w0,
1497     bool useLogWeights) {
1498   if (RH_ptrs_fracs.empty())
1499     return Cluster2ndMoments();
1500 
1501   double mid_eta(0), mid_phi(0), mid_x(0), mid_y(0);
1502 
1503   double Etot = EcalClusterToolsT<noZS>::getSumEnergy(RH_ptrs_fracs);
1504 
1505   double max_phi = -10.;
1506   double min_phi = 100.;
1507 
1508   std::vector<double> etaDetId;
1509   std::vector<double> phiDetId;
1510   std::vector<double> xDetId;
1511   std::vector<double> yDetId;
1512   std::vector<double> wiDetId;
1513 
1514   unsigned int nCry = 0;
1515   double denominator = 0.;
1516   bool isBarrel(true);
1517 
1518   // loop over rechits and compute weights:
1519   for (std::vector<std::pair<const EcalRecHit *, float>>::const_iterator rhf_ptr = RH_ptrs_fracs.begin();
1520        rhf_ptr != RH_ptrs_fracs.end();
1521        rhf_ptr++) {
1522     const EcalRecHit *rh_ptr = rhf_ptr->first;
1523 
1524     //get iEta, iPhi
1525     double temp_eta(0), temp_phi(0), temp_x(0), temp_y(0);
1526     isBarrel = rh_ptr->detid().subdetId() == EcalBarrel;
1527 
1528     if (isBarrel) {
1529       temp_eta = (getIEta(rh_ptr->detid()) > 0. ? getIEta(rh_ptr->detid()) + 84.5 : getIEta(rh_ptr->detid()) + 85.5);
1530       temp_phi = getIPhi(rh_ptr->detid()) - 0.5;
1531     } else {
1532       temp_eta = getIEta(rh_ptr->detid());
1533       temp_x = getNormedIX(rh_ptr->detid());
1534       temp_y = getNormedIY(rh_ptr->detid());
1535     }
1536 
1537     double temp_ene = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
1538 
1539     double temp_wi = ((useLogWeights) ? std::max(0.0, w0 + std::log(std::abs(temp_ene) / Etot)) : temp_ene);
1540 
1541     if (temp_phi > max_phi)
1542       max_phi = temp_phi;
1543     if (temp_phi < min_phi)
1544       min_phi = temp_phi;
1545     etaDetId.push_back(temp_eta);
1546     phiDetId.push_back(temp_phi);
1547     xDetId.push_back(temp_x);
1548     yDetId.push_back(temp_y);
1549     wiDetId.push_back(temp_wi);
1550     denominator += temp_wi;
1551     nCry++;
1552   }
1553 
1554   if (isBarrel) {
1555     // correct phi wrap-around:
1556     if (max_phi == 359.5 && min_phi == 0.5) {
1557       for (unsigned int i = 0; i < nCry; i++) {
1558         if (phiDetId[i] - 179. > 0.)
1559           phiDetId[i] -= 360.;
1560         mid_phi += phiDetId[i] * wiDetId[i];
1561         mid_eta += etaDetId[i] * wiDetId[i];
1562       }
1563     } else {
1564       for (unsigned int i = 0; i < nCry; i++) {
1565         mid_phi += phiDetId[i] * wiDetId[i];
1566         mid_eta += etaDetId[i] * wiDetId[i];
1567       }
1568     }
1569   } else {
1570     for (unsigned int i = 0; i < nCry; i++) {
1571       mid_eta += etaDetId[i] * wiDetId[i];
1572       mid_x += xDetId[i] * wiDetId[i];
1573       mid_y += yDetId[i] * wiDetId[i];
1574     }
1575   }
1576 
1577   mid_eta /= denominator;
1578   mid_phi /= denominator;
1579   mid_x /= denominator;
1580   mid_y /= denominator;
1581 
1582   // See = sigma eta eta
1583   // Spp = (B field corrected) sigma phi phi
1584   // See = (B field corrected) sigma eta phi
1585   double See = 0.;
1586   double Spp = 0.;
1587   double Sep = 0.;
1588   double deta(0), dphi(0);
1589   // compute (phi-corrected) covariance matrix:
1590   for (unsigned int i = 0; i < nCry; i++) {
1591     if (isBarrel) {
1592       deta = etaDetId[i] - mid_eta;
1593       dphi = phiDetId[i] - mid_phi;
1594     } else {
1595       deta = etaDetId[i] - mid_eta;
1596       float hitLocalR2 = (xDetId[i] - mid_x) * (xDetId[i] - mid_x) + (yDetId[i] - mid_y) * (yDetId[i] - mid_y);
1597       float hitR2 = xDetId[i] * xDetId[i] + yDetId[i] * yDetId[i];
1598       float meanR2 = mid_x * mid_x + mid_y * mid_y;
1599       float hitR = sqrt(hitR2);
1600       float meanR = sqrt(meanR2);
1601       float phi = acos((hitR2 + meanR2 - hitLocalR2) / (2 * hitR * meanR));
1602       dphi = hitR * phi;
1603     }
1604     See += (wiDetId[i] * deta * deta) / denominator;
1605     Spp += phiCorrectionFactor * (wiDetId[i] * dphi * dphi) / denominator;
1606     Sep += sqrt(phiCorrectionFactor) * (wiDetId[i] * deta * dphi) / denominator;
1607   }
1608 
1609   Cluster2ndMoments returnMoments;
1610 
1611   // compute matrix eigenvalues:
1612   returnMoments.sMaj = ((See + Spp) + sqrt((See - Spp) * (See - Spp) + 4. * Sep * Sep)) / 2.;
1613   returnMoments.sMin = ((See + Spp) - sqrt((See - Spp) * (See - Spp) + 4. * Sep * Sep)) / 2.;
1614 
1615   returnMoments.alpha = atan((See - Spp + sqrt((Spp - See) * (Spp - See) + 4. * Sep * Sep)) / (2. * Sep));
1616 
1617   return returnMoments;
1618 }
1619 
1620 //compute shower shapes: roundness and angle in a vector. Roundness is 0th element, Angle is 1st element.
1621 //description: uses classical mechanics inertia tensor.
1622 //             roundness is smaller_eValue/larger_eValue
1623 //             angle is the angle from the iEta axis to the smallest eVector (a.k.a. the shower's elongated axis)
1624 // this function uses only recHits belonging to a SC above energyThreshold (default 0)
1625 // you can select linear weighting = energy_recHit/total_energy         (weightedPositionMethod=0) default
1626 // or log weighting = max( 0.0, 4.2 + log(energy_recHit/total_energy) ) (weightedPositionMethod=1)
1627 template <bool noZS>
1628 std::vector<float> EcalClusterToolsT<noZS>::roundnessBarrelSuperClusters(
1629     const reco::SuperCluster &superCluster,
1630     const EcalRecHitCollection &recHits,
1631     int weightedPositionMethod,
1632     float energyThreshold) {  //int positionWeightingMethod=0){
1633   std::vector<std::pair<const EcalRecHit *, float>> RH_ptrs_fracs;
1634   const std::vector<std::pair<DetId, float>> &myHitsPair = superCluster.hitsAndFractions();
1635   for (unsigned int i = 0; i < myHitsPair.size(); ++i) {
1636     //get pointer to recHit object
1637     EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
1638     if (myRH != recHits.end() && myRH->energy() * (noZS ? 1.0 : myHitsPair[i].second) > energyThreshold) {
1639       //require rec hit to have positive energy
1640       RH_ptrs_fracs.push_back(std::make_pair(&(*myRH), myHitsPair[i].second));
1641     }
1642   }
1643   std::vector<float> temp =
1644       EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(RH_ptrs_fracs, weightedPositionMethod);
1645   return temp;
1646 }
1647 
1648 // this function uses all recHits within specified window ( with default values ieta_delta=2, iphi_delta=5) around SC's highest recHit.
1649 // recHits must pass an energy threshold "energyRHThresh" (default 0)
1650 // you can select linear weighting = energy_recHit/total_energy         (weightedPositionMethod=0)
1651 // or log weighting = max( 0.0, 4.2 + log(energy_recHit/total_energy) ) (weightedPositionMethod=1)
1652 template <bool noZS>
1653 std::vector<float> EcalClusterToolsT<noZS>::roundnessBarrelSuperClustersUserExtended(
1654     const reco::SuperCluster &superCluster,
1655     const EcalRecHitCollection &recHits,
1656     int ieta_delta,
1657     int iphi_delta,
1658     float energyRHThresh,
1659     int weightedPositionMethod) {
1660   std::vector<std::pair<const EcalRecHit *, float>> RH_ptrs_fracs;
1661   const std::vector<std::pair<DetId, float>> &myHitsPair = superCluster.hitsAndFractions();
1662   for (unsigned int i = 0; i < myHitsPair.size(); ++i) {
1663     //get pointer to recHit object
1664     EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
1665     if (myRH != recHits.end() && myRH->energy() * (noZS ? 1.0 : myHitsPair[i].second) > energyRHThresh)
1666       RH_ptrs_fracs.push_back(std::make_pair(&(*myRH), myHitsPair[i].second));
1667   }
1668 
1669   std::vector<int> seedPosition = EcalClusterToolsT<noZS>::getSeedPosition(RH_ptrs_fracs);
1670 
1671   for (EcalRecHitCollection::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++) {
1672     EBDetId EBdetIdi(rh->detid());
1673     float the_fraction = 0;
1674     //if(rh != recHits.end())
1675     bool inEtaWindow = (abs(deltaIEta(seedPosition[0], EBdetIdi.ieta())) <= ieta_delta);
1676     bool inPhiWindow = (abs(deltaIPhi(seedPosition[1], EBdetIdi.iphi())) <= iphi_delta);
1677     bool passEThresh = (rh->energy() > energyRHThresh);
1678     bool alreadyCounted = false;
1679 
1680     // figure out if the rechit considered now is already inside the SC
1681     bool is_SCrh_inside_recHits = false;
1682     for (unsigned int i = 0; i < myHitsPair.size(); i++) {
1683       EcalRecHitCollection::const_iterator SCrh = recHits.find(myHitsPair[i].first);
1684       if (SCrh != recHits.end()) {
1685         the_fraction = myHitsPair[i].second;
1686         is_SCrh_inside_recHits = true;
1687         if (rh->detid() == SCrh->detid())
1688           alreadyCounted = true;
1689       }
1690     }  //for loop over SC's recHits
1691 
1692     if (is_SCrh_inside_recHits && !alreadyCounted && passEThresh && inEtaWindow && inPhiWindow) {
1693       RH_ptrs_fracs.push_back(std::make_pair(&(*rh), the_fraction));
1694     }
1695 
1696   }  //for loop over rh
1697   return EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(RH_ptrs_fracs, weightedPositionMethod);
1698 }
1699 
1700 // this function computes the roundness and angle variables for vector of pointers to recHits you pass it
1701 // you can select linear weighting = energy_recHit/total_energy         (weightedPositionMethod=0)
1702 // or log weighting = max( 0.0, 4.2 + log(energy_recHit/total_energy) ) (weightedPositionMethod=1)
1703 template <bool noZS>
1704 std::vector<float> EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(
1705     const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs,
1706     int weightedPositionMethod) {  //int weightedPositionMethod = 0){
1707   //positionWeightingMethod = 0 linear weighting, 1 log energy weighting
1708 
1709   std::vector<float> shapes;  // this is the returning vector
1710 
1711   //make sure photon has more than one crystal; else roundness and angle suck
1712   if (RH_ptrs_fracs.size() < 2) {
1713     shapes.push_back(-3);
1714     shapes.push_back(-3);
1715     return shapes;
1716   }
1717 
1718   //Find highest E RH (Seed) and save info, compute sum total energy used
1719   std::vector<int> seedPosition = EcalClusterToolsT<noZS>::getSeedPosition(RH_ptrs_fracs);  // *recHits);
1720   int tempInt = seedPosition[0];
1721   if (tempInt < 0)
1722     tempInt++;
1723   float energyTotal = EcalClusterToolsT<noZS>::getSumEnergy(RH_ptrs_fracs);
1724 
1725   //1st loop over rechits: compute new weighted center position in coordinates centered on seed
1726   float centerIEta = 0.;
1727   float centerIPhi = 0.;
1728   float denominator = 0.;
1729 
1730   for (std::vector<std::pair<const EcalRecHit *, float>>::const_iterator rhf_ptr = RH_ptrs_fracs.begin();
1731        rhf_ptr != RH_ptrs_fracs.end();
1732        rhf_ptr++) {
1733     const EcalRecHit *rh_ptr = rhf_ptr->first;
1734     //get iEta, iPhi
1735     EBDetId EBdetIdi(rh_ptr->detid());
1736     if (fabs(energyTotal) < 0.0001) {
1737       // don't /0, bad!
1738       shapes.push_back(-2);
1739       shapes.push_back(-2);
1740       return shapes;
1741     }
1742     float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
1743     float weight = 0;
1744     if (std::abs(weightedPositionMethod) < 0.0001) {  //linear
1745       weight = rh_energy / energyTotal;
1746     } else {  //logrithmic
1747       weight = std::max(0.0, 4.2 + log(rh_energy / energyTotal));
1748     }
1749     denominator += weight;
1750     centerIEta += weight * deltaIEta(seedPosition[0], EBdetIdi.ieta());
1751     centerIPhi += weight * deltaIPhi(seedPosition[1], EBdetIdi.iphi());
1752   }
1753   if (fabs(denominator) < 0.0001) {
1754     // don't /0, bad!
1755     shapes.push_back(-2);
1756     shapes.push_back(-2);
1757     return shapes;
1758   }
1759   centerIEta = centerIEta / denominator;
1760   centerIPhi = centerIPhi / denominator;
1761 
1762   //2nd loop over rechits: compute inertia tensor
1763   TMatrixDSym inertia(2);  //initialize 2d inertia tensor
1764   double inertia00 = 0.;
1765   double inertia01 = 0.;  // = inertia10 b/c matrix should be symmetric
1766   double inertia11 = 0.;
1767   int i = 0;
1768   for (std::vector<std::pair<const EcalRecHit *, float>>::const_iterator rhf_ptr = RH_ptrs_fracs.begin();
1769        rhf_ptr != RH_ptrs_fracs.end();
1770        rhf_ptr++) {
1771     const EcalRecHit *rh_ptr = rhf_ptr->first;
1772     //get iEta, iPhi
1773     EBDetId EBdetIdi(rh_ptr->detid());
1774 
1775     if (fabs(energyTotal) < 0.0001) {
1776       // don't /0, bad!
1777       shapes.push_back(-2);
1778       shapes.push_back(-2);
1779       return shapes;
1780     }
1781     float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
1782     float weight = 0;
1783     if (std::abs(weightedPositionMethod) < 0.0001) {  //linear
1784       weight = rh_energy / energyTotal;
1785     } else {  //logrithmic
1786       weight = std::max(0.0, 4.2 + log(rh_energy / energyTotal));
1787     }
1788 
1789     float ieta_rh_to_center = deltaIEta(seedPosition[0], EBdetIdi.ieta()) - centerIEta;
1790     float iphi_rh_to_center = deltaIPhi(seedPosition[1], EBdetIdi.iphi()) - centerIPhi;
1791 
1792     inertia00 += weight * iphi_rh_to_center * iphi_rh_to_center;
1793     inertia01 -= weight * iphi_rh_to_center * ieta_rh_to_center;
1794     inertia11 += weight * ieta_rh_to_center * ieta_rh_to_center;
1795     i++;
1796   }
1797 
1798   inertia[0][0] = inertia00;
1799   inertia[0][1] = inertia01;  // use same number here
1800   inertia[1][0] = inertia01;  // and here to insure symmetry
1801   inertia[1][1] = inertia11;
1802 
1803   //step 1 find principal axes of inertia
1804   TMatrixD eVectors(2, 2);
1805   TVectorD eValues(2);
1806   //std::cout<<"EcalClusterToolsT<noZS>::showerRoundness- about to compute eVectors"<<std::endl;
1807   eVectors = inertia.EigenVectors(eValues);  //ordered highest eV to lowest eV (I checked!)
1808   //and eVectors are in columns of matrix! I checked!
1809   //and they are normalized to 1
1810 
1811   //step 2 select eta component of smaller eVal's eVector
1812   TVectorD smallerAxis(2);          //easiest to spin SC on this axis (smallest eVal)
1813   smallerAxis[0] = eVectors[0][1];  //row,col  //eta component
1814   smallerAxis[1] = eVectors[1][1];  //phi component
1815 
1816   //step 3 compute interesting quatities
1817   Double_t temp = fabs(smallerAxis[0]);  // closer to 1 ->beamhalo, closer to 0 something else
1818   if (fabs(eValues[0]) < 0.0001) {
1819     // don't /0, bad!
1820     shapes.push_back(-2);
1821     shapes.push_back(-2);
1822     return shapes;
1823   }
1824 
1825   float Roundness = eValues[1] / eValues[0];
1826   float Angle = acos(temp);
1827 
1828   if (-0.00001 < Roundness && Roundness < 0)
1829     Roundness = 0.;
1830   if (-0.00001 < Angle && Angle < 0)
1831     Angle = 0.;
1832 
1833   shapes.push_back(Roundness);
1834   shapes.push_back(Angle);
1835   return shapes;
1836 }
1837 
1838 template <bool noZS>
1839 int EcalClusterToolsT<noZS>::nrSaturatedCrysIn5x5(const DetId &id,
1840                                                   const EcalRecHitCollection *recHits,
1841                                                   const CaloTopology *topology) {
1842   int nrSat = 0;
1843   CaloRectangle rectangle{-2, 2, -2, 2};
1844   for (auto const &detId : rectangle(id, *topology)) {
1845     auto recHitIt = recHits->find(detId);
1846     if (recHitIt != recHits->end() && recHitIt->checkFlag(EcalRecHit::kSaturated)) {
1847       nrSat++;
1848     }
1849   }
1850   return nrSat;
1851 }
1852 
1853 //private functions useful for roundnessBarrelSuperClusters etc.
1854 //compute delta iphi between a seed and a particular recHit
1855 //iphi [1,360]
1856 //safe gaurds are put in to ensure the difference is between [-180,180]
1857 template <bool noZS>
1858 int EcalClusterToolsT<noZS>::deltaIPhi(int seed_iphi, int rh_iphi) {
1859   int rel_iphi = rh_iphi - seed_iphi;
1860   // take care of cyclic variable iphi [1,360]
1861   if (rel_iphi > 180)
1862     rel_iphi = rel_iphi - 360;
1863   if (rel_iphi < -180)
1864     rel_iphi = rel_iphi + 360;
1865   return rel_iphi;
1866 }
1867 
1868 //compute delta ieta between a seed and a particular recHit
1869 //ieta [-85,-1] and [1,85]
1870 //safe gaurds are put in to shift the negative ieta +1 to make an ieta=0 so differences are computed correctly
1871 template <bool noZS>
1872 int EcalClusterToolsT<noZS>::deltaIEta(int seed_ieta, int rh_ieta) {
1873   // get rid of the fact that there is no ieta=0
1874   if (seed_ieta < 0)
1875     seed_ieta++;
1876   if (rh_ieta < 0)
1877     rh_ieta++;
1878   int rel_ieta = rh_ieta - seed_ieta;
1879   return rel_ieta;
1880 }
1881 
1882 //return (ieta,iphi) of highest energy recHit of the recHits passed to this function
1883 template <bool noZS>
1884 std::vector<int> EcalClusterToolsT<noZS>::getSeedPosition(
1885     const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs) {
1886   std::vector<int> seedPosition;
1887   float eSeedRH = 0;
1888   int iEtaSeedRH = 0;
1889   int iPhiSeedRH = 0;
1890 
1891   for (auto const &rhf : RH_ptrs_fracs) {
1892     const EcalRecHit *rh_ptr = rhf.first;
1893     //get iEta, iPhi
1894     EBDetId EBdetIdi(rh_ptr->detid());
1895     float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf.second);
1896 
1897     if (eSeedRH < rh_energy) {
1898       eSeedRH = rh_energy;
1899       iEtaSeedRH = EBdetIdi.ieta();
1900       iPhiSeedRH = EBdetIdi.iphi();
1901     }
1902 
1903   }  // for loop
1904 
1905   seedPosition.push_back(iEtaSeedRH);
1906   seedPosition.push_back(iPhiSeedRH);
1907   return seedPosition;
1908 }
1909 
1910 // return the total energy of the recHits passed to this function
1911 template <bool noZS>
1912 float EcalClusterToolsT<noZS>::getSumEnergy(const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs) {
1913   float sumE = 0.;
1914   for (const auto &hAndF : RH_ptrs_fracs) {
1915     sumE += hAndF.first->energy() * (noZS ? 1.0 : hAndF.second);
1916   }
1917   return sumE;
1918 }
1919 
1920 typedef EcalClusterToolsT<false> EcalClusterTools;
1921 
1922 namespace noZS {
1923   typedef EcalClusterToolsT<true> EcalClusterTools;
1924 }
1925 
1926 #endif