CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EcalClusterTools.h
Go to the documentation of this file.
1 #ifndef RecoEcal_EgammaCoreTools_EcalClusterTools_h
2 #define RecoEcal_EgammaCoreTools_EcalClusterTools_h
3 
30 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
31 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
32 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
33 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
34 //#include "DataFormats/Math/interface/Point3D.h"
35 #include "DataFormats/Math/interface/Vector3D.h"
36 //includes for ShowerShape function to work
37 #include <vector>
38 #include <math.h>
39 #include <TMath.h>
40 #include <TMatrixT.h>
41 #include <TMatrixD.h>
42 #include <TVectorT.h>
43 #include <TVectorD.h>
44 
45 #include "DataFormats/DetId/interface/DetId.h"
46 #include "DataFormats/EcalDetId/interface/EBDetId.h"
47 #include "DataFormats/EcalDetId/interface/EEDetId.h"
48 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
49 
50 #include "FWCore/MessageLogger/interface/MessageLogger.h"
51 
52 
53 #include "CLHEP/Geometry/Transform3D.h"
54 
55 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
56 #include "Geometry/CaloTopology/interface/CaloTopology.h"
57 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
58 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
59 
60 
61 class DetId;
62 class CaloTopology;
63 class CaloGeometry;
64 
65 struct Cluster2ndMoments {
66 
67  // major and minor cluster moments wrt principale axes:
68  float sMaj;
69  float sMin;
70  // angle between sMaj and phi:
71  float alpha;
72 
73 };
74 
75 template<bool noZS>
76 class EcalClusterToolsT {
77  public:
78  EcalClusterToolsT() {};
79  ~EcalClusterToolsT() {};
80 
81  // various energies in the matrix nxn surrounding the maximum energy crystal of the input cluster
82  //we use an eta/phi coordinate system rather than phi/eta
83  //note e3x2 does not have a definate eta/phi geometry, it takes the maximum 3x2 block containing the
84  //seed regardless of whether that 3 in eta or phi
85  static float e1x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
86 
87 static float e1x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
88 
89  static float e3x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
90 
91 static float e3x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
92 
93  static float e1x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
94  //AA
95  //Take into account severities
96  static float e1x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
97 
98  static float e5x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
99 static float e5x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
100 
101  static float e2x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
102 static float e2x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
103 
104  static float e3x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
105 static float e3x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
106 
107  static float e3x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
108 static float e3x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
109 
110  static float e4x4( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology);
111 static float e4x4( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
112 
113  static float e5x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
114 static float e5x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
115 
116  // energy in the 2x5 strip right of the max crystal (does not contain max crystal)
117  // 2 crystals wide in eta, 5 wide in phi.
118  static float e2x5Right( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
119 static float e2x5Right( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
120  // energy in the 2x5 strip left of the max crystal (does not contain max crystal)
121 
122  static float e2x5Left( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
123 static float e2x5Left( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
124  // energy in the 5x2 strip above the max crystal (does not contain max crystal)
125  // 5 crystals wide in eta, 2 wide in phi.
126 
127  static float e2x5Top( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
128 static float e2x5Top( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
129  // energy in the 5x2 strip below the max crystal (does not contain max crystal)
130 
131  static float e2x5Bottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
132 static float e2x5Bottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
133  // energy in a 2x5 strip containing the seed (max) crystal.
134  // 2 crystals wide in eta, 5 wide in phi.
135  // it is the maximum of either (1x5left + 1x5center) or (1x5right + 1x5center)
136  static float e2x5Max( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
137 static float e2x5Max( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
138 
139  // energies in the crystal left, right, top, bottom w.r.t. to the most energetic crystal
140  static float eLeft( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
141 static float eLeft( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
142 
143  static float eRight( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
144 static float eRight( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
145 
146  static float eTop( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
147 static float eTop( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
148 
149  static float eBottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
150 static float eBottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
151  // the energy of the most energetic crystal in the cluster
152 
153  static float eMax( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits );
154 static float eMax( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
155 
156  // the energy of the second most energetic crystal in the cluster
157  static float e2nd( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits );
158 static float e2nd( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
159 
160  // get the DetId and the energy of the maximum energy crystal of the input cluster
161  static std::pair<DetId, float> getMaximum( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
162 static std::pair<DetId, float> getMaximum( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
163 
164  static std::vector<float> energyBasketFractionEta( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits );
165 static std::vector<float> energyBasketFractionEta( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
166 
167  static std::vector<float> energyBasketFractionPhi( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
168 static std::vector<float> energyBasketFractionPhi( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
169 
170  // return a vector v with v[0] = etaLat, v[1] = phiLat, v[2] = lat
171  static std::vector<float> lat( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW = true, float w0 = 4.7 );
172 static std::vector<float> lat( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv, bool logW = true,float w0 = 4.7);
173 
174  // return a vector v with v[0] = covEtaEta, v[1] = covEtaPhi, v[2] = covPhiPhi
175 
176  static std::vector<float> covariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, const CaloGeometry* geometry, float w0 = 4.7);
177 static std::vector<float> covariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, const CaloGeometry* geometry,const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv,float w0 = 4.7);
178  // return a vector v with v[0] = covIEtaIEta, v[1] = covIEtaIPhi, v[2] = covIPhiIPhi
179  //this function calculates differences in eta/phi in units of crystals not global eta/phi
180  //this is gives better performance in the crack regions of the calorimeter but gives otherwise identical results to covariances function
181  // except that it doesnt need an eta based correction funtion in the endcap
182  //it is multipled by an approprate crystal size to ensure it gives similar values to covariances(...)
183  //
184  //Warning: covIEtaIEta has been studied by egamma, but so far covIPhiIPhi hasnt been studied extensively so there could be a bug in
185  // the covIPhiIEta or covIPhiIPhi calculations. I dont think there is but as it hasnt been heavily used, there might be one
186  static std::vector<float> localCovariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, float w0 = 4.7);
187 static std::vector<float> localCovariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv,float w0 = 4.7);
188 
189  static std::vector<float> scLocalCovariances(const reco::SuperCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, float w0 = 4.7);
190 static std::vector<float> scLocalCovariances(const reco::SuperCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv, float w0 = 4.7);
191 
192  // cluster second moments with respect to principal axes:
193  static Cluster2ndMoments cluster2ndMoments( const reco::BasicCluster &basicCluster, const EcalRecHitCollection &recHits, double phiCorrectionFactor=0.8, double w0=4.7, bool useLogWeights=true);
194 
195  static Cluster2ndMoments cluster2ndMoments( const reco::SuperCluster &superCluster, const EcalRecHitCollection &recHits, double phiCorrectionFactor=0.8, double w0=4.7, bool useLogWeights=true);
196  static Cluster2ndMoments cluster2ndMoments( const std::vector<std::pair<const EcalRecHit*, float> >& RH_ptrs_fracs, double phiCorrectionFactor=0.8, double w0=4.7, bool useLogWeights=true);
197 
198  static double zernike20( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, double R0 = 6.6, bool logW = true, float w0 = 4.7 );
199  static double zernike42( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, double R0 = 6.6, bool logW = true, float w0 = 4.7 );
200 
201  // get the detId's of a matrix centered in the maximum energy crystal = (0,0)
202  // the size is specified by ixMin, ixMax, iyMin, iyMax in unit of crystals
203  static std::vector<DetId> matrixDetId( const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax );
204 static std::vector<DetId> matrixDetId( const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
205  // get the energy deposited in a matrix centered in the maximum energy crystal = (0,0)
206  // the size is specified by ixMin, ixMax, iyMin, iyMax in unit of crystals
207  static float matrixEnergy( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax );
208 
209  //AA
210  //Take into account the severities and flags
211  static float matrixEnergy( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
212  static float getFraction( const std::vector< std::pair<DetId, float> > &v_id, DetId id);
213  // get the DetId and the energy of the maximum energy crystal in a vector of DetId
214  static std::pair<DetId, float> getMaximum( const std::vector< std::pair<DetId, float> > &v_id, const EcalRecHitCollection *recHits);
215 static std::pair<DetId, float> getMaximum( const std::vector< std::pair<DetId, float> > &v_id, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
216  // get the energy of a DetId, return 0 if the DetId is not in the collection
217  static float recHitEnergy(DetId id, const EcalRecHitCollection *recHits);
218 
219  //AA
220  //Take into account severities and flags
221  static float recHitEnergy(DetId id, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
222 
223  //Shower shape variables return vector <Roundness, Angle> of a photon
224  static std::vector<float> roundnessBarrelSuperClusters( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int weightedPositionMethod = 0, float energyThreshold = 0.0);
225  static std::vector<float> roundnessBarrelSuperClustersUserExtended( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int ieta_delta=0, int iphi_delta=0, float energyRHThresh=0.00000, int weightedPositionMethod=0);
226  static std::vector<float> roundnessSelectedBarrelRecHits(const std::vector<std::pair<const EcalRecHit*,float> >&rhVector, int weightedPositionMethod = 0);
227  private:
228  struct EcalClusterEnergyDeposition
229  {
230  double deposited_energy;
231  double r;
232  double phi;
233  };
234 
235  static std::vector<EcalClusterEnergyDeposition> getEnergyDepTopology( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW, float w0 );
236 
237  static math::XYZVector meanClusterPosition( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, const CaloGeometry *geometry );
238 static math::XYZVector meanClusterPosition( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, const CaloGeometry *geometry, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv );
239  //return energy weighted mean distance from the seed crystal in number of crystals
240  //<iEta,iPhi>, iPhi is not defined for endcap and is returned as zero
241  static std::pair<float,float> mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology);
242 static std::pair<float,float> mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
243 
244  static std::pair<float,float> mean5x5PositionInXY(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology);
245 static std::pair<float,float> mean5x5PositionInXY(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv);
246 
247  static double f00(double r) { return 1; }
248  static double f11(double r) { return r; }
249  static double f20(double r) { return 2.0*r*r-1.0; }
250  static double f22(double r) { return r*r; }
251  static double f31(double r) { return 3.0*r*r*r - 2.0*r; }
252  static double f33(double r) { return r*r*r; }
253  static double f40(double r) { return 6.0*r*r*r*r-6.0*r*r+1.0; }
254  static double f42(double r) { return 4.0*r*r*r*r-3.0*r*r; }
255  static double f44(double r) { return r*r*r*r; }
256  static double f51(double r) { return 10.0*pow(r,5)-12.0*pow(r,3)+3.0*r; }
257  static double f53(double r) { return 5.0*pow(r,5) - 4.0*pow(r,3); }
258  static double f55(double r) { return pow(r,5); }
259 
260  static double absZernikeMoment( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 );
261  static double fast_AbsZernikeMoment(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 );
262  static double calc_AbsZernikeMoment(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 );
263 
264  static double factorial(int n) {
265  double res = 1.;
266  for (int i = 2; i <= n; ++i) res *= i;
267  return res;
268  }
269 
270  //useful functions for the localCovariances function
271  static float getIEta(const DetId& id);
272  static float getIPhi(const DetId& id);
273  static float getNormedIX(const DetId& id);
274  static float getNormedIY(const DetId& id);
275  static float getDPhiEndcap(const DetId& crysId,float meanX,float meanY);
276  static float getNrCrysDiffInEta(const DetId& crysId,const DetId& orginId);
277  static float getNrCrysDiffInPhi(const DetId& crysId,const DetId& orginId);
278 
279  //useful functions for showerRoundnessBarrel function
280  static int deltaIEta(int seed_ieta, int rh_ieta);
281  static int deltaIPhi(int seed_iphi, int rh_iphi);
282  static std::vector<int> getSeedPosition(const std::vector<std::pair<const EcalRecHit*,float> >&RH_ptrs);
283  static float getSumEnergy(const std::vector<std::pair<const EcalRecHit*,float> >&RH_ptrs_fracs);
284  static float computeWeight(float eRH, float energyTotal, int weightedPositionMethod);
285 
286 };
287 
288 // implementation
289 template<bool noZS>
290 float EcalClusterToolsT<noZS>::getFraction( const std::vector< std::pair<DetId, float> > &v_id, DetId id
291  ){
292  if(noZS) return 1.0;
293  float frac = 0.0;
294  for ( size_t i = 0; i < v_id.size(); ++i ) {
295  if(v_id[i].first.rawId()==id.rawId()){
296  frac= v_id[i].second;
297  break;
298  }
299  }
300  return frac;
301 }
302 
303 template<bool noZS>
304 std::pair<DetId, float> EcalClusterToolsT<noZS>::getMaximum( const std::vector< std::pair<DetId, float> > &v_id, const EcalRecHitCollection *recHits)
305 {
306  float max = 0;
307  DetId id(0);
308  for ( size_t i = 0; i < v_id.size(); ++i ) {
309  float energy = recHitEnergy( v_id[i].first, recHits ) * (noZS ? 1.0 : v_id[i].second);
310  if ( energy > max ) {
311  max = energy;
312  id = v_id[i].first;
313  }
314  }
315  return std::pair<DetId, float>(id, max);
316 }
317 
318 template<bool noZS>
319 std::pair<DetId, float> EcalClusterToolsT<noZS>::getMaximum( const std::vector< std::pair<DetId, float> > &v_id, const EcalRecHitCollection *recHits,const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv)
320 {
321  float max = 0;
322  DetId id(0);
323  for ( size_t i = 0; i < v_id.size(); ++i ) {
324  float energy = recHitEnergy( v_id[i].first, recHits,flagsexcl, severitiesexcl, sevLv ) * (noZS ? 1.0 : v_id[i].second);
325  if ( energy > max ) {
326  max = energy;
327  id = v_id[i].first;
328  }
329  }
330  return std::pair<DetId, float>(id, max);
331 }
332 
333 template<bool noZS>
334 std::pair<DetId, float> EcalClusterToolsT<noZS>::getMaximum( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits)
335 {
336  return getMaximum( cluster.hitsAndFractions(), recHits );
337 }
338 
339 template<bool noZS>
340 std::pair<DetId, float> EcalClusterToolsT<noZS>::getMaximum( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits,const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
341 {
342  return getMaximum( cluster.hitsAndFractions(), recHits, flagsexcl, severitiesexcl, sevLv );
343 }
344 
345 
346 template<bool noZS>
347 float EcalClusterToolsT<noZS>::recHitEnergy(DetId id, const EcalRecHitCollection *recHits)
348 {
349  if ( id == DetId(0) ) {
350  return 0;
351  } else {
352  EcalRecHitCollection::const_iterator it = recHits->find( id );
353  if ( it != recHits->end() ) {
354  if( noZS && ( it->checkFlag(EcalRecHit::kTowerRecovered) ||
355  it->checkFlag(EcalRecHit::kWeird) ||
356  (it->detid().subdetId() == EcalBarrel &&
357  it->checkFlag(EcalRecHit::kDiWeird) )
358  )
359  ) {
360  return 0.0;
361  } else {
362  return (*it).energy();
363  }
364  } else {
365  //throw cms::Exception("EcalRecHitNotFound") << "The recHit corresponding to the DetId" << id.rawId() << " not found in the EcalRecHitCollection";
366  // the recHit is not in the collection (hopefully zero suppressed)
367  return 0;
368  }
369  }
370  return 0;
371 }
372 
373 template<bool noZS>
374 float EcalClusterToolsT<noZS>::recHitEnergy(DetId id, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv)
375 {
376  if ( id == DetId(0) ) {
377  return 0;
378  } else {
379  EcalRecHitCollection::const_iterator it = recHits->find( id );
380  if ( it != recHits->end() ) {
381  // avoid anomalous channels (recoFlag based)
382  uint32_t rhFlag = (*it).recoFlag();
383  std::vector<int>::const_iterator vit = std::find( flagsexcl.begin(), flagsexcl.end(), rhFlag );
384  //if your flag was found to be one which is excluded, zero out
385  //this energy.
386  if ( vit != flagsexcl.end() ) return 0;
387 
388  int severityFlag = sevLv->severityLevel( it->id(), *recHits);
389  std::vector<int>::const_iterator sit = std::find(severitiesexcl.begin(), severitiesexcl.end(), severityFlag);
390  //if you were flagged by some condition (kWeird etc.)
391  //zero out this energy.
392  if (sit!= severitiesexcl.end())
393  return 0;
394  //If we make it here, you're a found, clean hit.
395  return (*it).energy();
396  } else {
397  //throw cms::Exception("EcalRecHitNotFound") << "The recHit corresponding to the DetId" << id.rawId() << " not found in the EcalRecHitCollection";
398  // the recHit is not in the collection (hopefully zero suppressed)
399  return 0;
400  }
401  }
402  return 0;
403 }
404 
405 
406 // Returns the energy in a rectangle of crystals
407 // specified in eta by ixMin and ixMax
408 // and in phi by iyMin and iyMax
409 //
410 // Reference picture (X=seed crystal)
411 // iy ___________
412 // 2 |_|_|_|_|_|
413 // 1 |_|_|_|_|_|
414 // 0 |_|_|X|_|_|
415 // -1 |_|_|_|_|_|
416 // -2 |_|_|_|_|_|
417 // -2 -1 0 1 2 ix
418 template<bool noZS>
419 float EcalClusterToolsT<noZS>::matrixEnergy( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax )
420 {
421  //take into account fractions
422  // fast version
423  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( id, topology->getSubdetectorTopology( id ) );
424  float energy = 0;
425  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
426  for ( int i = ixMin; i <= ixMax; ++i ) {
427  for ( int j = iyMin; j <= iyMax; ++j ) {
428  cursor.home();
429  cursor.offsetBy( i, j );
430  float frac=getFraction(v_id,*cursor);
431  energy += recHitEnergy( *cursor, recHits )*frac;
432  }
433  }
434  // slow elegant version
435  //float energy = 0;
436  //std::vector<DetId> v_id = matrixDetId( topology, id, ixMin, ixMax, iyMin, iyMax );
437  //for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
438  // energy += recHitEnergy( *it, recHits );
439  //}
440  return energy;
441 }
442 
443 template<bool noZS>
444 float EcalClusterToolsT<noZS>::matrixEnergy( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
445 {
446  // fast version
447  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( id, topology->getSubdetectorTopology( id ) );
448  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
449  float energy = 0;
450  for ( int i = ixMin; i <= ixMax; ++i ) {
451  for ( int j = iyMin; j <= iyMax; ++j ) {
452  cursor.home();
453  cursor.offsetBy( i, j );
454  float frac=getFraction(v_id,*cursor);
455  energy += recHitEnergy( *cursor, recHits, flagsexcl, severitiesexcl, sevLv )*frac;
456  }
457  }
458  return energy;
459 }
460 
461 
462 template<bool noZS>
463 std::vector<DetId> EcalClusterToolsT<noZS>::matrixDetId( const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax )
464 {
465  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( id, topology->getSubdetectorTopology( id ) );
466  std::vector<DetId> v;
467  for ( int i = ixMin; i <= ixMax; ++i ) {
468  for ( int j = iyMin; j <= iyMax; ++j ) {
469  cursor.home();
470  cursor.offsetBy( i, j );
471  if ( *cursor != DetId(0) ) v.push_back( *cursor );
472  }
473  }
474  return v;
475 }
476 
477 
478 template<bool noZS>
479 float EcalClusterToolsT<noZS>::e2x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
480 {
481  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
482  std::list<float> energies;
483  float max_E = matrixEnergy( cluster, recHits, topology, id, -1, 0, -1, 0 );
484  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -1, 0, 0, 1 ) );
485  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, 0, 1, 0, 1 ) );
486  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, 0, 1, -1, 0 ) );
487  return max_E;
488 }
489 
490 template<bool noZS>
491 float EcalClusterToolsT<noZS>::e2x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
492 {
493  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
494  float max_E = matrixEnergy( cluster, recHits, topology, id, -1, 0, -1, 0,flagsexcl, severitiesexcl, sevLv );
495  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -1, 0, 0, 1,flagsexcl, severitiesexcl, sevLv ) );
496  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, 0, 1, 0, 1,flagsexcl, severitiesexcl, sevLv ) );
497  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, 0, 1, -1, 0,flagsexcl, severitiesexcl, sevLv ) );
498  return max_E;
499 }
500 
501 template<bool noZS>
502 float EcalClusterToolsT<noZS>::e3x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
503 {
504  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
505  float max_E = matrixEnergy( cluster, recHits, topology, id, -1, 1, -1, 0 );
506  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, 0, 1, -1, 1 ) );
507  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -1, 1, 0, 1 ) );
508  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -1, 0, -1, 1 ) );
509  return max_E;
510 }
511 
512 template<bool noZS>
513 float EcalClusterToolsT<noZS>::e3x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
514 {
515  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
516  std::list<float> energies;
517  float max_E = matrixEnergy( cluster, recHits, topology, id, -1, 1, -1, 0,flagsexcl, severitiesexcl, sevLv );
518  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, 0, 1, -1, 1,flagsexcl, severitiesexcl, sevLv ) );
519  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -1, 1, 0, 1,flagsexcl, severitiesexcl, sevLv ) );
520  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -1, 0, -1, 1,flagsexcl, severitiesexcl, sevLv ) );
521  return max_E;
522 }
523 
524 template<bool noZS>
525 float EcalClusterToolsT<noZS>::e3x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
526 {
527  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
528  return matrixEnergy( cluster, recHits, topology, id, -1, 1, -1, 1 );
529 }
530 
531 template<bool noZS>
532 float EcalClusterToolsT<noZS>::e3x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
533 {
534  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
535  return matrixEnergy( cluster, recHits, topology, id, -1, 1, -1, 1,flagsexcl, severitiesexcl, sevLv );
536 }
537 
538 template<bool noZS>
539 float EcalClusterToolsT<noZS>::e4x4( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
540 {
541  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
542  float max_E = matrixEnergy( cluster, recHits, topology, id, -1, 2, -2, 1 );
543  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -2, 1, -2, 1 ) );
544  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -2, 1, -1, 2 ) );
545  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -1, 2, -1, 2 ) );
546  return max_E;
547 }
548 
549 template<bool noZS>
550 float EcalClusterToolsT<noZS>::e4x4( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
551 {
552  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
553  std::list<float> energies;
554  float max_E = matrixEnergy( cluster, recHits, topology, id, -1, 2, -2, 1,flagsexcl, severitiesexcl, sevLv );
555  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -2, 1, -2, 1,flagsexcl, severitiesexcl, sevLv ) );
556  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -2, 1, -1, 2,flagsexcl, severitiesexcl, sevLv ) );
557  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -1, 2, -1, 2,flagsexcl, severitiesexcl, sevLv ) );
558  return max_E;
559 }
560 
561 
562 template<bool noZS>
563 float EcalClusterToolsT<noZS>::e5x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
564 {
565  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
566  return matrixEnergy( cluster, recHits, topology, id, -2, 2, -2, 2 );
567 }
568 
569 template<bool noZS>
570 float EcalClusterToolsT<noZS>::e5x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
571 {
572  DetId id = getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).first;
573  return matrixEnergy( cluster, recHits, topology, id, -2, 2, -2, 2,flagsexcl, severitiesexcl, sevLv );
574 }
575 
576 template<bool noZS>
577 float EcalClusterToolsT<noZS>::eMax( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits )
578 {
579  return getMaximum( cluster.hitsAndFractions(), recHits ).second;
580 }
581 
582 template<bool noZS>
583 float EcalClusterToolsT<noZS>::eMax( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
584 {
585  return getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).second;
586 }
587 
588 template<bool noZS>
589 float EcalClusterToolsT<noZS>::e2nd( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits )
590 {
591  std::vector<float> energies;
592  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
593  energies.reserve( v_id.size() );
594  if ( v_id.size() < 2 ) return 0;
595  for ( size_t i = 0; i < v_id.size(); ++i ) {
596  energies.push_back( recHitEnergy( v_id[i].first, recHits ) * (noZS ? 1.0 : v_id[i].second) );
597  }
598  std::partial_sort( energies.begin(), energies.begin()+2, energies.end(), std::greater<float>() );
599  return energies[1];
600 
601 
602 }
603 
604 template<bool noZS>
605 float EcalClusterToolsT<noZS>::e2nd( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
606 {
607  std::vector<float> energies;
608  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
609  energies.reserve( v_id.size() );
610  if ( v_id.size() < 2 ) return 0;
611  for ( size_t i = 0; i < v_id.size(); ++i ) {
612  energies.push_back( recHitEnergy( v_id[i].first, recHits,flagsexcl, severitiesexcl, sevLv ) * (noZS ? 1.0 : v_id[i].second) );
613  }
614  std::partial_sort( energies.begin(), energies.begin()+2, energies.end(), std::greater<float>() );
615  return energies[1];
616 
617 
618 }
619 
620 template<bool noZS>
621 float EcalClusterToolsT<noZS>::e2x5Right( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
622 {
623  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
624  return matrixEnergy( cluster, recHits, topology, id, 1, 2, -2, 2 );
625 }
626 
627 template<bool noZS>
628 float EcalClusterToolsT<noZS>::e2x5Right( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
629 {
630  DetId id = getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).first;
631  return matrixEnergy( cluster, recHits, topology, id, 1, 2, -2, 2,flagsexcl, severitiesexcl, sevLv );
632 }
633 
634 template<bool noZS>
635 float EcalClusterToolsT<noZS>::e2x5Left( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
636 {
637  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
638  return matrixEnergy( cluster, recHits, topology, id, -2, -1, -2, 2 );
639 }
640 
641 template<bool noZS>
642 float EcalClusterToolsT<noZS>::e2x5Left( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
643 {
644  DetId id = getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).first;
645  return matrixEnergy( cluster, recHits, topology, id, -2, -1, -2, 2,flagsexcl, severitiesexcl, sevLv );
646 }
647 
648 
649 template<bool noZS>
650 float EcalClusterToolsT<noZS>::e2x5Top( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
651 {
652  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
653  return matrixEnergy( cluster, recHits, topology, id, -2, 2, 1, 2 );
654 }
655 
656 template<bool noZS>
657 float EcalClusterToolsT<noZS>::e2x5Top( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
658 {
659  DetId id = getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).first;
660  return matrixEnergy( cluster, recHits, topology, id, -2, 2, 1, 2,flagsexcl, severitiesexcl, sevLv );
661 }
662 
663 template<bool noZS>
664 float EcalClusterToolsT<noZS>::e2x5Bottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
665 {
666  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
667  return matrixEnergy( cluster, recHits, topology, id, -2, 2, -2, -1 );
668 }
669 
670 template<bool noZS>
671 float EcalClusterToolsT<noZS>::e2x5Bottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
672 {
673  DetId id = getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).first;
674  return matrixEnergy( cluster, recHits, topology, id, -2, 2, -2, -1,flagsexcl, severitiesexcl, sevLv );
675 }
676 
677 // Energy in 2x5 strip containing the max crystal.
678 // Adapted from code by Sam Harper
679 template<bool noZS>
680 float EcalClusterToolsT<noZS>::e2x5Max( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
681 {
682  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
683 
684  // 1x5 strip left of seed
685  float left = matrixEnergy( cluster, recHits, topology, id, -1, -1, -2, 2 );
686  // 1x5 strip right of seed
687  float right = matrixEnergy( cluster, recHits, topology, id, 1, 1, -2, 2 );
688  // 1x5 strip containing seed
689  float centre = matrixEnergy( cluster, recHits, topology, id, 0, 0, -2, 2 );
690 
691  // Return the maximum of (left+center) or (right+center) strip
692  return left > right ? left+centre : right+centre;
693 }
694 
695 template<bool noZS>
696 float EcalClusterToolsT<noZS>::e2x5Max( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
697 {
698  DetId id = getMaximum( cluster.hitsAndFractions(), recHits,flagsexcl, severitiesexcl, sevLv ).first;
699 
700  // 1x5 strip left of seed
701  float left = matrixEnergy( cluster, recHits, topology, id, -1, -1, -2, 2,flagsexcl, severitiesexcl, sevLv );
702  // 1x5 strip right of seed
703  float right = matrixEnergy( cluster, recHits, topology, id, 1, 1, -2, 2,flagsexcl, severitiesexcl, sevLv );
704  // 1x5 strip containing seed
705  float centre = matrixEnergy( cluster, recHits, topology, id, 0, 0, -2, 2,flagsexcl, severitiesexcl, sevLv );
706 
707  // Return the maximum of (left+center) or (right+center) strip
708  return left > right ? left+centre : right+centre;
709 }
710 
711 template<bool noZS>
712 float EcalClusterToolsT<noZS>::e1x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
713 {
714  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
715  return matrixEnergy( cluster, recHits, topology, id, 0, 0, -2, 2 );
716 }
717 
718 template<bool noZS>
719 float EcalClusterToolsT<noZS>::e1x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv)
720 {
721  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
722  return matrixEnergy( cluster, recHits, topology, id, 0, 0, -2, 2 ,
723  flagsexcl, severitiesexcl, sevLv);
724 }
725 
726 template<bool noZS>
727 float EcalClusterToolsT<noZS>::e5x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
728 {
729  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
730  return matrixEnergy( cluster, recHits, topology, id, -2, 2, 0, 0 );
731 }
732 
733 template<bool noZS>
734 float EcalClusterToolsT<noZS>::e5x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
735 {
736  DetId id = getMaximum( cluster.hitsAndFractions(), recHits).first;
737  return matrixEnergy( cluster, recHits, topology, id, -2, 2, 0, 0,flagsexcl, severitiesexcl, sevLv );
738 }
739 
740 template<bool noZS>
741 float EcalClusterToolsT<noZS>::e1x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
742 {
743  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
744  return matrixEnergy( cluster, recHits, topology, id, 0, 0, -1, 1 );
745 }
746 
747 template<bool noZS>
748 float EcalClusterToolsT<noZS>::e1x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
749 {
750  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
751  return matrixEnergy( cluster, recHits, topology, id, 0, 0, -1, 1,flagsexcl, severitiesexcl, sevLv );
752 }
753 
754 template<bool noZS>
755 float EcalClusterToolsT<noZS>::e3x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
756 {
757  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
758  return matrixEnergy( cluster, recHits, topology, id, -1, 1, 0, 0 );
759 }
760 
761 template<bool noZS>
762 float EcalClusterToolsT<noZS>::e3x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
763 {
764  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
765  return matrixEnergy( cluster, recHits, topology, id, -1, 1, 0, 0,flagsexcl, severitiesexcl, sevLv );
766 }
767 
768 template<bool noZS>
769 float EcalClusterToolsT<noZS>::eLeft( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
770 {
771  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
772  return matrixEnergy( cluster, recHits, topology, id, -1, -1, 0, 0 );
773 }
774 
775 template<bool noZS>
776 float EcalClusterToolsT<noZS>::eLeft( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
777 {
778  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
779  return matrixEnergy( cluster, recHits, topology, id, -1, -1, 0, 0,flagsexcl, severitiesexcl, sevLv );
780 }
781 
782 template<bool noZS>
783 float EcalClusterToolsT<noZS>::eRight( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
784 {
785  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
786  return matrixEnergy( cluster, recHits, topology, id, 1, 1, 0, 0 );
787 }
788 
789 template<bool noZS>
790 float EcalClusterToolsT<noZS>::eRight( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
791 {
792  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
793  return matrixEnergy( cluster, recHits, topology, id, 1, 1, 0, 0,flagsexcl, severitiesexcl, sevLv );
794 }
795 
796 template<bool noZS>
797 float EcalClusterToolsT<noZS>::eTop( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
798 {
799  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
800  return matrixEnergy( cluster, recHits, topology, id, 0, 0, 1, 1 );
801 }
802 
803 template<bool noZS>
804 float EcalClusterToolsT<noZS>::eTop( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
805 {
806  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
807  return matrixEnergy( cluster, recHits, topology, id, 0, 0, 1, 1,flagsexcl, severitiesexcl, sevLv );
808 }
809 
810 template<bool noZS>
811 float EcalClusterToolsT<noZS>::eBottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
812 {
813  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
814  return matrixEnergy( cluster, recHits, topology, id, 0, 0, -1, -1 );
815 }
816 
817 template<bool noZS>
818 float EcalClusterToolsT<noZS>::eBottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
819 {
820  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
821  return matrixEnergy( cluster, recHits, topology, id, 0, 0, -1, -1,flagsexcl, severitiesexcl, sevLv );
822 }
823 
824 template<bool noZS>
825 std::vector<float> EcalClusterToolsT<noZS>::energyBasketFractionEta( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits )
826 {
827  std::vector<float> basketFraction( 2 * EBDetId::kModulesPerSM );
828  float clusterEnergy = cluster.energy();
829  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
830  if ( v_id[0].first.subdetId() != EcalBarrel ) {
831  edm::LogWarning("EcalClusterToolsT<noZS>::energyBasketFractionEta") << "Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel basic-clusters. Returning empty vector.";
832  return basketFraction;
833  }
834  for ( size_t i = 0; i < v_id.size(); ++i ) {
835  basketFraction[ EBDetId(v_id[i].first).im()-1 + EBDetId(v_id[i].first).positiveZ()*EBDetId::kModulesPerSM ] += recHitEnergy( v_id[i].first, recHits ) * v_id[i].second / clusterEnergy;
836  }
837  std::sort( basketFraction.rbegin(), basketFraction.rend() );
838  return basketFraction;
839 }
840 
841 template<bool noZS>
842 std::vector<float> EcalClusterToolsT<noZS>::energyBasketFractionEta( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
843 {
844  std::vector<float> basketFraction( 2 * EBDetId::kModulesPerSM );
845  float clusterEnergy = cluster.energy();
846  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
847  if ( v_id[0].first.subdetId() != EcalBarrel ) {
848  edm::LogWarning("EcalClusterToolsT<noZS>::energyBasketFractionEta") << "Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel basic-clusters. Returning empty vector.";
849  return basketFraction;
850  }
851  for ( size_t i = 0; i < v_id.size(); ++i ) {
852  basketFraction[ EBDetId(v_id[i].first).im()-1 + EBDetId(v_id[i].first).positiveZ()*EBDetId::kModulesPerSM ] += recHitEnergy( v_id[i].first, recHits,flagsexcl, severitiesexcl, sevLv ) * (noZS ? 1.0 : v_id[i].second) / clusterEnergy;
853  }
854  std::sort( basketFraction.rbegin(), basketFraction.rend() );
855  return basketFraction;
856 }
857 
858 template<bool noZS>
859 std::vector<float> EcalClusterToolsT<noZS>::energyBasketFractionPhi( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits )
860 {
861  std::vector<float> basketFraction( 2 * (EBDetId::MAX_IPHI / EBDetId::kCrystalsInPhi) );
862  float clusterEnergy = cluster.energy();
863  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
864  if ( v_id[0].first.subdetId() != EcalBarrel ) {
865  edm::LogWarning("EcalClusterToolsT<noZS>::energyBasketFractionPhi") << "Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel basic-clusters. Returning empty vector.";
866  return basketFraction;
867  }
868  for ( size_t i = 0; i < v_id.size(); ++i ) {
869  basketFraction[ (EBDetId(v_id[i].first).iphi()-1)/EBDetId::kCrystalsInPhi + EBDetId(v_id[i].first).positiveZ()*EBDetId::kTowersInPhi] += recHitEnergy( v_id[i].first, recHits ) * (noZS ? 1.0 : v_id[i].second) / clusterEnergy;
870  }
871  std::sort( basketFraction.rbegin(), basketFraction.rend() );
872  return basketFraction;
873 }
874 
875 template<bool noZS>
876 std::vector<float> EcalClusterToolsT<noZS>::energyBasketFractionPhi( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
877 {
878  std::vector<float> basketFraction( 2 * (EBDetId::MAX_IPHI / EBDetId::kCrystalsInPhi) );
879  float clusterEnergy = cluster.energy();
880  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
881  if ( v_id[0].first.subdetId() != EcalBarrel ) {
882  edm::LogWarning("EcalClusterToolsT<noZS>::energyBasketFractionPhi") << "Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel basic-clusters. Returning empty vector.";
883  return basketFraction;
884  }
885  for ( size_t i = 0; i < v_id.size(); ++i ) {
886  basketFraction[ (EBDetId(v_id[i].first).iphi()-1)/EBDetId::kCrystalsInPhi + EBDetId(v_id[i].first).positiveZ()*EBDetId::kTowersInPhi] += recHitEnergy( v_id[i].first, recHits,flagsexcl, severitiesexcl, sevLv ) * (noZS ? 1.0 : v_id[i].second) / clusterEnergy;
887  }
888  std::sort( basketFraction.rbegin(), basketFraction.rend() );
889  return basketFraction;
890 }
891 
892 template<bool noZS>
893 std::vector<typename EcalClusterToolsT<noZS>::EcalClusterEnergyDeposition> EcalClusterToolsT<noZS>::getEnergyDepTopology( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW, float w0 )
894 {
895  std::vector<typename EcalClusterToolsT<noZS>::EcalClusterEnergyDeposition> energyDistribution;
896  // init a map of the energy deposition centered on the
897  // cluster centroid. This is for momenta calculation only.
898  CLHEP::Hep3Vector clVect(cluster.position().x(), cluster.position().y(), cluster.position().z());
899  CLHEP::Hep3Vector clDir(clVect);
900  clDir*=1.0/clDir.mag();
901  // in the transverse plane, axis perpendicular to clusterDir
902  CLHEP::Hep3Vector theta_axis(clDir.y(),-clDir.x(),0.0);
903  theta_axis *= 1.0/theta_axis.mag();
904  CLHEP::Hep3Vector phi_axis = theta_axis.cross(clDir);
905 
906  const std::vector< std::pair<DetId, float> >& clusterDetIds = cluster.hitsAndFractions();
907 
908  EcalClusterEnergyDeposition clEdep;
909  EcalRecHit testEcalRecHit;
910  std::vector< std::pair<DetId, float> >::const_iterator posCurrent;
911  // loop over crystals
912  for(posCurrent=clusterDetIds.begin(); posCurrent!=clusterDetIds.end(); ++posCurrent) {
913  EcalRecHitCollection::const_iterator itt = recHits->find( (*posCurrent).first );
914  testEcalRecHit=*itt;
915 
916  if(( (*posCurrent).first != DetId(0)) && (recHits->find( (*posCurrent).first ) != recHits->end())) {
917  clEdep.deposited_energy = testEcalRecHit.energy() * (noZS ? 1.0 : (*posCurrent).second);
918  // if logarithmic weight is requested, apply cut on minimum energy of the recHit
919  if(logW) {
920  //double w0 = parameterMap_.find("W0")->second;
921 
922  double weight = std::max(0.0, w0 + log(std::abs(clEdep.deposited_energy)/cluster.energy()) );
923  if(weight==0) {
924  LogDebug("ClusterShapeAlgo") << "Crystal has insufficient energy: E = "
925  << clEdep.deposited_energy << " GeV; skipping... ";
926  continue;
927  }
928  else LogDebug("ClusterShapeAlgo") << "===> got crystal. Energy = " << clEdep.deposited_energy << " GeV. ";
929  }
930  DetId id_ = (*posCurrent).first;
931  const CaloCellGeometry *this_cell = geometry->getSubdetectorGeometry(id_)->getGeometry(id_);
932  GlobalPoint cellPos = this_cell->getPosition();
933  CLHEP::Hep3Vector gblPos (cellPos.x(),cellPos.y(),cellPos.z()); //surface position?
934  // Evaluate the distance from the cluster centroid
935  CLHEP::Hep3Vector diff = gblPos - clVect;
936  // Important: for the moment calculation, only the "lateral distance" is important
937  // "lateral distance" r_i = distance of the digi position from the axis Origin-Cluster Center
938  // ---> subtract the projection on clDir
939  CLHEP::Hep3Vector DigiVect = diff - diff.dot(clDir)*clDir;
940  clEdep.r = DigiVect.mag();
941  LogDebug("ClusterShapeAlgo") << "E = " << clEdep.deposited_energy
942  << "\tdiff = " << diff.mag()
943  << "\tr = " << clEdep.r;
944  clEdep.phi = DigiVect.angle(theta_axis);
945  if(DigiVect.dot(phi_axis)<0) clEdep.phi = 2 * M_PI - clEdep.phi;
946  energyDistribution.push_back(clEdep);
947  }
948  }
949  return energyDistribution;
950 }
951 
952 template<bool noZS>
953 std::vector<float> EcalClusterToolsT<noZS>::lat( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW, float w0 )
954 {
955  std::vector<EcalClusterToolsT::EcalClusterEnergyDeposition> energyDistribution = getEnergyDepTopology( cluster, recHits, geometry, logW, w0 );
956 
957  std::vector<float> lat;
958  double r, redmoment=0;
959  double phiRedmoment = 0 ;
960  double etaRedmoment = 0 ;
961  int n,n1,n2,tmp;
962  int clusterSize=energyDistribution.size();
963  float etaLat_, phiLat_, lat_;
964  if (clusterSize<3) {
965  etaLat_ = 0.0 ;
966  lat_ = 0.0;
967  lat.push_back(0.);
968  lat.push_back(0.);
969  lat.push_back(0.);
970  return lat;
971  }
972 
973  n1=0; n2=1;
974  if (energyDistribution[1].deposited_energy >
975  energyDistribution[0].deposited_energy)
976  {
977  tmp=n2; n2=n1; n1=tmp;
978  }
979  for (int i=2; i<clusterSize; i++) {
980  n=i;
981  if (energyDistribution[i].deposited_energy >
982  energyDistribution[n1].deposited_energy)
983  {
984  tmp = n2;
985  n2 = n1; n1 = i; n=tmp;
986  } else {
987  if (energyDistribution[i].deposited_energy >
988  energyDistribution[n2].deposited_energy)
989  {
990  tmp=n2; n2=i; n=tmp;
991  }
992  }
993 
994  r = energyDistribution[n].r;
995  redmoment += r*r* energyDistribution[n].deposited_energy;
996  double rphi = r * cos (energyDistribution[n].phi) ;
997  phiRedmoment += rphi * rphi * energyDistribution[n].deposited_energy;
998  double reta = r * sin (energyDistribution[n].phi) ;
999  etaRedmoment += reta * reta * energyDistribution[n].deposited_energy;
1000  }
1001  double e1 = energyDistribution[n1].deposited_energy;
1002  double e2 = energyDistribution[n2].deposited_energy;
1003 
1004  lat_ = redmoment/(redmoment+2.19*2.19*(e1+e2));
1005  phiLat_ = phiRedmoment/(phiRedmoment+2.19*2.19*(e1+e2));
1006  etaLat_ = etaRedmoment/(etaRedmoment+2.19*2.19*(e1+e2));
1007 
1008  lat.push_back(etaLat_);
1009  lat.push_back(phiLat_);
1010  lat.push_back(lat_);
1011  return lat;
1012 }
1013 
1014 template<bool noZS>
1015 math::XYZVector EcalClusterToolsT<noZS>::meanClusterPosition( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, const CaloGeometry *geometry )
1016 {
1017  // find mean energy position of a 5x5 cluster around the maximum
1018  math::XYZVector meanPosition(0.0, 0.0, 0.0);
1019  const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
1020  std::vector<DetId> v_id = matrixDetId( topology, getMaximum( cluster, recHits ).first, -2, 2, -2, 2 );
1021  for( const std::pair<DetId,float>& hitAndFrac : hsAndFs ) {
1022  for( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
1023  if( hitAndFrac.first != *it && !noZS) continue;
1024  GlobalPoint positionGP = geometry->getSubdetectorGeometry( *it )->getGeometry( *it )->getPosition();
1025  math::XYZVector position(positionGP.x(),positionGP.y(),positionGP.z());
1026  meanPosition = meanPosition + recHitEnergy( *it, recHits ) * position * hitAndFrac.second;
1027  }
1028  if(noZS) break;
1029  }
1030  return meanPosition / e5x5( cluster, recHits, topology );
1031 }
1032 
1033 //================================================= meanClusterPosition===================================================================================
1034 template<bool noZS>
1035 math::XYZVector EcalClusterToolsT<noZS>::meanClusterPosition( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, const CaloGeometry *geometry, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv )
1036 {
1037  // find mean energy position of a 5x5 cluster around the maximum
1038  math::XYZVector meanPosition(0.0, 0.0, 0.0);
1039  const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
1040  std::vector<DetId> v_id = matrixDetId( topology, getMaximum( cluster, recHits ).first, -2, 2, -2, 2 );
1041  for( const std::pair<DetId,float>& hitAndFrac : hsAndFs ) {
1042  for( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
1043  if( hitAndFrac.first != *it && !noZS ) continue;
1044  GlobalPoint positionGP = geometry->getSubdetectorGeometry( *it )->getGeometry( *it )->getPosition();
1045  math::XYZVector position(positionGP.x(),positionGP.y(),positionGP.z());
1046  meanPosition = meanPosition + recHitEnergy( *it, recHits,flagsexcl, severitiesexcl, sevLv ) * position * hitAndFrac.second;
1047  }
1048  if(noZS) break;
1049  }
1050  return meanPosition / e5x5( cluster, recHits, topology,flagsexcl, severitiesexcl, sevLv );
1051 }
1052 
1053 //returns mean energy weighted eta/phi in crystals from the seed
1054 //iPhi is not defined for endcap and is returned as zero
1055 //return <eta,phi>
1056 //we have an issue in working out what to do for negative energies
1057 //I (Sam Harper) think it makes sense to ignore crystals with E<0 in the calculation as they are ignored
1058 //in the sigmaIEtaIEta calculation (well they arent fully ignored, they do still contribute to the e5x5 sum
1059 //in the sigmaIEtaIEta calculation but not here)
1060 template<bool noZS>
1061 std::pair<float,float> EcalClusterToolsT<noZS>::mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology)
1062 {
1063  DetId seedId = getMaximum( cluster, recHits ).first;
1064  float meanDEta=0.;
1065  float meanDPhi=0.;
1066  float energySum=0.;
1067 
1068  const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
1069  std::vector<DetId> v_id = matrixDetId( topology,seedId, -2, 2, -2, 2 );
1070  for( const std::pair<DetId,float>& hAndF : hsAndFs ) {
1071  for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
1072  if( hAndF.first != *it && !noZS ) continue;
1073  float energy = recHitEnergy(*it,recHits) * hAndF.second;
1074  if(energy<0.) continue;//skipping negative energy crystals
1075  meanDEta += energy * getNrCrysDiffInEta(*it,seedId);
1076  meanDPhi += energy * getNrCrysDiffInPhi(*it,seedId);
1077  energySum +=energy;
1078  }
1079  if(noZS) break;
1080  }
1081  meanDEta /=energySum;
1082  meanDPhi /=energySum;
1083  return std::pair<float,float>(meanDEta,meanDPhi);
1084 }
1085 
1086 template<bool noZS>
1087 std::pair<float,float> EcalClusterToolsT<noZS>::mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv)
1088 {
1089  DetId seedId = getMaximum( cluster, recHits ).first;
1090  float meanDEta=0.;
1091  float meanDPhi=0.;
1092  float energySum=0.;
1093 
1094  const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
1095  std::vector<DetId> v_id = matrixDetId( topology,seedId, -2, 2, -2, 2 );
1096  for( const std::pair<DetId,float>& hAndF : hsAndFs ) {
1097  for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
1098  if( hAndF.first != *it && !noZS ) continue;
1099  float energy = recHitEnergy(*it,recHits,flagsexcl, severitiesexcl, sevLv) * hAndF.second;
1100  if(energy<0.) continue;//skipping negative energy crystals
1101  meanDEta += energy * getNrCrysDiffInEta(*it,seedId);
1102  meanDPhi += energy * getNrCrysDiffInPhi(*it,seedId);
1103  energySum +=energy;
1104  }
1105  if(noZS) break;
1106  }
1107  meanDEta /=energySum;
1108  meanDPhi /=energySum;
1109  return std::pair<float,float>(meanDEta,meanDPhi);
1110 }
1111 
1112 //returns mean energy weighted x/y in normalised crystal coordinates
1113 //only valid for endcap, returns 0,0 for barrel
1114 //we have an issue in working out what to do for negative energies
1115 //I (Sam Harper) think it makes sense to ignore crystals with E<0 in the calculation as they are ignored
1116 //in the sigmaIEtaIEta calculation (well they arent fully ignored, they do still contribute to the e5x5 sum
1117 //in the sigmaIEtaIEta calculation but not here)
1118 template<bool noZS>
1119 std::pair<float,float> EcalClusterToolsT<noZS>::mean5x5PositionInXY(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology)
1120 {
1121  DetId seedId = getMaximum( cluster, recHits ).first;
1122 
1123  std::pair<float,float> meanXY(0.,0.);
1124  if(seedId.subdetId()==EcalBarrel) return meanXY;
1125 
1126  float energySum=0.;
1127 
1128  const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
1129  std::vector<DetId> v_id = matrixDetId( topology,seedId, -2, 2, -2, 2 );
1130  for( const std::pair<DetId,float>& hAndF : hsAndFs ) {
1131  for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
1132  if( hAndF.first != *it && !noZS) continue;
1133  float energy = recHitEnergy(*it,recHits) * hAndF.second;
1134  if(energy<0.) continue;//skipping negative energy crystals
1135  meanXY.first += energy * getNormedIX(*it);
1136  meanXY.second += energy * getNormedIY(*it);
1137  energySum +=energy;
1138  }
1139  if(noZS) break;
1140  }
1141  meanXY.first/=energySum;
1142  meanXY.second/=energySum;
1143  return meanXY;
1144 }
1145 
1146 template<bool noZS>
1147 std::pair<float,float> EcalClusterToolsT<noZS>::mean5x5PositionInXY(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv)
1148 {
1149  DetId seedId = getMaximum( cluster, recHits ).first;
1150 
1151  std::pair<float,float> meanXY(0.,0.);
1152  if(seedId.subdetId()==EcalBarrel) return meanXY;
1153 
1154  float energySum=0.;
1155 
1156  const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
1157  std::vector<DetId> v_id = matrixDetId( topology,seedId, -2, 2, -2, 2 );
1158  for( const std::pair<DetId,float>& hAndF : hsAndFs ) {
1159  for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
1160  if( hAndF.first != *it && !noZS) continue;
1161  float energy = recHitEnergy(*it,recHits,flagsexcl, severitiesexcl, sevLv) * (noZS ? 1.0 : hAndF.second);
1162  if(energy<0.) continue;//skipping negative energy crystals
1163  meanXY.first += energy * getNormedIX(*it);
1164  meanXY.second += energy * getNormedIY(*it);
1165  energySum +=energy;
1166  }
1167  if(noZS) break;
1168  }
1169  meanXY.first/=energySum;
1170  meanXY.second/=energySum;
1171  return meanXY;
1172 }
1173 
1174 template<bool noZS>
1175 std::vector<float> EcalClusterToolsT<noZS>::covariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, const CaloGeometry* geometry, float w0)
1176 {
1177  float e_5x5 = e5x5( cluster, recHits, topology );
1178  float covEtaEta, covEtaPhi, covPhiPhi;
1179  if (e_5x5 >= 0.) {
1180  //double w0_ = parameterMap_.find("W0")->second;
1181  const std::vector< std::pair<DetId, float>>& v_id =cluster.hitsAndFractions();
1182  math::XYZVector meanPosition = meanClusterPosition( cluster, recHits, topology, geometry );
1183 
1184  // now we can calculate the covariances
1185  double numeratorEtaEta = 0;
1186  double numeratorEtaPhi = 0;
1187  double numeratorPhiPhi = 0;
1188  double denominator = 0;
1189 
1190  DetId id = getMaximum( v_id, recHits ).first;
1191  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( id, topology->getSubdetectorTopology( id ) );
1192  for ( int i = -2; i <= 2; ++i ) {
1193  for ( int j = -2; j <= 2; ++j ) {
1194  cursor.home();
1195  cursor.offsetBy( i, j );
1196  float frac=getFraction(v_id,*cursor);
1197  float energy = recHitEnergy( *cursor, recHits )*frac;
1198 
1199  if ( energy <= 0 ) continue;
1200 
1201  GlobalPoint position = geometry->getSubdetectorGeometry(*cursor)->getGeometry(*cursor)->getPosition();
1202 
1203  double dPhi = position.phi() - meanPosition.phi();
1204  if (dPhi > + Geom::pi()) { dPhi = Geom::twoPi() - dPhi; }
1205  if (dPhi < - Geom::pi()) { dPhi = Geom::twoPi() + dPhi; }
1206 
1207  double dEta = position.eta() - meanPosition.eta();
1208  double w = 0.;
1209  w = std::max(0.0f, w0 + std::log( energy / e_5x5 ));
1210 
1211  denominator += w;
1212  numeratorEtaEta += w * dEta * dEta;
1213  numeratorEtaPhi += w * dEta * dPhi;
1214  numeratorPhiPhi += w * dPhi * dPhi;
1215  }
1216  }
1217 
1218  if (denominator != 0.0) {
1219  covEtaEta = numeratorEtaEta / denominator;
1220  covEtaPhi = numeratorEtaPhi / denominator;
1221  covPhiPhi = numeratorPhiPhi / denominator;
1222  } else {
1223  covEtaEta = 999.9;
1224  covEtaPhi = 999.9;
1225  covPhiPhi = 999.9;
1226  }
1227 
1228  } else {
1229  // Warn the user if there was no energy in the cells and return zeroes.
1230  // std::cout << "\ClusterShapeAlgo::Calculate_Covariances: no energy in supplied cells.\n";
1231  covEtaEta = 0;
1232  covEtaPhi = 0;
1233  covPhiPhi = 0;
1234  }
1235  std::vector<float> v;
1236  v.push_back( covEtaEta );
1237  v.push_back( covEtaPhi );
1238  v.push_back( covPhiPhi );
1239  return v;
1240 }
1241 
1242 //==================================================== Covariances===========================================================================
1243 template<bool noZS>
1244 std::vector<float> EcalClusterToolsT<noZS>::covariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, const CaloGeometry* geometry,const std::vector<int>& flagsexcl,const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv, float w0)
1245 {
1246  float e_5x5 = e5x5( cluster, recHits, topology,flagsexcl, severitiesexcl, sevLv );
1247  float covEtaEta, covEtaPhi, covPhiPhi;
1248  if (e_5x5 >= 0.) {
1249  //double w0_ = parameterMap_.find("W0")->second;
1250  const std::vector<std::pair<DetId, float>>& v_id= cluster.hitsAndFractions();
1251  math::XYZVector meanPosition = meanClusterPosition( cluster, recHits, topology, geometry,flagsexcl, severitiesexcl, sevLv );
1252 
1253  // now we can calculate the covariances
1254  double numeratorEtaEta = 0;
1255  double numeratorEtaPhi = 0;
1256  double numeratorPhiPhi = 0;
1257  double denominator = 0;
1258 
1259 
1260  DetId id = getMaximum( v_id, recHits ).first;
1261  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( id, topology->getSubdetectorTopology( id ) );
1262  for ( int i = -2; i <= 2; ++i ) {
1263  for ( int j = -2; j <= 2; ++j ) {
1264  cursor.home();
1265  cursor.offsetBy( i, j );
1266  float frac = getFraction(v_id,*cursor);
1267  float energy = recHitEnergy( *cursor, recHits,flagsexcl, severitiesexcl, sevLv )*frac;
1268 
1269  if ( energy <= 0 ) continue;
1270 
1271  GlobalPoint position = geometry->getSubdetectorGeometry(*cursor)->getGeometry(*cursor)->getPosition();
1272 
1273  double dPhi = position.phi() - meanPosition.phi();
1274  if (dPhi > + Geom::pi()) { dPhi = Geom::twoPi() - dPhi; }
1275  if (dPhi < - Geom::pi()) { dPhi = Geom::twoPi() + dPhi; }
1276 
1277  double dEta = position.eta() - meanPosition.eta();
1278  double w = 0.;
1279  w = std::max(0.0f, w0 + std::log( energy / e_5x5 ));
1280 
1281  denominator += w;
1282  numeratorEtaEta += w * dEta * dEta;
1283  numeratorEtaPhi += w * dEta * dPhi;
1284  numeratorPhiPhi += w * dPhi * dPhi;
1285  }
1286  }
1287 
1288  if (denominator != 0.0) {
1289  covEtaEta = numeratorEtaEta / denominator;
1290  covEtaPhi = numeratorEtaPhi / denominator;
1291  covPhiPhi = numeratorPhiPhi / denominator;
1292  } else {
1293  covEtaEta = 999.9;
1294  covEtaPhi = 999.9;
1295  covPhiPhi = 999.9;
1296  }
1297 
1298  } else {
1299  // Warn the user if there was no energy in the cells and return zeroes.
1300  // std::cout << "\ClusterShapeAlgo::Calculate_Covariances: no energy in supplied cells.\n";
1301  covEtaEta = 0;
1302  covEtaPhi = 0;
1303  covPhiPhi = 0;
1304  }
1305  std::vector<float> v;
1306  v.push_back( covEtaEta );
1307  v.push_back( covEtaPhi );
1308  v.push_back( covPhiPhi );
1309  return v;
1310 }
1311 
1312 //for covIEtaIEta,covIEtaIPhi and covIPhiIPhi are defined but only covIEtaIEta has been actively studied
1313 //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
1314 //it also does not require any eta correction function in the endcap
1315 //it is multipled by an approprate crystal size to ensure it gives similar values to covariances(...)
1316 template<bool noZS>
1317 std::vector<float> EcalClusterToolsT<noZS>::localCovariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology,float w0)
1318 {
1319 
1320  float e_5x5 = e5x5( cluster, recHits, topology );
1321  float covEtaEta, covEtaPhi, covPhiPhi;
1322 
1323  if (e_5x5 >= 0.) {
1324  //double w0_ = parameterMap_.find("W0")->second;
1325  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
1326  std::pair<float,float> mean5x5PosInNrCrysFromSeed = mean5x5PositionInLocalCrysCoord( cluster, recHits, topology );
1327  std::pair<float,float> mean5x5XYPos = mean5x5PositionInXY(cluster,recHits,topology);
1328 
1329  // now we can calculate the covariances
1330  double numeratorEtaEta = 0;
1331  double numeratorEtaPhi = 0;
1332  double numeratorPhiPhi = 0;
1333  double denominator = 0;
1334 
1335  //these allow us to scale the localCov by the crystal size
1336  //so that the localCovs have the same average value as the normal covs
1337  const double barrelCrysSize = 0.01745; //approximate size of crystal in eta,phi in barrel
1338  const double endcapCrysSize = 0.0447; //the approximate crystal size sigmaEtaEta was corrected to in the endcap
1339 
1340  DetId seedId = getMaximum( v_id, recHits ).first;
1341 
1342  bool isBarrel=seedId.subdetId()==EcalBarrel;
1343  const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1344 
1345  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( seedId, topology->getSubdetectorTopology( seedId ) );
1346 
1347  for ( int eastNr = -2; eastNr <= 2; ++eastNr ) { //east is eta in barrel
1348  for ( int northNr = -2; northNr <= 2; ++northNr ) { //north is phi in barrel
1349  cursor.home();
1350  cursor.offsetBy( eastNr, northNr);
1351  float frac = getFraction(v_id,*cursor);
1352  float energy = recHitEnergy( *cursor, recHits )*frac;
1353  if ( energy <= 0 ) continue;
1354 
1355  float dEta = getNrCrysDiffInEta(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.first;
1356  float dPhi = 0;
1357 
1358  if(isBarrel) dPhi = getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1359  else dPhi = getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1360 
1361 
1362  double w = std::max(0.0f,w0 + std::log( energy / e_5x5 ));
1363 
1364  denominator += w;
1365  numeratorEtaEta += w * dEta * dEta;
1366  numeratorEtaPhi += w * dEta * dPhi;
1367  numeratorPhiPhi += w * dPhi * dPhi;
1368  } //end east loop
1369  }//end north loop
1370 
1371 
1372  //multiplying by crysSize to make the values compariable to normal covariances
1373  if (denominator != 0.0) {
1374  covEtaEta = crysSize*crysSize* numeratorEtaEta / denominator;
1375  covEtaPhi = crysSize*crysSize* numeratorEtaPhi / denominator;
1376  covPhiPhi = crysSize*crysSize* numeratorPhiPhi / denominator;
1377  } else {
1378  covEtaEta = 999.9;
1379  covEtaPhi = 999.9;
1380  covPhiPhi = 999.9;
1381  }
1382 
1383 
1384  } else {
1385  // Warn the user if there was no energy in the cells and return zeroes.
1386  // std::cout << "\ClusterShapeAlgo::Calculate_Covariances: no energy in supplied cells.\n";
1387  covEtaEta = 0;
1388  covEtaPhi = 0;
1389  covPhiPhi = 0;
1390  }
1391  std::vector<float> v;
1392  v.push_back( covEtaEta );
1393  v.push_back( covEtaPhi );
1394  v.push_back( covPhiPhi );
1395  return v;
1396 }
1397 
1398 //==================================================================localCovariances======================================================================
1399 template<bool noZS>
1400 std::vector<float> EcalClusterToolsT<noZS>::localCovariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology,const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv,float w0)
1401 {
1402 
1403  float e_5x5 = e5x5( cluster, recHits, topology,flagsexcl, severitiesexcl, sevLv );
1404  float covEtaEta, covEtaPhi, covPhiPhi;
1405 
1406  if (e_5x5 >= 0.) {
1407  //double w0_ = parameterMap_.find("W0")->second;
1408  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
1409  std::pair<float,float> mean5x5PosInNrCrysFromSeed = mean5x5PositionInLocalCrysCoord( cluster, recHits, topology,flagsexcl, severitiesexcl, sevLv );
1410  std::pair<float,float> mean5x5XYPos = mean5x5PositionInXY(cluster,recHits,topology,flagsexcl, severitiesexcl, sevLv);
1411 
1412  // now we can calculate the covariances
1413  double numeratorEtaEta = 0;
1414  double numeratorEtaPhi = 0;
1415  double numeratorPhiPhi = 0;
1416  double denominator = 0;
1417 
1418  //these allow us to scale the localCov by the crystal size
1419  //so that the localCovs have the same average value as the normal covs
1420  const double barrelCrysSize = 0.01745; //approximate size of crystal in eta,phi in barrel
1421  const double endcapCrysSize = 0.0447; //the approximate crystal size sigmaEtaEta was corrected to in the endcap
1422 
1423  DetId seedId = getMaximum( v_id, recHits ).first;
1424 
1425  bool isBarrel=seedId.subdetId()==EcalBarrel;
1426  const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1427 
1428  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( seedId, topology->getSubdetectorTopology( seedId ) );
1429 
1430  for ( int eastNr = -2; eastNr <= 2; ++eastNr ) { //east is eta in barrel
1431  for ( int northNr = -2; northNr <= 2; ++northNr ) { //north is phi in barrel
1432  cursor.home();
1433  cursor.offsetBy( eastNr, northNr);
1434  float frac = getFraction(v_id,*cursor);
1435  float energy = recHitEnergy( *cursor, recHits,flagsexcl, severitiesexcl, sevLv)*frac;
1436  if ( energy <= 0 ) continue;
1437 
1438  float dEta = getNrCrysDiffInEta(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.first;
1439  float dPhi = 0;
1440 
1441  if(isBarrel) dPhi = getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1442  else dPhi = getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1443 
1444 
1445  double w = std::max(0.0f, w0 + std::log( energy / e_5x5 ));
1446 
1447  denominator += w;
1448  numeratorEtaEta += w * dEta * dEta;
1449  numeratorEtaPhi += w * dEta * dPhi;
1450  numeratorPhiPhi += w * dPhi * dPhi;
1451  } //end east loop
1452  }//end north loop
1453 
1454 
1455  //multiplying by crysSize to make the values compariable to normal covariances
1456  if (denominator != 0.0) {
1457  covEtaEta = crysSize*crysSize* numeratorEtaEta / denominator;
1458  covEtaPhi = crysSize*crysSize* numeratorEtaPhi / denominator;
1459  covPhiPhi = crysSize*crysSize* numeratorPhiPhi / denominator;
1460  } else {
1461  covEtaEta = 999.9;
1462  covEtaPhi = 999.9;
1463  covPhiPhi = 999.9;
1464  }
1465 
1466 
1467  } else {
1468  // Warn the user if there was no energy in the cells and return zeroes.
1469  // std::cout << "\ClusterShapeAlgo::Calculate_Covariances: no energy in supplied cells.\n";
1470  covEtaEta = 0;
1471  covEtaPhi = 0;
1472  covPhiPhi = 0;
1473  }
1474  std::vector<float> v;
1475  v.push_back( covEtaEta );
1476  v.push_back( covEtaPhi );
1477  v.push_back( covPhiPhi );
1478  return v;
1479 }
1480 
1481 template<bool noZS>
1482 double EcalClusterToolsT<noZS>::zernike20( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, double R0, bool logW, float w0 )
1483 {
1484  return absZernikeMoment( cluster, recHits, geometry, 2, 0, R0, logW, w0 );
1485 }
1486 
1487 template<bool noZS>
1488 double EcalClusterToolsT<noZS>::zernike42( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, double R0, bool logW, float w0 )
1489 {
1490  return absZernikeMoment( cluster, recHits, geometry, 4, 2, R0, logW, w0 );
1491 }
1492 
1493 template<bool noZS>
1494 double EcalClusterToolsT<noZS>::absZernikeMoment( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 )
1495 {
1496  // 1. Check if n,m are correctly
1497  if ((m>n) || ((n-m)%2 != 0) || (n<0) || (m<0)) return -1;
1498 
1499  // 2. Check if n,R0 are within validity Range :
1500  // n>20 or R0<2.19cm just makes no sense !
1501  if ((n>20) || (R0<=2.19)) return -1;
1502  if (n<=5) return fast_AbsZernikeMoment(cluster, recHits, geometry, n, m, R0, logW, w0 );
1503  else return calc_AbsZernikeMoment(cluster, recHits, geometry, n, m, R0, logW, w0 );
1504 }
1505 
1506 template<bool noZS>
1507 double EcalClusterToolsT<noZS>::fast_AbsZernikeMoment(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 )
1508 {
1509  double r,ph,e,Re=0,Im=0;
1510  double TotalEnergy = cluster.energy();
1511  int index = (n/2)*(n/2)+(n/2)+m;
1512  std::vector<EcalClusterEnergyDeposition> energyDistribution = getEnergyDepTopology( cluster, recHits, geometry, logW, w0 );
1513  int clusterSize = energyDistribution.size();
1514  if(clusterSize < 3) return 0.0;
1515 
1516  for (int i=0; i<clusterSize; i++)
1517  {
1518  r = energyDistribution[i].r / R0;
1519  if (r<1) {
1520  std::vector<double> pol;
1521  pol.push_back( f00(r) );
1522  pol.push_back( f11(r) );
1523  pol.push_back( f20(r) );
1524  pol.push_back( f22(r) );
1525  pol.push_back( f31(r) );
1526  pol.push_back( f33(r) );
1527  pol.push_back( f40(r) );
1528  pol.push_back( f42(r) );
1529  pol.push_back( f44(r) );
1530  pol.push_back( f51(r) );
1531  pol.push_back( f53(r) );
1532  pol.push_back( f55(r) );
1533  ph = (energyDistribution[i]).phi;
1534  e = energyDistribution[i].deposited_energy;
1535  Re = Re + e/TotalEnergy * pol[index] * cos( (double) m * ph);
1536  Im = Im - e/TotalEnergy * pol[index] * sin( (double) m * ph);
1537  }
1538  }
1539  return sqrt(Re*Re+Im*Im);
1540 }
1541 
1542 template<bool noZS>
1543 double EcalClusterToolsT<noZS>::calc_AbsZernikeMoment(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 )
1544 {
1545  double r, ph, e, Re=0, Im=0, f_nm;
1546  double TotalEnergy = cluster.energy();
1547  std::vector<EcalClusterEnergyDeposition> energyDistribution = getEnergyDepTopology( cluster, recHits, geometry, logW, w0 );
1548  int clusterSize=energyDistribution.size();
1549  if(clusterSize<3) return 0.0;
1550 
1551  for (int i = 0; i < clusterSize; ++i)
1552  {
1553  r = energyDistribution[i].r / R0;
1554  if (r < 1) {
1555  ph = energyDistribution[i].phi;
1556  e = energyDistribution[i].deposited_energy;
1557  f_nm = 0;
1558  for (int s=0; s<=(n-m)/2; s++) {
1559  if (s%2==0) {
1560  f_nm = f_nm + factorial(n-s)*pow(r,(double) (n-2*s))/(factorial(s)*factorial((n+m)/2-s)*factorial((n-m)/2-s));
1561  } else {
1562  f_nm = f_nm - factorial(n-s)*pow(r,(double) (n-2*s))/(factorial(s)*factorial((n+m)/2-s)*factorial((n-m)/2-s));
1563  }
1564  }
1565  Re = Re + e/TotalEnergy * f_nm * cos( (double) m*ph);
1566  Im = Im - e/TotalEnergy * f_nm * sin( (double) m*ph);
1567  }
1568  }
1569  return sqrt(Re*Re+Im*Im);
1570 }
1571 
1572 //returns the crystal 'eta' from the det id
1573 //it is defined as the number of crystals from the centre in the eta direction
1574 //for the barrel with its eta/phi geometry it is always integer
1575 //for the endcap it is fractional due to the x/y geometry
1576 template<bool noZS>
1577 float EcalClusterToolsT<noZS>::getIEta(const DetId& id)
1578 {
1579  if(id.det()==DetId::Ecal){
1580  if(id.subdetId()==EcalBarrel){
1581  EBDetId ebId(id);
1582  return ebId.ieta();
1583  }else if(id.subdetId()==EcalEndcap){
1584  float iXNorm = getNormedIX(id);
1585  float iYNorm = getNormedIY(id);
1586 
1587  return std::sqrt(iXNorm*iXNorm+iYNorm*iYNorm);
1588  }
1589  }
1590  return 0.;
1591 }
1592 
1593 
1594 //returns the crystal 'phi' from the det id
1595 //it is defined as the number of crystals from the centre in the phi direction
1596 //for the barrel with its eta/phi geometry it is always integer
1597 //for the endcap it is not defined
1598 template<bool noZS>
1599 float EcalClusterToolsT<noZS>::getIPhi(const DetId& id)
1600 {
1601  if(id.det()==DetId::Ecal){
1602  if(id.subdetId()==EcalBarrel){
1603  EBDetId ebId(id);
1604  return ebId.iphi();
1605  }
1606  }
1607  return 0.;
1608 }
1609 
1610 //want to map 1=-50,50=-1,51=1 and 100 to 50 so sub off one if zero or neg
1611 template<bool noZS>
1612 float EcalClusterToolsT<noZS>::getNormedIX(const DetId& id)
1613 {
1614  if(id.det()==DetId::Ecal && id.subdetId()==EcalEndcap){
1615  EEDetId eeId(id);
1616  int iXNorm = eeId.ix()-50;
1617  if(iXNorm<=0) iXNorm--;
1618  return iXNorm;
1619  }
1620  return 0;
1621 }
1622 
1623 //want to map 1=-50,50=-1,51=1 and 100 to 50 so sub off one if zero or neg
1624 template<bool noZS>
1625 float EcalClusterToolsT<noZS>::getNormedIY(const DetId& id)
1626 {
1627  if(id.det()==DetId::Ecal && id.subdetId()==EcalEndcap){
1628  EEDetId eeId(id);
1629  int iYNorm = eeId.iy()-50;
1630  if(iYNorm<=0) iYNorm--;
1631  return iYNorm;
1632  }
1633  return 0;
1634 }
1635 
1636 //nr crystals crysId is away from orgin id in eta
1637 template<bool noZS>
1638 float EcalClusterToolsT<noZS>::getNrCrysDiffInEta(const DetId& crysId,const DetId& orginId)
1639 {
1640  float crysIEta = getIEta(crysId);
1641  float orginIEta = getIEta(orginId);
1642  bool isBarrel = orginId.subdetId()==EcalBarrel;
1643 
1644  float nrCrysDiff = crysIEta-orginIEta;
1645 
1646  //no iEta=0 in barrel, so if go from positive to negative
1647  //need to reduce abs(detEta) by 1
1648  if(isBarrel){
1649  if(crysIEta*orginIEta<0){ // -1 to 1 transition
1650  if(crysIEta>0) nrCrysDiff--;
1651  else nrCrysDiff++;
1652  }
1653  }
1654  return nrCrysDiff;
1655 }
1656 
1657 //nr crystals crysId is away from orgin id in phi
1658 template<bool noZS>
1659 float EcalClusterToolsT<noZS>::getNrCrysDiffInPhi(const DetId& crysId,const DetId& orginId)
1660 {
1661  float crysIPhi = getIPhi(crysId);
1662  float orginIPhi = getIPhi(orginId);
1663  bool isBarrel = orginId.subdetId()==EcalBarrel;
1664 
1665  float nrCrysDiff = crysIPhi-orginIPhi;
1666 
1667  if(isBarrel){ //if barrel, need to map into 0-180
1668  if (nrCrysDiff > + 180) { nrCrysDiff = nrCrysDiff - 360; }
1669  if (nrCrysDiff < - 180) { nrCrysDiff = nrCrysDiff + 360; }
1670  }
1671  return nrCrysDiff;
1672 }
1673 
1674 //nr crystals crysId is away from mean phi in 5x5 in phi
1675 template<bool noZS>
1676 float EcalClusterToolsT<noZS>::getDPhiEndcap(const DetId& crysId,float meanX,float meanY)
1677 {
1678  float iXNorm = getNormedIX(crysId);
1679  float iYNorm = getNormedIY(crysId);
1680 
1681  float hitLocalR2 = (iXNorm-meanX)*(iXNorm-meanX)+(iYNorm-meanY)*(iYNorm-meanY);
1682  float hitR2 = iXNorm*iXNorm+iYNorm*iYNorm;
1683  float meanR2 = meanX*meanX+meanY*meanY;
1684  float hitR = sqrt(hitR2);
1685  float meanR = sqrt(meanR2);
1686 
1687  float tmp = (hitR2+meanR2-hitLocalR2)/(2*hitR*meanR);
1688  if (tmp<-1) tmp =-1;
1689  if (tmp>1) tmp=1;
1690  float phi = acos(tmp);
1691  float dPhi = hitR*phi;
1692 
1693  return dPhi;
1694 }
1695 
1696 template<bool noZS>
1697 std::vector<float> EcalClusterToolsT<noZS>::scLocalCovariances(const reco::SuperCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, float w0)
1698 {
1699  const reco::BasicCluster bcluster = *(cluster.seed());
1700 
1701  float e_5x5 = e5x5(bcluster, recHits, topology);
1702  float covEtaEta, covEtaPhi, covPhiPhi;
1703 
1704  if (e_5x5 >= 0.) {
1705  const std::vector<std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
1706  std::pair<float,float> mean5x5PosInNrCrysFromSeed = mean5x5PositionInLocalCrysCoord(bcluster, recHits, topology);
1707  std::pair<float,float> mean5x5XYPos = mean5x5PositionInXY(cluster,recHits,topology);
1708  // now we can calculate the covariances
1709  double numeratorEtaEta = 0;
1710  double numeratorEtaPhi = 0;
1711  double numeratorPhiPhi = 0;
1712  double denominator = 0;
1713 
1714  const double barrelCrysSize = 0.01745; //approximate size of crystal in eta,phi in barrel
1715  const double endcapCrysSize = 0.0447; //the approximate crystal size sigmaEtaEta was corrected to in the endcap
1716 
1717  DetId seedId = getMaximum(v_id, recHits).first;
1718  bool isBarrel=seedId.subdetId()==EcalBarrel;
1719 
1720  const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1721 
1722  for (size_t i = 0; i < v_id.size(); ++i) {
1723  CaloNavigator<DetId> cursor = CaloNavigator<DetId>(v_id[i].first, topology->getSubdetectorTopology(v_id[i].first));
1724  float frac = getFraction(v_id,*cursor);
1725  float energy = recHitEnergy(*cursor, recHits)*frac;
1726 
1727  if (energy <= 0) continue;
1728 
1729  float dEta = getNrCrysDiffInEta(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.first;
1730  float dPhi = 0;
1731  if(isBarrel) dPhi = getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1732  else dPhi = getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1733 
1734 
1735 
1736  double w = 0.;
1737  w = std::max(0.0f, w0 + std::log( energy / e_5x5 ));
1738 
1739  denominator += w;
1740  numeratorEtaEta += w * dEta * dEta;
1741  numeratorEtaPhi += w * dEta * dPhi;
1742  numeratorPhiPhi += w * dPhi * dPhi;
1743  }
1744 
1745  //multiplying by crysSize to make the values compariable to normal covariances
1746  if (denominator != 0.0) {
1747  covEtaEta = crysSize*crysSize* numeratorEtaEta / denominator;
1748  covEtaPhi = crysSize*crysSize* numeratorEtaPhi / denominator;
1749  covPhiPhi = crysSize*crysSize* numeratorPhiPhi / denominator;
1750  } else {
1751  covEtaEta = 999.9;
1752  covEtaPhi = 999.9;
1753  covPhiPhi = 999.9;
1754  }
1755 
1756  } else {
1757  // Warn the user if there was no energy in the cells and return zeroes.
1758  // std::cout << "\ClusterShapeAlgo::Calculate_Covariances: no energy in supplied cells.\n";
1759  covEtaEta = 0;
1760  covEtaPhi = 0;
1761  covPhiPhi = 0;
1762  }
1763 
1764  std::vector<float> v;
1765  v.push_back( covEtaEta );
1766  v.push_back( covEtaPhi );
1767  v.push_back( covPhiPhi );
1768 
1769  return v;
1770 }
1771 
1772 
1773 //================================================================== scLocalCovariances==============================================================
1774 template<bool noZS>
1775 std::vector<float> EcalClusterToolsT<noZS>::scLocalCovariances(const reco::SuperCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology,const std::vector<int>& flagsexcl, const std::vector<int>& severitiesexcl, const EcalSeverityLevelAlgo *sevLv, float w0)
1776 {
1777  const reco::BasicCluster bcluster = *(cluster.seed());
1778 
1779  float e_5x5 = e5x5(bcluster, recHits, topology);
1780  float covEtaEta, covEtaPhi, covPhiPhi;
1781 
1782  if (e_5x5 >= 0.) {
1783  const std::vector<std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
1784  std::pair<float,float> mean5x5PosInNrCrysFromSeed = mean5x5PositionInLocalCrysCoord(bcluster, recHits, topology,flagsexcl, severitiesexcl, sevLv);
1785  std::pair<float,float> mean5x5XYPos = mean5x5PositionInXY(cluster,recHits,topology,flagsexcl, severitiesexcl, sevLv);
1786  // now we can calculate the covariances
1787  double numeratorEtaEta = 0;
1788  double numeratorEtaPhi = 0;
1789  double numeratorPhiPhi = 0;
1790  double denominator = 0;
1791 
1792  const double barrelCrysSize = 0.01745; //approximate size of crystal in eta,phi in barrel
1793  const double endcapCrysSize = 0.0447; //the approximate crystal size sigmaEtaEta was corrected to in the endcap
1794 
1795  DetId seedId = getMaximum(v_id, recHits).first;
1796  bool isBarrel=seedId.subdetId()==EcalBarrel;
1797 
1798  const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1799 
1800  for (size_t i = 0; i < v_id.size(); ++i) {
1801  CaloNavigator<DetId> cursor = CaloNavigator<DetId>(v_id[i].first, topology->getSubdetectorTopology(v_id[i].first));
1802  float frac = getFraction(v_id,*cursor);
1803  float energy = recHitEnergy(*cursor, recHits,flagsexcl, severitiesexcl, sevLv)*frac;
1804 
1805  if (energy <= 0) continue;
1806 
1807  float dEta = getNrCrysDiffInEta(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.first;
1808  float dPhi = 0;
1809  if(isBarrel) dPhi = getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1810  else dPhi = getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1811 
1812 
1813 
1814  double w = 0.;
1815  w = std::max(0.0f, w0 + std::log( energy / e_5x5 ));
1816 
1817  denominator += w;
1818  numeratorEtaEta += w * dEta * dEta;
1819  numeratorEtaPhi += w * dEta * dPhi;
1820  numeratorPhiPhi += w * dPhi * dPhi;
1821  }
1822 
1823  //multiplying by crysSize to make the values compariable to normal covariances
1824  if (denominator != 0.0) {
1825  covEtaEta = crysSize*crysSize* numeratorEtaEta / denominator;
1826  covEtaPhi = crysSize*crysSize* numeratorEtaPhi / denominator;
1827  covPhiPhi = crysSize*crysSize* numeratorPhiPhi / denominator;
1828  } else {
1829  covEtaEta = 999.9;
1830  covEtaPhi = 999.9;
1831  covPhiPhi = 999.9;
1832  }
1833 
1834  } else {
1835  // Warn the user if there was no energy in the cells and return zeroes.
1836  // std::cout << "\ClusterShapeAlgo::Calculate_Covariances: no energy in supplied cells.\n";
1837  covEtaEta = 0;
1838  covEtaPhi = 0;
1839  covPhiPhi = 0;
1840  }
1841 
1842  std::vector<float> v;
1843  v.push_back( covEtaEta );
1844  v.push_back( covEtaPhi );
1845  v.push_back( covPhiPhi );
1846 
1847  return v;
1848 }
1849 
1850 // compute cluster second moments with respect to principal axes (eigenvectors of sEtaEta, sPhiPhi, sEtaPhi matrix)
1851 // store also angle alpha between major axis and phi.
1852 // takes into account shower elongation in phi direction due to magnetic field effect:
1853 // default value of 0.8 ensures sMaj = sMin for unconverted photons
1854 // (if phiCorrectionFactor=1 sMaj > sMin and alpha=0 also for unconverted photons)
1855 template<bool noZS>
1856 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments( const reco::BasicCluster &basicCluster, const EcalRecHitCollection &recHits, double phiCorrectionFactor, double w0, bool useLogWeights) {
1857 
1858  // for now implemented only for EB:
1859  // if( fabs( basicCluster.eta() ) < 1.479 ) {
1860 
1861  std::vector<std::pair<const EcalRecHit*, float> > RH_ptrs_fracs;
1862 
1863  const std::vector< std::pair<DetId, float> >& myHitsPair = basicCluster.hitsAndFractions();
1864 
1865  for(unsigned int i=0; i<myHitsPair.size(); i++){
1866  //get pointer to recHit object
1867  EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
1868  RH_ptrs_fracs.push_back( std::make_pair(&(*myRH) , myHitsPair[i].second) );
1869  }
1870 
1871  return EcalClusterToolsT<noZS>::cluster2ndMoments(RH_ptrs_fracs, phiCorrectionFactor, w0, useLogWeights);
1872 }
1873 
1874 template<bool noZS>
1875 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments( const reco::SuperCluster &superCluster, const EcalRecHitCollection &recHits, double phiCorrectionFactor, double w0, bool useLogWeights) {
1876 
1877  // for now returns second moments of supercluster seed cluster:
1878  Cluster2ndMoments returnMoments;
1879  returnMoments.sMaj = -1.;
1880  returnMoments.sMin = -1.;
1881  returnMoments.alpha = 0.;
1882 
1883  // for now implemented only for EB:
1884  // if( fabs( superCluster.eta() ) < 1.479 ) {
1885  returnMoments = EcalClusterToolsT<noZS>::cluster2ndMoments( *(superCluster.seed()), recHits, phiCorrectionFactor, w0, useLogWeights);
1886  // }
1887 
1888  return returnMoments;
1889 
1890 }
1891 
1892 template<bool noZS>
1893 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments( const std::vector<std::pair<const EcalRecHit*, float> >& RH_ptrs_fracs, double phiCorrectionFactor, double w0, bool useLogWeights) {
1894 
1895  double mid_eta(0),mid_phi(0),mid_x(0),mid_y(0);
1896 
1897  double Etot = EcalClusterToolsT<noZS>::getSumEnergy( RH_ptrs_fracs );
1898 
1899  double max_phi=-10.;
1900  double min_phi=100.;
1901 
1902 
1903  std::vector<double> etaDetId;
1904  std::vector<double> phiDetId;
1905  std::vector<double> xDetId;
1906  std::vector<double> yDetId;
1907  std::vector<double> wiDetId;
1908 
1909  unsigned int nCry=0;
1910  double denominator=0.;
1911  bool isBarrel(1);
1912 
1913  // loop over rechits and compute weights:
1914  for(std::vector<std::pair<const EcalRecHit*, float> >::const_iterator rhf_ptr = RH_ptrs_fracs.begin(); rhf_ptr != RH_ptrs_fracs.end(); rhf_ptr++){
1915  const EcalRecHit* rh_ptr = rhf_ptr->first;
1916 
1917 
1918  //get iEta, iPhi
1919  double temp_eta(0),temp_phi(0),temp_x(0),temp_y(0);
1920  isBarrel = rh_ptr->detid().subdetId()==EcalBarrel;
1921 
1922  if(isBarrel) {
1923  temp_eta = (getIEta(rh_ptr->detid()) > 0. ? getIEta(rh_ptr->detid()) + 84.5 : getIEta(rh_ptr->detid()) + 85.5);
1924  temp_phi= getIPhi(rh_ptr->detid()) - 0.5;
1925  }
1926  else {
1927  temp_eta = getIEta(rh_ptr->detid());
1928  temp_x = getNormedIX(rh_ptr->detid());
1929  temp_y = getNormedIY(rh_ptr->detid());
1930  }
1931 
1932  double temp_ene=rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
1933 
1934  double temp_wi=((useLogWeights) ?
1935  std::max(0.0, w0 + std::log( std::abs(temp_ene)/Etot ))
1936  : temp_ene);
1937 
1938 
1939  if(temp_phi>max_phi) max_phi=temp_phi;
1940  if(temp_phi<min_phi) min_phi=temp_phi;
1941  etaDetId.push_back(temp_eta);
1942  phiDetId.push_back(temp_phi);
1943  xDetId.push_back(temp_x);
1944  yDetId.push_back(temp_y);
1945  wiDetId.push_back(temp_wi);
1946  denominator+=temp_wi;
1947  nCry++;
1948  }
1949 
1950  if(isBarrel){
1951  // correct phi wrap-around:
1952  if(max_phi==359.5 && min_phi==0.5){
1953  for(unsigned int i=0; i<nCry; i++){
1954  if(phiDetId[i] - 179. > 0.) phiDetId[i]-=360.;
1955  mid_phi+=phiDetId[i]*wiDetId[i];
1956  mid_eta+=etaDetId[i]*wiDetId[i];
1957  }
1958  } else{
1959  for(unsigned int i=0; i<nCry; i++){
1960  mid_phi+=phiDetId[i]*wiDetId[i];
1961  mid_eta+=etaDetId[i]*wiDetId[i];
1962  }
1963  }
1964  }else{
1965  for(unsigned int i=0; i<nCry; i++){
1966  mid_eta+=etaDetId[i]*wiDetId[i];
1967  mid_x+=xDetId[i]*wiDetId[i];
1968  mid_y+=yDetId[i]*wiDetId[i];
1969  }
1970  }
1971 
1972  mid_eta/=denominator;
1973  mid_phi/=denominator;
1974  mid_x/=denominator;
1975  mid_y/=denominator;
1976 
1977 
1978  // See = sigma eta eta
1979  // Spp = (B field corrected) sigma phi phi
1980  // See = (B field corrected) sigma eta phi
1981  double See=0.;
1982  double Spp=0.;
1983  double Sep=0.;
1984  double deta(0),dphi(0);
1985  // compute (phi-corrected) covariance matrix:
1986  for(unsigned int i=0; i<nCry; i++) {
1987  if(isBarrel) {
1988  deta = etaDetId[i]-mid_eta;
1989  dphi = phiDetId[i]-mid_phi;
1990  } else {
1991  deta = etaDetId[i]-mid_eta;
1992  float hitLocalR2 = (xDetId[i]-mid_x)*(xDetId[i]-mid_x)+(yDetId[i]-mid_y)*(yDetId[i]-mid_y);
1993  float hitR2 = xDetId[i]*xDetId[i]+yDetId[i]*yDetId[i];
1994  float meanR2 = mid_x*mid_x+mid_y*mid_y;
1995  float hitR = sqrt(hitR2);
1996  float meanR = sqrt(meanR2);
1997  float phi = acos((hitR2+meanR2-hitLocalR2)/(2*hitR*meanR));
1998  dphi = hitR*phi;
1999 
2000  }
2001  See += (wiDetId[i]* deta * deta) / denominator;
2002  Spp += phiCorrectionFactor*(wiDetId[i]* dphi * dphi) / denominator;
2003  Sep += sqrt(phiCorrectionFactor)*(wiDetId[i]*deta*dphi) / denominator;
2004  }
2005 
2006  Cluster2ndMoments returnMoments;
2007 
2008  // compute matrix eigenvalues:
2009  returnMoments.sMaj = ((See + Spp) + sqrt((See - Spp)*(See - Spp) + 4.*Sep*Sep)) / 2.;
2010  returnMoments.sMin = ((See + Spp) - sqrt((See - Spp)*(See - Spp) + 4.*Sep*Sep)) / 2.;
2011 
2012  returnMoments.alpha = atan( (See - Spp + sqrt( (Spp - See)*(Spp - See) + 4.*Sep*Sep )) / (2.*Sep));
2013 
2014  return returnMoments;
2015 
2016 }
2017 
2018 //compute shower shapes: roundness and angle in a vector. Roundness is 0th element, Angle is 1st element.
2019 //description: uses classical mechanics inertia tensor.
2020 // roundness is smaller_eValue/larger_eValue
2021 // angle is the angle from the iEta axis to the smallest eVector (a.k.a. the shower's elongated axis)
2022 // this function uses only recHits belonging to a SC above energyThreshold (default 0)
2023 // you can select linear weighting = energy_recHit/total_energy (weightedPositionMethod=0) default
2024 // or log weighting = max( 0.0, 4.2 + log(energy_recHit/total_energy) ) (weightedPositionMethod=1)
2025 template<bool noZS>
2026 std::vector<float> EcalClusterToolsT<noZS>::roundnessBarrelSuperClusters( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int weightedPositionMethod, float energyThreshold){//int positionWeightingMethod=0){
2027  std::vector<std::pair<const EcalRecHit*, float> > RH_ptrs_fracs;
2028  const std::vector< std::pair<DetId, float> >& myHitsPair = superCluster.hitsAndFractions();
2029  for(unsigned int i=0; i< myHitsPair.size(); ++i){
2030  //get pointer to recHit object
2031  EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
2032  if( myRH != recHits.end() && myRH->energy()*(noZS ? 1.0 : myHitsPair[i].second) > energyThreshold){
2033  //require rec hit to have positive energy
2034  RH_ptrs_fracs.push_back( std::make_pair(&(*myRH) , myHitsPair[i].second) );
2035  }
2036  }
2037  std::vector<float> temp = EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(RH_ptrs_fracs,weightedPositionMethod);
2038  return temp;
2039 }
2040 
2041 // this function uses all recHits within specified window ( with default values ieta_delta=2, iphi_delta=5) around SC's highest recHit.
2042 // recHits must pass an energy threshold "energyRHThresh" (default 0)
2043 // you can select linear weighting = energy_recHit/total_energy (weightedPositionMethod=0)
2044 // or log weighting = max( 0.0, 4.2 + log(energy_recHit/total_energy) ) (weightedPositionMethod=1)
2045 template<bool noZS>
2046 std::vector<float> EcalClusterToolsT<noZS>::roundnessBarrelSuperClustersUserExtended( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int ieta_delta, int iphi_delta, float energyRHThresh, int weightedPositionMethod){
2047 
2048  std::vector<std::pair<const EcalRecHit*, float> > RH_ptrs_fracs;
2049  const std::vector< std::pair<DetId, float> >& myHitsPair = superCluster.hitsAndFractions();
2050  for(unsigned int i=0; i<myHitsPair.size(); ++i){
2051  //get pointer to recHit object
2052  EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
2053  if(myRH != recHits.end() && myRH->energy()*(noZS ? 1.0 : myHitsPair[i].second) > energyRHThresh)
2054  RH_ptrs_fracs.push_back( std::make_pair(&(*myRH) , myHitsPair[i].second) );
2055  }
2056 
2057 
2058  std::vector<int> seedPosition = EcalClusterToolsT<noZS>::getSeedPosition( RH_ptrs_fracs );
2059 
2060  for(EcalRecHitCollection::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++){
2061  EBDetId EBdetIdi( rh->detid() );
2062  float the_fraction = 0;
2063  //if(rh != recHits.end())
2064  bool inEtaWindow = ( abs( deltaIEta(seedPosition[0],EBdetIdi.ieta()) ) <= ieta_delta );
2065  bool inPhiWindow = ( abs( deltaIPhi(seedPosition[1],EBdetIdi.iphi()) ) <= iphi_delta );
2066  bool passEThresh = ( rh->energy() > energyRHThresh );
2067  bool alreadyCounted = false;
2068 
2069  // figure out if the rechit considered now is already inside the SC
2070  bool is_SCrh_inside_recHits = false;
2071  for(unsigned int i=0; i<myHitsPair.size(); i++){
2072  EcalRecHitCollection::const_iterator SCrh = recHits.find(myHitsPair[i].first);
2073  if(SCrh != recHits.end()){
2074  the_fraction = myHitsPair[i].second;
2075  is_SCrh_inside_recHits = true;
2076  if( rh->detid() == SCrh->detid() ) alreadyCounted = true;
2077  }
2078  }//for loop over SC's recHits
2079 
2080  if( is_SCrh_inside_recHits && !alreadyCounted && passEThresh && inEtaWindow && inPhiWindow){
2081  RH_ptrs_fracs.push_back( std::make_pair(&(*rh),the_fraction) );
2082  }
2083 
2084  }//for loop over rh
2085  return EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(RH_ptrs_fracs,weightedPositionMethod);
2086 }
2087 
2088 // this function computes the roundness and angle variables for vector of pointers to recHits you pass it
2089 // you can select linear weighting = energy_recHit/total_energy (weightedPositionMethod=0)
2090 // or log weighting = max( 0.0, 4.2 + log(energy_recHit/total_energy) ) (weightedPositionMethod=1)
2091 template<bool noZS>
2092 std::vector<float> EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits( const std::vector<std::pair<const EcalRecHit*,float> >& RH_ptrs_fracs, int weightedPositionMethod){//int weightedPositionMethod = 0){
2093  //positionWeightingMethod = 0 linear weighting, 1 log energy weighting
2094 
2095  std::vector<float> shapes; // this is the returning vector
2096 
2097  //make sure photon has more than one crystal; else roundness and angle suck
2098  if(RH_ptrs_fracs.size()<2){
2099  shapes.push_back( -3 );
2100  shapes.push_back( -3 );
2101  return shapes;
2102  }
2103 
2104  //Find highest E RH (Seed) and save info, compute sum total energy used
2105  std::vector<int> seedPosition = EcalClusterToolsT<noZS>::getSeedPosition( RH_ptrs_fracs );// *recHits);
2106  int tempInt = seedPosition[0];
2107  if(tempInt <0) tempInt++;
2108  float energyTotal = EcalClusterToolsT<noZS>::getSumEnergy( RH_ptrs_fracs );
2109 
2110  //1st loop over rechits: compute new weighted center position in coordinates centered on seed
2111  float centerIEta = 0.;
2112  float centerIPhi = 0.;
2113  float denominator = 0.;
2114 
2115  for(std::vector<std::pair<const EcalRecHit*,float> >::const_iterator rhf_ptr = RH_ptrs_fracs.begin(); rhf_ptr != RH_ptrs_fracs.end(); rhf_ptr++){
2116  const EcalRecHit* rh_ptr = rhf_ptr->first;
2117  //get iEta, iPhi
2118  EBDetId EBdetIdi( rh_ptr->detid() );
2119  if(fabs(energyTotal) < 0.0001){
2120  // don't /0, bad!
2121  shapes.push_back( -2 );
2122  shapes.push_back( -2 );
2123  return shapes;
2124  }
2125  float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
2126  float weight = 0;
2127  if(fabs(weightedPositionMethod)<0.0001){ //linear
2128  weight = rh_energy/energyTotal;
2129  }else{ //logrithmic
2130  weight = std::max(0.0, 4.2 + log(rh_energy/energyTotal));
2131  }
2132  denominator += weight;
2133  centerIEta += weight*deltaIEta(seedPosition[0],EBdetIdi.ieta());
2134  centerIPhi += weight*deltaIPhi(seedPosition[1],EBdetIdi.iphi());
2135  }
2136  if(fabs(denominator) < 0.0001){
2137  // don't /0, bad!
2138  shapes.push_back( -2 );
2139  shapes.push_back( -2 );
2140  return shapes;
2141  }
2142  centerIEta = centerIEta / denominator;
2143  centerIPhi = centerIPhi / denominator;
2144 
2145 
2146  //2nd loop over rechits: compute inertia tensor
2147  TMatrixDSym inertia(2); //initialize 2d inertia tensor
2148  double inertia00 = 0.;
2149  double inertia01 = 0.;// = inertia10 b/c matrix should be symmetric
2150  double inertia11 = 0.;
2151  int i = 0;
2152  for(std::vector<std::pair<const EcalRecHit*,float> >::const_iterator rhf_ptr = RH_ptrs_fracs.begin(); rhf_ptr != RH_ptrs_fracs.end(); rhf_ptr++){
2153  const EcalRecHit* rh_ptr = rhf_ptr->first;
2154  //get iEta, iPhi
2155  EBDetId EBdetIdi( rh_ptr->detid() );
2156 
2157  if(fabs(energyTotal) < 0.0001){
2158  // don't /0, bad!
2159  shapes.push_back( -2 );
2160  shapes.push_back( -2 );
2161  return shapes;
2162  }
2163  float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
2164  float weight = 0;
2165  if(fabs(weightedPositionMethod) < 0.0001){ //linear
2166  weight = rh_energy/energyTotal;
2167  }else{ //logrithmic
2168  weight = std::max(0.0, 4.2 + log(rh_energy/energyTotal));
2169  }
2170 
2171  float ieta_rh_to_center = deltaIEta(seedPosition[0],EBdetIdi.ieta()) - centerIEta;
2172  float iphi_rh_to_center = deltaIPhi(seedPosition[1],EBdetIdi.iphi()) - centerIPhi;
2173 
2174  inertia00 += weight*iphi_rh_to_center*iphi_rh_to_center;
2175  inertia01 -= weight*iphi_rh_to_center*ieta_rh_to_center;
2176  inertia11 += weight*ieta_rh_to_center*ieta_rh_to_center;
2177  i++;
2178  }
2179 
2180  inertia[0][0] = inertia00;
2181  inertia[0][1] = inertia01; // use same number here
2182  inertia[1][0] = inertia01; // and here to insure symmetry
2183  inertia[1][1] = inertia11;
2184 
2185 
2186  //step 1 find principal axes of inertia
2187  TMatrixD eVectors(2,2);
2188  TVectorD eValues(2);
2189  //std::cout<<"EcalClusterToolsT<noZS>::showerRoundness- about to compute eVectors"<<std::endl;
2190  eVectors=inertia.EigenVectors(eValues); //ordered highest eV to lowest eV (I checked!)
2191  //and eVectors are in columns of matrix! I checked!
2192  //and they are normalized to 1
2193 
2194 
2195 
2196  //step 2 select eta component of smaller eVal's eVector
2197  TVectorD smallerAxis(2);//easiest to spin SC on this axis (smallest eVal)
2198  smallerAxis[0]=eVectors[0][1];//row,col //eta component
2199  smallerAxis[1]=eVectors[1][1]; //phi component
2200 
2201  //step 3 compute interesting quatities
2202  Double_t temp = fabs(smallerAxis[0]);// closer to 1 ->beamhalo, closer to 0 something else
2203  if(fabs(eValues[0]) < 0.0001){
2204  // don't /0, bad!
2205  shapes.push_back( -2 );
2206  shapes.push_back( -2 );
2207  return shapes;
2208  }
2209 
2210  float Roundness = eValues[1]/eValues[0];
2211  float Angle=acos(temp);
2212 
2213  if( -0.00001 < Roundness && Roundness < 0) Roundness = 0.;
2214  if( -0.00001 < Angle && Angle < 0 ) Angle = 0.;
2215 
2216  shapes.push_back( Roundness );
2217  shapes.push_back( Angle );
2218  return shapes;
2219 
2220 }
2221 //private functions useful for roundnessBarrelSuperClusters etc.
2222 //compute delta iphi between a seed and a particular recHit
2223 //iphi [1,360]
2224 //safe gaurds are put in to ensure the difference is between [-180,180]
2225 template<bool noZS>
2226 int EcalClusterToolsT<noZS>::deltaIPhi(int seed_iphi, int rh_iphi){
2227  int rel_iphi = rh_iphi - seed_iphi;
2228  // take care of cyclic variable iphi [1,360]
2229  if(rel_iphi > 180) rel_iphi = rel_iphi - 360;
2230  if(rel_iphi < -180) rel_iphi = rel_iphi + 360;
2231  return rel_iphi;
2232 }
2233 
2234 //compute delta ieta between a seed and a particular recHit
2235 //ieta [-85,-1] and [1,85]
2236 //safe gaurds are put in to shift the negative ieta +1 to make an ieta=0 so differences are computed correctly
2237 template<bool noZS>
2238 int EcalClusterToolsT<noZS>::deltaIEta(int seed_ieta, int rh_ieta){
2239  // get rid of the fact that there is no ieta=0
2240  if(seed_ieta < 0) seed_ieta++;
2241  if(rh_ieta < 0) rh_ieta++;
2242  int rel_ieta = rh_ieta - seed_ieta;
2243  return rel_ieta;
2244 }
2245 
2246 //return (ieta,iphi) of highest energy recHit of the recHits passed to this function
2247 template<bool noZS>
2248 std::vector<int> EcalClusterToolsT<noZS>::getSeedPosition(const std::vector<std::pair<const EcalRecHit*, float> >& RH_ptrs_fracs){
2249  std::vector<int> seedPosition;
2250  float eSeedRH = 0;
2251  int iEtaSeedRH = 0;
2252  int iPhiSeedRH = 0;
2253 
2254  for(std::vector<std::pair<const EcalRecHit*,float> >::const_iterator rhf_ptr = RH_ptrs_fracs.begin(); rhf_ptr != RH_ptrs_fracs.end(); rhf_ptr++){
2255  const EcalRecHit* rh_ptr = rhf_ptr->first;
2256  //get iEta, iPhi
2257  EBDetId EBdetIdi( rh_ptr->detid() );
2258  float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
2259 
2260  if(eSeedRH < rh_energy){
2261  eSeedRH = rh_energy;
2262  iEtaSeedRH = EBdetIdi.ieta();
2263  iPhiSeedRH = EBdetIdi.iphi();
2264  }
2265 
2266  }// for loop
2267 
2268  seedPosition.push_back(iEtaSeedRH);
2269  seedPosition.push_back(iPhiSeedRH);
2270  return seedPosition;
2271 }
2272 
2273 // return the total energy of the recHits passed to this function
2274 template<bool noZS>
2275 float EcalClusterToolsT<noZS>::getSumEnergy(const std::vector<std::pair<const EcalRecHit*, float> >& RH_ptrs_fracs){
2276  float sumE = 0.;
2277  for( const auto& hAndF : RH_ptrs_fracs ) {
2278  sumE += hAndF.first->energy() * (noZS ? 1.0 : hAndF.second);
2279  }
2280  return sumE;
2281 }
2282 
2283 typedef EcalClusterToolsT<false> EcalClusterTools;
2284 
2285 namespace noZS {
2286  typedef EcalClusterToolsT<true> EcalClusterTools;
2287 }
2288 
2289 #endif
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
EcalClusterToolsT::zernike20
static double zernike20(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, double R0=6.6, bool logW=true, float w0=4.7)
Definition: EcalClusterTools.h:1482
CaloGeometry::getSubdetectorGeometry
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:43
i
int i
Definition: DBlmapReader.cc:9
EEDetId::ix
int ix() const
Definition: EEDetId.h:76
EcalClusterToolsT::e2x5Bottom
static float e2x5Bottom(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:664
EcalClusterToolsT::meanClusterPosition
static math::XYZVector meanClusterPosition(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, const CaloGeometry *geometry)
Definition: EcalClusterTools.h:1015
EcalClusterToolsT::covariances
static std::vector< float > covariances(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, const CaloGeometry *geometry, float w0=4.7)
Definition: EcalClusterTools.h:1175
EcalClusterTools
EcalClusterToolsT< false > EcalClusterTools
Definition: EcalClusterTools.h:2283
CaloRecHit::detid
const DetId & detid() const
Definition: CaloRecHit.h:20
EcalClusterToolsT::matrixDetId
static std::vector< DetId > matrixDetId(const CaloTopology *topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax)
Definition: EcalClusterTools.h:463
EcalClusterToolsT::eMax
static float eMax(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits)
Definition: EcalClusterTools.h:577
EcalClusterToolsT::roundnessBarrelSuperClustersUserExtended
static std::vector< float > roundnessBarrelSuperClustersUserExtended(const reco::SuperCluster &superCluster, const EcalRecHitCollection &recHits, int ieta_delta=0, int iphi_delta=0, float energyRHThresh=0.00000, int weightedPositionMethod=0)
Definition: EcalClusterTools.h:2046
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
EcalClusterToolsT::f33
static double f33(double r)
Definition: EcalClusterTools.h:252
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
create_public_lumi_plots.log
log
Definition: create_public_lumi_plots.py:1102
edm::SortedCollection< EcalRecHit >::const_iterator
std::vector< EcalRecHit >::const_iterator const_iterator
Definition: SortedCollection.h:80
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
EcalClusterToolsT::f11
static double f11(double r)
Definition: EcalClusterTools.h:248
reco::CaloCluster::hitsAndFractions
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
Definition: CaloCluster.h:190
EBDetId::kTowersInPhi
static const int kTowersInPhi
Definition: EBDetId.h:146
EcalClusterToolsT::deltaIEta
static int deltaIEta(int seed_ieta, int rh_ieta)
Definition: EcalClusterTools.h:2238
EcalClusterToolsT::getSeedPosition
static std::vector< int > getSeedPosition(const std::vector< std::pair< const EcalRecHit *, float > > &RH_ptrs)
Definition: EcalClusterTools.h:2248
spr::find
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
EcalClusterToolsT::calc_AbsZernikeMoment
static double calc_AbsZernikeMoment(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0)
Definition: EcalClusterTools.h:1543
EcalClusterToolsT::deltaIPhi
static int deltaIPhi(int seed_iphi, int rh_iphi)
Definition: EcalClusterTools.h:2226
EcalClusterToolsT::getEnergyDepTopology
static std::vector< EcalClusterEnergyDeposition > getEnergyDepTopology(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW, float w0)
Definition: EcalClusterTools.h:893
position
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
EcalClusterToolsT::roundnessBarrelSuperClusters
static std::vector< float > roundnessBarrelSuperClusters(const reco::SuperCluster &superCluster, const EcalRecHitCollection &recHits, int weightedPositionMethod=0, float energyThreshold=0.0)
Definition: EcalClusterTools.h:2026
EcalClusterToolsT::e2x5Top
static float e2x5Top(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:650
CaloSubdetectorGeometry::getGeometry
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
Definition: CaloSubdetectorGeometry.cc:44
EBDetId::iphi
int iphi() const
get the crystal iphi
Definition: EBDetId.h:53
EcalClusterToolsT::fast_AbsZernikeMoment
static double fast_AbsZernikeMoment(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0)
Definition: EcalClusterTools.h:1507
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:250
EBDetId::im
int im() const
get the number of module inside the SM (1-4)
Definition: EBDetId.h:66
EBDetId::kCrystalsInPhi
static const int kCrystalsInPhi
Definition: EBDetId.h:149
cuy.denominator
list denominator
Definition: cuy.py:484
EcalClusterToolsT::e3x1
static float e3x1(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:755
dPhi
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
EcalClusterToolsT::getDPhiEndcap
static float getDPhiEndcap(const DetId &crysId, float meanX, float meanY)
Definition: EcalClusterTools.h:1676
EcalClusterToolsT::getIEta
static float getIEta(const DetId &id)
Definition: EcalClusterTools.h:1577
EcalClusterToolsT::f55
static double f55(double r)
Definition: EcalClusterTools.h:258
EcalClusterToolsT::e2x2
static float e2x2(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:479
EcalClusterToolsT::getIPhi
static float getIPhi(const DetId &id)
Definition: EcalClusterTools.h:1599
CaloRecHit::energy
float energy() const
Definition: CaloRecHit.h:17
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
Geom::Phi::phi
T phi() const
Definition: Phi.h:41
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
EcalClusterToolsT::getMaximum
static std::pair< DetId, float > getMaximum(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits)
Definition: EcalClusterTools.h:334
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
EcalClusterToolsT::e2nd
static float e2nd(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits)
Definition: EcalClusterTools.h:589
EBDetId::kModulesPerSM
static const int kModulesPerSM
Definition: EBDetId.h:147
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
j
int j
Definition: DBlmapReader.cc:9
EcalClusterToolsT::f40
static double f40(double r)
Definition: EcalClusterTools.h:253
EcalClusterToolsT::e2x5Max
static float e2x5Max(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:680
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
EcalClusterToolsT::eBottom
static float eBottom(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:811
EcalClusterToolsT::getNrCrysDiffInEta
static float getNrCrysDiffInEta(const DetId &crysId, const DetId &orginId)
Definition: EcalClusterTools.h:1638
EcalClusterToolsT::e2x5Right
static float e2x5Right(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:621
EEDetId::iy
int iy() const
Definition: EEDetId.h:82
EcalClusterToolsT::scLocalCovariances
static std::vector< float > scLocalCovariances(const reco::SuperCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, float w0=4.7)
Definition: EcalClusterTools.h:1697
EcalClusterToolsT::getFraction
static float getFraction(const std::vector< std::pair< DetId, float > > &v_id, DetId id)
Definition: EcalClusterTools.h:290
EBDetId::ieta
int ieta() const
get the crystal ieta
Definition: EBDetId.h:51
EcalClusterToolsT::e3x2
static float e3x2(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:502
EcalClusterToolsT::f31
static double f31(double r)
Definition: EcalClusterTools.h:251
EcalClusterToolsT::getNormedIY
static float getNormedIY(const DetId &id)
Definition: EcalClusterTools.h:1625
first
bool first
Definition: L1TdeRCT.cc:79
DetId::subdetId
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
EcalClusterToolsT::factorial
static double factorial(int n)
Definition: EcalClusterTools.h:264
EcalClusterToolsT::f44
static double f44(double r)
Definition: EcalClusterTools.h:255
EcalClusterToolsT::absZernikeMoment
static double absZernikeMoment(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0)
Definition: EcalClusterTools.h:1494
EcalClusterToolsT::getSumEnergy
static float getSumEnergy(const std::vector< std::pair< const EcalRecHit *, float > > &RH_ptrs_fracs)
Definition: EcalClusterTools.h:2275
EcalClusterToolsT::matrixEnergy
static float matrixEnergy(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax)
Definition: EcalClusterTools.h:419
EcalClusterToolsT::getNrCrysDiffInPhi
static float getNrCrysDiffInPhi(const DetId &crysId, const DetId &orginId)
Definition: EcalClusterTools.h:1659
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:301
EcalClusterToolsT::mean5x5PositionInXY
static std::pair< float, float > mean5x5PositionInXY(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:1119
EcalClusterToolsT::cluster2ndMoments
static Cluster2ndMoments cluster2ndMoments(const reco::BasicCluster &basicCluster, const EcalRecHitCollection &recHits, double phiCorrectionFactor=0.8, double w0=4.7, bool useLogWeights=true)
Definition: EcalClusterTools.h:1856
EcalClusterToolsT::f00
static double f00(double r)
Definition: EcalClusterTools.h:247
EcalClusterToolsT::roundnessSelectedBarrelRecHits
static std::vector< float > roundnessSelectedBarrelRecHits(const std::vector< std::pair< const EcalRecHit *, float > > &rhVector, int weightedPositionMethod=0)
Definition: EcalClusterTools.h:2092
reco::e1
Float e1
Definition: deltaR.h:22
EcalClusterToolsT::eTop
static float eTop(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:797
DetId
Definition: DetId.h:18
EcalClusterToolsT::e4x4
static float e4x4(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:539
EcalClusterToolsT::e5x1
static float e5x1(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:727
EcalClusterToolsT::f53
static double f53(double r)
Definition: EcalClusterTools.h:257
M_PI
#define M_PI
Definition: BFit3D.cc:3
EcalClusterToolsT::energyBasketFractionEta
static std::vector< float > energyBasketFractionEta(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits)
Definition: EcalClusterTools.h:825
EBDetId::positiveZ
bool positiveZ() const
Definition: EBDetId.h:78
EBDetId::MAX_IPHI
static const int MAX_IPHI
Definition: EBDetId.h:144
noZS::EcalClusterTools
EcalClusterToolsT< true > EcalClusterTools
Definition: EcalClusterTools.h:2286
math::XYZVector
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:30
EcalClusterToolsT::lat
static std::vector< float > lat(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW=true, float w0=4.7)
Definition: EcalClusterTools.h:953
EcalClusterToolsT::recHitEnergy
static float recHitEnergy(DetId id, const EcalRecHitCollection *recHits)
Definition: EcalClusterTools.h:347
EcalClusterToolsT::e2x5Left
static float e2x5Left(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:635
EcalClusterToolsT::energyBasketFractionPhi
static std::vector< float > energyBasketFractionPhi(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits)
Definition: EcalClusterTools.h:859
CaloTopology::getSubdetectorTopology
const CaloSubdetectorTopology * getSubdetectorTopology(const DetId &id) const
access the subdetector Topology for the given subdetector directly
Definition: CaloTopology.cc:26
EcalClusterToolsT::f42
static double f42(double r)
Definition: EcalClusterTools.h:254
reco::e2
Float e2
Definition: deltaR.h:23
EcalClusterToolsT::f51
static double f51(double r)
Definition: EcalClusterTools.h:256
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:76
EcalClusterToolsT::f20
static double f20(double r)
Definition: EcalClusterTools.h:249
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
Angle
Definition: Angle.h:17
EcalClusterToolsT::f22
static double f22(double r)
Definition: EcalClusterTools.h:250
geometry
ESHandle< TrackerGeometry > geometry
Definition: TkLasBeamFitter.cc:176
edm::SortedCollection::find
iterator find(key_type k)
Definition: SortedCollection.h:262
Geom::pi
double pi()
Definition: Pi.h:31
EcalClusterToolsT::mean5x5PositionInLocalCrysCoord
static std::pair< float, float > mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:1061
Geom::twoPi
double twoPi()
Definition: Pi.h:32
factorial
int factorial(int n)
factorial function
Definition: L1GlobalTriggerFunctions.cc:13
EcalClusterToolsT::eRight
static float eRight(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:783
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
EcalClusterToolsT::e3x3
static float e3x3(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:525
reco::SuperCluster::seed
const CaloClusterPtr & seed() const
seed BasicCluster
Definition: SuperCluster.h:66
EcalClusterToolsT::zernike42
static double zernike42(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, double R0=6.6, bool logW=true, float w0=4.7)
Definition: EcalClusterTools.h:1488
w
T w() const
Definition: extBasic3DVector.h:245
EcalClusterToolsT::e1x3
static float e1x3(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:741
CastorDataFrameFilter_impl::energySum
double energySum(const DataFrame &df, int fs, int ls)
Definition: CastorDataFrameFilter.cc:22
histoStyle.weight
int weight
Definition: histoStyle.py:50
EcalClusterToolsT::getNormedIX
static float getNormedIX(const DetId &id)
Definition: EcalClusterTools.h:1612
CaloCellGeometry::getPosition
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
Definition: CaloCellGeometry.h:76
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
EcalClusterToolsT::eLeft
static float eLeft(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:769
EcalClusterToolsT::localCovariances
static std::vector< float > localCovariances(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, float w0=4.7)
Definition: EcalClusterTools.h:1317
EcalSeverityLevelAlgo::severityLevel
EcalSeverityLevel::SeverityLevel severityLevel(const DetId &id, const EcalRecHitCollection &rhs) const
Evaluate status from id.
Definition: EcalSeverityLevelAlgo.cc:82
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
EcalClusterToolsT::e1x5
static float e1x5(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:712
EcalClusterToolsT::computeWeight
static float computeWeight(float eRH, float energyTotal, int weightedPositionMethod)
edm::SortedCollection::begin
const_iterator begin() const
Definition: SortedCollection.h:294
EcalClusterToolsT::e5x5
static float e5x5(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:563
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5