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