CMS 3D CMS Logo

PFEnergyCalibration.cc
Go to the documentation of this file.
3 
5 
6 #include <TMath.h>
7 #include <math.h>
8 #include <vector>
9 #include <TF1.h>
10 #include <map>
11 #include <algorithm>
12 #include <numeric>
13 
14 using namespace std;
15 
16 PFEnergyCalibration::PFEnergyCalibration() : pfCalibrations(0), esEEInterCalib_(0)
17 {
19 }
20 
22 {
23 
24 }
25 
26 void
28 
29  threshE = 3.5;
30  threshH = 2.5;
31 
32  //calibChrisClean.C calibration parameters shubham May 01, 2017
33  faBarrel = std::make_unique<TF1>("faBarrel","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
34  faBarrel->SetParameter(0,-13.9219);
35  faBarrel->SetParameter(1,14.9124);
36  faBarrel->SetParameter(2,5.38578);
37  faBarrel->SetParameter(3,0.861981);
38  faBarrel->SetParameter(4,-0.00759275);
39  faBarrel->SetParameter(5,0.00373563);
40  faBarrel->SetParameter(6,-1.17946);
41  faBarrel->SetParameter(7,-1.69561);
42  fbBarrel = std::make_unique<TF1>("fbBarrel","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
43  fbBarrel->SetParameter(0,2.18802);
44  fbBarrel->SetParameter(1,0.522053);
45  fbBarrel->SetParameter(2,-4.67354);
46  fbBarrel->SetParameter(3,12.109);
47  fbBarrel->SetParameter(4,1.75318);
48  fbBarrel->SetParameter(5,0.187919);
49  fbBarrel->SetParameter(6,-6.26234);
50  fbBarrel->SetParameter(7,-0.607392);
51  fcBarrel = std::make_unique<TF1>("fcBarrel","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
52  fcBarrel->SetParameter(0,1.46854);
53  fcBarrel->SetParameter(1,0.830152);
54  fcBarrel->SetParameter(2,-5.86601);
55  fcBarrel->SetParameter(3,2.08229);
56  fcBarrel->SetParameter(4,0.575016);
57  fcBarrel->SetParameter(5,0.0291232);
58  fcBarrel->SetParameter(6,0.364802);
59  fcBarrel->SetParameter(7,-1.50142);
60  faEtaBarrelEH = std::make_unique<TF1>("faEtaBarrelEH","[0]+[1]*exp(-x/[2])",1.,1000.);
61  faEtaBarrelEH->SetParameter(0,0.0185555);
62  faEtaBarrelEH->SetParameter(1,-0.0470674);
63  faEtaBarrelEH->SetParameter(2,396.959);
64  fbEtaBarrelEH = std::make_unique<TF1>("fbEtaBarrelEH","[0]+[1]*exp(-x/[2])",1.,1000.);
65  fbEtaBarrelEH->SetParameter(0,0.0396458);
66  fbEtaBarrelEH->SetParameter(1,0.114128);
67  fbEtaBarrelEH->SetParameter(2,251.405);
68  faEtaBarrelH = std::make_unique<TF1>("faEtaBarrelH","[0]+[1]*x",1.,1000.);
69  faEtaBarrelH->SetParameter(0,0.00434994);
70  faEtaBarrelH->SetParameter(1,-5.16564e-06);
71  fbEtaBarrelH = std::make_unique<TF1>("fbEtaBarrelH","[0]+[1]*exp(-x/[2])",1.,1000.);
72  fbEtaBarrelH->SetParameter(0,-0.0232604);
73  fbEtaBarrelH->SetParameter(1,0.0937525);
74  fbEtaBarrelH->SetParameter(2,34.9935);
75 
76  faEndcap = std::make_unique<TF1>("faEndcap","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
77  faEndcap->SetParameter(0,0.962468);
78  faEndcap->SetParameter(1,11.9536);
79  faEndcap->SetParameter(2,-27.7088);
80  faEndcap->SetParameter(3,0.755474);
81  faEndcap->SetParameter(4,0.0791012);
82  faEndcap->SetParameter(5,0.0011082);
83  faEndcap->SetParameter(6,0.158734);
84  faEndcap->SetParameter(7,-2.1);
85  fbEndcap = std::make_unique<TF1>("fbEndcap","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
86  fbEndcap->SetParameter(0,-0.43671);
87  fbEndcap->SetParameter(1,2.90096);
88  fbEndcap->SetParameter(2,-5.10099);
89  fbEndcap->SetParameter(3,1.20771);
90  fbEndcap->SetParameter(4,-1.30656);
91  fbEndcap->SetParameter(5,0.0189607);
92  fbEndcap->SetParameter(6,0.270027);
93  fbEndcap->SetParameter(7,-2.30372);
94  fcEndcap = std::make_unique<TF1>("fcEndcap","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",1.,1000.);
95  fcEndcap->SetParameter(0,1.11623);
96  fcEndcap->SetParameter(1,1.66795);
97  fcEndcap->SetParameter(2,-5.49393);
98  fcEndcap->SetParameter(3,7.73734);
99  fcEndcap->SetParameter(4,0.259307);
100  fcEndcap->SetParameter(5,0.0708707);
101  fcEndcap->SetParameter(6,0.950338);
102  fcEndcap->SetParameter(7,-0.922411);
103  faEtaEndcapEH = std::make_unique<TF1>("faEtaEndcapEH","[0]+[1]*exp(-x/[2])",1.,1000.);
104  faEtaEndcapEH->SetParameter(0,384.307);
105  faEtaEndcapEH->SetParameter(1,-384.305);
106  faEtaEndcapEH->SetParameter(2,1.14535e+07);
107  fbEtaEndcapEH = std::make_unique<TF1>("fbEtaEndcapEH","[0]+[1]*exp(-x/[2])",1.,1000.);
108  fbEtaEndcapEH->SetParameter(0,0.0282454);
109  fbEtaEndcapEH->SetParameter(1,-0.107485);
110  fbEtaEndcapEH->SetParameter(2,96.1833);
111  faEtaEndcapH = std::make_unique<TF1>("faEtaEndcapH","[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",1.,1000.);
112  faEtaEndcapH->SetParameter(0,-0.0106029);
113  faEtaEndcapH->SetParameter(1,-0.692207);
114  faEtaEndcapH->SetParameter(2,0.0542991);
115  faEtaEndcapH->SetParameter(3,-0.171435);
116  faEtaEndcapH->SetParameter(4,-61.2277);
117  fbEtaEndcapH = std::make_unique<TF1>("fbEtaEndcapH","[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",1.,1000.);
118  fbEtaEndcapH->SetParameter(0,0.0514927);
119  fbEtaEndcapH->SetParameter(1,-0.291726);
120  fbEtaEndcapH->SetParameter(2,5.32989);
121  fbEtaEndcapH->SetParameter(3,0.309057);
122  fbEtaEndcapH->SetParameter(4,-103.359);
123 }
124 
125 void
126 PFEnergyCalibration::energyEmHad(double t, double& e, double&h, double eta, double phi) const {
127 
128  // Use calorimetric energy as true energy for neutral particles
129  double tt = t;
130  double ee = e;
131  double hh = h;
132  double a = 1.;
133  double b = 1.;
134  double etaCorrE = 1.;
135  double etaCorrH = 1.;
136  auto absEta = std::abs(eta);
137  t = min(999.9,max(tt,e+h));
138  if ( t < 1. ) return;
139 
140  // Barrel calibration
141  if ( absEta < 1.48 ) {
142  // The energy correction
143  a = e>0. ? aBarrel(t) : 1.;
144  b = e>0. ? bBarrel(t) : cBarrel(t);
145  double thresh = e > 0. ? threshE : threshH;
146 
147  // Protection against negative calibration
148  if ( a < -0.25 || b < -0.25 ) {
149  a = 1.;
150  b = 1.;
151  thresh = 0.;
152  }
153 
154  // The new estimate of the true energy
155  t = min(999.9,max(tt, thresh+a*e+b*h));
156 
157  // The angular correction
158  if ( e > 0. && thresh > 0. ) {
159  etaCorrE = 1.0 + aEtaBarrelEH(t) + 1.3*bEtaBarrelEH(t)*absEta*absEta;
160  etaCorrH = 1.0;
161  } else {
162  etaCorrE = 1.0 + aEtaBarrelH(t) + 1.3*bEtaBarrelH(t)*absEta*absEta;
163  etaCorrH = 1.0 + aEtaBarrelH(t) + bEtaBarrelH(t)*absEta*absEta;
164  }
165  if ( e > 0. && thresh > 0. )
166  e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
167  if ( h > 0. && thresh > 0. ) {
168  h = threshH + etaCorrH * b * h;
169  }
170 
171  // Endcap calibration
172  } else {
173  // The energy correction
174  a = e>0. ? aEndcap(t) : 1.;
175  b = e>0. ? bEndcap(t) : cEndcap(t);
176  double thresh = e > 0. ? threshE : threshH;
177 
178  // Protection against negative calibration
179  if ( a < -0.25 || b < -0.25 ) {
180  a = 1.;
181  b = 1.;
182  thresh = 0.;
183  }
184 
185  // The new estimate of the true energy
186  t = min(999.9,max(tt, thresh+a*e+b*h));
187 
188  // The angular correction
189  double dEta = std::abs( absEta - 1.5 );
190  double etaPow = dEta * dEta * dEta * dEta;
191 
192  if ( e > 0. && thresh > 0. ) {
193  etaCorrE = 1. + aEtaEndcapEH(t) + 1.3*bEtaEndcapEH(t)*(0.04 + etaPow);
194  etaCorrH = 1. + aEtaEndcapEH(t) + bEtaEndcapEH(t)*(0.04 + etaPow);
195  } else {
196  etaCorrE = 1. + aEtaEndcapH(t) + 1.3*bEtaEndcapH(t)*(0.04 + etaPow);
197  if(absEta<2.5) {
198  etaCorrH = 1. + aEtaEndcapH(t) + 0.05*bEtaEndcapH(t);
199  }
200  else {
201  etaCorrH = 1. + aEtaEndcapH(t) + bEtaEndcapH(t)*(0.04 + etaPow);
202  }
203  }
204 
205  //t = min(999.9,max(tt, thresh + etaCorrE*a*e + etaCorrH*b*h));
206 
207  if ( e > 0. && thresh > 0. )
208  e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
209  if ( h > 0. && thresh > 0. ) {
210  h = threshH + etaCorrH * b * h;
211  }
212  }
213 
214  // Protection
215  if ( e < 0. || h < 0. ) {
216 
217  // Some protection against crazy calibration
218  if ( e < 0. ) e = ee;
219  if ( h < 0. ) h = hh;
220  }
221 
222  // And that's it !
223 
224 
225 }
226 
227 // The calibration functions
228 double
230 
231  if ( pfCalibrations ) {
235 
236  } else {
237 
238  return faBarrel->Eval(x);
239 
240  }
241 }
242 
243 double
245 
246  if ( pfCalibrations ) {
247 
251 
252  } else {
253 
254  return fbBarrel->Eval(x);
255 
256  }
257 }
258 
259 double
261 
262  if ( pfCalibrations ) {
263 
267 
268  } else {
269 
270  return fcBarrel->Eval(x);
271 
272  }
273 }
274 
275 double
277 
278  if ( pfCalibrations ) {
279 
283 
284  } else {
285 
286  return faEtaBarrelEH->Eval(x);
287 
288  }
289 }
290 
291 double
293 
294  if ( pfCalibrations ) {
295 
299 
300  } else {
301 
302  return fbEtaBarrelEH->Eval(x);
303 
304  }
305 }
306 
307 double
309 
310  if ( pfCalibrations ) {
311 
315 
316  } else {
317 
318  return faEtaBarrelH->Eval(x);
319 
320  }
321 }
322 
323 double
325 
326  if ( pfCalibrations ) {
327 
331 
332  } else {
333 
334  return fbEtaBarrelH->Eval(x);
335 
336  }
337 }
338 
339 double
341 
342  if ( pfCalibrations ) {
343 
347 
348  } else {
349 
350  return faEndcap->Eval(x);
351 
352  }
353 }
354 
355 double
357 
358  if ( pfCalibrations ) {
359 
363 
364  } else {
365 
366  return fbEndcap->Eval(x);
367 
368  }
369 }
370 
371 double
373 
374  if ( pfCalibrations ) {
375 
379 
380  } else {
381 
382  return fcEndcap->Eval(x);
383 
384  }
385 }
386 
387 double
389 
390  if ( pfCalibrations ) {
391 
395 
396  } else {
397 
398  return faEtaEndcapEH->Eval(x);
399 
400  }
401 }
402 
403 double
405 
406  if ( pfCalibrations ) {
407 
411 
412  } else {
413 
414  return fbEtaEndcapEH->Eval(x);
415 
416  }
417 }
418 
419 double
421 
422  if ( pfCalibrations ) {
423 
427 
428  } else {
429 
430  return faEtaEndcapH->Eval(x);
431 
432  }
433 }
434 
435 double
437 
438  if ( pfCalibrations ) {
439 
443 
444  } else {
445 
446  return fbEtaEndcapH->Eval(x);
447 
448  }
449 }
450 
451 
452 double
454  std::vector<double> &EclustersPS1,
455  std::vector<double> &EclustersPS2,
456  bool crackCorrection ) const {
457  double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
458  double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
459  return energyEm(clusterEcal, ePS1, ePS2, crackCorrection);
460 }
461 
462 double
464  double ePS1,
465  double ePS2,
466  bool crackCorrection ) const {
467  double eEcal = clusterEcal.energy();
468  //temporaty ugly fix
469  reco::PFCluster myPFCluster=clusterEcal;
470  myPFCluster.calculatePositionREP();
471  double eta = myPFCluster.positionREP().eta();
472  double phi = myPFCluster.positionREP().phi();
473 
474  double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi, crackCorrection);
475  // if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
476  return calibrated;
477 }
478 
480  std::vector<double> &EclustersPS1,
481  std::vector<double> &EclustersPS2,
482  double& ps1,double& ps2,
483  bool crackCorrection) const {
484  double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
485  double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
486  return energyEm(clusterEcal, ePS1, ePS2, ps1, ps2, crackCorrection);
487 }
489  double ePS1, double ePS2,
490  double& ps1,double& ps2,
491  bool crackCorrection) const {
492  double eEcal = clusterEcal.energy();
493  //temporaty ugly fix
494  reco::PFCluster myPFCluster=clusterEcal;
495  myPFCluster.calculatePositionREP();
496  double eta = myPFCluster.positionREP().eta();
497  double phi = myPFCluster.positionREP().phi();
498 
499  double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi,ps1,ps2,crackCorrection);
500  // if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
501  return calibrated;
502 }
503 
504 
505 std::ostream& operator<<(std::ostream& out,
506  const PFEnergyCalibration& calib) {
507 
508  if(!out ) return out;
509 
510  out<<"PFEnergyCalibration -- "<<endl;
511 
512  if ( calib.pfCalibrations ) {
513 
514  static std::map<std::string, PerformanceResult::ResultType> functType;
515 
516  functType["PFfa_BARREL"] = PerformanceResult::PFfa_BARREL;
517  functType["PFfa_ENDCAP"] = PerformanceResult::PFfa_ENDCAP;
518  functType["PFfb_BARREL"] = PerformanceResult::PFfb_BARREL;
519  functType["PFfb_ENDCAP"] = PerformanceResult::PFfb_ENDCAP;
520  functType["PFfc_BARREL"] = PerformanceResult::PFfc_BARREL;
521  functType["PFfc_ENDCAP"] = PerformanceResult::PFfc_ENDCAP;
522  functType["PFfaEta_BARRELH"] = PerformanceResult::PFfaEta_BARRELH;
523  functType["PFfaEta_ENDCAPH"] = PerformanceResult::PFfaEta_ENDCAPH;
524  functType["PFfbEta_BARRELH"] = PerformanceResult::PFfbEta_BARRELH;
525  functType["PFfbEta_ENDCAPH"] = PerformanceResult::PFfbEta_ENDCAPH;
526  functType["PFfaEta_BARRELEH"] = PerformanceResult::PFfaEta_BARRELEH;
527  functType["PFfaEta_ENDCAPEH"] = PerformanceResult::PFfaEta_ENDCAPEH;
528  functType["PFfbEta_BARRELEH"] = PerformanceResult::PFfbEta_BARRELEH;
529  functType["PFfbEta_ENDCAPEH"] = PerformanceResult::PFfbEta_ENDCAPEH;
530 
531  for(std::map<std::string,PerformanceResult::ResultType>::const_iterator
532  func = functType.begin();
533  func != functType.end();
534  ++func) {
535 
536  cout << "Function: " << func->first << endl;
537  PerformanceResult::ResultType fType = func->second;
538  calib.pfCalibrations->printFormula(fType);
539  }
540 
541  } else {
542 
543  std::cout << "Default calibration functions : " << std::endl;
544 
545  calib.faBarrel->Print();
546  calib.fbBarrel->Print();
547  calib.fcBarrel->Print();
548  calib.faEtaBarrelEH->Print();
549  calib.fbEtaBarrelEH->Print();
550  calib.faEtaBarrelH->Print();
551  calib.fbEtaBarrelH->Print();
552  calib.faEndcap->Print();
553  calib.fbEndcap->Print();
554  calib.fcEndcap->Print();
555  calib.faEtaEndcapEH->Print();
556  calib.fbEtaEndcapEH->Print();
557  calib.faEtaEndcapH->Print();
558  calib.fbEtaEndcapH->Print();
559 
560  }
561 
562  return out;
563 }
564 
565 
566 
567 
580 
581 
582 
588 
589 
590 //useful to compute the signed distance to the closest crack in the barrel
591 double
592 PFEnergyCalibration::minimum(double a,double b) const {
593  if(TMath::Abs(b)<TMath::Abs(a)) a=b;
594  return a;
595 }
596 
597 namespace {
598  constexpr double pi= M_PI;// 3.14159265358979323846;
599 
600 
601  std::vector<double> fillcPhi() {
602  std::vector<double> retValue;
603  retValue.resize(18,0);
604  retValue[0]=2.97025;
605  for(unsigned i=1;i<=17;++i) retValue[i]=retValue[0]-2*i*pi/18;
606 
607  return retValue;
608  }
609 
610  //Location of the 18 phi-cracks
611  const std::vector<double> cPhi = fillcPhi();
612 }
613 
614 //compute the unsigned distance to the closest phi-crack in the barrel
615 double
616 PFEnergyCalibration::dCrackPhi(double phi, double eta) const {
617 
618 
619  //Shift of this location if eta<0
620  constexpr double delta_cPhi=0.00638;
621 
622  double m; //the result
623 
624  //the location is shifted
625  if(eta<0) phi +=delta_cPhi;
626 
627  if (phi>=-pi && phi<=pi){
628 
629  //the problem of the extrema
630  if (phi<cPhi[17] || phi>=cPhi[0]){
631  if (phi<0) phi+= 2*pi;
632  m = minimum(phi -cPhi[0],phi-cPhi[17]-2*pi);
633  }
634 
635  //between these extrema...
636  else{
637  bool OK = false;
638  unsigned i=16;
639  while(!OK){
640  if (phi<cPhi[i]){
641  m=minimum(phi-cPhi[i+1],phi-cPhi[i]);
642  OK=true;
643  }
644  else i-=1;
645  }
646  }
647  }
648  else{
649  m=0.; //if there is a problem, we assum that we are in a crack
650  std::cout<<"Problem in dminphi"<<std::endl;
651  }
652  if(eta<0) m=-m; //because of the disymetry
653  return m;
654 }
655 
656 // corrects the effect of phi-cracks
657 double
658 PFEnergyCalibration::CorrPhi(double phi, double eta) const {
659 
660  // we use 3 gaussians to correct the phi-cracks effect
661  constexpr double p1= 5.59379e-01;
662  constexpr double p2= -1.26607e-03;
663  constexpr double p3= 9.61133e-04;
664 
665  constexpr double p4= 1.81691e-01;
666  constexpr double p5= -4.97535e-03;
667  constexpr double p6= 1.31006e-03;
668 
669  constexpr double p7= 1.38498e-01;
670  constexpr double p8= 1.18599e-04;
671  constexpr double p9= 2.01858e-03;
672 
673 
674  double dminphi = dCrackPhi(phi,eta);
675 
676  double result = (1+p1*TMath::Gaus(dminphi,p2,p3)+p4*TMath::Gaus(dminphi,p5,p6)+p7*TMath::Gaus(dminphi,p8,p9));
677 
678  return result;
679 }
680 
681 
682 // corrects the effect of |eta|-cracks
683 double
685 
686  // we use a gaussian with a screwness for each of the 5 |eta|-cracks
687  constexpr double a[] = {6.13349e-01, 5.08146e-01, 4.44480e-01, 3.3487e-01, 7.65627e-01}; // amplitude
688  constexpr double m[] = {-1.79514e-02, 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00}; // mean
689  constexpr double s[] = {7.92382e-03, 3.06028e-03, 3.36139e-03, 3.94521e-03, 8.63950e-04}; // sigma
690  constexpr double sa[] = {1.27228e+01, 3.81517e-02, 1.63507e-01, -6.56480e-02, 1.87160e-01}; // screwness amplitude
691  constexpr double ss[] = {5.48753e-02, -1.00223e-02, 2.22866e-03, 4.26288e-04, 2.67937e-03}; // screwness sigma
692  double result = 1;
693 
694  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]));
695 
696  return result;
697 }
698 
699 
700 //corrects the global behaviour in the barrel
701 double
702 PFEnergyCalibration::CorrBarrel(double E, double eta) const {
703 
704  //Energy dependency
705  /*
706  //YM Parameters 52XX:
707  constexpr double p0=1.00000e+00;
708  constexpr double p1=3.27753e+01;
709  constexpr double p2=2.28552e-02;
710  constexpr double p3=3.06139e+00;
711  constexpr double p4=2.25135e-01;
712  constexpr double p5=1.47824e+00;
713  constexpr double p6=1.09e-02;
714  constexpr double p7=4.19343e+01;
715  */
716  constexpr double p0 = 0.9944;
717  constexpr double p1 = 9.827;
718  constexpr double p2 = 1.503;
719  constexpr double p3 = 1.196;
720  constexpr double p4 = 0.3349;
721  constexpr double p5 = 0.89;
722  constexpr double p6 = 0.004361;
723  constexpr double p7 = 51.51;
724  //Eta dependency
725  constexpr double p8=2.705593e-03;
726 
727  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);
728 
729  return result;
730 }
731 
732 
733 
742 
743 
744 //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
745 // Etot = Beta*eEcal + Gamma*(ePS1 + Alpha*ePS2)
746 
747 double
749 
750  //Energy dependency
751  constexpr double p0 = 5.97621e-01;
752 
753  //Eta dependency
754  constexpr double p1 =-1.86407e-01;
755  constexpr double p2 = 3.85197e-01;
756 
757  //so that <feta()> = 1
758  constexpr double norm = (p1+p2*(2.6+1.656)/2);
759 
760  double result = p0*(p1+p2*eta)/norm;
761 
762  return result;
763 }
764 
765 double
766 PFEnergyCalibration::Beta(double E, double eta) const {
767 
768  //Energy dependency
769  constexpr double p0 = 0.032;
770  constexpr double p1 = 9.70394e-02;
771  constexpr double p2 = 2.23072e+01;
772  constexpr double p3 = 100;
773 
774  //Eta dependency
775  constexpr double p4 = 1.02496e+00 ;
776  constexpr double p5 = -4.40176e-03 ;
777 
778  //so that <feta()> = 1
779  constexpr double norm = (p4+p5*(2.6+1.656)/2);
780 
781  double result = (1.0012+p0*TMath::Exp(-E/p3)+p1*TMath::Exp(-E/p2))*(p4+p5*eta)/norm;
782  return result;
783 }
784 
785 
786 double
787 PFEnergyCalibration::Gamma(double etaEcal) const {
788 
789  //Energy dependency
790  constexpr double p0 = 2.49752e-02;
791 
792  //Eta dependency
793  constexpr double p1 = 6.48816e-02;
794  constexpr double p2 = -1.59517e-02;
795 
796  //so that <feta()> = 1
797  constexpr double norm = (p1+p2*(2.6+1.656)/2);
798 
799  double result = p0*(p1+p2*etaEcal)/norm;
800 
801  return result;
802 }
803 
804 
805 
811 
812 
813 // returns the corrected energy in the barrel (0,1.48)
814 // Global Behaviour, phi and eta cracks are taken into account
815 double
816 PFEnergyCalibration::EcorrBarrel(double E, double eta, double phi,
817  bool crackCorrection ) const {
818 
819  // double result = E*CorrBarrel(E,eta)*CorrEta(eta)*CorrPhi(phi,eta);
820  double correction = crackCorrection ? std::max(CorrEta(eta),CorrPhi(phi,eta)) : 1.;
821  double result = E * CorrBarrel(E,eta) * correction;
822 
823  return result;
824 }
825 
826 
827 // returns the corrected energy in the area between the barrel and the PS (1.48,1.65)
828 double
830 
831  //Energy dependency
832  constexpr double p0 =1;
833  constexpr double p1 =0.18;
834  constexpr double p2 =8.;
835 
836  //Eta dependency
837  constexpr double p3 =0.3;
838  constexpr double p4 =1.11;
839  constexpr double p5 =0.025;
840  constexpr double p6 =1.49;
841  constexpr double p7 =0.6;
842 
843  //so that <feta()> = 1
844  constexpr double norm = 1.21;
845 
846  double result = E*(p0+p1*TMath::Exp(-E/p2))*(p3+p4*TMath::Gaus(eta,p6,p5)+p7*eta)/norm;
847 
848  return result;
849 }
850 
851 
852 // returns the corrected energy in the PS (1.65,2.6)
853 // only when (ePS1>0)||(ePS2>0)
854 double
855 PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal) const {
856 
857  // gives the good weights to each subdetector
858  double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+Gamma(etaEcal)*(ePS1+Alpha(etaEcal)*ePS2)/9e-5 ;
859 
860  //Correction of the residual energy dependency
861  constexpr double p0 = 1.00;
862  constexpr double p1 = 2.18;
863  constexpr double p2 =1.94;
864  constexpr double p3 =4.13;
865  constexpr double p4 =1.127;
866 
867  double result = E*(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
868 
869  return result;
870 }
871 
872 // returns the corrected energy in the PS (1.65,2.6)
873 // only when (ePS1>0)||(ePS2>0)
874 double
875 PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal,double & outputPS1, double & outputPS2) const {
876 
877  // gives the good weights to each subdetector
878  double gammaprime=Gamma(etaEcal)/9e-5;
879 
880  if(outputPS1 == 0 && outputPS2 == 0 && esEEInterCalib_ != 0){
881  // both ES planes working
882  // scaling factor accounting for data-mc
883  outputPS1=gammaprime*ePS1 * esEEInterCalib_->getGammaLow0();
884  outputPS2=gammaprime*Alpha(etaEcal)*ePS2 * esEEInterCalib_->getGammaLow3();
885  }
886  else if(outputPS1 == 0 && outputPS2 == -1 && esEEInterCalib_ != 0){
887  // ESP1 only working
888  double corrTotES = gammaprime*ePS1 * esEEInterCalib_->getGammaLow0() * esEEInterCalib_->getGammaLow1();
889  outputPS1 = gammaprime*ePS1 * esEEInterCalib_->getGammaLow0();
890  outputPS2 = corrTotES - outputPS1;
891  }
892  else if(outputPS1 == -1 && outputPS2 == 0 && esEEInterCalib_ != 0){
893  // ESP2 only working
894  double corrTotES = gammaprime*Alpha(etaEcal)*ePS2 * esEEInterCalib_->getGammaLow3() * esEEInterCalib_->getGammaLow2();
895  outputPS2 = gammaprime*Alpha(etaEcal)*ePS2 * esEEInterCalib_->getGammaLow3();
896  outputPS1 = corrTotES - outputPS2;
897  }
898  else{
899  // none working
900  outputPS1 = gammaprime*ePS1;
901  outputPS2 = gammaprime*Alpha(etaEcal)*ePS2;
902  }
903 
904  double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+outputPS1+outputPS2;
905 
906  //Correction of the residual energy dependency
907  constexpr double p0 = 1.00;
908  constexpr double p1 = 2.18;
909  constexpr double p2 =1.94;
910  constexpr double p3 =4.13;
911  constexpr double p4 =1.127;
912 
913  double corrfac=(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
914  outputPS1*=corrfac;
915  outputPS2*=corrfac;
916  double result = E*corrfac;
917 
918  return result;
919 }
920 
921 
922 // returns the corrected energy in the PS (1.65,2.6)
923 // only when (ePS1=0)&&(ePS2=0)
924 double
925 PFEnergyCalibration::EcorrPS_ePSNil(double eEcal,double eta) const {
926 
927  //Energy dependency
928  constexpr double p0= 1.02;
929  constexpr double p1= 0.165;
930  constexpr double p2= 6.5 ;
931  constexpr double p3= 2.1 ;
932 
933  //Eta dependency
934  constexpr double p4 = 1.02496e+00 ;
935  constexpr double p5 = -4.40176e-03 ;
936 
937  //so that <feta()> = 1
938  constexpr double norm = (p4+p5*(2.6+1.656)/2);
939 
940  double result = eEcal*(p0+p1*TMath::Exp(-TMath::Abs(eEcal-p3)/p2))*(p4+p5*eta)/norm;
941 
942  return result;
943 }
944 
945 
946 // returns the corrected energy in the area between the end of the PS and the end of the endcap (2.6,2.98)
947 double
949 
950  //Energy dependency
951  constexpr double p0 =1;
952  constexpr double p1 = 0.058;
953  constexpr double p2 =12.5;
954  constexpr double p3 =-1.05444e+00;
955  constexpr double p4 =-5.39557e+00;
956  constexpr double p5 =8.38444e+00;
957  constexpr double p6 = 6.10998e-01 ;
958 
959  //Eta dependency
960  constexpr double p7 =1.06161e+00;
961  constexpr double p8 = 0.41;
962  constexpr double p9 =2.918;
963  constexpr double p10 =0.0181;
964  constexpr double p11= 2.05;
965  constexpr double p12 =2.99;
966  constexpr double p13=0.0287;
967 
968  //so that <feta()> = 1
969  constexpr double norm=1.045;
970 
971  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;
972  return result;
973 }
974 
975 
976 
977 
978 // returns the corrected energy everywhere
979 // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
980 double
981 PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,
982  double eta,double phi,
983  bool crackCorrection ) const {
984 
985  constexpr double endBarrel=1.48;
986  constexpr double beginingPS=1.65;
987  constexpr double endPS=2.6;
988  constexpr double endEndCap=2.98;
989 
990  double result=0;
991 
992  eta=TMath::Abs(eta);
993 
994  if(eEcal>0){
995  if(eta <= endBarrel) result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
996  else if(eta <= beginingPS) result = EcorrZoneBeforePS(eEcal,eta);
997  else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
998  else if(eta < endPS) result = EcorrPS(eEcal,ePS1,ePS2,eta);
999  else if(eta < endEndCap) result = EcorrZoneAfterPS(eEcal,eta);
1000  else result =eEcal;
1001  }
1002  else result = eEcal;// useful if eEcal=0 or eta>2.98
1003  //protection
1004  if(result<eEcal) result=eEcal;
1005  return result;
1006 }
1007 
1008 // returns the corrected energy everywhere
1009 // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
1010 double
1011 PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,double eta,double phi,double& ps1,double&ps2,bool crackCorrection) const {
1012 
1013  constexpr double endBarrel=1.48;
1014  constexpr double beginingPS=1.65;
1015  constexpr double endPS=2.6;
1016  constexpr double endEndCap=2.98;
1017 
1018  double result=0;
1019 
1020  eta=TMath::Abs(eta);
1021 
1022  if(eEcal>0){
1023  if(eta <= endBarrel) result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
1024  else if(eta <= beginingPS) result = EcorrZoneBeforePS(eEcal,eta);
1025  else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
1026  else if(eta < endPS) result = EcorrPS(eEcal,ePS1,ePS2,eta,ps1,ps2);
1027  else if(eta < endEndCap) result = EcorrZoneAfterPS(eEcal,eta);
1028  else result =eEcal;
1029  }
1030  else result = eEcal;// useful if eEcal=0 or eta>2.98
1031  // protection
1032  if(result<eEcal) result=eEcal;
1033  return result;
1034 }
std::unique_ptr< TF1 > faEtaBarrelH
const PerformancePayloadFromTFormula * pfCalibrations
double aEtaEndcapEH(double x) const
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
std::unique_ptr< TF1 > faEtaBarrelEH
double bBarrel(double x) const
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:47
T Sign(T A, T B)
Definition: MathUtil.h:54
double aEtaBarrelEH(double x) const
double aEndcap(double x) const
double bEndcap(double x) const
const ESEEIntercalibConstants * esEEInterCalib_
double bEtaEndcapH(double x) const
double Ecorr(double eEcal, double ePS1, double ePS2, double eta, double phi, bool crackCorrection=true) const
std::unique_ptr< TF1 > fbEtaBarrelEH
double aEtaBarrelH(double x) const
double EcorrZoneBeforePS(double E, double eta) const
#define constexpr
double Alpha(double eta) const
void energyEmHad(double t, double &e, double &h, double eta, double phi) const
const Double_t pi
void calculatePositionREP()
computes posrep_ once and for all
Definition: PFCluster.h:100
double bEtaBarrelEH(double x) const
double p4[4]
Definition: TauolaWrapper.h:92
std::unique_ptr< TF1 > fbBarrel
std::unique_ptr< TF1 > faEtaEndcapEH
friend std::ostream & operator<<(std::ostream &out, const PFEnergyCalibration &calib)
double cEndcap(double x) const
const REPPoint & positionREP() const
cluster position: rho, eta, phi
Definition: PFCluster.h:97
T Abs(T a)
Definition: MathUtil.h:49
double Beta(double E, double eta) const
std::unique_ptr< TF1 > faBarrel
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::unique_ptr< TF1 > faEndcap
T min(T a, T b)
Definition: MathUtil.h:58
double p2[4]
Definition: TauolaWrapper.h:90
double energy() const
cluster energy
Definition: PFCluster.h:82
float getResult(PerformanceResult::ResultType, const BinningPointByMap &) const
#define M_PI
double bEtaEndcapEH(double x) const
double EcorrZoneAfterPS(double E, double eta) const
std::pair< int, edm::FunctionWithDict > OK
Definition: findMethod.cc:136
bool insert(BinningVariables::BinningVariablesType, float)
double bEtaBarrelH(double x) const
std::unique_ptr< TF1 > fbEtaBarrelH
double b
Definition: hdecay.h:120
std::unique_ptr< TF1 > faEtaEndcapH
double CorrBarrel(double E, double eta) const
double Gamma(double etaEcal) const
double energyEm(const reco::PFCluster &clusterEcal, std::vector< double > &EclustersPS1, std::vector< double > &EclustersPS2, bool crackCorrection=true) const
double dCrackPhi(double phi, double eta) const
double p1[4]
Definition: TauolaWrapper.h:89
double a
Definition: hdecay.h:121
std::unique_ptr< TF1 > fbEndcap
std::unique_ptr< TF1 > fcBarrel
void printFormula(PerformanceResult::ResultType res) const
double CorrEta(double eta) const
double aEtaEndcapH(double x) const
double aBarrel(double x) const
double CorrPhi(double phi, double eta) const
double EcorrPS_ePSNil(double eEcal, double eta) const
std::unique_ptr< TF1 > fbEtaEndcapEH
double cBarrel(double x) const
double EcorrBarrel(double E, double eta, double phi, bool crackCorrection=true) const
std::unique_ptr< TF1 > fcEndcap
*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
std::unique_ptr< TF1 > fbEtaEndcapH
double EcorrPS(double eEcal, double ePS1, double ePS2, double etaEcal) const
double p3[4]
Definition: TauolaWrapper.h:91
double minimum(double a, double b) const