CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PFEnergyCalibration.cc
Go to the documentation of this file.
1 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
2 #include "CondFormats/PhysicsToolsObjects/interface/PerformancePayloadFromTFormula.h"
3 #include <TMath.h>
4 #include <math.h>
5 #include <vector>
6 #include <TF1.h>
7 #include <map>
8 
9 using namespace std;
10 
11 PFEnergyCalibration::PFEnergyCalibration() : pfCalibrations(0)
12 {
13  initializeCalibrationFunctions();
14 }
15 
16 PFEnergyCalibration::~PFEnergyCalibration()
17 {
18 
19  delete faBarrel;
20  delete fbBarrel;
21  delete fcBarrel;
22  delete faEtaBarrel;
23  delete fbEtaBarrel;
24  delete faEndcap;
25  delete fbEndcap;
26  delete fcEndcap;
27  delete faEtaEndcap;
28  delete fbEtaEndcap;
29 
30 }
31 
32 void
33 PFEnergyCalibration::initializeCalibrationFunctions() {
34 
35  // NEW NEW with HCAL pre-calibration
36 
37  threshE = 3.5;
38  threshH = 2.5;
39 
40  // Barrel (fit made with |eta| < 1.2)
41  faBarrel = new TF1("faBarrel","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
42  fbBarrel = new TF1("fbBarrel","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
43  fcBarrel = new TF1("fcBarrel","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
44  faEtaBarrel = new TF1("faEtaBarrel","[0]+[1]*exp(-x/[2])",1.,1000.);
45  fbEtaBarrel = new TF1("fbEtaBarrel","[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",1.,1000.);
46  faBarrel->SetParameter(0,1.15665);
47  fbBarrel->SetParameter(0,0.994603);
48  fcBarrel->SetParameter(0,0.956544);
49  faEtaBarrel->SetParameter(0,0.014664);
50  fbEtaBarrel->SetParameter(0,0.00975451);
51  faBarrel->SetParameter(1,0.165627);
52  fbBarrel->SetParameter(1,0.13632);
53  fcBarrel->SetParameter(1,0.0857207);
54  faEtaBarrel->SetParameter(1,-0.0426776);
55  fbEtaBarrel->SetParameter(1,0.102247);
56  faBarrel->SetParameter(2,0.827718);
57  fbBarrel->SetParameter(2,-0.758013);
58  fcBarrel->SetParameter(2,-0.44347);
59  faEtaBarrel->SetParameter(2,431.054);
60  fbEtaBarrel->SetParameter(2,436.21);
61  faBarrel->SetParameter(3,231.339);
62  fbBarrel->SetParameter(3,183.627);
63  fcBarrel->SetParameter(3,63.3479);
64  faBarrel->SetParameter(4,2.45332);
65  fbBarrel->SetParameter(4,1);
66  fcBarrel->SetParameter(4,1.24174);
67  faBarrel->SetParameter(5,29.6603);
68  fbBarrel->SetParameter(5,39.6784);
69  fcBarrel->SetParameter(5,12.322);
70 
71  // End-caps (fit made with eta
72  faEndcap = new TF1("faEndcap","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
73  fbEndcap = new TF1("fbEndcap","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
74  fcEndcap = new TF1("fcEndcap","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
75  faEtaEndcap = new TF1("faEtaEndcap","[0]+[1]*exp(-x/[2])",1.,1000.);
76  fbEtaEndcap = new TF1("fbEtaEndcap","[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",1.,1000.);
77  faEndcap->SetParameter(0,1.1272);
78  fbEndcap->SetParameter(0,0.982824);
79  fcEndcap->SetParameter(0,0.950244);
80  faEtaEndcap->SetParameter(0,-0.000582903);
81  fbEtaEndcap->SetParameter(0,0.0267319);
82  faEndcap->SetParameter(1,0.258536);
83  fbEndcap->SetParameter(1,0.0977533);
84  fcEndcap->SetParameter(1,0.00564779);
85  faEtaEndcap->SetParameter(1,-0.000482148);
86  fbEtaEndcap->SetParameter(1,-0.554552);
87  faEndcap->SetParameter(2,0.808071);
88  fbEndcap->SetParameter(2,0.155416);
89  fcEndcap->SetParameter(2,0.227162);
90  faEtaEndcap->SetParameter(2,209.466);
91  fbEtaEndcap->SetParameter(2,1.71188);
92  faEndcap->SetParameter(3,214.039);
93  fbEndcap->SetParameter(3,240.379);
94  fcEndcap->SetParameter(3,207.786);
95  fbEtaEndcap->SetParameter(3,0.235834);
96  faEndcap->SetParameter(4,2);
97  fbEndcap->SetParameter(4,1.2);
98  fcEndcap->SetParameter(4,1.32824);
99  fbEtaEndcap->SetParameter(4,-135.431);
100  faEndcap->SetParameter(5,47.2602);
101  fbEndcap->SetParameter(5,78.3083);
102  fcEndcap->SetParameter(5,22.1825);
103 
104 }
105 
106 void
107 PFEnergyCalibration::energyEmHad(double t, double& e, double&h, double eta, double phi) const {
108 
109 
110  // Use calorimetric energy as true energy for neutral particles
111  double tt = t;
112  double ee = e;
113  double hh = h;
114  double a = 1.;
115  double b = 1.;
116  double etaCorrE = 1.;
117  double etaCorrH = 1.;
118  t = min(999.9,max(tt,e+h));
119  if ( t < 1. ) return;
120 
121  // Barrel calibration
122  if ( fabs(eta) < 1.48 ) {
123 
124  // The energy correction
125  a = e>0. ? aBarrel(t) : 1.;
126  b = e>0. ? bBarrel(t) : cBarrel(t);
127  double thresh = e > 0. ? threshE : threshH;
128 
129  // Protection against negative calibration - to be tuned
130  if ( a < -0.25 || b < -0.25 ) {
131  a = 1.;
132  b = 1.;
133  thresh = 0.;
134  }
135 
136  // The new estimate of the true energy
137  t = min(999.9,max(tt, thresh+a*e+b*h));
138 
139  // The angular correction for ECAL hadronic deposits
140  etaCorrE = 1. + aEtaBarrel(t) + bEtaBarrel(t)*fabs(eta)*fabs(eta);
141  etaCorrH = 1.;
142  // etaCorr = 1.;
143  t = max(tt, thresh+etaCorrE*a*e+etaCorrH*b*h);
144 
145  if ( e > 0. && thresh > 0. )
146  e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
147  if ( h > 0. && thresh > 0. )
148  h = threshH + etaCorrH * b * h;
149 
150  /*
151  if ( e < 0. || h < 0. ) {
152  std::cout << "Warning : Energy correction ! " << std::endl
153  << "eta,tt,e,h,a,b = " << eta << " " << tt << " "
154  << ee << "/" << e << " " << hh << "/" << h << " " << a << " " << b << std::endl;
155  }
156 
157  if ( etaCorrE > 2. || etaCorrE < 0.5 ||
158  etaCorrH > 2. || etaCorrH < 0.5 )
159  std::cout << "Warning : Angular correction ! " << std::endl
160  << "etaCorrE,etaCorrH,eta,t = "
161  << etaCorrE << " " << etaCorrH << " " << eta << " " << t << std::endl;
162  */
163 
164  // Endcap calibration
165  } else {
166 
167  // The energy correction
168  a = e>0. ? aEndcap(t) : 1.;
169  b = e>0. ? bEndcap(t) : cEndcap(t);
170  double thresh = e > 0. ? threshE : threshH;
171 
172  if ( a < -0.25 || b < -0.25 ) {
173  a = 1.;
174  b = 1.;
175  thresh = 0.;
176  }
177 
178  // The new estimate of the true energy
179  t = min(999.9,max(tt, thresh+a*e+b*h));
180 
181  // The angular correction
182  double dEta = fabs ( fabs(eta) - 1.5 );
183  double etaPow = dEta * dEta * dEta * dEta;
184  //etaCorrE = 1. + aEtaEndcap(t) + 0.5*bEtaEndcap(t)*etaPow;
185  etaCorrE = 1. + aEtaEndcap(t) + bEtaEndcap(t)*etaPow;
186  etaCorrH = 1. + aEtaEndcap(t) + bEtaEndcap(t)*etaPow;
187  /*
188  if ( etaCorr > 2. || etaCorr < 0.5 )
189  std::cout << "Warning : Angular correction ! " << std::endl
190  << "etaCorr,eta,t = " << etaCorr << " " << eta << " " << tt
191  << " ee,hh,e,h = " << e << " " << h << " " << a*e << " " << b*h
192  << std::endl;
193  */
194 
195  t = min(999.9,max(tt, thresh + etaCorrE*a*e + etaCorrH*b*h));
196 
197  if ( e > 0. && thresh > 0. )
198  e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
199  if ( h > 0. && thresh > 0. )
200  h = threshH + b * etaCorrH * h;
201 
202 
203  }
204 
205  // Protection
206  if ( e < 0. || h < 0. ) {
207  /*
208  std::cout << "Warning : Energy correction ! " << std::endl
209  << "eta,tt,e,h,a,b = " << eta << " " << tt << " "
210  << ee << "/" << e << " " << hh << "/" << h << " " << a << " " << b << std::endl;
211  */
212  // Some protection against crazy calibration
213  if ( e < 0. ) e = ee;
214  if ( h < 0. ) h = hh;
215  }
216 
217  // And that's it !
218 
219 
220 }
221 
222 // The calibration functions
223 double
224 PFEnergyCalibration::aBarrel(double x) const {
225 
226  if ( pfCalibrations ) {
227 
228  BinningPointByMap point;
229  point.insert(BinningVariables::JetEt, x);
230  return pfCalibrations->getResult(PerformanceResult::PFfa_BARREL,point);
231 
232  } else {
233 
234  return faBarrel->Eval(x);
235 
236  }
237 }
238 
239 double
240 PFEnergyCalibration::bBarrel(double x) const {
241 
242  if ( pfCalibrations ) {
243 
244  BinningPointByMap point;
245  point.insert(BinningVariables::JetEt, x);
246  return pfCalibrations->getResult(PerformanceResult::PFfb_BARREL,point);
247 
248  } else {
249 
250  return fbBarrel->Eval(x);
251 
252  }
253 }
254 
255 double
256 PFEnergyCalibration::cBarrel(double x) const {
257 
258  if ( pfCalibrations ) {
259 
260  BinningPointByMap point;
261  point.insert(BinningVariables::JetEt, x);
262  return pfCalibrations->getResult(PerformanceResult::PFfc_BARREL,point);
263 
264  } else {
265 
266  return fcBarrel->Eval(x);
267 
268  }
269 }
270 
271 double
272 PFEnergyCalibration::aEtaBarrel(double x) const {
273 
274  if ( pfCalibrations ) {
275 
276  BinningPointByMap point;
277  point.insert(BinningVariables::JetEt, x);
278  return pfCalibrations->getResult(PerformanceResult::PFfaEta_BARREL,point);
279 
280  } else {
281 
282  return faEtaBarrel->Eval(x);
283 
284  }
285 }
286 
287 double
288 PFEnergyCalibration::bEtaBarrel(double x) const {
289 
290  if ( pfCalibrations ) {
291 
292  BinningPointByMap point;
293  point.insert(BinningVariables::JetEt, x);
294  return pfCalibrations->getResult(PerformanceResult::PFfbEta_BARREL,point);
295 
296  } else {
297 
298  return fbEtaBarrel->Eval(x);
299 
300  }
301 }
302 
303 double
304 PFEnergyCalibration::aEndcap(double x) const {
305 
306  if ( pfCalibrations ) {
307 
308  BinningPointByMap point;
309  point.insert(BinningVariables::JetEt, x);
310  return pfCalibrations->getResult(PerformanceResult::PFfa_ENDCAP,point);
311 
312  } else {
313 
314  return faEndcap->Eval(x);
315 
316  }
317 }
318 
319 double
320 PFEnergyCalibration::bEndcap(double x) const {
321 
322  if ( pfCalibrations ) {
323 
324  BinningPointByMap point;
325  point.insert(BinningVariables::JetEt, x);
326  return pfCalibrations->getResult(PerformanceResult::PFfb_ENDCAP,point);
327 
328  } else {
329 
330  return fbEndcap->Eval(x);
331 
332  }
333 }
334 
335 double
336 PFEnergyCalibration::cEndcap(double x) const {
337 
338  if ( pfCalibrations ) {
339 
340  BinningPointByMap point;
341  point.insert(BinningVariables::JetEt, x);
342  return pfCalibrations->getResult(PerformanceResult::PFfc_ENDCAP,point);
343 
344  } else {
345 
346  return fcEndcap->Eval(x);
347 
348  }
349 }
350 
351 double
352 PFEnergyCalibration::aEtaEndcap(double x) const {
353 
354  if ( pfCalibrations ) {
355 
356  BinningPointByMap point;
357  point.insert(BinningVariables::JetEt, x);
358  return pfCalibrations->getResult(PerformanceResult::PFfaEta_ENDCAP,point);
359 
360  } else {
361 
362  return faEtaEndcap->Eval(x);
363 
364  }
365 }
366 
367 double
368 PFEnergyCalibration::bEtaEndcap(double x) const {
369 
370  if ( pfCalibrations ) {
371 
372  BinningPointByMap point;
373  point.insert(BinningVariables::JetEt, x);
374  return pfCalibrations->getResult(PerformanceResult::PFfbEta_ENDCAP,point);
375 
376  } else {
377 
378  return fbEtaEndcap->Eval(x);
379 
380  }
381 }
382 
383 
384 double
385 PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
386  std::vector<double> &EclustersPS1,
387  std::vector<double> &EclustersPS2,
388  bool crackCorrection ){
389  double eEcal = clusterEcal.energy();
390  //temporaty ugly fix
391  reco::PFCluster myPFCluster=clusterEcal;
392  myPFCluster.calculatePositionREP();
393  double eta = myPFCluster.positionREP().eta();
394  double phi = myPFCluster.positionREP().phi();
395 
396  double ePS1 = 0;
397  double ePS2 = 0;
398 
399  for(unsigned i=0;i<EclustersPS1.size();i++) ePS1 += EclustersPS1[i];
400  for(unsigned i=0;i<EclustersPS2.size();i++) ePS2 += EclustersPS2[i];
401 
402  double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi, crackCorrection);
403  if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
404  return calibrated;
405 }
406 
407 double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
408  std::vector<double> &EclustersPS1,
409  std::vector<double> &EclustersPS2,
410  double& ps1,double& ps2,
411  bool crackCorrection){
412  double eEcal = clusterEcal.energy();
413  //temporaty ugly fix
414  reco::PFCluster myPFCluster=clusterEcal;
415  myPFCluster.calculatePositionREP();
416  double eta = myPFCluster.positionREP().eta();
417  double phi = myPFCluster.positionREP().phi();
418 
419  double ePS1 = 0;
420  double ePS2 = 0;
421 
422  for(unsigned i=0;i<EclustersPS1.size();i++) ePS1 += EclustersPS1[i];
423  for(unsigned i=0;i<EclustersPS2.size();i++) ePS2 += EclustersPS2[i];
424 
425  double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi,ps1,ps2,crackCorrection);
426  if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
427  return calibrated;
428 }
429 
430 
431 std::ostream& operator<<(std::ostream& out,
432  const PFEnergyCalibration& calib) {
433 
434  if(!out ) return out;
435 
436  out<<"PFEnergyCalibration -- "<<endl;
437 
438  if ( calib.pfCalibrations ) {
439 
440  std::cout << "Functions taken from the global tags : " << std::endl;
441 
442  static std::map<std::string, PerformanceResult::ResultType> functType;
443 
444  functType["PFfa_BARREL"] = PerformanceResult::PFfa_BARREL;
445  functType["PFfa_ENDCAP"] = PerformanceResult::PFfa_ENDCAP;
446  functType["PFfb_BARREL"] = PerformanceResult::PFfb_BARREL;
447  functType["PFfb_ENDCAP"] = PerformanceResult::PFfb_ENDCAP;
448  functType["PFfc_BARREL"] = PerformanceResult::PFfc_BARREL;
449  functType["PFfc_ENDCAP"] = PerformanceResult::PFfc_ENDCAP;
450  functType["PFfaEta_BARREL"] = PerformanceResult::PFfaEta_BARREL;
451  functType["PFfaEta_ENDCAP"] = PerformanceResult::PFfaEta_ENDCAP;
452  functType["PFfbEta_BARREL"] = PerformanceResult::PFfbEta_BARREL;
453  functType["PFfbEta_ENDCAP"] = PerformanceResult::PFfbEta_ENDCAP;
454 
455  for(std::map<std::string,PerformanceResult::ResultType>::const_iterator
456  func = functType.begin();
457  func != functType.end();
458  ++func) {
459 
460  cout << "Function: " << func->first << endl;
461  PerformanceResult::ResultType fType = func->second;
462  calib.pfCalibrations->printFormula(fType);
463  }
464 
465  } else {
466 
467  std::cout << "Default calibration functions : " << std::endl;
468 
469  calib.faBarrel->Print();
470  calib.fbBarrel->Print();
471  calib.fcBarrel->Print();
472  calib.faEtaBarrel->Print();
473  calib.fbEtaBarrel->Print();
474  calib.faEndcap->Print();
475  calib.fbEndcap->Print();
476  calib.fcEndcap->Print();
477  calib.faEtaEndcap->Print();
478  calib.fbEtaEndcap->Print();
479  }
480 
481  return out;
482 }
483 
484 
485 
486 
499 
500 
501 
507 
508 
509 //useful to compute the signed distance to the closest crack in the barrel
510 double
511 PFEnergyCalibration::minimum(double a,double b){
512  if(TMath::Abs(b)<TMath::Abs(a)) a=b;
513  return a;
514 }
515 
516 
517 //compute the unsigned distance to the closest phi-crack in the barrel
518 double
519 PFEnergyCalibration::dCrackPhi(double phi, double eta){
520 
521  static double pi= M_PI;// 3.14159265358979323846;
522 
523  //Location of the 18 phi-cracks
524  static std::vector<double> cPhi;
525  if(cPhi.size()==0)
526  {
527  cPhi.resize(18,0);
528  cPhi[0]=2.97025;
529  for(unsigned i=1;i<=17;++i) cPhi[i]=cPhi[0]-2*i*pi/18;
530  }
531 
532  //Shift of this location if eta<0
533  static double delta_cPhi=0.00638;
534 
535  double m; //the result
536 
537  //the location is shifted
538  if(eta<0) phi +=delta_cPhi;
539 
540  if (phi>=-pi && phi<=pi){
541 
542  //the problem of the extrema
543  if (phi<cPhi[17] || phi>=cPhi[0]){
544  if (phi<0) phi+= 2*pi;
545  m = minimum(phi -cPhi[0],phi-cPhi[17]-2*pi);
546  }
547 
548  //between these extrema...
549  else{
550  bool OK = false;
551  unsigned i=16;
552  while(!OK){
553  if (phi<cPhi[i]){
554  m=minimum(phi-cPhi[i+1],phi-cPhi[i]);
555  OK=true;
556  }
557  else i-=1;
558  }
559  }
560  }
561  else{
562  m=0.; //if there is a problem, we assum that we are in a crack
563  std::cout<<"Problem in dminphi"<<std::endl;
564  }
565  if(eta<0) m=-m; //because of the disymetry
566  return m;
567 }
568 
569 // corrects the effect of phi-cracks
570 double
571 PFEnergyCalibration::CorrPhi(double phi, double eta) {
572 
573  // we use 3 gaussians to correct the phi-cracks effect
574  static double p1= 5.59379e-01;
575  static double p2= -1.26607e-03;
576  static double p3= 9.61133e-04;
577 
578  static double p4= 1.81691e-01;
579  static double p5= -4.97535e-03;
580  static double p6= 1.31006e-03;
581 
582  static double p7= 1.38498e-01;
583  static double p8= 1.18599e-04;
584  static double p9= 2.01858e-03;
585 
586 
587  double dminphi = dCrackPhi(phi,eta);
588 
589  double result = (1+p1*TMath::Gaus(dminphi,p2,p3)+p4*TMath::Gaus(dminphi,p5,p6)+p7*TMath::Gaus(dminphi,p8,p9));
590 
591  return result;
592 }
593 
594 
595 // corrects the effect of |eta|-cracks
596 double
597 PFEnergyCalibration::CorrEta(double eta){
598 
599  // we use a gaussian with a screwness for each of the 5 |eta|-cracks
600  static std::vector<double> a; //amplitude
601  static std::vector<double> m; //mean
602  static std::vector<double> s; //sigma
603  static std::vector<double> sa; // screwness amplitude
604  static std::vector<double> ss; // screwness sigma
605 
606  if(a.size()==0)
607  {
608  a.push_back(6.13349e-01) ;a.push_back(5.08146e-01) ;a.push_back(4.44480e-01) ;a.push_back(3.3487e-01) ;a.push_back(7.65627e-01) ;
609  m.push_back(-1.79514e-02);m.push_back(4.44747e-01) ;m.push_back(7.92824e-01) ;m.push_back(1.14090e+00) ;m.push_back(1.47464e+00) ;
610  s.push_back(7.92382e-03) ;s.push_back(3.06028e-03) ;s.push_back(3.36139e-03) ;s.push_back(3.94521e-03) ;s.push_back(8.63950e-04) ;
611  sa.push_back(1.27228e+01);sa.push_back(3.81517e-02) ;sa.push_back(1.63507e-01);sa.push_back(-6.56480e-02);sa.push_back(1.87160e-01);
612  ss.push_back(5.48753e-02);ss.push_back(-1.00223e-02);ss.push_back(2.22866e-03);ss.push_back(4.26288e-04) ;ss.push_back(2.67937e-03);
613  }
614  double result = 1;
615 
616  for(unsigned i=0;i<=4;i++) result+=a[i]*TMath::Gaus(eta,m[i],s[i])*(1+sa[i]*TMath::Sign(1.,eta-m[i])*TMath::Exp(-TMath::Abs(eta-m[i])/ss[i]));
617 
618  return result;
619 }
620 
621 
622 //corrects the global behaviour in the barrel
623 double
624 PFEnergyCalibration::CorrBarrel(double E, double eta) {
625 
626  //Energy dependency
627  static double p0=1.00000e+00;
628  static double p1=3.27753e+01;
629  static double p2=2.28552e-02;
630  static double p3=3.06139e+00;
631  static double p4=2.25135e-01;
632  static double p5=1.47824e+00;
633  static double p6=1.09e-02;
634  static double p7=4.19343e+01;
635 
636  //Eta dependency
637  static double p8=2.705593e-03;
638 
639  double result = (p0+1/(p1+p2*TMath::Power(E,p3))+p4*TMath::Exp(-E/p5)+p6*TMath::Exp(-E*E/(p7*p7)))*(1+p8*eta*eta);
640 
641  return result;
642 }
643 
644 
645 
654 
655 
656 //Alpha, Beta, Gamma give the weight of each sub-detector (PS layer1, PS layer2 and Ecal) in the areas of the endcaps where there is a PS
657 // Etot = Beta*eEcal + Gamma*(ePS1 + Alpha*ePS2)
658 
659 double
660 PFEnergyCalibration::Alpha(double eta) {
661 
662  //Energy dependency
663  static double p0 = 5.97621e-01;
664 
665  //Eta dependency
666  static double p1 =-1.86407e-01;
667  static double p2 = 3.85197e-01;
668 
669  //so that <feta()> = 1
670  static double norm = (p1+p2*(2.6+1.656)/2);
671 
672  double result = p0*(p1+p2*eta)/norm;
673 
674  return result;
675 }
676 
677 double
678 PFEnergyCalibration::Beta(double E, double eta) {
679 
680  //Energy dependency
681  static double p0 = 0.032;
682  static double p1 = 9.70394e-02;
683  static double p2 = 2.23072e+01;
684  static double p3 = 100;
685 
686  //Eta dependency
687  static double p4 = 1.02496e+00 ;
688  static double p5 = -4.40176e-03 ;
689 
690  //so that <feta()> = 1
691  static double norm = (p4+p5*(2.6+1.656)/2);
692 
693  double result = (1.0012+p0*TMath::Exp(-E/p3)+p1*TMath::Exp(-E/p2))*(p4+p5*eta)/norm;
694  return result;
695 }
696 
697 
698 double
699 PFEnergyCalibration::Gamma(double etaEcal) {
700 
701  //Energy dependency
702  static double p0 = 2.49752e-02;
703 
704  //Eta dependency
705  static double p1 = 6.48816e-02;
706  static double p2 = -1.59517e-02;
707 
708  //so that <feta()> = 1
709  static double norm = (p1+p2*(2.6+1.656)/2);
710 
711  double result = p0*(p1+p2*etaEcal)/norm;
712 
713  return result;
714 }
715 
716 
717 
723 
724 
725 // returns the corrected energy in the barrel (0,1.48)
726 // Global Behaviour, phi and eta cracks are taken into account
727 double
728 PFEnergyCalibration::EcorrBarrel(double E, double eta, double phi,
729  bool crackCorrection ){
730 
731  // double result = E*CorrBarrel(E,eta)*CorrEta(eta)*CorrPhi(phi,eta);
732  double correction = crackCorrection ? std::max(CorrEta(eta),CorrPhi(phi,eta)) : 1.;
733  double result = E * CorrBarrel(E,eta) * correction;
734 
735  return result;
736 }
737 
738 
739 // returns the corrected energy in the area between the barrel and the PS (1.48,1.65)
740 double
741 PFEnergyCalibration::EcorrZoneBeforePS(double E, double eta){
742 
743  //Energy dependency
744  static double p0 =1;
745  static double p1 =0.18;
746  static double p2 =8.;
747 
748  //Eta dependency
749  static double p3 =0.3;
750  static double p4 =1.11;
751  static double p5 =0.025;
752  static double p6 =1.49;
753  static double p7 =0.6;
754 
755  //so that <feta()> = 1
756  static double norm = 1.21;
757 
758  double result = E*(p0+p1*TMath::Exp(-E/p2))*(p3+p4*TMath::Gaus(eta,p6,p5)+p7*eta)/norm;
759 
760  return result;
761 }
762 
763 
764 // returns the corrected energy in the PS (1.65,2.6)
765 // only when (ePS1>0)||(ePS2>0)
766 double
767 PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal) {
768 
769  // gives the good weights to each subdetector
770  double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+Gamma(etaEcal)*(ePS1+Alpha(etaEcal)*ePS2)/9e-5 ;
771 
772  //Correction of the residual energy dependency
773  static double p0 = 1.00;
774  static double p1 = 2.18;
775  static double p2 =1.94;
776  static double p3 =4.13;
777  static double p4 =1.127;
778 
779  double result = E*(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
780 
781  return result;
782 }
783 
784 // returns the corrected energy in the PS (1.65,2.6)
785 // only when (ePS1>0)||(ePS2>0)
786 double
787 PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal,double & outputPS1, double & outputPS2) {
788 
789  // gives the good weights to each subdetector
790  double gammaprime=Gamma(etaEcal)/9e-5;
791  outputPS1=gammaprime*ePS1;
792  outputPS2=gammaprime*Alpha(etaEcal)*ePS2;
793  double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+outputPS1+outputPS2;
794 
795  //Correction of the residual energy dependency
796  static double p0 = 1.00;
797  static double p1 = 2.18;
798  static double p2 =1.94;
799  static double p3 =4.13;
800  static double p4 =1.127;
801 
802  double corrfac=(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
803  outputPS1*=corrfac;
804  outputPS2*=corrfac;
805  double result = E*corrfac;
806 
807  return result;
808 }
809 
810 
811 // returns the corrected energy in the PS (1.65,2.6)
812 // only when (ePS1=0)&&(ePS2=0)
813 double
814 PFEnergyCalibration::EcorrPS_ePSNil(double eEcal,double eta){
815 
816  //Energy dependency
817  static double p0= 1.02;
818  static double p1= 0.165;
819  static double p2= 6.5 ;
820  static double p3= 2.1 ;
821 
822  //Eta dependency
823  static double p4 = 1.02496e+00 ;
824  static double p5 = -4.40176e-03 ;
825 
826  //so that <feta()> = 1
827  static double norm = (p4+p5*(2.6+1.656)/2);
828 
829  double result = eEcal*(p0+p1*TMath::Exp(-TMath::Abs(eEcal-p3)/p2))*(p4+p5*eta)/norm;
830 
831  return result;
832 }
833 
834 
835 // returns the corrected energy in the area between the end of the PS and the end of the endcap (2.6,2.98)
836 double
837 PFEnergyCalibration::EcorrZoneAfterPS(double E, double eta){
838 
839  //Energy dependency
840  static double p0 =1;
841  static double p1 = 0.058;
842  static double p2 =12.5;
843  static double p3 =-1.05444e+00;
844  static double p4 =-5.39557e+00;
845  static double p5 =8.38444e+00;
846  static double p6 = 6.10998e-01 ;
847 
848  //Eta dependency
849  static double p7 =1.06161e+00;
850  static double p8 = 0.41;
851  static double p9 =2.918;
852  static double p10 =0.0181;
853  static double p11= 2.05;
854  static double p12 =2.99;
855  static double p13=0.0287;
856 
857  //so that <feta()> = 1
858  static double norm=1.045;
859 
860  double result = E*(p0+p1*TMath::Exp(-(E-p3)/p2)+1/(p4+p5*TMath::Power(E,p6)))*(p7+p8*TMath::Gaus(eta,p9,p10)+p11*TMath::Gaus(eta,p12,p13))/norm;
861  return result;
862 }
863 
864 
865 
866 
867 // returns the corrected energy everywhere
868 // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
869 double
870 PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,
871  double eta,double phi,
872  bool crackCorrection ) {
873 
874  static double endBarrel=1.48;
875  static double beginingPS=1.65;
876  static double endPS=2.6;
877  static double endEndCap=2.98;
878 
879  double result=0;
880 
881  eta=TMath::Abs(eta);
882 
883  if(eEcal>0){
884  if(eta <= endBarrel) result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
885  else if(eta <= beginingPS) result = EcorrZoneBeforePS(eEcal,eta);
886  else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
887  else if(eta < endPS) result = EcorrPS(eEcal,ePS1,ePS2,eta);
888  else if(eta < endEndCap) result = EcorrZoneAfterPS(eEcal,eta);
889  else result =eEcal;
890  }
891  else result = eEcal;// useful if eEcal=0 or eta>2.98
892  //protection
893  if(result<eEcal) result=eEcal;
894  return result;
895 }
896 
897 // returns the corrected energy everywhere
898 // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
899 double
900 PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,double eta,double phi,double& ps1,double&ps2,bool crackCorrection) {
901 
902  static double endBarrel=1.48;
903  static double beginingPS=1.65;
904  static double endPS=2.6;
905  static double endEndCap=2.98;
906 
907  double result=0;
908 
909  eta=TMath::Abs(eta);
910 
911  if(eEcal>0){
912  if(eta <= endBarrel) result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
913  else if(eta <= beginingPS) result = EcorrZoneBeforePS(eEcal,eta);
914  else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
915  else if(eta < endPS) result = EcorrPS(eEcal,ePS1,ePS2,eta,ps1,ps2);
916  else if(eta < endEndCap) result = EcorrZoneAfterPS(eEcal,eta);
917  else result =eEcal;
918  }
919  else result = eEcal;// useful if eEcal=0 or eta>2.98
920  // protection
921  if(result<eEcal) result=eEcal;
922  return result;
923 }
PerformancePayloadFromTFormula::getResult
float getResult(PerformanceResult::ResultType, BinningPointByMap) const
Definition: PerformancePayloadFromTFormula.cc:8
i
int i
Definition: DBlmapReader.cc:9
PFEnergyCalibration::Beta
double Beta(double E, double eta)
Definition: PFEnergyCalibration.cc:678
PFEnergyCalibration::EcorrZoneAfterPS
double EcorrZoneAfterPS(double E, double eta)
Definition: PFEnergyCalibration.cc:837
PFEnergyCalibration::Alpha
double Alpha(double eta)
Definition: PFEnergyCalibration.cc:660
PFEnergyCalibration::pfCalibrations
const PerformancePayloadFromTFormula * pfCalibrations
Definition: PFEnergyCalibration.h:76
PFEnergyCalibration::bBarrel
double bBarrel(double x) const
Definition: PFEnergyCalibration.cc:240
reco::PFCluster
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:42
PFEnergyCalibration::aEndcap
double aEndcap(double x) const
Definition: PFEnergyCalibration.cc:304
PFEnergyCalibration::bEndcap
double bEndcap(double x) const
Definition: PFEnergyCalibration.cc:320
PFEnergyCalibration::EcorrBarrel
double EcorrBarrel(double E, double eta, double phi, bool crackCorrection=true)
Definition: PFEnergyCalibration.cc:728
PFEnergyCalibration::EcorrZoneBeforePS
double EcorrZoneBeforePS(double E, double eta)
Definition: PFEnergyCalibration.cc:741
PFEnergyCalibration::EcorrPS_ePSNil
double EcorrPS_ePSNil(double eEcal, double eta)
Definition: PFEnergyCalibration.cc:814
PFEnergyCalibration::CorrEta
double CorrEta(double eta)
Definition: PFEnergyCalibration.cc:597
min
#define min(a, b)
Definition: mlp_lapack.h:161
eta
T eta() const
Definition: Basic3DVectorLD.h:168
operator<<
std::ostream & operator<<(std::ostream &out, const ALILine &li)
Definition: ALILine.cc:187
PFEnergyCalibration::energyEmHad
void energyEmHad(double t, double &e, double &h, double eta, double phi) const
Definition: PFEnergyCalibration.cc:107
PFEnergyCalibration::Gamma
double Gamma(double etaEcal)
Definition: PFEnergyCalibration.cc:699
calib
MVATrainerComputer * calib
Definition: MVATrainer.cc:64
PFEnergyCalibration::CorrPhi
double CorrPhi(double phi, double eta)
Definition: PFEnergyCalibration.cc:571
PFEnergyCalibration::CorrBarrel
double CorrBarrel(double E, double eta)
Definition: PFEnergyCalibration.cc:624
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
reco::PFCluster::calculatePositionREP
void calculatePositionREP()
computes posrep_ once and for all
Definition: PFCluster.h:79
PFEnergyCalibration::bEtaEndcap
double bEtaEndcap(double x) const
Definition: PFEnergyCalibration.cc:368
p4
double p4[4]
Definition: TauolaWrapper.h:92
lumiNorm.norm
tuple norm
Definition: lumiNorm.py:78
PFEnergyCalibration::cEndcap
double cEndcap(double x) const
Definition: PFEnergyCalibration.cc:336
query.result
tuple result
Definition: query.py:137
reco::PFCluster::positionREP
const REPPoint & positionREP() const
cluster position: rho, eta, phi
Definition: PFCluster.h:76
p2
double p2[4]
Definition: TauolaWrapper.h:90
PFEnergyCalibration::dCrackPhi
double dCrackPhi(double phi, double eta)
Definition: PFEnergyCalibration.cc:519
reco::PFCluster::energy
double energy() const
cluster energy
Definition: PFCluster.h:73
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
M_PI
#define M_PI
Definition: BFit3D.cc:3
BinningPointByMap::insert
bool insert(BinningVariables::BinningVariablesType, float)
Definition: BinningPointByMap.cc:3
b
double b
Definition: hdecay.h:120
PFEnergyCalibration::EcorrPS
double EcorrPS(double eEcal, double ePS1, double ePS2, double etaEcal)
Definition: PFEnergyCalibration.cc:767
PFEnergyCalibration::aEtaBarrel
double aEtaBarrel(double x) const
Definition: PFEnergyCalibration.cc:272
p1
double p1[4]
Definition: TauolaWrapper.h:89
PFEnergyCalibration::Ecorr
double Ecorr(double eEcal, double ePS1, double ePS2, double eta, double phi, bool crackCorrection=true)
Definition: PFEnergyCalibration.cc:870
a
double a
Definition: hdecay.h:121
GOODCOLL_filter_cfg.thresh
tuple thresh
Definition: GOODCOLL_filter_cfg.py:74
PFEnergyCalibration::aEtaEndcap
double aEtaEndcap(double x) const
Definition: PFEnergyCalibration.cc:352
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:41
PerformancePayloadFromTFormula::printFormula
void printFormula(PerformanceResult::ResultType res) const
Definition: PerformancePayloadFromTFormula.cc:69
PFEnergyCalibration::bEtaBarrel
double bEtaBarrel(double x) const
Definition: PFEnergyCalibration.cc:288
asciidump.s
string s
Definition: asciidump.py:422
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
x
Definition: DDAxes.h:10
PFEnergyCalibration::minimum
double minimum(double a, double b)
Definition: PFEnergyCalibration.cc:511
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
PFEnergyCalibration::aBarrel
double aBarrel(double x) const
Definition: PFEnergyCalibration.cc:224
PFEnergyCalibration::cBarrel
double cBarrel(double x) const
Definition: PFEnergyCalibration.cc:256
point
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
PFEnergyCalibration::energyEm
double energyEm(const reco::PFCluster &clusterEcal, std::vector< double > &EclustersPS1, std::vector< double > &EclustersPS2, bool crackCorrection=true)
Definition: PFEnergyCalibration.cc:385
p3
double p3[4]
Definition: TauolaWrapper.h:91
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5