CMS 3D CMS Logo

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