CMS 3D CMS Logo
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 
88  static float e3x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
89 
90 
91  static float e1x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
92 
93  static float e5x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
94 
95  static float e2x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
96 
97  static float e3x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
98 
99  static float e3x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
100 
101  static float e4x4( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology);
102 
103  static float e5x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
104  static int n5x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
105 
106  // energy in the 2x5 strip right of the max crystal (does not contain max crystal)
107  // 2 crystals wide in eta, 5 wide in phi.
108  static float e2x5Right( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
109  // energy in the 2x5 strip left of the max crystal (does not contain max crystal)
110 
111  static float e2x5Left( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
112  // energy in the 5x2 strip above the max crystal (does not contain max crystal)
113  // 5 crystals wide in eta, 2 wide in phi.
114 
115  static float e2x5Top( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
116  // energy in the 5x2 strip below the max crystal (does not contain max crystal)
117 
118  static float e2x5Bottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
119  // energy in a 2x5 strip containing the seed (max) crystal.
120  // 2 crystals wide in eta, 5 wide in phi.
121  // it is the maximum of either (1x5left + 1x5center) or (1x5right + 1x5center)
122  static float e2x5Max( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
123 
124  // energies in the crystal left, right, top, bottom w.r.t. to the most energetic crystal
125  static float eLeft( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
126 
127  static float eRight( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
128 
129  static float eTop( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
130 
131  static float eBottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology );
132  // the energy of the most energetic crystal in the cluster
133 
134  static float eMax( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits );
135 
136  // the energy of the second most energetic crystal in the cluster
137  static float e2nd( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits );
138 
139  // get the DetId and the energy of the maximum energy crystal of the input cluster
140  static std::pair<DetId, float> getMaximum( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
141 
142  static std::vector<float> energyBasketFractionEta( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits );
143 
144  static std::vector<float> energyBasketFractionPhi( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
145 
146  // return a vector v with v[0] = etaLat, v[1] = phiLat, v[2] = lat
147  static std::vector<float> lat( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW = true, float w0 = 4.7 );
148 
149  // return a vector v with v[0] = covEtaEta, v[1] = covEtaPhi, v[2] = covPhiPhi
150 
151  static std::vector<float> covariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, const CaloGeometry* geometry, float w0 = 4.7);
152 
153  // return a vector v with v[0] = covIEtaIEta, v[1] = covIEtaIPhi, v[2] = covIPhiIPhi
154  //this function calculates differences in eta/phi in units of crystals not global eta/phi
155  //this is gives better performance in the crack regions of the calorimeter but gives otherwise identical results to covariances function
156  // except that it doesnt need an eta based correction funtion in the endcap
157  //it is multipled by an approprate crystal size to ensure it gives similar values to covariances(...)
158  //
159  //Warning: covIEtaIEta has been studied by egamma, but so far covIPhiIPhi hasnt been studied extensively so there could be a bug in
160  // the covIPhiIEta or covIPhiIPhi calculations. I dont think there is but as it hasnt been heavily used, there might be one
161  static std::vector<float> localCovariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, float w0 = 4.7);
162 
163  static std::vector<float> scLocalCovariances(const reco::SuperCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, float w0 = 4.7);
164 
165  // cluster second moments with respect to principal axes:
166  static Cluster2ndMoments cluster2ndMoments( const reco::BasicCluster &basicCluster, const EcalRecHitCollection &recHits, double phiCorrectionFactor=0.8, double w0=4.7, bool useLogWeights=true);
167 
168  static Cluster2ndMoments cluster2ndMoments( const reco::SuperCluster &superCluster, const EcalRecHitCollection &recHits, double phiCorrectionFactor=0.8, double w0=4.7, bool useLogWeights=true);
169  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);
170 
171  static double zernike20( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, double R0 = 6.6, bool logW = true, float w0 = 4.7 );
172  static double zernike42( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, double R0 = 6.6, bool logW = true, float w0 = 4.7 );
173 
174  // get the detId's of a matrix centered in the maximum energy crystal = (0,0)
175  // the size is specified by ixMin, ixMax, iyMin, iyMax in unit of crystals
176  static std::vector<DetId> matrixDetId( const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax );
177 
178  // get the energy deposited in a matrix centered in the maximum energy crystal = (0,0)
179  // the size is specified by ixMin, ixMax, iyMin, iyMax in unit of crystals
180  static float matrixEnergy( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax );
181  static int matrixSize( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax );
182 
183  static float getFraction( const std::vector< std::pair<DetId, float> > &v_id, DetId id);
184  // get the DetId and the energy of the maximum energy crystal in a vector of DetId
185  static std::pair<DetId, float> getMaximum( const std::vector< std::pair<DetId, float> > &v_id, const EcalRecHitCollection *recHits);
186 
187  // get the energy of a DetId, return 0 if the DetId is not in the collection
188  static float recHitEnergy(DetId id, const EcalRecHitCollection *recHits);
189 
190  //Shower shape variables return vector <Roundness, Angle> of a photon
191  static std::vector<float> roundnessBarrelSuperClusters( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int weightedPositionMethod = 0, float energyThreshold = 0.0);
192  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);
193  static std::vector<float> roundnessSelectedBarrelRecHits(const std::vector<std::pair<const EcalRecHit*,float> >&rhVector, int weightedPositionMethod = 0);
194 
195  //works out the number of staturated crystals in 5x5
196  static int nrSaturatedCrysIn5x5(const DetId& id,const EcalRecHitCollection* recHits,const CaloTopology *topology);
197  private:
198  struct EcalClusterEnergyDeposition
199  {
200  double deposited_energy;
201  double r;
202  double phi;
203  };
204 
205  static std::vector<EcalClusterEnergyDeposition> getEnergyDepTopology( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW, float w0 );
206 
207  static math::XYZVector meanClusterPosition( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, const CaloGeometry *geometry );
208 
209  //return energy weighted mean distance from the seed crystal in number of crystals
210  //<iEta,iPhi>, iPhi is not defined for endcap and is returned as zero
211  static std::pair<float,float> mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology);
212 
213  static std::pair<float,float> mean5x5PositionInXY(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology);
214 
215  static double f00(double r) { return 1; }
216  static double f11(double r) { return r; }
217  static double f20(double r) { return 2.0*r*r-1.0; }
218  static double f22(double r) { return r*r; }
219  static double f31(double r) { return 3.0*r*r*r - 2.0*r; }
220  static double f33(double r) { return r*r*r; }
221  static double f40(double r) { return 6.0*r*r*r*r-6.0*r*r+1.0; }
222  static double f42(double r) { return 4.0*r*r*r*r-3.0*r*r; }
223  static double f44(double r) { return r*r*r*r; }
224  static double f51(double r) { return 10.0*pow(r,5)-12.0*pow(r,3)+3.0*r; }
225  static double f53(double r) { return 5.0*pow(r,5) - 4.0*pow(r,3); }
226  static double f55(double r) { return pow(r,5); }
227 
228  static double absZernikeMoment( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 );
229  static double fast_AbsZernikeMoment(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 );
230  static double calc_AbsZernikeMoment(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 );
231 
232  static double factorial(int n) {
233  double res = 1.;
234  for (int i = 2; i <= n; ++i) res *= i;
235  return res;
236  }
237 
238  //useful functions for the localCovariances function
239  static float getIEta(const DetId& id);
240  static float getIPhi(const DetId& id);
241  static float getNormedIX(const DetId& id);
242  static float getNormedIY(const DetId& id);
243  static float getDPhiEndcap(const DetId& crysId,float meanX,float meanY);
244  static float getNrCrysDiffInEta(const DetId& crysId,const DetId& orginId);
245  static float getNrCrysDiffInPhi(const DetId& crysId,const DetId& orginId);
246 
247  //useful functions for showerRoundnessBarrel function
248  static int deltaIEta(int seed_ieta, int rh_ieta);
249  static int deltaIPhi(int seed_iphi, int rh_iphi);
250  static std::vector<int> getSeedPosition(const std::vector<std::pair<const EcalRecHit*,float> >&RH_ptrs);
251  static float getSumEnergy(const std::vector<std::pair<const EcalRecHit*,float> >&RH_ptrs_fracs);
252  static float computeWeight(float eRH, float energyTotal, int weightedPositionMethod);
253 
254 
255 };
256 
257 // implementation
258 template<bool noZS>
259 float EcalClusterToolsT<noZS>::getFraction( const std::vector< std::pair<DetId, float> > &v_id, DetId id
260  ){
261  if(noZS) return 1.0;
262  float frac = 0.0;
263  for ( size_t i = 0; i < v_id.size(); ++i ) {
264  if(v_id[i].first.rawId()==id.rawId()){
265  frac= v_id[i].second;
266  break;
267  }
268  }
269  return frac;
270 }
271 
272 template<bool noZS>
273 std::pair<DetId, float> EcalClusterToolsT<noZS>::getMaximum( const std::vector< std::pair<DetId, float> > &v_id, const EcalRecHitCollection *recHits)
274 {
275  float max = 0;
276  DetId id(0);
277  for ( size_t i = 0; i < v_id.size(); ++i ) {
278  float energy = recHitEnergy( v_id[i].first, recHits ) * (noZS ? 1.0 : v_id[i].second);
279  if ( energy > max ) {
280  max = energy;
281  id = v_id[i].first;
282  }
283  }
284  return std::pair<DetId, float>(id, max);
285 }
286 
287 template<bool noZS>
288 std::pair<DetId, float> EcalClusterToolsT<noZS>::getMaximum( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits)
289 {
290  return getMaximum( cluster.hitsAndFractions(), recHits );
291 }
292 
293 
294 template<bool noZS>
295 float EcalClusterToolsT<noZS>::recHitEnergy(DetId id, const EcalRecHitCollection *recHits)
296 {
297  if ( id == DetId(0) ) {
298  return 0;
299  } else {
300  EcalRecHitCollection::const_iterator it = recHits->find( id );
301  if ( it != recHits->end() ) {
302  if( noZS && ( it->checkFlag(EcalRecHit::kTowerRecovered) ||
303  it->checkFlag(EcalRecHit::kWeird) ||
304  (it->detid().subdetId() == EcalBarrel &&
305  it->checkFlag(EcalRecHit::kDiWeird) )
306  )
307  ) {
308  return 0.0;
309  } else {
310  return (*it).energy();
311  }
312  } else {
313  //throw cms::Exception("EcalRecHitNotFound") << "The recHit corresponding to the DetId" << id.rawId() << " not found in the EcalRecHitCollection";
314  // the recHit is not in the collection (hopefully zero suppressed)
315  return 0;
316  }
317  }
318  return 0;
319 }
320 
321 
322 // Returns the energy in a rectangle of crystals
323 // specified in eta by ixMin and ixMax
324 // and in phi by iyMin and iyMax
325 //
326 // Reference picture (X=seed crystal)
327 // iy ___________
328 // 2 |_|_|_|_|_|
329 // 1 |_|_|_|_|_|
330 // 0 |_|_|X|_|_|
331 // -1 |_|_|_|_|_|
332 // -2 |_|_|_|_|_|
333 // -2 -1 0 1 2 ix
334 template<bool noZS>
335 float EcalClusterToolsT<noZS>::matrixEnergy( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax )
336 {
337  //take into account fractions
338  // fast version
339  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( id, topology->getSubdetectorTopology( id ) );
340  float energy = 0;
341  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
342  for ( int i = ixMin; i <= ixMax; ++i ) {
343  for ( int j = iyMin; j <= iyMax; ++j ) {
344  cursor.home();
345  cursor.offsetBy( i, j );
346  float frac=getFraction(v_id,*cursor);
347  energy += recHitEnergy( *cursor, recHits )*frac;
348  }
349  }
350  // slow elegant version
351  //float energy = 0;
352  //std::vector<DetId> v_id = matrixDetId( topology, id, ixMin, ixMax, iyMin, iyMax );
353  //for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
354  // energy += recHitEnergy( *it, recHits );
355  //}
356  return energy;
357 }
358 
359 template<bool noZS>
360 int EcalClusterToolsT<noZS>::matrixSize( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax )
361 {
362  // fast version
363  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( id, topology->getSubdetectorTopology( id ) );
364  int result = 0;
365  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
366  for ( int i = ixMin; i <= ixMax; ++i ) {
367  for ( int j = iyMin; j <= iyMax; ++j ) {
368  cursor.home();
369  cursor.offsetBy( i, j );
370  float frac=getFraction(v_id,*cursor);
371  float energy = recHitEnergy( *cursor, recHits )*frac;
372  if (energy > 0) result++;
373  }
374  }
375  return result;
376 }
377 
378 
379 template<bool noZS>
380 std::vector<DetId> EcalClusterToolsT<noZS>::matrixDetId( const CaloTopology* topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax )
381 {
382  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( id, topology->getSubdetectorTopology( id ) );
383  std::vector<DetId> v;
384  for ( int i = ixMin; i <= ixMax; ++i ) {
385  for ( int j = iyMin; j <= iyMax; ++j ) {
386  cursor.home();
387  cursor.offsetBy( i, j );
388  if ( *cursor != DetId(0) ) v.push_back( *cursor );
389  }
390  }
391  return v;
392 }
393 
394 
395 template<bool noZS>
396 float EcalClusterToolsT<noZS>::e2x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
397 {
398  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
399  std::list<float> energies;
400  float max_E = matrixEnergy( cluster, recHits, topology, id, -1, 0, -1, 0 );
401  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -1, 0, 0, 1 ) );
402  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, 0, 1, 0, 1 ) );
403  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, 0, 1, -1, 0 ) );
404  return max_E;
405 }
406 
407 template<bool noZS>
408 float EcalClusterToolsT<noZS>::e3x2( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
409 {
410  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
411  float max_E = matrixEnergy( cluster, recHits, topology, id, -1, 1, -1, 0 );
412  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, 0, 1, -1, 1 ) );
413  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -1, 1, 0, 1 ) );
414  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -1, 0, -1, 1 ) );
415  return max_E;
416 }
417 
418 template<bool noZS>
419 float EcalClusterToolsT<noZS>::e3x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
420 {
421  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
422  return matrixEnergy( cluster, recHits, topology, id, -1, 1, -1, 1 );
423 }
424 
425 template<bool noZS>
426 float EcalClusterToolsT<noZS>::e4x4( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
427 {
428  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
429  float max_E = matrixEnergy( cluster, recHits, topology, id, -1, 2, -2, 1 );
430  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -2, 1, -2, 1 ) );
431  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -2, 1, -1, 2 ) );
432  max_E = std::max( max_E, matrixEnergy( cluster, recHits, topology, id, -1, 2, -1, 2 ) );
433  return max_E;
434 }
435 
436 template<bool noZS>
437 float EcalClusterToolsT<noZS>::e5x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
438 {
439  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
440  return matrixEnergy( cluster, recHits, topology, id, -2, 2, -2, 2 );
441 }
442 
443 template<bool noZS>
444 int EcalClusterToolsT<noZS>::n5x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
445 {
446  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
447  return matrixSize( cluster, recHits, topology, id, -2, 2, -2, 2 );
448 }
449 
450 template<bool noZS>
451 float EcalClusterToolsT<noZS>::eMax( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits )
452 {
453  return getMaximum( cluster.hitsAndFractions(), recHits ).second;
454 }
455 
456 template<bool noZS>
457 float EcalClusterToolsT<noZS>::e2nd( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits )
458 {
459  std::vector<float> energies;
460  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
461  energies.reserve( v_id.size() );
462  if ( v_id.size() < 2 ) return 0;
463  for ( size_t i = 0; i < v_id.size(); ++i ) {
464  energies.push_back( recHitEnergy( v_id[i].first, recHits ) * (noZS ? 1.0 : v_id[i].second) );
465  }
466  std::partial_sort( energies.begin(), energies.begin()+2, energies.end(), std::greater<float>() );
467  return energies[1];
468 
469 
470 }
471 
472 template<bool noZS>
473 float EcalClusterToolsT<noZS>::e2x5Right( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
474 {
475  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
476  return matrixEnergy( cluster, recHits, topology, id, 1, 2, -2, 2 );
477 }
478 
479 template<bool noZS>
480 float EcalClusterToolsT<noZS>::e2x5Left( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
481 {
482  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
483  return matrixEnergy( cluster, recHits, topology, id, -2, -1, -2, 2 );
484 }
485 
486 template<bool noZS>
487 float EcalClusterToolsT<noZS>::e2x5Top( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
488 {
489  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
490  return matrixEnergy( cluster, recHits, topology, id, -2, 2, 1, 2 );
491 }
492 
493 template<bool noZS>
494 float EcalClusterToolsT<noZS>::e2x5Bottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
495 {
496  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
497  return matrixEnergy( cluster, recHits, topology, id, -2, 2, -2, -1 );
498 }
499 
500 // Energy in 2x5 strip containing the max crystal.
501 // Adapted from code by Sam Harper
502 template<bool noZS>
503 float EcalClusterToolsT<noZS>::e2x5Max( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
504 {
505  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
506 
507  // 1x5 strip left of seed
508  float left = matrixEnergy( cluster, recHits, topology, id, -1, -1, -2, 2 );
509  // 1x5 strip right of seed
510  float right = matrixEnergy( cluster, recHits, topology, id, 1, 1, -2, 2 );
511  // 1x5 strip containing seed
512  float centre = matrixEnergy( cluster, recHits, topology, id, 0, 0, -2, 2 );
513 
514  // Return the maximum of (left+center) or (right+center) strip
515  return left > right ? left+centre : right+centre;
516 }
517 
518 template<bool noZS>
519 float EcalClusterToolsT<noZS>::e1x5( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
520 {
521  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
522  return matrixEnergy( cluster, recHits, topology, id, 0, 0, -2, 2 );
523 }
524 
525 template<bool noZS>
526 float EcalClusterToolsT<noZS>::e5x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
527 {
528  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
529  return matrixEnergy( cluster, recHits, topology, id, -2, 2, 0, 0 );
530 }
531 
532 template<bool noZS>
533 float EcalClusterToolsT<noZS>::e1x3( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
534 {
535  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
536  return matrixEnergy( cluster, recHits, topology, id, 0, 0, -1, 1 );
537 }
538 
539 template<bool noZS>
540 float EcalClusterToolsT<noZS>::e3x1( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
541 {
542  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
543  return matrixEnergy( cluster, recHits, topology, id, -1, 1, 0, 0 );
544 }
545 
546 template<bool noZS>
547 float EcalClusterToolsT<noZS>::eLeft( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
548 {
549  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
550  return matrixEnergy( cluster, recHits, topology, id, -1, -1, 0, 0 );
551 }
552 
553 template<bool noZS>
554 float EcalClusterToolsT<noZS>::eRight( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
555 {
556  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
557  return matrixEnergy( cluster, recHits, topology, id, 1, 1, 0, 0 );
558 }
559 
560 template<bool noZS>
561 float EcalClusterToolsT<noZS>::eTop( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
562 {
563  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
564  return matrixEnergy( cluster, recHits, topology, id, 0, 0, 1, 1 );
565 }
566 
567 template<bool noZS>
568 float EcalClusterToolsT<noZS>::eBottom( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology* topology )
569 {
570  DetId id = getMaximum( cluster.hitsAndFractions(), recHits ).first;
571  return matrixEnergy( cluster, recHits, topology, id, 0, 0, -1, -1 );
572 }
573 
574 template<bool noZS>
575 std::vector<float> EcalClusterToolsT<noZS>::energyBasketFractionEta( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits )
576 {
577  std::vector<float> basketFraction( 2 * EBDetId::kModulesPerSM );
578  float clusterEnergy = cluster.energy();
579  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
580  if ( v_id[0].first.subdetId() != EcalBarrel ) {
581  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.";
582  return basketFraction;
583  }
584  for ( size_t i = 0; i < v_id.size(); ++i ) {
585  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;
586  }
587  std::sort( basketFraction.rbegin(), basketFraction.rend() );
588  return basketFraction;
589 }
590 
591 template<bool noZS>
592 std::vector<float> EcalClusterToolsT<noZS>::energyBasketFractionPhi( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits )
593 {
594  std::vector<float> basketFraction( 2 * (EBDetId::MAX_IPHI / EBDetId::kCrystalsInPhi) );
595  float clusterEnergy = cluster.energy();
596  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
597  if ( v_id[0].first.subdetId() != EcalBarrel ) {
598  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.";
599  return basketFraction;
600  }
601  for ( size_t i = 0; i < v_id.size(); ++i ) {
602  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;
603  }
604  std::sort( basketFraction.rbegin(), basketFraction.rend() );
605  return basketFraction;
606 }
607 
608 template<bool noZS>
609 std::vector<typename EcalClusterToolsT<noZS>::EcalClusterEnergyDeposition> EcalClusterToolsT<noZS>::getEnergyDepTopology( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW, float w0 )
610 {
611  std::vector<typename EcalClusterToolsT<noZS>::EcalClusterEnergyDeposition> energyDistribution;
612  // init a map of the energy deposition centered on the
613  // cluster centroid. This is for momenta calculation only.
614  CLHEP::Hep3Vector clVect(cluster.position().x(), cluster.position().y(), cluster.position().z());
615  CLHEP::Hep3Vector clDir(clVect);
616  clDir*=1.0/clDir.mag();
617  // in the transverse plane, axis perpendicular to clusterDir
618  CLHEP::Hep3Vector theta_axis(clDir.y(),-clDir.x(),0.0);
619  theta_axis *= 1.0/theta_axis.mag();
620  CLHEP::Hep3Vector phi_axis = theta_axis.cross(clDir);
621 
622  const std::vector< std::pair<DetId, float> >& clusterDetIds = cluster.hitsAndFractions();
623 
624  EcalClusterEnergyDeposition clEdep;
625  EcalRecHit testEcalRecHit;
626  std::vector< std::pair<DetId, float> >::const_iterator posCurrent;
627  // loop over crystals
628  for(posCurrent=clusterDetIds.begin(); posCurrent!=clusterDetIds.end(); ++posCurrent) {
629  EcalRecHitCollection::const_iterator itt = recHits->find( (*posCurrent).first );
630  testEcalRecHit=*itt;
631 
632  if(( (*posCurrent).first != DetId(0)) && (recHits->find( (*posCurrent).first ) != recHits->end())) {
633  clEdep.deposited_energy = testEcalRecHit.energy() * (noZS ? 1.0 : (*posCurrent).second);
634  // if logarithmic weight is requested, apply cut on minimum energy of the recHit
635  if(logW) {
636  //double w0 = parameterMap_.find("W0")->second;
637 
638  double weight = std::max(0.0, w0 + log(std::abs(clEdep.deposited_energy)/cluster.energy()) );
639  if(weight==0) {
640  LogDebug("ClusterShapeAlgo") << "Crystal has insufficient energy: E = "
641  << clEdep.deposited_energy << " GeV; skipping... ";
642  continue;
643  }
644  else LogDebug("ClusterShapeAlgo") << "===> got crystal. Energy = " << clEdep.deposited_energy << " GeV. ";
645  }
646  DetId id_ = (*posCurrent).first;
647  const CaloCellGeometry *this_cell = geometry->getSubdetectorGeometry(id_)->getGeometry(id_);
648  GlobalPoint cellPos = this_cell->getPosition();
649  CLHEP::Hep3Vector gblPos (cellPos.x(),cellPos.y(),cellPos.z()); //surface position?
650  // Evaluate the distance from the cluster centroid
651  CLHEP::Hep3Vector diff = gblPos - clVect;
652  // Important: for the moment calculation, only the "lateral distance" is important
653  // "lateral distance" r_i = distance of the digi position from the axis Origin-Cluster Center
654  // ---> subtract the projection on clDir
655  CLHEP::Hep3Vector DigiVect = diff - diff.dot(clDir)*clDir;
656  clEdep.r = DigiVect.mag();
657  LogDebug("ClusterShapeAlgo") << "E = " << clEdep.deposited_energy
658  << "\tdiff = " << diff.mag()
659  << "\tr = " << clEdep.r;
660  clEdep.phi = DigiVect.angle(theta_axis);
661  if(DigiVect.dot(phi_axis)<0) clEdep.phi = 2 * M_PI - clEdep.phi;
662  energyDistribution.push_back(clEdep);
663  }
664  }
665  return energyDistribution;
666 }
667 
668 template<bool noZS>
669 std::vector<float> EcalClusterToolsT<noZS>::lat( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW, float w0 )
670 {
671  std::vector<EcalClusterToolsT::EcalClusterEnergyDeposition> energyDistribution = getEnergyDepTopology( cluster, recHits, geometry, logW, w0 );
672 
673  std::vector<float> lat;
674  double r, redmoment=0;
675  double phiRedmoment = 0 ;
676  double etaRedmoment = 0 ;
677  int n,n1,n2,tmp;
678  int clusterSize=energyDistribution.size();
679  float etaLat_, phiLat_, lat_;
680  if (clusterSize<3) {
681  etaLat_ = 0.0 ;
682  lat_ = 0.0;
683  lat.push_back(0.);
684  lat.push_back(0.);
685  lat.push_back(0.);
686  return lat;
687  }
688 
689  n1=0; n2=1;
690  if (energyDistribution[1].deposited_energy >
691  energyDistribution[0].deposited_energy)
692  {
693  tmp=n2; n2=n1; n1=tmp;
694  }
695  for (int i=2; i<clusterSize; i++) {
696  n=i;
697  if (energyDistribution[i].deposited_energy >
698  energyDistribution[n1].deposited_energy)
699  {
700  tmp = n2;
701  n2 = n1; n1 = i; n=tmp;
702  } else {
703  if (energyDistribution[i].deposited_energy >
704  energyDistribution[n2].deposited_energy)
705  {
706  tmp=n2; n2=i; n=tmp;
707  }
708  }
709 
710  r = energyDistribution[n].r;
711  redmoment += r*r* energyDistribution[n].deposited_energy;
712  double rphi = r * cos (energyDistribution[n].phi) ;
713  phiRedmoment += rphi * rphi * energyDistribution[n].deposited_energy;
714  double reta = r * sin (energyDistribution[n].phi) ;
715  etaRedmoment += reta * reta * energyDistribution[n].deposited_energy;
716  }
717  double e1 = energyDistribution[n1].deposited_energy;
718  double e2 = energyDistribution[n2].deposited_energy;
719 
720  lat_ = redmoment/(redmoment+2.19*2.19*(e1+e2));
721  phiLat_ = phiRedmoment/(phiRedmoment+2.19*2.19*(e1+e2));
722  etaLat_ = etaRedmoment/(etaRedmoment+2.19*2.19*(e1+e2));
723 
724  lat.push_back(etaLat_);
725  lat.push_back(phiLat_);
726  lat.push_back(lat_);
727  return lat;
728 }
729 
730 template<bool noZS>
731 math::XYZVector EcalClusterToolsT<noZS>::meanClusterPosition( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, const CaloGeometry *geometry )
732 {
733  // find mean energy position of a 5x5 cluster around the maximum
734  math::XYZVector meanPosition(0.0, 0.0, 0.0);
735  const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
736  std::vector<DetId> v_id = matrixDetId( topology, getMaximum( cluster, recHits ).first, -2, 2, -2, 2 );
737  for( const std::pair<DetId,float>& hitAndFrac : hsAndFs ) {
738  for( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
739  if( hitAndFrac.first != *it && !noZS) continue;
740  GlobalPoint positionGP = geometry->getSubdetectorGeometry( *it )->getGeometry( *it )->getPosition();
741  math::XYZVector position(positionGP.x(),positionGP.y(),positionGP.z());
742  meanPosition = meanPosition + recHitEnergy( *it, recHits ) * position * hitAndFrac.second;
743  }
744  if(noZS) break;
745  }
746  return meanPosition / e5x5( cluster, recHits, topology );
747 }
748 
749 //returns mean energy weighted eta/phi in crystals from the seed
750 //iPhi is not defined for endcap and is returned as zero
751 //return <eta,phi>
752 //we have an issue in working out what to do for negative energies
753 //I (Sam Harper) think it makes sense to ignore crystals with E<0 in the calculation as they are ignored
754 //in the sigmaIEtaIEta calculation (well they arent fully ignored, they do still contribute to the e5x5 sum
755 //in the sigmaIEtaIEta calculation but not here)
756 template<bool noZS>
757 std::pair<float,float> EcalClusterToolsT<noZS>::mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology)
758 {
759  DetId seedId = getMaximum( cluster, recHits ).first;
760  float meanDEta=0.;
761  float meanDPhi=0.;
762  float energySum=0.;
763 
764  const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
765  std::vector<DetId> v_id = matrixDetId( topology,seedId, -2, 2, -2, 2 );
766  for( const std::pair<DetId,float>& hAndF : hsAndFs ) {
767  for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
768  if( hAndF.first != *it && !noZS ) continue;
769  float energy = recHitEnergy(*it,recHits) * hAndF.second;
770  if(energy<0.) continue;//skipping negative energy crystals
771  meanDEta += energy * getNrCrysDiffInEta(*it,seedId);
772  meanDPhi += energy * getNrCrysDiffInPhi(*it,seedId);
773  energySum +=energy;
774  }
775  if(noZS) break;
776  }
777  meanDEta /=energySum;
778  meanDPhi /=energySum;
779  return std::pair<float,float>(meanDEta,meanDPhi);
780 }
781 
782 //returns mean energy weighted x/y in normalised crystal coordinates
783 //only valid for endcap, returns 0,0 for barrel
784 //we have an issue in working out what to do for negative energies
785 //I (Sam Harper) think it makes sense to ignore crystals with E<0 in the calculation as they are ignored
786 //in the sigmaIEtaIEta calculation (well they arent fully ignored, they do still contribute to the e5x5 sum
787 //in the sigmaIEtaIEta calculation but not here)
788 template<bool noZS>
789 std::pair<float,float> EcalClusterToolsT<noZS>::mean5x5PositionInXY(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology)
790 {
791  DetId seedId = getMaximum( cluster, recHits ).first;
792 
793  std::pair<float,float> meanXY(0.,0.);
794  if(seedId.subdetId()==EcalBarrel) return meanXY;
795 
796  float energySum=0.;
797 
798  const std::vector<std::pair<DetId,float> >& hsAndFs = cluster.hitsAndFractions();
799  std::vector<DetId> v_id = matrixDetId( topology,seedId, -2, 2, -2, 2 );
800  for( const std::pair<DetId,float>& hAndF : hsAndFs ) {
801  for ( std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it ) {
802  if( hAndF.first != *it && !noZS) continue;
803  float energy = recHitEnergy(*it,recHits) * hAndF.second;
804  if(energy<0.) continue;//skipping negative energy crystals
805  meanXY.first += energy * getNormedIX(*it);
806  meanXY.second += energy * getNormedIY(*it);
807  energySum +=energy;
808  }
809  if(noZS) break;
810  }
811  meanXY.first/=energySum;
812  meanXY.second/=energySum;
813  return meanXY;
814 }
815 
816 template<bool noZS>
817 std::vector<float> EcalClusterToolsT<noZS>::covariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits, const CaloTopology *topology, const CaloGeometry* geometry, float w0)
818 {
819  float e_5x5 = e5x5( cluster, recHits, topology );
820  float covEtaEta, covEtaPhi, covPhiPhi;
821  if (e_5x5 >= 0.) {
822  //double w0_ = parameterMap_.find("W0")->second;
823  const std::vector< std::pair<DetId, float>>& v_id =cluster.hitsAndFractions();
824  math::XYZVector meanPosition = meanClusterPosition( cluster, recHits, topology, geometry );
825 
826  // now we can calculate the covariances
827  double numeratorEtaEta = 0;
828  double numeratorEtaPhi = 0;
829  double numeratorPhiPhi = 0;
830  double denominator = 0;
831 
832  DetId id = getMaximum( v_id, recHits ).first;
833  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( id, topology->getSubdetectorTopology( id ) );
834  for ( int i = -2; i <= 2; ++i ) {
835  for ( int j = -2; j <= 2; ++j ) {
836  cursor.home();
837  cursor.offsetBy( i, j );
838  float frac=getFraction(v_id,*cursor);
839  float energy = recHitEnergy( *cursor, recHits )*frac;
840 
841  if ( energy <= 0 ) continue;
842 
843  GlobalPoint position = geometry->getSubdetectorGeometry(*cursor)->getGeometry(*cursor)->getPosition();
844 
845  double dPhi = position.phi() - meanPosition.phi();
846  if (dPhi > + Geom::pi()) { dPhi = Geom::twoPi() - dPhi; }
847  if (dPhi < - Geom::pi()) { dPhi = Geom::twoPi() + dPhi; }
848 
849  double dEta = position.eta() - meanPosition.eta();
850  double w = 0.;
851  w = std::max(0.0f, w0 + std::log( energy / e_5x5 ));
852 
853  denominator += w;
854  numeratorEtaEta += w * dEta * dEta;
855  numeratorEtaPhi += w * dEta * dPhi;
856  numeratorPhiPhi += w * dPhi * dPhi;
857  }
858  }
859 
860  if (denominator != 0.0) {
861  covEtaEta = numeratorEtaEta / denominator;
862  covEtaPhi = numeratorEtaPhi / denominator;
863  covPhiPhi = numeratorPhiPhi / denominator;
864  } else {
865  covEtaEta = 999.9;
866  covEtaPhi = 999.9;
867  covPhiPhi = 999.9;
868  }
869 
870  } else {
871  // Warn the user if there was no energy in the cells and return zeroes.
872  // std::cout << "\ClusterShapeAlgo::Calculate_Covariances: no energy in supplied cells.\n";
873  covEtaEta = 0;
874  covEtaPhi = 0;
875  covPhiPhi = 0;
876  }
877  std::vector<float> v;
878  v.push_back( covEtaEta );
879  v.push_back( covEtaPhi );
880  v.push_back( covPhiPhi );
881  return v;
882 }
883 
884 //for covIEtaIEta,covIEtaIPhi and covIPhiIPhi are defined but only covIEtaIEta has been actively studied
885 //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
886 //it also does not require any eta correction function in the endcap
887 //it is multipled by an approprate crystal size to ensure it gives similar values to covariances(...)
888 template<bool noZS>
889 std::vector<float> EcalClusterToolsT<noZS>::localCovariances(const reco::BasicCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology,float w0)
890 {
891 
892  float e_5x5 = e5x5( cluster, recHits, topology );
893  float covEtaEta, covEtaPhi, covPhiPhi;
894 
895  if (e_5x5 >= 0.) {
896  //double w0_ = parameterMap_.find("W0")->second;
897  const std::vector< std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
898  std::pair<float,float> mean5x5PosInNrCrysFromSeed = mean5x5PositionInLocalCrysCoord( cluster, recHits, topology );
899  std::pair<float,float> mean5x5XYPos = mean5x5PositionInXY(cluster,recHits,topology);
900 
901  // now we can calculate the covariances
902  double numeratorEtaEta = 0;
903  double numeratorEtaPhi = 0;
904  double numeratorPhiPhi = 0;
905  double denominator = 0;
906 
907  //these allow us to scale the localCov by the crystal size
908  //so that the localCovs have the same average value as the normal covs
909  const double barrelCrysSize = 0.01745; //approximate size of crystal in eta,phi in barrel
910  const double endcapCrysSize = 0.0447; //the approximate crystal size sigmaEtaEta was corrected to in the endcap
911 
912  DetId seedId = getMaximum( v_id, recHits ).first;
913 
914  bool isBarrel=seedId.subdetId()==EcalBarrel;
915  const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
916 
917  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( seedId, topology->getSubdetectorTopology( seedId ) );
918 
919  for ( int eastNr = -2; eastNr <= 2; ++eastNr ) { //east is eta in barrel
920  for ( int northNr = -2; northNr <= 2; ++northNr ) { //north is phi in barrel
921  cursor.home();
922  cursor.offsetBy( eastNr, northNr);
923  float frac = getFraction(v_id,*cursor);
924  float energy = recHitEnergy( *cursor, recHits )*frac;
925  if ( energy <= 0 ) continue;
926 
927  float dEta = getNrCrysDiffInEta(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.first;
928  float dPhi = 0;
929 
930  if(isBarrel) dPhi = getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
931  else dPhi = getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
932 
933 
934  double w = std::max(0.0f,w0 + std::log( energy / e_5x5 ));
935 
936  denominator += w;
937  numeratorEtaEta += w * dEta * dEta;
938  numeratorEtaPhi += w * dEta * dPhi;
939  numeratorPhiPhi += w * dPhi * dPhi;
940  } //end east loop
941  }//end north loop
942 
943 
944  //multiplying by crysSize to make the values compariable to normal covariances
945  if (denominator != 0.0) {
946  covEtaEta = crysSize*crysSize* numeratorEtaEta / denominator;
947  covEtaPhi = crysSize*crysSize* numeratorEtaPhi / denominator;
948  covPhiPhi = crysSize*crysSize* numeratorPhiPhi / denominator;
949  } else {
950  covEtaEta = 999.9;
951  covEtaPhi = 999.9;
952  covPhiPhi = 999.9;
953  }
954 
955 
956  } else {
957  // Warn the user if there was no energy in the cells and return zeroes.
958  // std::cout << "\ClusterShapeAlgo::Calculate_Covariances: no energy in supplied cells.\n";
959  covEtaEta = 0;
960  covEtaPhi = 0;
961  covPhiPhi = 0;
962  }
963  std::vector<float> v;
964  v.push_back( covEtaEta );
965  v.push_back( covEtaPhi );
966  v.push_back( covPhiPhi );
967  return v;
968 }
969 
970 template<bool noZS>
971 double EcalClusterToolsT<noZS>::zernike20( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, double R0, bool logW, float w0 )
972 {
973  return absZernikeMoment( cluster, recHits, geometry, 2, 0, R0, logW, w0 );
974 }
975 
976 template<bool noZS>
977 double EcalClusterToolsT<noZS>::zernike42( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, double R0, bool logW, float w0 )
978 {
979  return absZernikeMoment( cluster, recHits, geometry, 4, 2, R0, logW, w0 );
980 }
981 
982 template<bool noZS>
983 double EcalClusterToolsT<noZS>::absZernikeMoment( const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, int n, int m, double R0, bool logW, float w0 )
984 {
985  // 1. Check if n,m are correctly
986  if ((m>n) || ((n-m)%2 != 0) || (n<0) || (m<0)) return -1;
987 
988  // 2. Check if n,R0 are within validity Range :
989  // n>20 or R0<2.19cm just makes no sense !
990  if ((n>20) || (R0<=2.19)) return -1;
991  if (n<=5) return fast_AbsZernikeMoment(cluster, recHits, geometry, n, m, R0, logW, w0 );
992  else return calc_AbsZernikeMoment(cluster, recHits, geometry, n, m, R0, logW, w0 );
993 }
994 
995 template<bool noZS>
996 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 )
997 {
998  double r,ph,e,Re=0,Im=0;
999  double TotalEnergy = cluster.energy();
1000  int index = (n/2)*(n/2)+(n/2)+m;
1001  std::vector<EcalClusterEnergyDeposition> energyDistribution = getEnergyDepTopology( cluster, recHits, geometry, logW, w0 );
1002  int clusterSize = energyDistribution.size();
1003  if(clusterSize < 3) return 0.0;
1004 
1005  for (int i=0; i<clusterSize; i++)
1006  {
1007  r = energyDistribution[i].r / R0;
1008  if (r<1) {
1009  std::vector<double> pol;
1010  pol.push_back( f00(r) );
1011  pol.push_back( f11(r) );
1012  pol.push_back( f20(r) );
1013  pol.push_back( f22(r) );
1014  pol.push_back( f31(r) );
1015  pol.push_back( f33(r) );
1016  pol.push_back( f40(r) );
1017  pol.push_back( f42(r) );
1018  pol.push_back( f44(r) );
1019  pol.push_back( f51(r) );
1020  pol.push_back( f53(r) );
1021  pol.push_back( f55(r) );
1022  ph = (energyDistribution[i]).phi;
1023  e = energyDistribution[i].deposited_energy;
1024  Re = Re + e/TotalEnergy * pol[index] * cos( (double) m * ph);
1025  Im = Im - e/TotalEnergy * pol[index] * sin( (double) m * ph);
1026  }
1027  }
1028  return sqrt(Re*Re+Im*Im);
1029 }
1030 
1031 template<bool noZS>
1032 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 )
1033 {
1034  double r, ph, e, Re=0, Im=0, f_nm;
1035  double TotalEnergy = cluster.energy();
1036  std::vector<EcalClusterEnergyDeposition> energyDistribution = getEnergyDepTopology( cluster, recHits, geometry, logW, w0 );
1037  int clusterSize=energyDistribution.size();
1038  if(clusterSize<3) return 0.0;
1039 
1040  for (int i = 0; i < clusterSize; ++i)
1041  {
1042  r = energyDistribution[i].r / R0;
1043  if (r < 1) {
1044  ph = energyDistribution[i].phi;
1045  e = energyDistribution[i].deposited_energy;
1046  f_nm = 0;
1047  for (int s=0; s<=(n-m)/2; s++) {
1048  if (s%2==0) {
1049  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));
1050  } else {
1051  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));
1052  }
1053  }
1054  Re = Re + e/TotalEnergy * f_nm * cos( (double) m*ph);
1055  Im = Im - e/TotalEnergy * f_nm * sin( (double) m*ph);
1056  }
1057  }
1058  return sqrt(Re*Re+Im*Im);
1059 }
1060 
1061 //returns the crystal 'eta' from the det id
1062 //it is defined as the number of crystals from the centre in the eta direction
1063 //for the barrel with its eta/phi geometry it is always integer
1064 //for the endcap it is fractional due to the x/y geometry
1065 template<bool noZS>
1066 float EcalClusterToolsT<noZS>::getIEta(const DetId& id)
1067 {
1068  if(id.det()==DetId::Ecal){
1069  if(id.subdetId()==EcalBarrel){
1070  EBDetId ebId(id);
1071  return ebId.ieta();
1072  }else if(id.subdetId()==EcalEndcap){
1073  float iXNorm = getNormedIX(id);
1074  float iYNorm = getNormedIY(id);
1075 
1076  return std::sqrt(iXNorm*iXNorm+iYNorm*iYNorm);
1077  }
1078  }
1079  return 0.;
1080 }
1081 
1082 
1083 //returns the crystal 'phi' from the det id
1084 //it is defined as the number of crystals from the centre in the phi direction
1085 //for the barrel with its eta/phi geometry it is always integer
1086 //for the endcap it is not defined
1087 template<bool noZS>
1088 float EcalClusterToolsT<noZS>::getIPhi(const DetId& id)
1089 {
1090  if(id.det()==DetId::Ecal){
1091  if(id.subdetId()==EcalBarrel){
1092  EBDetId ebId(id);
1093  return ebId.iphi();
1094  }
1095  }
1096  return 0.;
1097 }
1098 
1099 //want to map 1=-50,50=-1,51=1 and 100 to 50 so sub off one if zero or neg
1100 template<bool noZS>
1101 float EcalClusterToolsT<noZS>::getNormedIX(const DetId& id)
1102 {
1103  if(id.det()==DetId::Ecal && id.subdetId()==EcalEndcap){
1104  EEDetId eeId(id);
1105  int iXNorm = eeId.ix()-50;
1106  if(iXNorm<=0) iXNorm--;
1107  return iXNorm;
1108  }
1109  return 0;
1110 }
1111 
1112 //want to map 1=-50,50=-1,51=1 and 100 to 50 so sub off one if zero or neg
1113 template<bool noZS>
1114 float EcalClusterToolsT<noZS>::getNormedIY(const DetId& id)
1115 {
1116  if(id.det()==DetId::Ecal && id.subdetId()==EcalEndcap){
1117  EEDetId eeId(id);
1118  int iYNorm = eeId.iy()-50;
1119  if(iYNorm<=0) iYNorm--;
1120  return iYNorm;
1121  }
1122  return 0;
1123 }
1124 
1125 //nr crystals crysId is away from orgin id in eta
1126 template<bool noZS>
1127 float EcalClusterToolsT<noZS>::getNrCrysDiffInEta(const DetId& crysId,const DetId& orginId)
1128 {
1129  float crysIEta = getIEta(crysId);
1130  float orginIEta = getIEta(orginId);
1131  bool isBarrel = orginId.subdetId()==EcalBarrel;
1132 
1133  float nrCrysDiff = crysIEta-orginIEta;
1134 
1135  //no iEta=0 in barrel, so if go from positive to negative
1136  //need to reduce abs(detEta) by 1
1137  if(isBarrel){
1138  if(crysIEta*orginIEta<0){ // -1 to 1 transition
1139  if(crysIEta>0) nrCrysDiff--;
1140  else nrCrysDiff++;
1141  }
1142  }
1143  return nrCrysDiff;
1144 }
1145 
1146 //nr crystals crysId is away from orgin id in phi
1147 template<bool noZS>
1148 float EcalClusterToolsT<noZS>::getNrCrysDiffInPhi(const DetId& crysId,const DetId& orginId)
1149 {
1150  float crysIPhi = getIPhi(crysId);
1151  float orginIPhi = getIPhi(orginId);
1152  bool isBarrel = orginId.subdetId()==EcalBarrel;
1153 
1154  float nrCrysDiff = crysIPhi-orginIPhi;
1155 
1156  if(isBarrel){ //if barrel, need to map into 0-180
1157  if (nrCrysDiff > + 180) { nrCrysDiff = nrCrysDiff - 360; }
1158  if (nrCrysDiff < - 180) { nrCrysDiff = nrCrysDiff + 360; }
1159  }
1160  return nrCrysDiff;
1161 }
1162 
1163 //nr crystals crysId is away from mean phi in 5x5 in phi
1164 template<bool noZS>
1165 float EcalClusterToolsT<noZS>::getDPhiEndcap(const DetId& crysId,float meanX,float meanY)
1166 {
1167  float iXNorm = getNormedIX(crysId);
1168  float iYNorm = getNormedIY(crysId);
1169 
1170  float hitLocalR2 = (iXNorm-meanX)*(iXNorm-meanX)+(iYNorm-meanY)*(iYNorm-meanY);
1171  float hitR2 = iXNorm*iXNorm+iYNorm*iYNorm;
1172  float meanR2 = meanX*meanX+meanY*meanY;
1173  float hitR = sqrt(hitR2);
1174  float meanR = sqrt(meanR2);
1175 
1176  float tmp = (hitR2+meanR2-hitLocalR2)/(2*hitR*meanR);
1177  if (tmp<-1) tmp =-1;
1178  if (tmp>1) tmp=1;
1179  float phi = acos(tmp);
1180  float dPhi = hitR*phi;
1181 
1182  return dPhi;
1183 }
1184 
1185 template<bool noZS>
1186 std::vector<float> EcalClusterToolsT<noZS>::scLocalCovariances(const reco::SuperCluster &cluster, const EcalRecHitCollection* recHits,const CaloTopology *topology, float w0)
1187 {
1188  const reco::BasicCluster bcluster = *(cluster.seed());
1189 
1190  float e_5x5 = e5x5(bcluster, recHits, topology);
1191  float covEtaEta, covEtaPhi, covPhiPhi;
1192 
1193  if (e_5x5 >= 0.) {
1194  const std::vector<std::pair<DetId, float> >& v_id = cluster.hitsAndFractions();
1195  std::pair<float,float> mean5x5PosInNrCrysFromSeed = mean5x5PositionInLocalCrysCoord(bcluster, recHits, topology);
1196  std::pair<float,float> mean5x5XYPos = mean5x5PositionInXY(cluster,recHits,topology);
1197  // now we can calculate the covariances
1198  double numeratorEtaEta = 0;
1199  double numeratorEtaPhi = 0;
1200  double numeratorPhiPhi = 0;
1201  double denominator = 0;
1202 
1203  const double barrelCrysSize = 0.01745; //approximate size of crystal in eta,phi in barrel
1204  const double endcapCrysSize = 0.0447; //the approximate crystal size sigmaEtaEta was corrected to in the endcap
1205 
1206  DetId seedId = getMaximum(v_id, recHits).first;
1207  bool isBarrel=seedId.subdetId()==EcalBarrel;
1208 
1209  const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1210 
1211  for (size_t i = 0; i < v_id.size(); ++i) {
1212  CaloNavigator<DetId> cursor = CaloNavigator<DetId>(v_id[i].first, topology->getSubdetectorTopology(v_id[i].first));
1213  float frac = getFraction(v_id,*cursor);
1214  float energy = recHitEnergy(*cursor, recHits)*frac;
1215 
1216  if (energy <= 0) continue;
1217 
1218  float dEta = getNrCrysDiffInEta(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.first;
1219  float dPhi = 0;
1220  if(isBarrel) dPhi = getNrCrysDiffInPhi(*cursor,seedId) - mean5x5PosInNrCrysFromSeed.second;
1221  else dPhi = getDPhiEndcap(*cursor,mean5x5XYPos.first,mean5x5XYPos.second);
1222 
1223 
1224 
1225  double w = 0.;
1226  w = std::max(0.0f, w0 + std::log( energy / e_5x5 ));
1227 
1228  denominator += w;
1229  numeratorEtaEta += w * dEta * dEta;
1230  numeratorEtaPhi += w * dEta * dPhi;
1231  numeratorPhiPhi += w * dPhi * dPhi;
1232  }
1233 
1234  //multiplying by crysSize to make the values compariable to normal covariances
1235  if (denominator != 0.0) {
1236  covEtaEta = crysSize*crysSize* numeratorEtaEta / denominator;
1237  covEtaPhi = crysSize*crysSize* numeratorEtaPhi / denominator;
1238  covPhiPhi = crysSize*crysSize* numeratorPhiPhi / denominator;
1239  } else {
1240  covEtaEta = 999.9;
1241  covEtaPhi = 999.9;
1242  covPhiPhi = 999.9;
1243  }
1244 
1245  } else {
1246  // Warn the user if there was no energy in the cells and return zeroes.
1247  // std::cout << "\ClusterShapeAlgo::Calculate_Covariances: no energy in supplied cells.\n";
1248  covEtaEta = 0;
1249  covEtaPhi = 0;
1250  covPhiPhi = 0;
1251  }
1252 
1253  std::vector<float> v;
1254  v.push_back( covEtaEta );
1255  v.push_back( covEtaPhi );
1256  v.push_back( covPhiPhi );
1257 
1258  return v;
1259 }
1260 
1261 
1262 // compute cluster second moments with respect to principal axes (eigenvectors of sEtaEta, sPhiPhi, sEtaPhi matrix)
1263 // store also angle alpha between major axis and phi.
1264 // takes into account shower elongation in phi direction due to magnetic field effect:
1265 // default value of 0.8 ensures sMaj = sMin for unconverted photons
1266 // (if phiCorrectionFactor=1 sMaj > sMin and alpha=0 also for unconverted photons)
1267 template<bool noZS>
1268 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments( const reco::BasicCluster &basicCluster, const EcalRecHitCollection &recHits, double phiCorrectionFactor, double w0, bool useLogWeights) {
1269 
1270  // for now implemented only for EB:
1271  // if( fabs( basicCluster.eta() ) < 1.479 ) {
1272 
1273  std::vector<std::pair<const EcalRecHit*, float> > RH_ptrs_fracs;
1274 
1275  const std::vector< std::pair<DetId, float> >& myHitsPair = basicCluster.hitsAndFractions();
1276 
1277  for(unsigned int i=0; i<myHitsPair.size(); i++){
1278  //get pointer to recHit object
1279  EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
1280  RH_ptrs_fracs.push_back( std::make_pair(&(*myRH) , myHitsPair[i].second) );
1281  }
1282 
1283  return EcalClusterToolsT<noZS>::cluster2ndMoments(RH_ptrs_fracs, phiCorrectionFactor, w0, useLogWeights);
1284 }
1285 
1286 template<bool noZS>
1287 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments( const reco::SuperCluster &superCluster, const EcalRecHitCollection &recHits, double phiCorrectionFactor, double w0, bool useLogWeights) {
1288 
1289  // for now returns second moments of supercluster seed cluster:
1290  Cluster2ndMoments returnMoments;
1291  returnMoments.sMaj = -1.;
1292  returnMoments.sMin = -1.;
1293  returnMoments.alpha = 0.;
1294 
1295  // for now implemented only for EB:
1296  // if( fabs( superCluster.eta() ) < 1.479 ) {
1297  returnMoments = EcalClusterToolsT<noZS>::cluster2ndMoments( *(superCluster.seed()), recHits, phiCorrectionFactor, w0, useLogWeights);
1298  // }
1299 
1300  return returnMoments;
1301 
1302 }
1303 
1304 template<bool noZS>
1305 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments( const std::vector<std::pair<const EcalRecHit*, float> >& RH_ptrs_fracs, double phiCorrectionFactor, double w0, bool useLogWeights) {
1306 
1307  double mid_eta(0),mid_phi(0),mid_x(0),mid_y(0);
1308 
1309  double Etot = EcalClusterToolsT<noZS>::getSumEnergy( RH_ptrs_fracs );
1310 
1311  double max_phi=-10.;
1312  double min_phi=100.;
1313 
1314 
1315  std::vector<double> etaDetId;
1316  std::vector<double> phiDetId;
1317  std::vector<double> xDetId;
1318  std::vector<double> yDetId;
1319  std::vector<double> wiDetId;
1320 
1321  unsigned int nCry=0;
1322  double denominator=0.;
1323  bool isBarrel(1);
1324 
1325  // loop over rechits and compute weights:
1326  for(std::vector<std::pair<const EcalRecHit*, float> >::const_iterator rhf_ptr = RH_ptrs_fracs.begin(); rhf_ptr != RH_ptrs_fracs.end(); rhf_ptr++){
1327  const EcalRecHit* rh_ptr = rhf_ptr->first;
1328 
1329 
1330  //get iEta, iPhi
1331  double temp_eta(0),temp_phi(0),temp_x(0),temp_y(0);
1332  isBarrel = rh_ptr->detid().subdetId()==EcalBarrel;
1333 
1334  if(isBarrel) {
1335  temp_eta = (getIEta(rh_ptr->detid()) > 0. ? getIEta(rh_ptr->detid()) + 84.5 : getIEta(rh_ptr->detid()) + 85.5);
1336  temp_phi= getIPhi(rh_ptr->detid()) - 0.5;
1337  }
1338  else {
1339  temp_eta = getIEta(rh_ptr->detid());
1340  temp_x = getNormedIX(rh_ptr->detid());
1341  temp_y = getNormedIY(rh_ptr->detid());
1342  }
1343 
1344  double temp_ene=rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
1345 
1346  double temp_wi=((useLogWeights) ?
1347  std::max(0.0, w0 + std::log( std::abs(temp_ene)/Etot ))
1348  : temp_ene);
1349 
1350 
1351  if(temp_phi>max_phi) max_phi=temp_phi;
1352  if(temp_phi<min_phi) min_phi=temp_phi;
1353  etaDetId.push_back(temp_eta);
1354  phiDetId.push_back(temp_phi);
1355  xDetId.push_back(temp_x);
1356  yDetId.push_back(temp_y);
1357  wiDetId.push_back(temp_wi);
1358  denominator+=temp_wi;
1359  nCry++;
1360  }
1361 
1362  if(isBarrel){
1363  // correct phi wrap-around:
1364  if(max_phi==359.5 && min_phi==0.5){
1365  for(unsigned int i=0; i<nCry; i++){
1366  if(phiDetId[i] - 179. > 0.) phiDetId[i]-=360.;
1367  mid_phi+=phiDetId[i]*wiDetId[i];
1368  mid_eta+=etaDetId[i]*wiDetId[i];
1369  }
1370  } else{
1371  for(unsigned int i=0; i<nCry; i++){
1372  mid_phi+=phiDetId[i]*wiDetId[i];
1373  mid_eta+=etaDetId[i]*wiDetId[i];
1374  }
1375  }
1376  }else{
1377  for(unsigned int i=0; i<nCry; i++){
1378  mid_eta+=etaDetId[i]*wiDetId[i];
1379  mid_x+=xDetId[i]*wiDetId[i];
1380  mid_y+=yDetId[i]*wiDetId[i];
1381  }
1382  }
1383 
1384  mid_eta/=denominator;
1385  mid_phi/=denominator;
1386  mid_x/=denominator;
1387  mid_y/=denominator;
1388 
1389 
1390  // See = sigma eta eta
1391  // Spp = (B field corrected) sigma phi phi
1392  // See = (B field corrected) sigma eta phi
1393  double See=0.;
1394  double Spp=0.;
1395  double Sep=0.;
1396  double deta(0),dphi(0);
1397  // compute (phi-corrected) covariance matrix:
1398  for(unsigned int i=0; i<nCry; i++) {
1399  if(isBarrel) {
1400  deta = etaDetId[i]-mid_eta;
1401  dphi = phiDetId[i]-mid_phi;
1402  } else {
1403  deta = etaDetId[i]-mid_eta;
1404  float hitLocalR2 = (xDetId[i]-mid_x)*(xDetId[i]-mid_x)+(yDetId[i]-mid_y)*(yDetId[i]-mid_y);
1405  float hitR2 = xDetId[i]*xDetId[i]+yDetId[i]*yDetId[i];
1406  float meanR2 = mid_x*mid_x+mid_y*mid_y;
1407  float hitR = sqrt(hitR2);
1408  float meanR = sqrt(meanR2);
1409  float phi = acos((hitR2+meanR2-hitLocalR2)/(2*hitR*meanR));
1410  dphi = hitR*phi;
1411 
1412  }
1413  See += (wiDetId[i]* deta * deta) / denominator;
1414  Spp += phiCorrectionFactor*(wiDetId[i]* dphi * dphi) / denominator;
1415  Sep += sqrt(phiCorrectionFactor)*(wiDetId[i]*deta*dphi) / denominator;
1416  }
1417 
1418  Cluster2ndMoments returnMoments;
1419 
1420  // compute matrix eigenvalues:
1421  returnMoments.sMaj = ((See + Spp) + sqrt((See - Spp)*(See - Spp) + 4.*Sep*Sep)) / 2.;
1422  returnMoments.sMin = ((See + Spp) - sqrt((See - Spp)*(See - Spp) + 4.*Sep*Sep)) / 2.;
1423 
1424  returnMoments.alpha = atan( (See - Spp + sqrt( (Spp - See)*(Spp - See) + 4.*Sep*Sep )) / (2.*Sep));
1425 
1426  return returnMoments;
1427 
1428 }
1429 
1430 //compute shower shapes: roundness and angle in a vector. Roundness is 0th element, Angle is 1st element.
1431 //description: uses classical mechanics inertia tensor.
1432 // roundness is smaller_eValue/larger_eValue
1433 // angle is the angle from the iEta axis to the smallest eVector (a.k.a. the shower's elongated axis)
1434 // this function uses only recHits belonging to a SC above energyThreshold (default 0)
1435 // you can select linear weighting = energy_recHit/total_energy (weightedPositionMethod=0) default
1436 // or log weighting = max( 0.0, 4.2 + log(energy_recHit/total_energy) ) (weightedPositionMethod=1)
1437 template<bool noZS>
1438 std::vector<float> EcalClusterToolsT<noZS>::roundnessBarrelSuperClusters( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int weightedPositionMethod, float energyThreshold){//int positionWeightingMethod=0){
1439  std::vector<std::pair<const EcalRecHit*, float> > RH_ptrs_fracs;
1440  const std::vector< std::pair<DetId, float> >& myHitsPair = superCluster.hitsAndFractions();
1441  for(unsigned int i=0; i< myHitsPair.size(); ++i){
1442  //get pointer to recHit object
1443  EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
1444  if( myRH != recHits.end() && myRH->energy()*(noZS ? 1.0 : myHitsPair[i].second) > energyThreshold){
1445  //require rec hit to have positive energy
1446  RH_ptrs_fracs.push_back( std::make_pair(&(*myRH) , myHitsPair[i].second) );
1447  }
1448  }
1449  std::vector<float> temp = EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(RH_ptrs_fracs,weightedPositionMethod);
1450  return temp;
1451 }
1452 
1453 // this function uses all recHits within specified window ( with default values ieta_delta=2, iphi_delta=5) around SC's highest recHit.
1454 // recHits must pass an energy threshold "energyRHThresh" (default 0)
1455 // you can select linear weighting = energy_recHit/total_energy (weightedPositionMethod=0)
1456 // or log weighting = max( 0.0, 4.2 + log(energy_recHit/total_energy) ) (weightedPositionMethod=1)
1457 template<bool noZS>
1458 std::vector<float> EcalClusterToolsT<noZS>::roundnessBarrelSuperClustersUserExtended( const reco::SuperCluster &superCluster ,const EcalRecHitCollection &recHits, int ieta_delta, int iphi_delta, float energyRHThresh, int weightedPositionMethod){
1459 
1460  std::vector<std::pair<const EcalRecHit*, float> > RH_ptrs_fracs;
1461  const std::vector< std::pair<DetId, float> >& myHitsPair = superCluster.hitsAndFractions();
1462  for(unsigned int i=0; i<myHitsPair.size(); ++i){
1463  //get pointer to recHit object
1464  EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
1465  if(myRH != recHits.end() && myRH->energy()*(noZS ? 1.0 : myHitsPair[i].second) > energyRHThresh)
1466  RH_ptrs_fracs.push_back( std::make_pair(&(*myRH) , myHitsPair[i].second) );
1467  }
1468 
1469 
1470  std::vector<int> seedPosition = EcalClusterToolsT<noZS>::getSeedPosition( RH_ptrs_fracs );
1471 
1472  for(EcalRecHitCollection::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++){
1473  EBDetId EBdetIdi( rh->detid() );
1474  float the_fraction = 0;
1475  //if(rh != recHits.end())
1476  bool inEtaWindow = ( abs( deltaIEta(seedPosition[0],EBdetIdi.ieta()) ) <= ieta_delta );
1477  bool inPhiWindow = ( abs( deltaIPhi(seedPosition[1],EBdetIdi.iphi()) ) <= iphi_delta );
1478  bool passEThresh = ( rh->energy() > energyRHThresh );
1479  bool alreadyCounted = false;
1480 
1481  // figure out if the rechit considered now is already inside the SC
1482  bool is_SCrh_inside_recHits = false;
1483  for(unsigned int i=0; i<myHitsPair.size(); i++){
1484  EcalRecHitCollection::const_iterator SCrh = recHits.find(myHitsPair[i].first);
1485  if(SCrh != recHits.end()){
1486  the_fraction = myHitsPair[i].second;
1487  is_SCrh_inside_recHits = true;
1488  if( rh->detid() == SCrh->detid() ) alreadyCounted = true;
1489  }
1490  }//for loop over SC's recHits
1491 
1492  if( is_SCrh_inside_recHits && !alreadyCounted && passEThresh && inEtaWindow && inPhiWindow){
1493  RH_ptrs_fracs.push_back( std::make_pair(&(*rh),the_fraction) );
1494  }
1495 
1496  }//for loop over rh
1497  return EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(RH_ptrs_fracs,weightedPositionMethod);
1498 }
1499 
1500 // this function computes the roundness and angle variables for vector of pointers to recHits you pass it
1501 // you can select linear weighting = energy_recHit/total_energy (weightedPositionMethod=0)
1502 // or log weighting = max( 0.0, 4.2 + log(energy_recHit/total_energy) ) (weightedPositionMethod=1)
1503 template<bool noZS>
1504 std::vector<float> EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits( const std::vector<std::pair<const EcalRecHit*,float> >& RH_ptrs_fracs, int weightedPositionMethod){//int weightedPositionMethod = 0){
1505  //positionWeightingMethod = 0 linear weighting, 1 log energy weighting
1506 
1507  std::vector<float> shapes; // this is the returning vector
1508 
1509  //make sure photon has more than one crystal; else roundness and angle suck
1510  if(RH_ptrs_fracs.size()<2){
1511  shapes.push_back( -3 );
1512  shapes.push_back( -3 );
1513  return shapes;
1514  }
1515 
1516  //Find highest E RH (Seed) and save info, compute sum total energy used
1517  std::vector<int> seedPosition = EcalClusterToolsT<noZS>::getSeedPosition( RH_ptrs_fracs );// *recHits);
1518  int tempInt = seedPosition[0];
1519  if(tempInt <0) tempInt++;
1520  float energyTotal = EcalClusterToolsT<noZS>::getSumEnergy( RH_ptrs_fracs );
1521 
1522  //1st loop over rechits: compute new weighted center position in coordinates centered on seed
1523  float centerIEta = 0.;
1524  float centerIPhi = 0.;
1525  float denominator = 0.;
1526 
1527  for(std::vector<std::pair<const EcalRecHit*,float> >::const_iterator rhf_ptr = RH_ptrs_fracs.begin(); rhf_ptr != RH_ptrs_fracs.end(); rhf_ptr++){
1528  const EcalRecHit* rh_ptr = rhf_ptr->first;
1529  //get iEta, iPhi
1530  EBDetId EBdetIdi( rh_ptr->detid() );
1531  if(fabs(energyTotal) < 0.0001){
1532  // don't /0, bad!
1533  shapes.push_back( -2 );
1534  shapes.push_back( -2 );
1535  return shapes;
1536  }
1537  float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
1538  float weight = 0;
1539  if(std::abs(weightedPositionMethod)<0.0001){ //linear
1540  weight = rh_energy/energyTotal;
1541  }else{ //logrithmic
1542  weight = std::max(0.0, 4.2 + log(rh_energy/energyTotal));
1543  }
1544  denominator += weight;
1545  centerIEta += weight*deltaIEta(seedPosition[0],EBdetIdi.ieta());
1546  centerIPhi += weight*deltaIPhi(seedPosition[1],EBdetIdi.iphi());
1547  }
1548  if(fabs(denominator) < 0.0001){
1549  // don't /0, bad!
1550  shapes.push_back( -2 );
1551  shapes.push_back( -2 );
1552  return shapes;
1553  }
1554  centerIEta = centerIEta / denominator;
1555  centerIPhi = centerIPhi / denominator;
1556 
1557 
1558  //2nd loop over rechits: compute inertia tensor
1559  TMatrixDSym inertia(2); //initialize 2d inertia tensor
1560  double inertia00 = 0.;
1561  double inertia01 = 0.;// = inertia10 b/c matrix should be symmetric
1562  double inertia11 = 0.;
1563  int i = 0;
1564  for(std::vector<std::pair<const EcalRecHit*,float> >::const_iterator rhf_ptr = RH_ptrs_fracs.begin(); rhf_ptr != RH_ptrs_fracs.end(); rhf_ptr++){
1565  const EcalRecHit* rh_ptr = rhf_ptr->first;
1566  //get iEta, iPhi
1567  EBDetId EBdetIdi( rh_ptr->detid() );
1568 
1569  if(fabs(energyTotal) < 0.0001){
1570  // don't /0, bad!
1571  shapes.push_back( -2 );
1572  shapes.push_back( -2 );
1573  return shapes;
1574  }
1575  float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
1576  float weight = 0;
1577  if(std::abs(weightedPositionMethod) < 0.0001){ //linear
1578  weight = rh_energy/energyTotal;
1579  }else{ //logrithmic
1580  weight = std::max(0.0, 4.2 + log(rh_energy/energyTotal));
1581  }
1582 
1583  float ieta_rh_to_center = deltaIEta(seedPosition[0],EBdetIdi.ieta()) - centerIEta;
1584  float iphi_rh_to_center = deltaIPhi(seedPosition[1],EBdetIdi.iphi()) - centerIPhi;
1585 
1586  inertia00 += weight*iphi_rh_to_center*iphi_rh_to_center;
1587  inertia01 -= weight*iphi_rh_to_center*ieta_rh_to_center;
1588  inertia11 += weight*ieta_rh_to_center*ieta_rh_to_center;
1589  i++;
1590  }
1591 
1592  inertia[0][0] = inertia00;
1593  inertia[0][1] = inertia01; // use same number here
1594  inertia[1][0] = inertia01; // and here to insure symmetry
1595  inertia[1][1] = inertia11;
1596 
1597 
1598  //step 1 find principal axes of inertia
1599  TMatrixD eVectors(2,2);
1600  TVectorD eValues(2);
1601  //std::cout<<"EcalClusterToolsT<noZS>::showerRoundness- about to compute eVectors"<<std::endl;
1602  eVectors=inertia.EigenVectors(eValues); //ordered highest eV to lowest eV (I checked!)
1603  //and eVectors are in columns of matrix! I checked!
1604  //and they are normalized to 1
1605 
1606 
1607 
1608  //step 2 select eta component of smaller eVal's eVector
1609  TVectorD smallerAxis(2);//easiest to spin SC on this axis (smallest eVal)
1610  smallerAxis[0]=eVectors[0][1];//row,col //eta component
1611  smallerAxis[1]=eVectors[1][1]; //phi component
1612 
1613  //step 3 compute interesting quatities
1614  Double_t temp = fabs(smallerAxis[0]);// closer to 1 ->beamhalo, closer to 0 something else
1615  if(fabs(eValues[0]) < 0.0001){
1616  // don't /0, bad!
1617  shapes.push_back( -2 );
1618  shapes.push_back( -2 );
1619  return shapes;
1620  }
1621 
1622  float Roundness = eValues[1]/eValues[0];
1623  float Angle=acos(temp);
1624 
1625  if( -0.00001 < Roundness && Roundness < 0) Roundness = 0.;
1626  if( -0.00001 < Angle && Angle < 0 ) Angle = 0.;
1627 
1628  shapes.push_back( Roundness );
1629  shapes.push_back( Angle );
1630  return shapes;
1631 
1632 }
1633 
1634 
1635 template<bool noZS>
1636 int EcalClusterToolsT<noZS>::nrSaturatedCrysIn5x5(const DetId& id,const EcalRecHitCollection* recHits,const CaloTopology *topology)
1637 {
1638  int nrSat=0;
1639  CaloNavigator<DetId> cursor = CaloNavigator<DetId>( id, topology->getSubdetectorTopology( id ) );
1640 
1641  for ( int eastNr = -2; eastNr <= 2; ++eastNr ) { //east is eta in barrel
1642  for ( int northNr = -2; northNr <= 2; ++northNr ) { //north is phi in barrel
1643  cursor.home();
1644  cursor.offsetBy( eastNr, northNr);
1645  DetId id = *cursor;
1646  auto recHitIt = recHits->find(id);
1647  if(recHitIt!=recHits->end() &&
1648  recHitIt->checkFlag(EcalRecHit::kSaturated)){
1649  nrSat++;
1650  }
1651 
1652  }
1653  }
1654  return nrSat;
1655 }
1656 
1657 
1658 //private functions useful for roundnessBarrelSuperClusters etc.
1659 //compute delta iphi between a seed and a particular recHit
1660 //iphi [1,360]
1661 //safe gaurds are put in to ensure the difference is between [-180,180]
1662 template<bool noZS>
1663 int EcalClusterToolsT<noZS>::deltaIPhi(int seed_iphi, int rh_iphi){
1664  int rel_iphi = rh_iphi - seed_iphi;
1665  // take care of cyclic variable iphi [1,360]
1666  if(rel_iphi > 180) rel_iphi = rel_iphi - 360;
1667  if(rel_iphi < -180) rel_iphi = rel_iphi + 360;
1668  return rel_iphi;
1669 }
1670 
1671 //compute delta ieta between a seed and a particular recHit
1672 //ieta [-85,-1] and [1,85]
1673 //safe gaurds are put in to shift the negative ieta +1 to make an ieta=0 so differences are computed correctly
1674 template<bool noZS>
1675 int EcalClusterToolsT<noZS>::deltaIEta(int seed_ieta, int rh_ieta){
1676  // get rid of the fact that there is no ieta=0
1677  if(seed_ieta < 0) seed_ieta++;
1678  if(rh_ieta < 0) rh_ieta++;
1679  int rel_ieta = rh_ieta - seed_ieta;
1680  return rel_ieta;
1681 }
1682 
1683 //return (ieta,iphi) of highest energy recHit of the recHits passed to this function
1684 template<bool noZS>
1685 std::vector<int> EcalClusterToolsT<noZS>::getSeedPosition(const std::vector<std::pair<const EcalRecHit*, float> >& RH_ptrs_fracs){
1686  std::vector<int> seedPosition;
1687  float eSeedRH = 0;
1688  int iEtaSeedRH = 0;
1689  int iPhiSeedRH = 0;
1690 
1691  for(std::vector<std::pair<const EcalRecHit*,float> >::const_iterator rhf_ptr = RH_ptrs_fracs.begin(); rhf_ptr != RH_ptrs_fracs.end(); rhf_ptr++){
1692  const EcalRecHit* rh_ptr = rhf_ptr->first;
1693  //get iEta, iPhi
1694  EBDetId EBdetIdi( rh_ptr->detid() );
1695  float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
1696 
1697  if(eSeedRH < rh_energy){
1698  eSeedRH = rh_energy;
1699  iEtaSeedRH = EBdetIdi.ieta();
1700  iPhiSeedRH = EBdetIdi.iphi();
1701  }
1702 
1703  }// for loop
1704 
1705  seedPosition.push_back(iEtaSeedRH);
1706  seedPosition.push_back(iPhiSeedRH);
1707  return seedPosition;
1708 }
1709 
1710 // return the total energy of the recHits passed to this function
1711 template<bool noZS>
1712 float EcalClusterToolsT<noZS>::getSumEnergy(const std::vector<std::pair<const EcalRecHit*, float> >& RH_ptrs_fracs){
1713  float sumE = 0.;
1714  for( const auto& hAndF : RH_ptrs_fracs ) {
1715  sumE += hAndF.first->energy() * (noZS ? 1.0 : hAndF.second);
1716  }
1717  return sumE;
1718 }
1719 
1720 
1721 typedef EcalClusterToolsT<false> EcalClusterTools;
1722 
1723 namespace noZS {
1724  typedef EcalClusterToolsT<true> EcalClusterTools;
1725 }
1726 
1727 #endif
LogDebug
#define LogDebug(id)
Definition: PFTrackAlgoTools.cc:503
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:971
EcalClusterToolsT::nrSaturatedCrysIn5x5
static int nrSaturatedCrysIn5x5(const DetId &id, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:1636
CaloGeometry::getSubdetectorGeometry
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:45
photonPostprocessing_cfi.eMax
eMax
Definition: photonPostprocessing_cfi.py:22
pfBoostedDoubleSVAK8TagInfos_cfi.R0
R0
Definition: pfBoostedDoubleSVAK8TagInfos_cfi.py:8
mps_fire.i
i
Definition: mps_fire.py:156
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:494
EcalClusterToolsT::meanClusterPosition
static math::XYZVector meanClusterPosition(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, const CaloGeometry *geometry)
Definition: EcalClusterTools.h:731
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:817
w
const double w
Definition: UKUtility.cc:23
EcalClusterTools
EcalClusterToolsT< false > EcalClusterTools
Definition: EcalClusterTools.h:1721
mps_fire.result
result
Definition: mps_fire.py:113
ecaldqm::topology
CaloTopology const * topology(0)
EcalClusterToolsT::matrixDetId
static std::vector< DetId > matrixDetId(const CaloTopology *topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax)
Definition: EcalClusterTools.h:380
EcalClusterToolsT::eMax
static float eMax(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits)
Definition: EcalClusterTools.h:451
hpstanc_transforms.max
max
Definition: hpstanc_transforms.py:5
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:1458
EcalClusterToolsT::n5x5
static int n5x5(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:444
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
EcalClusterToolsT::f33
static double f33(double r)
Definition: EcalClusterTools.h:220
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
EcalRecHit::detid
const DetId & detid() const
Definition: EcalRecHit.h:72
edm::SortedCollection< EcalRecHit >::const_iterator
std::vector< EcalRecHit >::const_iterator const_iterator
Definition: SortedCollection.h:82
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
EcalClusterToolsT::f11
static double f11(double r)
Definition: EcalClusterTools.h:216
weight
Definition: weight.py:1
reco::CaloCluster::hitsAndFractions
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
Definition: CaloCluster.h:195
EBDetId::kTowersInPhi
static const int kTowersInPhi
Definition: EBDetId.h:146
EcalClusterToolsT::deltaIEta
static int deltaIEta(int seed_ieta, int rh_ieta)
Definition: EcalClusterTools.h:1675
EcalClusterToolsT::getSeedPosition
static std::vector< int > getSeedPosition(const std::vector< std::pair< const EcalRecHit *, float > > &RH_ptrs)
Definition: EcalClusterTools.h:1685
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
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:1032
EcalClusterToolsT::deltaIPhi
static int deltaIPhi(int seed_iphi, int rh_iphi)
Definition: EcalClusterTools.h:1663
lat
Definition: DQMStore.h:25
EcalClusterToolsT::getEnergyDepTopology
static std::vector< EcalClusterEnergyDeposition > getEnergyDepTopology(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloGeometry *geometry, bool logW, float w0)
Definition: EcalClusterTools.h:609
res
Definition: Electron.h:4
EcalClusterToolsT::roundnessBarrelSuperClusters
static std::vector< float > roundnessBarrelSuperClusters(const reco::SuperCluster &superCluster, const EcalRecHitCollection &recHits, int weightedPositionMethod=0, float energyThreshold=0.0)
Definition: EcalClusterTools.h:1438
EcalClusterToolsT::e2x5Top
static float e2x5Top(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:487
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:996
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
CaloNavigator::offsetBy
T offsetBy(int deltaX, int deltaY) const
Free movement of arbitray steps.
Definition: CaloNavigator.h:80
EBDetId::kCrystalsInPhi
static const int kCrystalsInPhi
Definition: EBDetId.h:149
EcalClusterToolsT::e3x1
static float e3x1(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:540
EcalClusterToolsT::getDPhiEndcap
static float getDPhiEndcap(const DetId &crysId, float meanX, float meanY)
Definition: EcalClusterTools.h:1165
EcalClusterToolsT::getIEta
static float getIEta(const DetId &id)
Definition: EcalClusterTools.h:1066
EcalClusterToolsT::f55
static double f55(double r)
Definition: EcalClusterTools.h:226
EcalClusterToolsT::e2x2
static float e2x2(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:396
EcalClusterToolsT::getIPhi
static float getIPhi(const DetId &id)
Definition: EcalClusterTools.h:1088
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
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:288
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:457
EBDetId::kModulesPerSM
static const int kModulesPerSM
Definition: EBDetId.h:147
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
EcalClusterToolsT::f40
static double f40(double r)
Definition: EcalClusterTools.h:221
EcalClusterToolsT::e2x5Max
static float e2x5Max(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:503
EcalRecHit::energy
float energy() const
Definition: EcalRecHit.h:68
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:568
EcalClusterToolsT::getNrCrysDiffInEta
static float getNrCrysDiffInEta(const DetId &crysId, const DetId &orginId)
Definition: EcalClusterTools.h:1127
EcalClusterToolsT::e2x5Right
static float e2x5Right(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:473
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:1186
EcalClusterToolsT::getFraction
static float getFraction(const std::vector< std::pair< DetId, float > > &v_id, DetId id)
Definition: EcalClusterTools.h:259
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:408
EcalClusterToolsT::f31
static double f31(double r)
Definition: EcalClusterTools.h:219
EcalClusterToolsT::getNormedIY
static float getNormedIY(const DetId &id)
Definition: EcalClusterTools.h:1114
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:232
EcalClusterToolsT::f44
static double f44(double r)
Definition: EcalClusterTools.h:223
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:983
EcalClusterToolsT::getSumEnergy
static float getSumEnergy(const std::vector< std::pair< const EcalRecHit *, float > > &RH_ptrs_fracs)
Definition: EcalClusterTools.h:1712
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:335
EcalClusterToolsT::getNrCrysDiffInPhi
static float getNrCrysDiffInPhi(const DetId &crysId, const DetId &orginId)
Definition: EcalClusterTools.h:1148
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:51
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:310
EcalClusterToolsT::mean5x5PositionInXY
static std::pair< float, float > mean5x5PositionInXY(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:789
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:1268
EcalClusterToolsT::f00
static double f00(double r)
Definition: EcalClusterTools.h:215
EcalClusterToolsT::roundnessSelectedBarrelRecHits
static std::vector< float > roundnessSelectedBarrelRecHits(const std::vector< std::pair< const EcalRecHit *, float > > &rhVector, int weightedPositionMethod=0)
Definition: EcalClusterTools.h:1504
reco::e1
Float e1
Definition: deltaR.h:20
CaloNavigator::home
void home() const
move the navigator back to the starting point
Definition: CaloNavigator.h:114
EcalClusterToolsT::eTop
static float eTop(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:561
DetId
Definition: DetId.h:18
EcalClusterToolsT::e4x4
static float e4x4(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:426
EcalClusterToolsT::e5x1
static float e5x1(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:526
EcalClusterToolsT::f53
static double f53(double r)
Definition: EcalClusterTools.h:225
EcalClusterToolsT::energyBasketFractionEta
static std::vector< float > energyBasketFractionEta(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits)
Definition: EcalClusterTools.h:575
EBDetId::positiveZ
bool positiveZ() const
Definition: EBDetId.h:78
EBDetId::MAX_IPHI
static const int MAX_IPHI
Definition: EBDetId.h:144
cuy.denominator
denominator
Definition: cuy.py:484
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:669
EcalClusterToolsT::recHitEnergy
static float recHitEnergy(DetId id, const EcalRecHitCollection *recHits)
Definition: EcalClusterTools.h:295
EcalClusterToolsT::e2x5Left
static float e2x5Left(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:480
EcalClusterToolsT::energyBasketFractionPhi
static std::vector< float > energyBasketFractionPhi(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits)
Definition: EcalClusterTools.h:592
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:222
reco::e2
Float e2
Definition: deltaR.h:21
gedGsfElectrons_cfi.isBarrel
isBarrel
Definition: gedGsfElectrons_cfi.py:72
EcalClusterToolsT::matrixSize
static int matrixSize(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, DetId id, int ixMin, int ixMax, int iyMin, int iyMax)
Definition: EcalClusterTools.h:360
EcalClusterToolsT::f51
static double f51(double r)
Definition: EcalClusterTools.h:224
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:76
EcalClusterToolsT::f20
static double f20(double r)
Definition: EcalClusterTools.h:217
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:218
mps_update.diff
diff
Definition: mps_update.py:85
edm::SortedCollection::find
iterator find(key_type k)
Definition: SortedCollection.h:271
heepElectronID_HEEPV50_CSA14_25ns_cff.e5x5
e5x5
Definition: heepElectronID_HEEPV50_CSA14_25ns_cff.py:55
EcalClusterToolsT::mean5x5PositionInLocalCrysCoord
static std::pair< float, float > mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:757
position
static int position[264][3]
Definition: ReadPGInfo.cc:509
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:554
EcalClusterToolsT::e3x3
static float e3x3(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:419
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:977
EcalClusterToolsT::e1x3
static float e1x3(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:533
CastorDataFrameFilter_impl::energySum
double energySum(const DataFrame &df, int fs, int ls)
Definition: CastorDataFrameFilter.cc:22
Geom::pi
constexpr double pi()
Definition: Pi.h:31
EcalClusterToolsT::getNormedIX
static float getNormedIX(const DetId &id)
Definition: EcalClusterTools.h:1101
mps_merge.weight
weight
Definition: mps_merge.py:86
CaloCellGeometry::getPosition
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
Definition: CaloCellGeometry.h:84
Geom::twoPi
constexpr double twoPi()
Definition: Pi.h:32
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:547
EcalClusterToolsT::localCovariances
static std::vector< float > localCovariances(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology, float w0=4.7)
Definition: EcalClusterTools.h:889
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:519
hcalTTPDigis_cfi.id
id
Definition: hcalTTPDigis_cfi.py:7
heepElectronID_HEEPV50_CSA14_25ns_cff.e1x5
e1x5
Definition: heepElectronID_HEEPV50_CSA14_25ns_cff.py:53
edm::SortedCollection::begin
const_iterator begin() const
Definition: SortedCollection.h:303
EcalClusterToolsT::e5x5
static float e5x5(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
Definition: EcalClusterTools.h:437
Generated for CMSSW Reference Manual by doxygen 1.8.11