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 #include <algorithm>
9 
10 using namespace std;
11 
12 PFEnergyCalibration::PFEnergyCalibration() : pfCalibrations(0)
13 {
14  initializeCalibrationFunctions();
15 }
16 
17 PFEnergyCalibration::~PFEnergyCalibration()
18 {
19 
20  delete faBarrel;
21  delete fbBarrel;
22  delete fcBarrel;
23  delete faEtaBarrel;
24  delete fbEtaBarrel;
25  delete faEndcap;
26  delete fbEndcap;
27  delete fcEndcap;
28  delete faEtaEndcap;
29  delete fbEtaEndcap;
30 
31 }
32 
33 void
34 PFEnergyCalibration::initializeCalibrationFunctions() {
35 
36  // NEW NEW with HCAL pre-calibration
37 
38  threshE = 3.5;
39  threshH = 2.5;
40 
41  // Barrel (fit made with |eta| < 1.2)
42  faBarrel = new TF1("faBarrel","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
43  faBarrel->SetParameter(0,-13.8672);
44  faBarrel->SetParameter(1,14.9688);
45  faBarrel->SetParameter(2,3.42084);
46  faBarrel->SetParameter(3,1.02649);
47  faBarrel->SetParameter(4,-0.0197432);
48  faBarrel->SetParameter(5,6.776e-16);
49  faBarrel->SetParameter(6,-1.3274);
50  faBarrel->SetParameter(7,-7.14684);
51  fbBarrel = new TF1("fbBarrel","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
52  fbBarrel->SetParameter(0,1.70114);
53  fbBarrel->SetParameter(1,0.404676);
54  fbBarrel->SetParameter(2,-3.88962);
55  fbBarrel->SetParameter(3,1.2109e+06);
56  fbBarrel->SetParameter(4,0.970741);
57  fbBarrel->SetParameter(5,0.0527482);
58  fbBarrel->SetParameter(6,2.60552);
59  fbBarrel->SetParameter(7,-0.8956);
60  fcBarrel = new TF1("fcBarrel","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
61  fcBarrel->SetParameter(0,1.91262);
62  fcBarrel->SetParameter(1,0.631347);
63  fcBarrel->SetParameter(2,-4.11001);
64  fcBarrel->SetParameter(3,58.2297);
65  fcBarrel->SetParameter(4,1.10332);
66  fcBarrel->SetParameter(5,0.0372895);
67  fcBarrel->SetParameter(6,0.86592);
68  fcBarrel->SetParameter(7,-1.33708);
69  faEtaBarrel = new TF1("faEtaBarrel","[0]+[1]*exp(-x/[2])",1.,1000.);
70  faEtaBarrel->SetParameter(0,0.0114108);
71  faEtaBarrel->SetParameter(1,-0.0394225);
72  faEtaBarrel->SetParameter(2,177.569);
73  fbEtaBarrel = new TF1("fbEtaBarrel","[0]+[1]*exp(-x/[2])",1.,1000.);
74  fbEtaBarrel->SetParameter(0,0.049778);
75  fbEtaBarrel->SetParameter(1,0.0971732);
76  fbEtaBarrel->SetParameter(2,136.546);
77 
78  // End-caps (fit made with eta
79  faEndcap = new TF1("faEndcap","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
80  faEndcap->SetParameter(0,0.930193);
81  faEndcap->SetParameter(1,11.9536);
82  faEndcap->SetParameter(2,-30.0337);
83  faEndcap->SetParameter(3,0.76133);
84  faEndcap->SetParameter(4,0.0776373);
85  faEndcap->SetParameter(5,7.3809e-10);
86  faEndcap->SetParameter(6,0.158734);
87  faEndcap->SetParameter(7,-6.92163);
88  fbEndcap = new TF1("fbEndcap","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
89  fbEndcap->SetParameter(0,-0.436687);
90  fbEndcap->SetParameter(1,2.73698);
91  fbEndcap->SetParameter(2,-3.1509);
92  fbEndcap->SetParameter(3,1.20536);
93  fbEndcap->SetParameter(4,-1.39685);
94  fbEndcap->SetParameter(5,0.0180331);
95  fbEndcap->SetParameter(6,0.270058);
96  fbEndcap->SetParameter(7,-2.30372);
97  fcEndcap = new TF1("fcEndcap","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
98  fcEndcap->SetParameter(0,1.13795);
99  fcEndcap->SetParameter(1,1.21698);
100  fcEndcap->SetParameter(2,-3.81192);
101  fcEndcap->SetParameter(3,115.409);
102  fcEndcap->SetParameter(4,0.673456);
103  fcEndcap->SetParameter(5,0.217077);
104  fcEndcap->SetParameter(6,1.95596);
105  fcEndcap->SetParameter(7,-0.252215);
106  faEtaEndcap = new TF1("faEtaEndcap","[0]+[1]*exp(-x/[2])",1.,1000.);
107  faEtaEndcap->SetParameter(0,0.00706648);
108  faEtaEndcap->SetParameter(1,-0.0279056);
109  faEtaEndcap->SetParameter(2,45.9808);
110  fbEtaEndcap = new TF1("fbEtaEndcap","[0]+[1]*exp(-x/[2])",1.,1000.);
111  fbEtaEndcap->SetParameter(0,0.000368865);
112  fbEtaEndcap->SetParameter(1,0.309993);
113  fbEtaEndcap->SetParameter(2,29.5954);
114 
115 
116 }
117 
118 void
119 PFEnergyCalibration::energyEmHad(double t, double& e, double&h, double eta, double phi) const {
120 
121 
122  // Use calorimetric energy as true energy for neutral particles
123  double tt = t;
124  double ee = e;
125  double hh = h;
126  double a = 1.;
127  double b = 1.;
128  double etaCorrE = 1.;
129  double etaCorrH = 1.;
130  t = min(999.9,max(tt,e+h));
131  if ( t < 1. ) return;
132 
133  // Barrel calibration
134  if ( std::abs(eta) < 1.48 ) {
135 
136  // The energy correction
137  a = e>0. ? aBarrel(t) : 1.;
138  b = e>0. ? bBarrel(t) : cBarrel(t);
139  double thresh = e > 0. ? threshE : threshH;
140 
141  // Protection against negative calibration - to be tuned
142  if ( a < -0.25 || b < -0.25 ) {
143  a = 1.;
144  b = 1.;
145  thresh = 0.;
146  }
147 
148  // The new estimate of the true energy
149  t = min(999.9,max(tt, thresh+a*e+b*h));
150 
151  // The angular correction for ECAL hadronic deposits
152  etaCorrE = 1. + aEtaBarrel(t) + 1.3*bEtaBarrel(t)*std::abs(eta)*std::abs(eta);
153  etaCorrH = 1.;
154  // etaCorr = 1.;
155  t = max(tt, thresh+etaCorrE*a*e+etaCorrH*b*h);
156 
157  if ( e > 0. && thresh > 0. )
158  e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
159  if ( h > 0. && thresh > 0. )
160  h = threshH + etaCorrH * b * h;
161 
162  // Endcap calibration
163  } else {
164 
165  // The energy correction
166  a = e>0. ? aEndcap(t) : 1.;
167  b = e>0. ? bEndcap(t) : cEndcap(t);
168  double thresh = e > 0. ? threshE : threshH;
169 
170  if ( a < -0.25 || b < -0.25 ) {
171  a = 1.;
172  b = 1.;
173  thresh = 0.;
174  }
175 
176  // The new estimate of the true energy
177  t = min(999.9,max(tt, thresh+a*e+b*h));
178 
179  // The angular correction
180  //MM: only necessary for EH-hadrons only. 1.3 factor just helps the parametrization
181  double dEta = std::abs( std::abs(eta) - 1.5 );
182  double etaPow = dEta * dEta * dEta * dEta;
183  etaCorrE = 1. + aEtaEndcap(t) + 1.3*bEtaEndcap(t)*etaPow;
184  etaCorrH = 1.;
185 
186  t = min(999.9,max(tt, thresh + etaCorrE*a*e + etaCorrH*b*h));
187 
188  if ( e > 0. && thresh > 0. )
189  e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
190  if ( h > 0. && thresh > 0. )
191  h = threshH + b * etaCorrH * h;
192 
193 
194  }
195 
196  // Protection
197  if ( e < 0. || h < 0. ) {
198 
199  // Some protection against crazy calibration
200  if ( e < 0. ) e = ee;
201  if ( h < 0. ) h = hh;
202  }
203 
204  // And that's it !
205 
206 
207 }
208 
209 // The calibration functions
210 double
211 PFEnergyCalibration::aBarrel(double x) const {
212 
213  if ( pfCalibrations ) {
214 
215  BinningPointByMap point;
216  point.insert(BinningVariables::JetEt, x);
217  return pfCalibrations->getResult(PerformanceResult::PFfa_BARREL,point);
218 
219  } else {
220 
221  return faBarrel->Eval(x);
222 
223  }
224 }
225 
226 double
227 PFEnergyCalibration::bBarrel(double x) const {
228 
229  if ( pfCalibrations ) {
230 
231  BinningPointByMap point;
232  point.insert(BinningVariables::JetEt, x);
233  return pfCalibrations->getResult(PerformanceResult::PFfb_BARREL,point);
234 
235  } else {
236 
237  return fbBarrel->Eval(x);
238 
239  }
240 }
241 
242 double
243 PFEnergyCalibration::cBarrel(double x) const {
244 
245  if ( pfCalibrations ) {
246 
247  BinningPointByMap point;
248  point.insert(BinningVariables::JetEt, x);
249  return pfCalibrations->getResult(PerformanceResult::PFfc_BARREL,point);
250 
251  } else {
252 
253  return fcBarrel->Eval(x);
254 
255  }
256 }
257 
258 double
259 PFEnergyCalibration::aEtaBarrel(double x) const {
260 
261  if ( pfCalibrations ) {
262 
263  BinningPointByMap point;
264  point.insert(BinningVariables::JetEt, x);
265  return pfCalibrations->getResult(PerformanceResult::PFfaEta_BARREL,point);
266 
267  } else {
268 
269  return faEtaBarrel->Eval(x);
270 
271  }
272 }
273 
274 double
275 PFEnergyCalibration::bEtaBarrel(double x) const {
276 
277  if ( pfCalibrations ) {
278 
279  BinningPointByMap point;
280  point.insert(BinningVariables::JetEt, x);
281  return pfCalibrations->getResult(PerformanceResult::PFfbEta_BARREL,point);
282 
283  } else {
284 
285  return fbEtaBarrel->Eval(x);
286 
287  }
288 }
289 
290 double
291 PFEnergyCalibration::aEndcap(double x) const {
292 
293  if ( pfCalibrations ) {
294 
295  BinningPointByMap point;
296  point.insert(BinningVariables::JetEt, x);
297  return pfCalibrations->getResult(PerformanceResult::PFfa_ENDCAP,point);
298 
299  } else {
300 
301  return faEndcap->Eval(x);
302 
303  }
304 }
305 
306 double
307 PFEnergyCalibration::bEndcap(double x) const {
308 
309  if ( pfCalibrations ) {
310 
311  BinningPointByMap point;
312  point.insert(BinningVariables::JetEt, x);
313  return pfCalibrations->getResult(PerformanceResult::PFfb_ENDCAP,point);
314 
315  } else {
316 
317  return fbEndcap->Eval(x);
318 
319  }
320 }
321 
322 double
323 PFEnergyCalibration::cEndcap(double x) const {
324 
325  if ( pfCalibrations ) {
326 
327  BinningPointByMap point;
328  point.insert(BinningVariables::JetEt, x);
329  return pfCalibrations->getResult(PerformanceResult::PFfc_ENDCAP,point);
330 
331  } else {
332 
333  return fcEndcap->Eval(x);
334 
335  }
336 }
337 
338 double
339 PFEnergyCalibration::aEtaEndcap(double x) const {
340 
341  if ( pfCalibrations ) {
342 
343  BinningPointByMap point;
344  point.insert(BinningVariables::JetEt, x);
345  return pfCalibrations->getResult(PerformanceResult::PFfaEta_ENDCAP,point);
346 
347  } else {
348 
349  return faEtaEndcap->Eval(x);
350 
351  }
352 }
353 
354 double
355 PFEnergyCalibration::bEtaEndcap(double x) const {
356 
357  if ( pfCalibrations ) {
358 
359  BinningPointByMap point;
360  point.insert(BinningVariables::JetEt, x);
361  return pfCalibrations->getResult(PerformanceResult::PFfbEta_ENDCAP,point);
362 
363  } else {
364 
365  return fbEtaEndcap->Eval(x);
366 
367  }
368 }
369 
370 
371 double
372 PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
373  std::vector<double> &EclustersPS1,
374  std::vector<double> &EclustersPS2,
375  bool crackCorrection ){
376  double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
377  double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
378  return energyEm(clusterEcal, ePS1, ePS2, crackCorrection);
379 }
380 
381 double
382 PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
383  double ePS1,
384  double ePS2,
385  bool crackCorrection ){
386  double eEcal = clusterEcal.energy();
387  //temporaty ugly fix
388  reco::PFCluster myPFCluster=clusterEcal;
389  myPFCluster.calculatePositionREP();
390  double eta = myPFCluster.positionREP().eta();
391  double phi = myPFCluster.positionREP().phi();
392 
393  double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi, crackCorrection);
394  // if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
395  return calibrated;
396 }
397 
398 double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
399  std::vector<double> &EclustersPS1,
400  std::vector<double> &EclustersPS2,
401  double& ps1,double& ps2,
402  bool crackCorrection){
403  double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
404  double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
405  return energyEm(clusterEcal, ePS1, ePS2, ps1, ps2, crackCorrection);
406 }
407 double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
408  double ePS1, double ePS2,
409  double& ps1,double& ps2,
410  bool crackCorrection){
411  double eEcal = clusterEcal.energy();
412  //temporaty ugly fix
413  reco::PFCluster myPFCluster=clusterEcal;
414  myPFCluster.calculatePositionREP();
415  double eta = myPFCluster.positionREP().eta();
416  double phi = myPFCluster.positionREP().phi();
417 
418  double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi,ps1,ps2,crackCorrection);
419  // if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
420  return calibrated;
421 }
422 
423 
424 std::ostream& operator<<(std::ostream& out,
425  const PFEnergyCalibration& calib) {
426 
427  if(!out ) return out;
428 
429  out<<"PFEnergyCalibration -- "<<endl;
430 
431  if ( calib.pfCalibrations ) {
432 
433  std::cout << "Functions taken from the global tags : " << std::endl;
434 
435  static std::map<std::string, PerformanceResult::ResultType> functType;
436 
437  functType["PFfa_BARREL"] = PerformanceResult::PFfa_BARREL;
438  functType["PFfa_ENDCAP"] = PerformanceResult::PFfa_ENDCAP;
439  functType["PFfb_BARREL"] = PerformanceResult::PFfb_BARREL;
440  functType["PFfb_ENDCAP"] = PerformanceResult::PFfb_ENDCAP;
441  functType["PFfc_BARREL"] = PerformanceResult::PFfc_BARREL;
442  functType["PFfc_ENDCAP"] = PerformanceResult::PFfc_ENDCAP;
443  functType["PFfaEta_BARREL"] = PerformanceResult::PFfaEta_BARREL;
444  functType["PFfaEta_ENDCAP"] = PerformanceResult::PFfaEta_ENDCAP;
445  functType["PFfbEta_BARREL"] = PerformanceResult::PFfbEta_BARREL;
446  functType["PFfbEta_ENDCAP"] = PerformanceResult::PFfbEta_ENDCAP;
447 
448  for(std::map<std::string,PerformanceResult::ResultType>::const_iterator
449  func = functType.begin();
450  func != functType.end();
451  ++func) {
452 
453  cout << "Function: " << func->first << endl;
454  PerformanceResult::ResultType fType = func->second;
455  calib.pfCalibrations->printFormula(fType);
456  }
457 
458  } else {
459 
460  std::cout << "Default calibration functions : " << std::endl;
461 
462  calib.faBarrel->Print();
463  calib.fbBarrel->Print();
464  calib.fcBarrel->Print();
465  calib.faEtaBarrel->Print();
466  calib.fbEtaBarrel->Print();
467  calib.faEndcap->Print();
468  calib.fbEndcap->Print();
469  calib.fcEndcap->Print();
470  calib.faEtaEndcap->Print();
471  calib.fbEtaEndcap->Print();
472  }
473 
474  return out;
475 }
476 
477 
478 
479 
492 
493 
494 
500 
501 
502 //useful to compute the signed distance to the closest crack in the barrel
503 double
504 PFEnergyCalibration::minimum(double a,double b){
505  if(TMath::Abs(b)<TMath::Abs(a)) a=b;
506  return a;
507 }
508 
509 namespace {
510  constexpr double pi= M_PI;// 3.14159265358979323846;
511 
512 
513  std::vector<double> fillcPhi() {
514  std::vector<double> retValue;
515  retValue.resize(18,0);
516  retValue[0]=2.97025;
517  for(unsigned i=1;i<=17;++i) retValue[i]=retValue[0]-2*i*pi/18;
518 
519  return retValue;
520  }
521 
522  //Location of the 18 phi-cracks
523  const std::vector<double> cPhi = fillcPhi();
524 }
525 
526 //compute the unsigned distance to the closest phi-crack in the barrel
527 double
528 PFEnergyCalibration::dCrackPhi(double phi, double eta){
529 
530 
531  //Shift of this location if eta<0
532  constexpr double delta_cPhi=0.00638;
533 
534  double m; //the result
535 
536  //the location is shifted
537  if(eta<0) phi +=delta_cPhi;
538 
539  if (phi>=-pi && phi<=pi){
540 
541  //the problem of the extrema
542  if (phi<cPhi[17] || phi>=cPhi[0]){
543  if (phi<0) phi+= 2*pi;
544  m = minimum(phi -cPhi[0],phi-cPhi[17]-2*pi);
545  }
546 
547  //between these extrema...
548  else{
549  bool OK = false;
550  unsigned i=16;
551  while(!OK){
552  if (phi<cPhi[i]){
553  m=minimum(phi-cPhi[i+1],phi-cPhi[i]);
554  OK=true;
555  }
556  else i-=1;
557  }
558  }
559  }
560  else{
561  m=0.; //if there is a problem, we assum that we are in a crack
562  std::cout<<"Problem in dminphi"<<std::endl;
563  }
564  if(eta<0) m=-m; //because of the disymetry
565  return m;
566 }
567 
568 // corrects the effect of phi-cracks
569 double
570 PFEnergyCalibration::CorrPhi(double phi, double eta) {
571 
572  // we use 3 gaussians to correct the phi-cracks effect
573  constexpr double p1= 5.59379e-01;
574  constexpr double p2= -1.26607e-03;
575  constexpr double p3= 9.61133e-04;
576 
577  constexpr double p4= 1.81691e-01;
578  constexpr double p5= -4.97535e-03;
579  constexpr double p6= 1.31006e-03;
580 
581  constexpr double p7= 1.38498e-01;
582  constexpr double p8= 1.18599e-04;
583  constexpr double p9= 2.01858e-03;
584 
585 
586  double dminphi = dCrackPhi(phi,eta);
587 
588  double result = (1+p1*TMath::Gaus(dminphi,p2,p3)+p4*TMath::Gaus(dminphi,p5,p6)+p7*TMath::Gaus(dminphi,p8,p9));
589 
590  return result;
591 }
592 
593 
594 // corrects the effect of |eta|-cracks
595 double
596 PFEnergyCalibration::CorrEta(double eta){
597 
598  // we use a gaussian with a screwness for each of the 5 |eta|-cracks
599  constexpr double a[] = {6.13349e-01, 5.08146e-01, 4.44480e-01, 3.3487e-01, 7.65627e-01}; // amplitude
600  constexpr double m[] = {-1.79514e-02, 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00}; // mean
601  constexpr double s[] = {7.92382e-03, 3.06028e-03, 3.36139e-03, 3.94521e-03, 8.63950e-04}; // sigma
602  constexpr double sa[] = {1.27228e+01, 3.81517e-02, 1.63507e-01, -6.56480e-02, 1.87160e-01}; // screwness amplitude
603  constexpr double ss[] = {5.48753e-02, -1.00223e-02, 2.22866e-03, 4.26288e-04, 2.67937e-03}; // screwness sigma
604  double result = 1;
605 
606  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]));
607 
608  return result;
609 }
610 
611 
612 //corrects the global behaviour in the barrel
613 double
614 PFEnergyCalibration::CorrBarrel(double E, double eta) {
615 
616  //Energy dependency
617  /*
618  //YM Parameters 52XX:
619  constexpr double p0=1.00000e+00;
620  constexpr double p1=3.27753e+01;
621  constexpr double p2=2.28552e-02;
622  constexpr double p3=3.06139e+00;
623  constexpr double p4=2.25135e-01;
624  constexpr double p5=1.47824e+00;
625  constexpr double p6=1.09e-02;
626  constexpr double p7=4.19343e+01;
627  */
628  constexpr double p0 = 0.9944;
629  constexpr double p1 = 9.827;
630  constexpr double p2 = 1.503;
631  constexpr double p3 = 1.196;
632  constexpr double p4 = 0.3349;
633  constexpr double p5 = 0.89;
634  constexpr double p6 = 0.004361;
635  constexpr double p7 = 51.51;
636  //Eta dependency
637  constexpr 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  constexpr double p0 = 5.97621e-01;
664 
665  //Eta dependency
666  constexpr double p1 =-1.86407e-01;
667  constexpr double p2 = 3.85197e-01;
668 
669  //so that <feta()> = 1
670  constexpr 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  constexpr double p0 = 0.032;
682  constexpr double p1 = 9.70394e-02;
683  constexpr double p2 = 2.23072e+01;
684  constexpr double p3 = 100;
685 
686  //Eta dependency
687  constexpr double p4 = 1.02496e+00 ;
688  constexpr double p5 = -4.40176e-03 ;
689 
690  //so that <feta()> = 1
691  constexpr 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  constexpr double p0 = 2.49752e-02;
703 
704  //Eta dependency
705  constexpr double p1 = 6.48816e-02;
706  constexpr double p2 = -1.59517e-02;
707 
708  //so that <feta()> = 1
709  constexpr 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  constexpr double p0 =1;
745  constexpr double p1 =0.18;
746  constexpr double p2 =8.;
747 
748  //Eta dependency
749  constexpr double p3 =0.3;
750  constexpr double p4 =1.11;
751  constexpr double p5 =0.025;
752  constexpr double p6 =1.49;
753  constexpr double p7 =0.6;
754 
755  //so that <feta()> = 1
756  constexpr 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  constexpr double p0 = 1.00;
774  constexpr double p1 = 2.18;
775  constexpr double p2 =1.94;
776  constexpr double p3 =4.13;
777  constexpr 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  constexpr double p0 = 1.00;
797  constexpr double p1 = 2.18;
798  constexpr double p2 =1.94;
799  constexpr double p3 =4.13;
800  constexpr 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  constexpr double p0= 1.02;
818  constexpr double p1= 0.165;
819  constexpr double p2= 6.5 ;
820  constexpr double p3= 2.1 ;
821 
822  //Eta dependency
823  constexpr double p4 = 1.02496e+00 ;
824  constexpr double p5 = -4.40176e-03 ;
825 
826  //so that <feta()> = 1
827  constexpr 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  constexpr double p0 =1;
841  constexpr double p1 = 0.058;
842  constexpr double p2 =12.5;
843  constexpr double p3 =-1.05444e+00;
844  constexpr double p4 =-5.39557e+00;
845  constexpr double p5 =8.38444e+00;
846  constexpr double p6 = 6.10998e-01 ;
847 
848  //Eta dependency
849  constexpr double p7 =1.06161e+00;
850  constexpr double p8 = 0.41;
851  constexpr double p9 =2.918;
852  constexpr double p10 =0.0181;
853  constexpr double p11= 2.05;
854  constexpr double p12 =2.99;
855  constexpr double p13=0.0287;
856 
857  //so that <feta()> = 1
858  constexpr 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  constexpr double endBarrel=1.48;
875  constexpr double beginingPS=1.65;
876  constexpr double endPS=2.6;
877  constexpr 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  constexpr double endBarrel=1.48;
903  constexpr double beginingPS=1.65;
904  constexpr double endPS=2.6;
905  constexpr 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 }
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:82
reco::OK
pair< int, edm::FunctionWithDict > OK
Definition: findMethod.cc:72
PFEnergyCalibration::bBarrel
double bBarrel(double x) const
Definition: PFEnergyCalibration.cc:227
reco::PFCluster
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:47
Sign
T Sign(T A, T B)
Definition: MathUtil.h:54
PFEnergyCalibration::aEndcap
double aEndcap(double x) const
Definition: PFEnergyCalibration.cc:291
PFEnergyCalibration::bEndcap
double bEndcap(double x) const
Definition: PFEnergyCalibration.cc:307
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:596
eta
T eta() const
Definition: Basic3DVectorLD.h:177
operator<<
std::ostream & operator<<(std::ostream &out, const ALILine &li)
Definition: ALILine.cc:187
constexpr
#define constexpr
PFEnergyCalibration::energyEmHad
void energyEmHad(double t, double &e, double &h, double eta, double phi) const
Definition: PFEnergyCalibration.cc:119
pi
const Double_t pi
Definition: trackSplitPlot.h:35
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:570
PFEnergyCalibration::CorrBarrel
double CorrBarrel(double E, double eta)
Definition: PFEnergyCalibration.cc:614
reco::PFCluster::calculatePositionREP
void calculatePositionREP()
computes posrep_ once and for all
Definition: PFCluster.h:93
PFEnergyCalibration::bEtaEndcap
double bEtaEndcap(double x) const
Definition: PFEnergyCalibration.cc:355
p4
double p4[4]
Definition: TauolaWrapper.h:92
PFEnergyCalibration::cEndcap
double cEndcap(double x) const
Definition: PFEnergyCalibration.cc:323
query.result
tuple result
Definition: query.py:137
reco::PFCluster::positionREP
const REPPoint & positionREP() const
cluster position: rho, eta, phi
Definition: PFCluster.h:90
Abs
T Abs(T a)
Definition: MathUtil.h:49
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
min
T min(T a, T b)
Definition: MathUtil.h:58
p2
double p2[4]
Definition: TauolaWrapper.h:90
PFEnergyCalibration::dCrackPhi
double dCrackPhi(double phi, double eta)
Definition: PFEnergyCalibration.cc:528
M_PI
#define M_PI
Definition: readTestVector.cc:25
reco::PFCluster::energy
double energy() const
cluster energy
Definition: PFCluster.h:79
PerformancePayloadFromTFormula::getResult
float getResult(PerformanceResult::ResultType, const BinningPointByMap &) const
Definition: PerformancePayloadFromTFormula.cc:18
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
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:259
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:339
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
PerformancePayloadFromTFormula::printFormula
void printFormula(PerformanceResult::ResultType res) const
Definition: PerformancePayloadFromTFormula.cc:69
PFEnergyCalibration::bEtaBarrel
double bEtaBarrel(double x) const
Definition: PFEnergyCalibration.cc:275
x
Definition: DDAxes.h:10
PFEnergyCalibration::minimum
double minimum(double a, double b)
Definition: PFEnergyCalibration.cc:504
PFEnergyCalibration::aBarrel
double aBarrel(double x) const
Definition: PFEnergyCalibration.cc:211
PFEnergyCalibration::cBarrel
double cBarrel(double x) const
Definition: PFEnergyCalibration.cc:243
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:372
p3
double p3[4]
Definition: TauolaWrapper.h:91
phi
Definition: DDAxes.h:10
Generated for CMSSW Reference Manual by doxygen 1.8.5