CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CxCalculator.cc
Go to the documentation of this file.
1 // ROOT includes
2 #include <Math/VectorUtil.h>
3 
4 #include "RecoHI/HiEgammaAlgos/interface/CxCalculator.h"
5 #include "FWCore/MessageLogger/interface/MessageLogger.h"
6 
7 #include "Geometry/Records/interface/CaloGeometryRecord.h"
8 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
9 
10 #include "DataFormats/Common/interface/Handle.h"
11 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
12 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
13 #include "DataFormats/Math/interface/Vector3D.h"
14 
15 
16 using namespace edm;
17 using namespace reco;
18 using namespace std;
19 using namespace ROOT::Math::VectorUtil;
20 
21 #define PI 3.141592653589793238462643383279502884197169399375105820974945
22 
23 
24 double CxCalculator::getJurassicArea( double r1, double r2, double width) {
25 
26  float theta1 = asin( width / r1);
27  float theta2 = asin( width / r2);
28  float theA = sqrt ( r1*r1 + r2*r2 - 2 * r1 * r2 * cos ( theta1 - theta2) );
29  float area1 = 0.5 * r1*r1 * ( 3.141592 - 2 * theta1 ) ;
30  float area2 = 0.5 * r2*r2 * ( 3.141592 - 2 * theta2 ) ;
31  float area3 = width * theA;
32  float finalArea = 2 * ( area1 - area2 - area3);
33  return finalArea;
34 }
35 
36 
37 
38 CxCalculator::CxCalculator (const edm::Event &iEvent, const edm::EventSetup &iSetup, edm::InputTag barrelLabel, edm::InputTag endcapLabel)
39 {
40 //InputTag("islandBasicClusters:islandBarrelBasicClusters")
41 //InputTag("islandBasicClusters:islandEndcapBasicClusters")
42  Handle<BasicClusterCollection> pEBclusters;
43  iEvent.getByLabel(barrelLabel, pEBclusters);
44  if(pEBclusters.isValid())
45  fEBclusters_ = pEBclusters.product();
46  else
47  fEBclusters_ = NULL;
48 
49  Handle<BasicClusterCollection> pEEclusters;
50  iEvent.getByLabel(endcapLabel, pEEclusters);
51  if(pEEclusters.isValid())
52  fEEclusters_ = pEEclusters.product();
53  else
54  fEEclusters_ = NULL;
55 
56  ESHandle<CaloGeometry> geometryHandle;
57  iSetup.get<CaloGeometryRecord>().get(geometryHandle);
58  if(geometryHandle.isValid())
59  geometry_ = geometryHandle.product();
60  else
61  geometry_ = NULL;
62 
63 }
64 
65 double CxCalculator::getCx(const reco::SuperClusterRef cluster, double x, double threshold)
66 {
67  using namespace edm;
68  using namespace reco;
69 
70  if(!fEBclusters_) {
71 // LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
72  return -100;
73  }
74 
75  if(!fEEclusters_) {
76 // LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
77  return -100;
78  }
79 
80  math::XYZVector SClusPoint(cluster->position().x(),
81  cluster->position().y(),
82  cluster->position().z());
83 
84  double TotalEt = 0;
85 
86  TotalEt = - cluster->rawEnergy()/cosh(cluster->eta());
87 
88  // Loop over barrel basic clusters
89  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
90  iclu != fEBclusters_->end(); ++iclu) {
91  const BasicCluster *clu = &(*iclu);
92  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
93  double eta = ClusPoint.eta();
94 
95  double dR = ROOT::Math::VectorUtil::DeltaR(ClusPoint,SClusPoint);
96 
97  if (dR<x*0.1) {
98  double et = clu->energy()/cosh(eta);
99  if (et<threshold) et=0;
100  TotalEt += et;
101  }
102  }
103 
104  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
105  iclu != fEEclusters_->end(); ++iclu) {
106  const BasicCluster *clu = &(*iclu);
107  const GlobalPoint clusPoint(clu->x(),clu->y(),clu->z());
108  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
109  double eta = ClusPoint.eta();
110 
111  double dR = ROOT::Math::VectorUtil::DeltaR(ClusPoint,SClusPoint);
112 
113  if (dR<x*0.1) {
114  double et = clu->energy()/cosh(eta);
115  if (et<threshold) et=0;
116  TotalEt += et;
117  }
118  }
119 
120  return TotalEt;
121 }
122 
123 double CxCalculator::getCxRemoveSC(const reco::SuperClusterRef cluster, double x, double threshold)
124 {
125  // Calculate Cx and remove the basicClusters used by superCluster
126 
127  using namespace edm;
128  using namespace reco;
129 
130  if(!fEBclusters_) {
131  LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
132  return -100;
133  }
134 
135  if(!fEEclusters_) {
136  LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
137  return -100;
138  }
139 
140  math::XYZVector SClusPoint(cluster->position().x(),
141  cluster->position().y(),
142  cluster->position().z());
143 
144  double TotalEt = 0;
145 
146  TotalEt = 0;
147 
148  // Loop over barrel basic clusters
149  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
150  iclu != fEBclusters_->end(); ++iclu) {
151  const BasicCluster *clu = &(*iclu);
152  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
153  double eta = ClusPoint.eta();
154 
155  double dR = ROOT::Math::VectorUtil::DeltaR(ClusPoint,SClusPoint);
156 
157  // check if this basic cluster is used in the target supercluster
158  bool inSuperCluster = checkUsed(cluster,clu);
159 
160  if (dR<x*0.1&&inSuperCluster==false) {
161  double et = clu->energy()/cosh(eta);
162  if (et<threshold) et=0;
163  TotalEt += et;
164  }
165 
166  }
167 
168  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
169  iclu != fEEclusters_->end(); ++iclu) {
170  const BasicCluster *clu = &(*iclu);
171  const GlobalPoint clusPoint(clu->x(),clu->y(),clu->z());
172  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
173  double eta = ClusPoint.eta();
174 
175  double dR = ROOT::Math::VectorUtil::DeltaR(ClusPoint,SClusPoint);
176 
177  // check if this basic cluster is used in the target supercluster
178  bool inSuperCluster = checkUsed(cluster,clu);
179 
180  if (dR<x*0.1&&inSuperCluster==false) {
181  double et = clu->energy()/cosh(eta);
182  if (et<threshold) et=0;
183  TotalEt += et;
184  }
185  }
186 
187  return TotalEt;
188 }
189 
190 double CxCalculator::getCCx(const reco::SuperClusterRef cluster, double x, double threshold)
191 {
192  using namespace edm;
193  using namespace reco;
194 
195 
196  if(!fEBclusters_) {
197  // LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
198  return -100;
199  }
200 
201  if(!fEEclusters_) {
202  // LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
203  return -100;
204  }
205 
206  double SClusterEta = cluster->eta();
207  double SClusterPhi = cluster->phi();
208  double TotalEt = 0;
209 
210  TotalEt = - cluster->rawEnergy()/cosh(cluster->eta());
211 
212  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
213  iclu != fEBclusters_->end(); ++iclu) {
214  const BasicCluster *clu = &(*iclu);
215  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
216  double eta = ClusPoint.eta();
217 
218  double dEta = fabs(eta-SClusterEta);
219 
220  if (dEta<x*0.1) {
221  double et = clu->energy()/cosh(eta);
222  if (et<threshold) et=0;
223  TotalEt += et;
224  }
225  }
226 
227  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
228  iclu != fEEclusters_->end(); ++iclu) {
229  const BasicCluster *clu = &(*iclu);
230  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
231  double eta = ClusPoint.eta();
232  double phi = ClusPoint.phi();
233 
234  double dEta = fabs(eta-SClusterEta);
235  double dPhi = fabs(phi-SClusterPhi);
236  while (dPhi>2*PI) dPhi-=2*PI;
237  if (dPhi>PI) dPhi=2*PI-dPhi;
238 
239  if (dEta<x*0.1) {
240  double et = clu->energy()/cosh(eta);
241  if (et<threshold) et=0;
242  TotalEt += et;
243  }
244  }
245 
246  double Cx = getCx(cluster,x,threshold);
247  double CCx = Cx - TotalEt / 40.0 * x;
248 
249  return CCx;
250 }
251 
252 
253 
254 
255 double CxCalculator::getJc(const reco::SuperClusterRef cluster, double r1, double r2, double jWidth, double threshold)
256 {
257  using namespace edm;
258  using namespace reco;
259  if(!fEBclusters_) {
260  // LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
261  return -100;
262  }
263  if(!fEEclusters_) {
264  return -100;
265  }
266  double SClusterEta = cluster->eta();
267  double SClusterPhi = cluster->phi();
268  double TotalEt = 0;
269 
270  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
271  iclu != fEBclusters_->end(); ++iclu) {
272  const BasicCluster *clu = &(*iclu);
273  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
274  double eta = ClusPoint.eta();
275  double phi = ClusPoint.phi();
276 
277  double dEta = fabs(eta-SClusterEta);
278  double dPhi = phi-SClusterPhi;
279  if ( dPhi < -PI ) dPhi = dPhi + 2*PI ;
280  if ( dPhi > PI ) dPhi = dPhi - 2*PI ;
281  if ( fabs(dPhi) > PI ) cout << " error!!! dphi > 2pi : " << dPhi << endl;
282  double dR = sqrt(dEta*dEta+dPhi*dPhi);
283 
284  // Jurassic Cone /////
285  if ( dR > r1 ) continue;
286  if ( dR < r2 ) continue;
287  if ( fabs(dEta) < jWidth) continue;
289  double theEt = clu->energy()/cosh(eta);
290  if (theEt<threshold) continue;
291  TotalEt += theEt;
292  }
293 
294  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
295  iclu != fEEclusters_->end(); ++iclu) {
296  const BasicCluster *clu = &(*iclu);
297  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
298  double eta = ClusPoint.eta();
299  double phi = ClusPoint.phi();
300  double dEta = fabs(eta-SClusterEta);
301  double dPhi = phi-SClusterPhi;
302  if ( dPhi < -PI ) dPhi = dPhi + 2*PI ;
303  if ( dPhi > PI ) dPhi = dPhi - 2*PI ;
304  if ( fabs(dPhi) >PI ) cout << " error!!! dphi > 2pi : " << dPhi << endl;
305  double dR = sqrt(dEta*dEta+dPhi*dPhi);
306  // Jurassic Cone /////
307  if ( dR > r1 ) continue;
308  if ( dR < r2 ) continue;
309  if ( fabs(dEta) < jWidth) continue;
311  double theEt = clu->energy()/cosh(eta);
312  if (theEt<threshold) continue;
313  TotalEt += theEt;
314  }
315  return TotalEt;
316 }
317 
318 
319 double CxCalculator::getJcc(const reco::SuperClusterRef cluster, double r1, double r2, double jWidth, double threshold)
320 {
321 
322  using namespace edm;
323  using namespace reco;
324  if(!fEBclusters_) {
325  // LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
326  return -100;
327  }
328  if(!fEEclusters_) {
329  return -100;
330  }
331  double SClusterEta = cluster->eta();
332  double SClusterPhi = cluster->phi();
333  double TotalEt = 0;
334 
335  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
336  iclu != fEBclusters_->end(); ++iclu) {
337  const BasicCluster *clu = &(*iclu);
338  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
339  double eta = ClusPoint.eta();
340  double phi = ClusPoint.phi();
341 
342  double dEta = fabs(eta-SClusterEta);
343  double dPhi = phi-SClusterPhi;
344  if ( dPhi < -PI ) dPhi = dPhi + 2*PI ;
345  if ( dPhi > PI ) dPhi = dPhi - 2*PI ;
346  // double dR = sqrt(dEta*dEta+dPhi*dPhi);
347 
349  if ( fabs(dEta) > r1 ) continue;
350  if ( fabs(dPhi) <r1 ) continue;
352 
353  double theEt = clu->energy()/cosh(eta);
354  if (theEt<threshold) continue;
355  TotalEt += theEt;
356  }
357  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
358  iclu != fEEclusters_->end(); ++iclu) {
359  const BasicCluster *clu = &(*iclu);
360  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
361  double eta = ClusPoint.eta();
362  double phi = ClusPoint.phi();
363  double dEta = fabs(eta-SClusterEta);
364  double dPhi = phi-SClusterPhi;
365  if ( dPhi < -PI ) dPhi = dPhi + 2*PI ;
366  if ( dPhi > PI ) dPhi = dPhi - 2*PI ;
367  // double dR = sqrt(dEta*dEta+dPhi*dPhi);
368 
370  if ( fabs(dEta) > r1 ) continue;
371  if ( fabs(dPhi) < r1 ) continue;
373 
374  double theEt = clu->energy()/cosh(eta);
375  if (theEt<threshold) continue;
376  TotalEt += theEt;
377  }
378 
379  double areaStrip = 4*PI*r1 - 4*r1*r1;
380  double areaJura = getJurassicArea(r1,r2, jWidth) ;
381  double theCJ = getJc(cluster,r1, r2, jWidth, threshold);
382  // cout << "areJura = " << areaJura << endl;
383  // cout << "areaStrip " << areaStrip << endl;
384  double theCCJ = theCJ - TotalEt * areaJura / areaStrip ;
385  return theCCJ;
386 }
387 
388 
389 
390 double CxCalculator::getCCxRemoveSC(const reco::SuperClusterRef cluster, double x, double threshold)
391 {
392  using namespace edm;
393  using namespace reco;
394 
395 
396  if(!fEBclusters_) {
397  LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
398  return -100;
399  }
400 
401  if(!fEEclusters_) {
402  LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
403  return -100;
404  }
405 
406  double SClusterEta = cluster->eta();
407  double SClusterPhi = cluster->phi();
408  double TotalEt = 0;
409 
410  TotalEt = 0;
411 
412  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
413  iclu != fEBclusters_->end(); ++iclu) {
414  const BasicCluster *clu = &(*iclu);
415  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
416  double eta = ClusPoint.eta();
417 
418  double dEta = fabs(eta-SClusterEta);
419 
420  // check if this basic cluster is used in the target supercluster
421  bool inSuperCluster = checkUsed(cluster,clu);
422 
423  if (dEta<x*0.1&&inSuperCluster==false) {
424  double et = clu->energy()/cosh(eta);
425  if (et<threshold) et=0;
426  TotalEt += et;
427  }
428  }
429 
430  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
431  iclu != fEEclusters_->end(); ++iclu) {
432  const BasicCluster *clu = &(*iclu);
433  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
434  double eta = ClusPoint.eta();
435  double phi = ClusPoint.phi();
436 
437  double dEta = fabs(eta-SClusterEta);
438  double dPhi = fabs(phi-SClusterPhi);
439  while (dPhi>2*PI) dPhi-=2*PI;
440  if (dPhi>PI) dPhi=2*PI-dPhi;
441 
442  // check if this basic cluster is used in the target supercluster
443  bool inSuperCluster = checkUsed(cluster,clu);
444 
445  if (dEta<x*0.1&&inSuperCluster==false) {
446  double et = clu->energy()/cosh(eta);
447  if (et<threshold) et=0;
448  TotalEt += et;
449  }
450  }
451 
452  double Cx = getCxRemoveSC(cluster,x,threshold);
453  double CCx = Cx - TotalEt / 40.0 * x;
454 
455  return CCx;
456 }
457 
458 
459 bool CxCalculator::checkUsed(const reco::SuperClusterRef sc, const reco::BasicCluster* bc)
460 {
461  reco::CaloCluster_iterator theEclust = sc->clustersBegin();
462 
463  // Loop over the basicClusters inside the target superCluster
464  for(;theEclust != sc->clustersEnd(); theEclust++) {
465  if ((**theEclust) == (*bc) ) return true; //matched, so it's used.
466  }
467  return false;
468 }
469 
470 double CxCalculator::getBCMax(const reco::SuperClusterRef cluster,int i)
471 {
472  reco::CaloCluster_iterator theEclust = cluster->clustersBegin();
473 
474  double energyMax=0,energySecond=0;
475  // Loop over the basicClusters inside the target superCluster
476  for(;theEclust != cluster->clustersEnd(); theEclust++) {
477  if ((*theEclust)->energy()>energyMax ) {
478  energySecond=energyMax;
479  energyMax=(*theEclust)->energy();
480  } else if ((*theEclust)->energy()>energySecond) {
481  energySecond=(*theEclust)->energy();
482  }
483  }
484  if (i==1) return energyMax;
485  return energySecond;
486 }
487 
488 
489 double CxCalculator::getCorrection(const reco::SuperClusterRef cluster, double x, double y,double threshold)
490 {
491  using namespace edm;
492  using namespace reco;
493 
494  // doesn't really work now ^^; (Yen-Jie)
495  if(!fEBclusters_) {
496  LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
497  return -100;
498  }
499 
500  if(!fEEclusters_) {
501  LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
502  return -100;
503  }
504 
505  double SClusterEta = cluster->eta();
506  double SClusterPhi = cluster->phi();
507  double TotalEnergy = 0;
508  double TotalBC = 0;
509 
510  TotalEnergy = 0;
511 
512  double Area = PI * (-x*x+y*y) / 100.0;
513  double nCrystal = Area / 0.0174 / 0.0174; // ignore the difference between endcap and barrel for the moment....
514 
515  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
516  iclu != fEBclusters_->end(); ++iclu) {
517  const BasicCluster *clu = &(*iclu);
518  const GlobalPoint clusPoint(clu->x(),clu->y(),clu->z());
519  double eta = clusPoint.eta();
520  double phi = clusPoint.phi();
521  double dEta = fabs(eta-SClusterEta);
522  double dPhi = fabs(phi-SClusterPhi);
523  while (dPhi>2*PI) dPhi-=2*PI;
524  if (dPhi>PI) dPhi=2*PI-dPhi;
525  double dR = sqrt(dEta*dEta+dPhi*dPhi);
526 
527  if (dR>x*0.1&&dR<y*0.1) {
528  double e = clu->energy();
529  if (e<threshold) e=0;
530  TotalEnergy += e;
531  if (e!=0) TotalBC+=clu->size(); // number of crystals
532 
533  }
534  }
535 
536  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
537  iclu != fEEclusters_->end(); ++iclu) {
538  const BasicCluster *clu = &(*iclu);
539  const GlobalPoint clusPoint(clu->x(),clu->y(),clu->z());
540  double eta = clusPoint.eta();
541  double phi = clusPoint.phi();
542  double dEta = fabs(eta-SClusterEta);
543  double dPhi = fabs(phi-SClusterPhi);
544  while (dPhi>2*PI) dPhi-=2*PI;
545  if (dPhi>PI) dPhi=2*PI-dPhi;
546  double dR = sqrt(dEta*dEta+dPhi*dPhi);
547 
548  if (dR>x*0.1&&dR<y*0.1) {
549  double e = clu->energy();
550  if (e<threshold) e=0;
551  TotalEnergy += e;
552  if (e!=0) TotalBC += clu->size(); // number of crystals
553  }
554  }
555 
556 
557  if (TotalBC==0) return 0;
558  return TotalEnergy/nCrystal;
559 }
560 
561 double CxCalculator::getAvgBCEt(const reco::SuperClusterRef cluster, double x,double phi1,double phi2, double threshold)
562 // x: eta cut, phi1: deltaPhiMin cut, phi2: deltaPhiMax
563 {
564  using namespace edm;
565  using namespace reco;
566 
567 
568  if(!fEBclusters_) {
569  LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
570  return -100;
571  }
572 
573  if(!fEEclusters_) {
574  LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
575  return -100;
576  }
577 
578  double SClusterEta = cluster->eta();
579  double SClusterPhi = cluster->phi();
580 
581  double TotalEt = 0; // Total E
582  double TotalN = 0; // Total N
583 
584  TotalEt = - cluster->rawEnergy()/cosh(cluster->eta());
585 
586  if (fabs(SClusterEta) < 1.479) {
587  //Barrel
588  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
589  iclu != fEBclusters_->end(); ++iclu) {
590  const BasicCluster *clu = &(*iclu);
591  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
592  double eta = ClusPoint.eta();
593  double phi = ClusPoint.phi();
594 
595  double dEta = fabs(eta-SClusterEta);
596  double dPhi = fabs(phi-SClusterPhi);
597  while (dPhi>2*PI) dPhi-=2*PI;
598 
599  bool inSuperCluster = checkUsed(cluster,clu);
600 
601  if (dEta<x*0.1&&inSuperCluster==false&&dPhi>phi1*0.1&&dPhi<phi2*0.1) {
602  double et = clu->energy()/cosh(eta);
603  if (et<threshold) et=0;
604  TotalEt += et;
605  TotalN ++;
606  }
607  }
608  } else {
609  //Endcap
610  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
611  iclu != fEEclusters_->end(); ++iclu) {
612  const BasicCluster *clu = &(*iclu);
613  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
614  double eta = ClusPoint.eta();
615  double phi = ClusPoint.phi();
616 
617  double dEta = fabs(eta-SClusterEta);
618  double dPhi = fabs(phi-SClusterPhi);
619  while (dPhi>2*PI) dPhi-=2*PI;
620 
621  bool inSuperCluster = checkUsed(cluster,clu);
622 
623  if (dEta<x*0.1&&inSuperCluster==false&&dPhi>phi1*0.1&&dPhi<phi2*0.1) {
624  double et = clu->energy()/cosh(eta);
625  if (et<threshold) et=0;
626  TotalEt += et;
627  TotalN ++;
628  }
629  }
630  }
631  return TotalEt / TotalN;
632 }
633 
634 double CxCalculator::getNBC(const reco::SuperClusterRef cluster, double x,double phi1,double phi2, double threshold)
635 // x: eta cut, phi1: deltaPhiMin cut, phi2: deltaPhiMax
636 {
637  using namespace edm;
638  using namespace reco;
639 
640 
641  if(!fEBclusters_) {
642  LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
643  return -100;
644  }
645 
646  if(!fEEclusters_) {
647  LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
648  return -100;
649  }
650 
651  double SClusterEta = cluster->eta();
652  double SClusterPhi = cluster->phi();
653 
654  double TotalEt = 0; // Total E
655  double TotalN = 0; // Total N
656 
657  TotalEt = - cluster->rawEnergy()/cosh(cluster->eta());
658 
659 
660 
661  if (fabs(SClusterEta) < 1.479) {
662  //Barrel
663  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
664  iclu != fEBclusters_->end(); ++iclu) {
665  const BasicCluster *clu = &(*iclu);
666  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
667  double eta = ClusPoint.eta();
668  double phi = ClusPoint.phi();
669 
670  double dEta = fabs(eta-SClusterEta);
671  double dPhi = fabs(phi-SClusterPhi);
672  while (dPhi>2*PI) dPhi-=2*PI;
673 
674  bool inSuperCluster = checkUsed(cluster,clu);
675 
676  if (dEta<x*0.1&&inSuperCluster==false&&dPhi>phi1*0.1&&dPhi<phi2*0.1) {
677  double et = clu->energy()/cosh(eta);
678  if (et<threshold) et=0;
679  TotalEt += et;
680  TotalN ++;
681  }
682  }
683  } else {
684  //Endcap
685  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
686  iclu != fEEclusters_->end(); ++iclu) {
687  const BasicCluster *clu = &(*iclu);
688  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
689  double eta = ClusPoint.eta();
690  double phi = ClusPoint.phi();
691 
692  double dEta = fabs(eta-SClusterEta);
693  double dPhi = fabs(phi-SClusterPhi);
694  while (dPhi>2*PI) dPhi-=2*PI;
695 
696  bool inSuperCluster = checkUsed(cluster,clu);
697 
698  if (dEta<x*0.1&&inSuperCluster==false&&dPhi>phi1*0.1&&dPhi<phi2*0.1) {
699  double et = clu->energy()/cosh(eta);
700  if (et<threshold) et=0;
701  TotalEt += et;
702  TotalN ++;
703  }
704  }
705  }
706  return TotalN;
707 }
i
int i
Definition: DBlmapReader.cc:9
PI
#define PI
Definition: HcalBeamMonitor.cc:14
CxCalculator::getJc
double getJc(const reco::SuperClusterRef cluster, double r1=0.4, double r2=0.06, double jWidth=0.04, double threshold=0)
Definition: CxCalculator.cc:255
CxCalculator::getAvgBCEt
double getAvgBCEt(const reco::SuperClusterRef clus, double eta, double phi1, double phi2, double threshold)
Definition: CxCalculator.cc:561
CxCalculator::getJcc
double getJcc(const reco::SuperClusterRef cluster, double r1=0.4, double r2=0.06, double jWidth=0.04, double threshold=0)
Definition: CxCalculator.cc:319
NULL
#define NULL
Definition: scimark2.h:8
CxCalculator::CxCalculator
CxCalculator(const edm::Event &iEvent, const edm::EventSetup &iSetup, edm::InputTag barrelLabel, edm::InputTag endcapLabel)
Definition: CxCalculator.cc:38
CxCalculator::getNBC
double getNBC(const reco::SuperClusterRef clus, double eta, double phi1, double phi2, double threshold)
Definition: CxCalculator.cc:634
eta
T eta() const
Definition: Basic3DVectorLD.h:172
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:105
CxCalculator::getCCxRemoveSC
double getCCxRemoveSC(const reco::SuperClusterRef clus, double i, double threshold)
Definition: CxCalculator.cc:390
iEvent
int iEvent
Definition: GenABIO.cc:243
dPhi
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
CxCalculator::checkUsed
bool checkUsed(const reco::SuperClusterRef clus, const reco::BasicCluster *clu)
Definition: CxCalculator.cc:459
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
dtDQMClient_cfg.threshold
tuple threshold
Definition: dtDQMClient_cfg.py:16
CxCalculator::getCCx
double getCCx(const reco::SuperClusterRef clus, double i, double threshold)
Definition: CxCalculator.cc:190
PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.dR
tuple dR
Definition: PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.py:19
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:76
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
math::XYZVector
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
CxCalculator::getCxRemoveSC
double getCxRemoveSC(const reco::SuperClusterRef clus, double i, double threshold)
Definition: CxCalculator.cc:123
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:62
CxCalculator::getCx
double getCx(const reco::SuperClusterRef clus, double i, double threshold)
Definition: CxCalculator.cc:65
edm::Handle::product
T const * product() const
Definition: Handle.h:74
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
x
Definition: DDAxes.h:10
CxCalculator::getJurassicArea
double getJurassicArea(double r1, double r2, double width)
Definition: CxCalculator.cc:24
edm::ESHandleBase::isValid
bool isValid() const
Definition: ESHandle.h:37
CxCalculator::getCorrection
double getCorrection(const reco::SuperClusterRef clus, double i, double j, double threshold)
Definition: CxCalculator.cc:489
CxCalculator::getBCMax
double getBCMax(const reco::SuperClusterRef clus, int i)
Definition: CxCalculator.cc:470
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5