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