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 CxCalculator::CxCalculator (const edm::Event &iEvent, const edm::EventSetup &iSetup, edm::InputTag barrelLabel, edm::InputTag endcapLabel)
24 {
25 //InputTag("islandBasicClusters:islandBarrelBasicClusters")
26 //InputTag("islandBasicClusters:islandEndcapBasicClusters")
27  Handle<BasicClusterCollection> pEBclusters;
28  iEvent.getByLabel(barrelLabel, pEBclusters);
29  fEBclusters_ = pEBclusters.product();
30 
31  Handle<BasicClusterCollection> pEEclusters;
32  iEvent.getByLabel(endcapLabel, pEEclusters);
33  fEEclusters_ = pEEclusters.product();
34 
35  ESHandle<CaloGeometry> geometryHandle;
36  iSetup.get<CaloGeometryRecord>().get(geometryHandle);
37 
38  geometry_ = geometryHandle.product();
39 }
40 
41 double CxCalculator::getCx(const reco::SuperClusterRef cluster, double x, double threshold)
42 {
43  using namespace edm;
44  using namespace reco;
45 
46  if(!fEBclusters_) {
47  LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
48  return -100;
49  }
50 
51  if(!fEEclusters_) {
52  LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
53  return -100;
54  }
55 
56  math::XYZVector SClusPoint(cluster->position().x(),
57  cluster->position().y(),
58  cluster->position().z());
59 
60  double TotalEt = 0;
61 
62  TotalEt = - cluster->rawEnergy()/cosh(cluster->eta());
63 
64  // Loop over barrel basic clusters
65  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
66  iclu != fEBclusters_->end(); ++iclu) {
67  const BasicCluster *clu = &(*iclu);
68  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
69  double eta = ClusPoint.eta();
70 
71  double dR = ROOT::Math::VectorUtil::DeltaR(ClusPoint,SClusPoint);
72 
73  if (dR<x*0.1) {
74  double et = clu->energy()/cosh(eta);
75  if (et<threshold) et=0;
76  TotalEt += et;
77  }
78  }
79 
80  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
81  iclu != fEEclusters_->end(); ++iclu) {
82  const BasicCluster *clu = &(*iclu);
83  const GlobalPoint clusPoint(clu->x(),clu->y(),clu->z());
84  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
85  double eta = ClusPoint.eta();
86 
87  double dR = ROOT::Math::VectorUtil::DeltaR(ClusPoint,SClusPoint);
88 
89  if (dR<x*0.1) {
90  double et = clu->energy()/cosh(eta);
91  if (et<threshold) et=0;
92  TotalEt += et;
93  }
94  }
95 
96  return TotalEt;
97 }
98 
99 double CxCalculator::getCxRemoveSC(const reco::SuperClusterRef cluster, double x, double threshold)
100 {
101  // Calculate Cx and remove the basicClusters used by superCluster
102 
103  using namespace edm;
104  using namespace reco;
105 
106  if(!fEBclusters_) {
107  LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
108  return -100;
109  }
110 
111  if(!fEEclusters_) {
112  LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
113  return -100;
114  }
115 
116  math::XYZVector SClusPoint(cluster->position().x(),
117  cluster->position().y(),
118  cluster->position().z());
119 
120  double TotalEt = 0;
121 
122  TotalEt = 0;
123 
124  // Loop over barrel basic clusters
125  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
126  iclu != fEBclusters_->end(); ++iclu) {
127  const BasicCluster *clu = &(*iclu);
128  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
129  double eta = ClusPoint.eta();
130 
131  double dR = ROOT::Math::VectorUtil::DeltaR(ClusPoint,SClusPoint);
132 
133  // check if this basic cluster is used in the target supercluster
134  bool inSuperCluster = checkUsed(cluster,clu);
135 
136  if (dR<x*0.1&&inSuperCluster==false) {
137  double et = clu->energy()/cosh(eta);
138  if (et<threshold) et=0;
139  TotalEt += et;
140  }
141 
142  }
143 
144  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
145  iclu != fEEclusters_->end(); ++iclu) {
146  const BasicCluster *clu = &(*iclu);
147  const GlobalPoint clusPoint(clu->x(),clu->y(),clu->z());
148  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
149  double eta = ClusPoint.eta();
150 
151  double dR = ROOT::Math::VectorUtil::DeltaR(ClusPoint,SClusPoint);
152 
153  // check if this basic cluster is used in the target supercluster
154  bool inSuperCluster = checkUsed(cluster,clu);
155 
156  if (dR<x*0.1&&inSuperCluster==false) {
157  double et = clu->energy()/cosh(eta);
158  if (et<threshold) et=0;
159  TotalEt += et;
160  }
161  }
162 
163  return TotalEt;
164 }
165 
166 double CxCalculator::getCCx(const reco::SuperClusterRef cluster, double x, double threshold)
167 {
168  using namespace edm;
169  using namespace reco;
170 
171 
172  if(!fEBclusters_) {
173  LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
174  return -100;
175  }
176 
177  if(!fEEclusters_) {
178  LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
179  return -100;
180  }
181 
182  double SClusterEta = cluster->eta();
183  double SClusterPhi = cluster->phi();
184  double TotalEt = 0;
185 
186  TotalEt = - cluster->rawEnergy()/cosh(cluster->eta());
187 
188  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
189  iclu != fEBclusters_->end(); ++iclu) {
190  const BasicCluster *clu = &(*iclu);
191  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
192  double eta = ClusPoint.eta();
193 
194  double dEta = fabs(eta-SClusterEta);
195 
196  if (dEta<x*0.1) {
197  double et = clu->energy()/cosh(eta);
198  if (et<threshold) et=0;
199  TotalEt += et;
200  }
201  }
202 
203  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
204  iclu != fEEclusters_->end(); ++iclu) {
205  const BasicCluster *clu = &(*iclu);
206  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
207  double eta = ClusPoint.eta();
208  double phi = ClusPoint.phi();
209 
210  double dEta = fabs(eta-SClusterEta);
211  double dPhi = fabs(phi-SClusterPhi);
212  while (dPhi>2*PI) dPhi-=2*PI;
213  if (dPhi>PI) dPhi=2*PI-dPhi;
214 
215  if (dEta<x*0.1) {
216  double et = clu->energy()/cosh(eta);
217  if (et<threshold) et=0;
218  TotalEt += et;
219  }
220  }
221 
222  double Cx = getCx(cluster,x,threshold);
223  double CCx = Cx - TotalEt / 40.0 * x;
224 
225  return CCx;
226 }
227 
228 
229 double CxCalculator::getCCxRemoveSC(const reco::SuperClusterRef cluster, double x, double threshold)
230 {
231  using namespace edm;
232  using namespace reco;
233 
234 
235  if(!fEBclusters_) {
236  LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
237  return -100;
238  }
239 
240  if(!fEEclusters_) {
241  LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
242  return -100;
243  }
244 
245  double SClusterEta = cluster->eta();
246  double SClusterPhi = cluster->phi();
247  double TotalEt = 0;
248 
249  TotalEt = 0;
250 
251  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
252  iclu != fEBclusters_->end(); ++iclu) {
253  const BasicCluster *clu = &(*iclu);
254  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
255  double eta = ClusPoint.eta();
256 
257  double dEta = fabs(eta-SClusterEta);
258 
259  // check if this basic cluster is used in the target supercluster
260  bool inSuperCluster = checkUsed(cluster,clu);
261 
262  if (dEta<x*0.1&&inSuperCluster==false) {
263  double et = clu->energy()/cosh(eta);
264  if (et<threshold) et=0;
265  TotalEt += et;
266  }
267  }
268 
269  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
270  iclu != fEEclusters_->end(); ++iclu) {
271  const BasicCluster *clu = &(*iclu);
272  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
273  double eta = ClusPoint.eta();
274  double phi = ClusPoint.phi();
275 
276  double dEta = fabs(eta-SClusterEta);
277  double dPhi = fabs(phi-SClusterPhi);
278  while (dPhi>2*PI) dPhi-=2*PI;
279  if (dPhi>PI) dPhi=2*PI-dPhi;
280 
281  // check if this basic cluster is used in the target supercluster
282  bool inSuperCluster = checkUsed(cluster,clu);
283 
284  if (dEta<x*0.1&&inSuperCluster==false) {
285  double et = clu->energy()/cosh(eta);
286  if (et<threshold) et=0;
287  TotalEt += et;
288  }
289  }
290 
291  double Cx = getCxRemoveSC(cluster,x,threshold);
292  double CCx = Cx - TotalEt / 40.0 * x;
293 
294  return CCx;
295 }
296 
297 
298 bool CxCalculator::checkUsed(const reco::SuperClusterRef sc, const reco::BasicCluster* bc)
299 {
300  reco::CaloCluster_iterator theEclust = sc->clustersBegin();
301 
302  // Loop over the basicClusters inside the target superCluster
303  for(;theEclust != sc->clustersEnd(); theEclust++) {
304  if ((**theEclust) == (*bc) ) return true; //matched, so it's used.
305  }
306  return false;
307 }
308 
309 double CxCalculator::getBCMax(const reco::SuperClusterRef cluster,int i)
310 {
311  reco::CaloCluster_iterator theEclust = cluster->clustersBegin();
312 
313  double energyMax=0,energySecond=0;
314  // Loop over the basicClusters inside the target superCluster
315  for(;theEclust != cluster->clustersEnd(); theEclust++) {
316  if ((*theEclust)->energy()>energyMax ) {
317  energySecond=energyMax;
318  energyMax=(*theEclust)->energy();
319  } else if ((*theEclust)->energy()>energySecond) {
320  energySecond=(*theEclust)->energy();
321  }
322  }
323  if (i==1) return energyMax;
324  return energySecond;
325 }
326 
327 
328 double CxCalculator::getCorrection(const reco::SuperClusterRef cluster, double x, double y,double threshold)
329 {
330  using namespace edm;
331  using namespace reco;
332 
333  // doesn't really work now ^^; (Yen-Jie)
334  if(!fEBclusters_) {
335  LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
336  return -100;
337  }
338 
339  if(!fEEclusters_) {
340  LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
341  return -100;
342  }
343 
344  double SClusterEta = cluster->eta();
345  double SClusterPhi = cluster->phi();
346  double TotalEnergy = 0;
347  double TotalBC = 0;
348 
349  TotalEnergy = 0;
350 
351  double Area = PI * (-x*x+y*y) / 100.0;
352  double nCrystal = Area / 0.0174 / 0.0174; // ignore the difference between endcap and barrel for the moment....
353 
354  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
355  iclu != fEBclusters_->end(); ++iclu) {
356  const BasicCluster *clu = &(*iclu);
357  const GlobalPoint clusPoint(clu->x(),clu->y(),clu->z());
358  double eta = clusPoint.eta();
359  double phi = clusPoint.phi();
360  double dEta = fabs(eta-SClusterEta);
361  double dPhi = fabs(phi-SClusterPhi);
362  while (dPhi>2*PI) dPhi-=2*PI;
363  if (dPhi>PI) dPhi=2*PI-dPhi;
364  double dR = sqrt(dEta*dEta+dPhi*dPhi);
365 
366  if (dR>x*0.1&&dR<y*0.1) {
367  double e = clu->energy();
368  if (e<threshold) e=0;
369  TotalEnergy += e;
370  if (e!=0) TotalBC+=clu->size(); // number of crystals
371 
372  }
373  }
374 
375  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
376  iclu != fEEclusters_->end(); ++iclu) {
377  const BasicCluster *clu = &(*iclu);
378  const GlobalPoint clusPoint(clu->x(),clu->y(),clu->z());
379  double eta = clusPoint.eta();
380  double phi = clusPoint.phi();
381  double dEta = fabs(eta-SClusterEta);
382  double dPhi = fabs(phi-SClusterPhi);
383  while (dPhi>2*PI) dPhi-=2*PI;
384  if (dPhi>PI) dPhi=2*PI-dPhi;
385  double dR = sqrt(dEta*dEta+dPhi*dPhi);
386 
387  if (dR>x*0.1&&dR<y*0.1) {
388  double e = clu->energy();
389  if (e<threshold) e=0;
390  TotalEnergy += e;
391  if (e!=0) TotalBC += clu->size(); // number of crystals
392  }
393  }
394 
395 
396  if (TotalBC==0) return 0;
397  return TotalEnergy/nCrystal;
398 }
399 
400 double CxCalculator::getAvgBCEt(const reco::SuperClusterRef cluster, double x,double phi1,double phi2, double threshold)
401 // x: eta cut, phi1: deltaPhiMin cut, phi2: deltaPhiMax
402 {
403  using namespace edm;
404  using namespace reco;
405 
406 
407  if(!fEBclusters_) {
408  LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
409  return -100;
410  }
411 
412  if(!fEEclusters_) {
413  LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
414  return -100;
415  }
416 
417  double SClusterEta = cluster->eta();
418  double SClusterPhi = cluster->phi();
419 
420  double TotalEt = 0; // Total E
421  double TotalN = 0; // Total N
422 
423  TotalEt = - cluster->rawEnergy()/cosh(cluster->eta());
424 
425  if (fabs(SClusterEta) < 1.479) {
426  //Barrel
427  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
428  iclu != fEBclusters_->end(); ++iclu) {
429  const BasicCluster *clu = &(*iclu);
430  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
431  double eta = ClusPoint.eta();
432  double phi = ClusPoint.phi();
433 
434  double dEta = fabs(eta-SClusterEta);
435  double dPhi = fabs(phi-SClusterPhi);
436  while (dPhi>2*PI) dPhi-=2*PI;
437 
438  bool inSuperCluster = checkUsed(cluster,clu);
439 
440  if (dEta<x*0.1&&inSuperCluster==false&&dPhi>phi1*0.1&&dPhi<phi2*0.1) {
441  double et = clu->energy()/cosh(eta);
442  if (et<threshold) et=0;
443  TotalEt += et;
444  TotalN ++;
445  }
446  }
447  } else {
448  //Endcap
449  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
450  iclu != fEEclusters_->end(); ++iclu) {
451  const BasicCluster *clu = &(*iclu);
452  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
453  double eta = ClusPoint.eta();
454  double phi = ClusPoint.phi();
455 
456  double dEta = fabs(eta-SClusterEta);
457  double dPhi = fabs(phi-SClusterPhi);
458  while (dPhi>2*PI) dPhi-=2*PI;
459 
460  bool inSuperCluster = checkUsed(cluster,clu);
461 
462  if (dEta<x*0.1&&inSuperCluster==false&&dPhi>phi1*0.1&&dPhi<phi2*0.1) {
463  double et = clu->energy()/cosh(eta);
464  if (et<threshold) et=0;
465  TotalEt += et;
466  TotalN ++;
467  }
468  }
469  }
470  return TotalEt / TotalN;
471 }
472 
473 double CxCalculator::getNBC(const reco::SuperClusterRef cluster, double x,double phi1,double phi2, double threshold)
474 // x: eta cut, phi1: deltaPhiMin cut, phi2: deltaPhiMax
475 {
476  using namespace edm;
477  using namespace reco;
478 
479 
480  if(!fEBclusters_) {
481  LogError("CxCalculator") << "Error! Can't get EBclusters for event.";
482  return -100;
483  }
484 
485  if(!fEEclusters_) {
486  LogError("CxCalculator") << "Error! Can't get EEclusters for event.";
487  return -100;
488  }
489 
490  double SClusterEta = cluster->eta();
491  double SClusterPhi = cluster->phi();
492 
493  double TotalEt = 0; // Total E
494  double TotalN = 0; // Total N
495 
496  TotalEt = - cluster->rawEnergy()/cosh(cluster->eta());
497 
498 
499 
500  if (fabs(SClusterEta) < 1.479) {
501  //Barrel
502  for(BasicClusterCollection::const_iterator iclu = fEBclusters_->begin();
503  iclu != fEBclusters_->end(); ++iclu) {
504  const BasicCluster *clu = &(*iclu);
505  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
506  double eta = ClusPoint.eta();
507  double phi = ClusPoint.phi();
508 
509  double dEta = fabs(eta-SClusterEta);
510  double dPhi = fabs(phi-SClusterPhi);
511  while (dPhi>2*PI) dPhi-=2*PI;
512 
513  bool inSuperCluster = checkUsed(cluster,clu);
514 
515  if (dEta<x*0.1&&inSuperCluster==false&&dPhi>phi1*0.1&&dPhi<phi2*0.1) {
516  double et = clu->energy()/cosh(eta);
517  if (et<threshold) et=0;
518  TotalEt += et;
519  TotalN ++;
520  }
521  }
522  } else {
523  //Endcap
524  for(BasicClusterCollection::const_iterator iclu = fEEclusters_->begin();
525  iclu != fEEclusters_->end(); ++iclu) {
526  const BasicCluster *clu = &(*iclu);
527  math::XYZVector ClusPoint(clu->x(),clu->y(),clu->z());
528  double eta = ClusPoint.eta();
529  double phi = ClusPoint.phi();
530 
531  double dEta = fabs(eta-SClusterEta);
532  double dPhi = fabs(phi-SClusterPhi);
533  while (dPhi>2*PI) dPhi-=2*PI;
534 
535  bool inSuperCluster = checkUsed(cluster,clu);
536 
537  if (dEta<x*0.1&&inSuperCluster==false&&dPhi>phi1*0.1&&dPhi<phi2*0.1) {
538  double et = clu->energy()/cosh(eta);
539  if (et<threshold) et=0;
540  TotalEt += et;
541  TotalN ++;
542  }
543  }
544  }
545  return TotalN;
546 }
i
int i
Definition: DBlmapReader.cc:9
PI
#define PI
Definition: HcalBeamMonitor.cc:14
CxCalculator::getAvgBCEt
double getAvgBCEt(const reco::SuperClusterRef clus, double eta, double phi1, double phi2, double threshold)
Definition: CxCalculator.cc:400
CxCalculator::CxCalculator
CxCalculator(const edm::Event &iEvent, const edm::EventSetup &iSetup, edm::InputTag barrelLabel, edm::InputTag endcapLabel)
Definition: CxCalculator.cc:23
CxCalculator::getNBC
double getNBC(const reco::SuperClusterRef clus, double eta, double phi1, double phi2, double threshold)
Definition: CxCalculator.cc:473
eta
T eta() const
Definition: Basic3DVectorLD.h:168
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:101
CxCalculator::getCCxRemoveSC
double getCCxRemoveSC(const reco::SuperClusterRef clus, double i, double threshold)
Definition: CxCalculator.cc:229
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:298
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
CxCalculator::getCCx
double getCCx(const reco::SuperClusterRef clus, double i, double threshold)
Definition: CxCalculator.cc:166
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:355
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:99
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:41
edm::Handle::product
T const * product() const
Definition: Handle.h:74
x
Definition: DDAxes.h:10
CxCalculator::getCorrection
double getCorrection(const reco::SuperClusterRef clus, double i, double j, double threshold)
Definition: CxCalculator.cc:328
CxCalculator::getBCMax
double getBCMax(const reco::SuperClusterRef clus, int i)
Definition: CxCalculator.cc:309
dtT0WireCalibration_cfg.threshold
tuple threshold
Definition: dtT0WireCalibration_cfg.py:10
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5