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