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