CMS 3D CMS Logo

PFEnergyCalibration.cc
Go to the documentation of this file.
3 
5 
6 #include <TMath.h>
7 #include <cmath>
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(nullptr), esEEInterCalib_(nullptr)
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 
192  if ( e > 0. && thresh > 0. ) {
193  if(absEta<2.5) {
194  etaCorrE = 1. + aEtaEndcapEH(t) ;
195  }
196  else {
197  etaPow = dEta * dEta;
198  etaCorrE = 1. + aEtaEndcapEH(t) + 1.3*bEtaEndcapEH(t)*(0.6 + etaPow);
199  }
200  etaPow = dEta * dEta * dEta * dEta;
201  etaCorrH = 1. + aEtaEndcapEH(t) + bEtaEndcapEH(t)*(0.04 + etaPow);
202  } else {
203  etaCorrE = 1. + aEtaEndcapH(t) + 1.3*bEtaEndcapH(t)*(0.04 + etaPow);
204  if(absEta<2.5) {
205  etaCorrH = 1. + aEtaEndcapH(t) + 0.05*bEtaEndcapH(t);
206  }
207  else {
208  etaCorrH = 1. + aEtaEndcapH(t) + bEtaEndcapH(t)*(etaPow - 1.1);
209  }
210  }
211 
212  //t = min(999.9,max(tt, thresh + etaCorrE*a*e + etaCorrH*b*h));
213 
214  if ( e > 0. && thresh > 0. )
215  e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
216  if ( h > 0. && thresh > 0. ) {
217  h = threshH + etaCorrH * b * h;
218  }
219  }
220 
221  // Protection
222  if ( e < 0. || h < 0. ) {
223 
224  // Some protection against crazy calibration
225  if ( e < 0. ) e = ee;
226  if ( h < 0. ) h = hh;
227  }
228 
229  // And that's it !
230 
231 
232 }
233 
234 // The calibration functions
235 double
237 
238  if ( pfCalibrations ) {
242 
243  } else {
244 
245  return faBarrel->Eval(x);
246 
247  }
248 }
249 
250 double
252 
253  if ( pfCalibrations ) {
254 
258 
259  } else {
260 
261  return fbBarrel->Eval(x);
262 
263  }
264 }
265 
266 double
268 
269  if ( pfCalibrations ) {
270 
274 
275  } else {
276 
277  return fcBarrel->Eval(x);
278 
279  }
280 }
281 
282 double
284 
285  if ( pfCalibrations ) {
286 
290 
291  } else {
292 
293  return faEtaBarrelEH->Eval(x);
294 
295  }
296 }
297 
298 double
300 
301  if ( pfCalibrations ) {
302 
306 
307  } else {
308 
309  return fbEtaBarrelEH->Eval(x);
310 
311  }
312 }
313 
314 double
316 
317  if ( pfCalibrations ) {
318 
322 
323  } else {
324 
325  return faEtaBarrelH->Eval(x);
326 
327  }
328 }
329 
330 double
332 
333  if ( pfCalibrations ) {
334 
338 
339  } else {
340 
341  return fbEtaBarrelH->Eval(x);
342 
343  }
344 }
345 
346 double
348 
349  if ( pfCalibrations ) {
350 
354 
355  } else {
356 
357  return faEndcap->Eval(x);
358 
359  }
360 }
361 
362 double
364 
365  if ( pfCalibrations ) {
366 
370 
371  } else {
372 
373  return fbEndcap->Eval(x);
374 
375  }
376 }
377 
378 double
380 
381  if ( pfCalibrations ) {
382 
386 
387  } else {
388 
389  return fcEndcap->Eval(x);
390 
391  }
392 }
393 
394 double
396 
397  if ( pfCalibrations ) {
398 
402 
403  } else {
404 
405  return faEtaEndcapEH->Eval(x);
406 
407  }
408 }
409 
410 double
412 
413  if ( pfCalibrations ) {
414 
418 
419  } else {
420 
421  return fbEtaEndcapEH->Eval(x);
422 
423  }
424 }
425 
426 double
428 
429  if ( pfCalibrations ) {
430 
434 
435  } else {
436 
437  return faEtaEndcapH->Eval(x);
438 
439  }
440 }
441 
442 double
444 
445  if ( pfCalibrations ) {
446 
450 
451  } else {
452 
453  return fbEtaEndcapH->Eval(x);
454 
455  }
456 }
457 
458 
459 double
461  std::vector<double> &EclustersPS1,
462  std::vector<double> &EclustersPS2,
463  bool crackCorrection ) const {
464  double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
465  double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
466  return energyEm(clusterEcal, ePS1, ePS2, crackCorrection);
467 }
468 
469 double
471  double ePS1,
472  double ePS2,
473  bool crackCorrection ) const {
474  double eEcal = clusterEcal.energy();
475  //temporaty ugly fix
476  reco::PFCluster myPFCluster=clusterEcal;
477  myPFCluster.calculatePositionREP();
478  double eta = myPFCluster.positionREP().eta();
479  double phi = myPFCluster.positionREP().phi();
480 
481  double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi, crackCorrection);
482  // if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
483  return calibrated;
484 }
485 
487  std::vector<double> &EclustersPS1,
488  std::vector<double> &EclustersPS2,
489  double& ps1,double& ps2,
490  bool crackCorrection) const {
491  double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
492  double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
493  return energyEm(clusterEcal, ePS1, ePS2, ps1, ps2, crackCorrection);
494 }
496  double ePS1, double ePS2,
497  double& ps1,double& ps2,
498  bool crackCorrection) const {
499  double eEcal = clusterEcal.energy();
500  //temporaty ugly fix
501  reco::PFCluster myPFCluster=clusterEcal;
502  myPFCluster.calculatePositionREP();
503  double eta = myPFCluster.positionREP().eta();
504  double phi = myPFCluster.positionREP().phi();
505 
506  double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi,ps1,ps2,crackCorrection);
507  // if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
508  return calibrated;
509 }
510 
511 
512 std::ostream& operator<<(std::ostream& out,
513  const PFEnergyCalibration& calib) {
514 
515  if(!out ) return out;
516 
517  out<<"PFEnergyCalibration -- "<<endl;
518 
519  if ( calib.pfCalibrations ) {
520 
521  static std::map<std::string, PerformanceResult::ResultType> functType;
522 
523  functType["PFfa_BARREL"] = PerformanceResult::PFfa_BARREL;
524  functType["PFfa_ENDCAP"] = PerformanceResult::PFfa_ENDCAP;
525  functType["PFfb_BARREL"] = PerformanceResult::PFfb_BARREL;
526  functType["PFfb_ENDCAP"] = PerformanceResult::PFfb_ENDCAP;
527  functType["PFfc_BARREL"] = PerformanceResult::PFfc_BARREL;
528  functType["PFfc_ENDCAP"] = PerformanceResult::PFfc_ENDCAP;
529  functType["PFfaEta_BARRELH"] = PerformanceResult::PFfaEta_BARRELH;
530  functType["PFfaEta_ENDCAPH"] = PerformanceResult::PFfaEta_ENDCAPH;
531  functType["PFfbEta_BARRELH"] = PerformanceResult::PFfbEta_BARRELH;
532  functType["PFfbEta_ENDCAPH"] = PerformanceResult::PFfbEta_ENDCAPH;
533  functType["PFfaEta_BARRELEH"] = PerformanceResult::PFfaEta_BARRELEH;
534  functType["PFfaEta_ENDCAPEH"] = PerformanceResult::PFfaEta_ENDCAPEH;
535  functType["PFfbEta_BARRELEH"] = PerformanceResult::PFfbEta_BARRELEH;
536  functType["PFfbEta_ENDCAPEH"] = PerformanceResult::PFfbEta_ENDCAPEH;
537 
538  for(std::map<std::string,PerformanceResult::ResultType>::const_iterator
539  func = functType.begin();
540  func != functType.end();
541  ++func) {
542 
543  cout << "Function: " << func->first << endl;
544  PerformanceResult::ResultType fType = func->second;
545  calib.pfCalibrations->printFormula(fType);
546  }
547 
548  } else {
549 
550  std::cout << "Default calibration functions : " << std::endl;
551 
552  calib.faBarrel->Print();
553  calib.fbBarrel->Print();
554  calib.fcBarrel->Print();
555  calib.faEtaBarrelEH->Print();
556  calib.fbEtaBarrelEH->Print();
557  calib.faEtaBarrelH->Print();
558  calib.fbEtaBarrelH->Print();
559  calib.faEndcap->Print();
560  calib.fbEndcap->Print();
561  calib.fcEndcap->Print();
562  calib.faEtaEndcapEH->Print();
563  calib.fbEtaEndcapEH->Print();
564  calib.faEtaEndcapH->Print();
565  calib.fbEtaEndcapH->Print();
566 
567  }
568 
569  return out;
570 }
571 
572 
573 
574 
587 
588 
589 
595 
596 
597 //useful to compute the signed distance to the closest crack in the barrel
598 double
599 PFEnergyCalibration::minimum(double a,double b) const {
600  if(TMath::Abs(b)<TMath::Abs(a)) a=b;
601  return a;
602 }
603 
604 namespace {
605  constexpr double pi= M_PI;// 3.14159265358979323846;
606 
607 
608  std::vector<double> fillcPhi() {
609  std::vector<double> retValue;
610  retValue.resize(18,0);
611  retValue[0]=2.97025;
612  for(unsigned i=1;i<=17;++i) retValue[i]=retValue[0]-2*i*pi/18;
613 
614  return retValue;
615  }
616 
617  //Location of the 18 phi-cracks
618  const std::vector<double> cPhi = fillcPhi();
619 }
620 
621 //compute the unsigned distance to the closest phi-crack in the barrel
622 double
623 PFEnergyCalibration::dCrackPhi(double phi, double eta) const {
624 
625 
626  //Shift of this location if eta<0
627  constexpr double delta_cPhi=0.00638;
628 
629  double m; //the result
630 
631  //the location is shifted
632  if(eta<0) phi +=delta_cPhi;
633 
634  if (phi>=-pi && phi<=pi){
635 
636  //the problem of the extrema
637  if (phi<cPhi[17] || phi>=cPhi[0]){
638  if (phi<0) phi+= 2*pi;
639  m = minimum(phi -cPhi[0],phi-cPhi[17]-2*pi);
640  }
641 
642  //between these extrema...
643  else{
644  bool OK = false;
645  unsigned i=16;
646  while(!OK){
647  if (phi<cPhi[i]){
648  m=minimum(phi-cPhi[i+1],phi-cPhi[i]);
649  OK=true;
650  }
651  else i-=1;
652  }
653  }
654  }
655  else{
656  m=0.; //if there is a problem, we assum that we are in a crack
657  std::cout<<"Problem in dminphi"<<std::endl;
658  }
659  if(eta<0) m=-m; //because of the disymetry
660  return m;
661 }
662 
663 // corrects the effect of phi-cracks
664 double
665 PFEnergyCalibration::CorrPhi(double phi, double eta) const {
666 
667  // we use 3 gaussians to correct the phi-cracks effect
668  constexpr double p1= 5.59379e-01;
669  constexpr double p2= -1.26607e-03;
670  constexpr double p3= 9.61133e-04;
671 
672  constexpr double p4= 1.81691e-01;
673  constexpr double p5= -4.97535e-03;
674  constexpr double p6= 1.31006e-03;
675 
676  constexpr double p7= 1.38498e-01;
677  constexpr double p8= 1.18599e-04;
678  constexpr double p9= 2.01858e-03;
679 
680 
681  double dminphi = dCrackPhi(phi,eta);
682 
683  double result = (1+p1*TMath::Gaus(dminphi,p2,p3)+p4*TMath::Gaus(dminphi,p5,p6)+p7*TMath::Gaus(dminphi,p8,p9));
684 
685  return result;
686 }
687 
688 
689 // corrects the effect of |eta|-cracks
690 double
692 
693  // we use a gaussian with a screwness for each of the 5 |eta|-cracks
694  constexpr double a[] = {6.13349e-01, 5.08146e-01, 4.44480e-01, 3.3487e-01, 7.65627e-01}; // amplitude
695  constexpr double m[] = {-1.79514e-02, 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00}; // mean
696  constexpr double s[] = {7.92382e-03, 3.06028e-03, 3.36139e-03, 3.94521e-03, 8.63950e-04}; // sigma
697  constexpr double sa[] = {1.27228e+01, 3.81517e-02, 1.63507e-01, -6.56480e-02, 1.87160e-01}; // screwness amplitude
698  constexpr double ss[] = {5.48753e-02, -1.00223e-02, 2.22866e-03, 4.26288e-04, 2.67937e-03}; // screwness sigma
699  double result = 1;
700 
701  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]));
702 
703  return result;
704 }
705 
706 
707 //corrects the global behaviour in the barrel
708 double
709 PFEnergyCalibration::CorrBarrel(double E, double eta) const {
710 
711  //Energy dependency
712  /*
713  //YM Parameters 52XX:
714  constexpr double p0=1.00000e+00;
715  constexpr double p1=3.27753e+01;
716  constexpr double p2=2.28552e-02;
717  constexpr double p3=3.06139e+00;
718  constexpr double p4=2.25135e-01;
719  constexpr double p5=1.47824e+00;
720  constexpr double p6=1.09e-02;
721  constexpr double p7=4.19343e+01;
722  */
723  constexpr double p0 = 0.9944;
724  constexpr double p1 = 9.827;
725  constexpr double p2 = 1.503;
726  constexpr double p3 = 1.196;
727  constexpr double p4 = 0.3349;
728  constexpr double p5 = 0.89;
729  constexpr double p6 = 0.004361;
730  constexpr double p7 = 51.51;
731  //Eta dependency
732  constexpr double p8=2.705593e-03;
733 
734  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);
735 
736  return result;
737 }
738 
739 
740 
749 
750 
751 //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
752 // Etot = Beta*eEcal + Gamma*(ePS1 + Alpha*ePS2)
753 
754 double
756 
757  //Energy dependency
758  constexpr double p0 = 5.97621e-01;
759 
760  //Eta dependency
761  constexpr double p1 =-1.86407e-01;
762  constexpr double p2 = 3.85197e-01;
763 
764  //so that <feta()> = 1
765  constexpr double norm = (p1+p2*(2.6+1.656)/2);
766 
767  double result = p0*(p1+p2*eta)/norm;
768 
769  return result;
770 }
771 
772 double
773 PFEnergyCalibration::Beta(double E, double eta) const {
774 
775  //Energy dependency
776  constexpr double p0 = 0.032;
777  constexpr double p1 = 9.70394e-02;
778  constexpr double p2 = 2.23072e+01;
779  constexpr double p3 = 100;
780 
781  //Eta dependency
782  constexpr double p4 = 1.02496e+00 ;
783  constexpr double p5 = -4.40176e-03 ;
784 
785  //so that <feta()> = 1
786  constexpr double norm = (p4+p5*(2.6+1.656)/2);
787 
788  double result = (1.0012+p0*TMath::Exp(-E/p3)+p1*TMath::Exp(-E/p2))*(p4+p5*eta)/norm;
789  return result;
790 }
791 
792 
793 double
794 PFEnergyCalibration::Gamma(double etaEcal) const {
795 
796  //Energy dependency
797  constexpr double p0 = 2.49752e-02;
798 
799  //Eta dependency
800  constexpr double p1 = 6.48816e-02;
801  constexpr double p2 = -1.59517e-02;
802 
803  //so that <feta()> = 1
804  constexpr double norm = (p1+p2*(2.6+1.656)/2);
805 
806  double result = p0*(p1+p2*etaEcal)/norm;
807 
808  return result;
809 }
810 
811 
812 
818 
819 
820 // returns the corrected energy in the barrel (0,1.48)
821 // Global Behaviour, phi and eta cracks are taken into account
822 double
823 PFEnergyCalibration::EcorrBarrel(double E, double eta, double phi,
824  bool crackCorrection ) const {
825 
826  // double result = E*CorrBarrel(E,eta)*CorrEta(eta)*CorrPhi(phi,eta);
827  double correction = crackCorrection ? std::max(CorrEta(eta),CorrPhi(phi,eta)) : 1.;
828  double result = E * CorrBarrel(E,eta) * correction;
829 
830  return result;
831 }
832 
833 
834 // returns the corrected energy in the area between the barrel and the PS (1.48,1.65)
835 double
837 
838  //Energy dependency
839  constexpr double p0 =1;
840  constexpr double p1 =0.18;
841  constexpr double p2 =8.;
842 
843  //Eta dependency
844  constexpr double p3 =0.3;
845  constexpr double p4 =1.11;
846  constexpr double p5 =0.025;
847  constexpr double p6 =1.49;
848  constexpr double p7 =0.6;
849 
850  //so that <feta()> = 1
851  constexpr double norm = 1.21;
852 
853  double result = E*(p0+p1*TMath::Exp(-E/p2))*(p3+p4*TMath::Gaus(eta,p6,p5)+p7*eta)/norm;
854 
855  return result;
856 }
857 
858 
859 // returns the corrected energy in the PS (1.65,2.6)
860 // only when (ePS1>0)||(ePS2>0)
861 double
862 PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal) const {
863 
864  // gives the good weights to each subdetector
865  double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+Gamma(etaEcal)*(ePS1+Alpha(etaEcal)*ePS2)/9e-5 ;
866 
867  //Correction of the residual energy dependency
868  constexpr double p0 = 1.00;
869  constexpr double p1 = 2.18;
870  constexpr double p2 =1.94;
871  constexpr double p3 =4.13;
872  constexpr double p4 =1.127;
873 
874  double result = E*(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
875 
876  return result;
877 }
878 
879 // returns the corrected energy in the PS (1.65,2.6)
880 // only when (ePS1>0)||(ePS2>0)
881 double
882 PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal,double & outputPS1, double & outputPS2) const {
883 
884  // gives the good weights to each subdetector
885  double gammaprime=Gamma(etaEcal)/9e-5;
886 
887  if(outputPS1 == 0 && outputPS2 == 0 && esEEInterCalib_ != nullptr){
888  // both ES planes working
889  // scaling factor accounting for data-mc
890  outputPS1=gammaprime*ePS1 * esEEInterCalib_->getGammaLow0();
891  outputPS2=gammaprime*Alpha(etaEcal)*ePS2 * esEEInterCalib_->getGammaLow3();
892  }
893  else if(outputPS1 == 0 && outputPS2 == -1 && esEEInterCalib_ != nullptr){
894  // ESP1 only working
895  double corrTotES = gammaprime*ePS1 * esEEInterCalib_->getGammaLow0() * esEEInterCalib_->getGammaLow1();
896  outputPS1 = gammaprime*ePS1 * esEEInterCalib_->getGammaLow0();
897  outputPS2 = corrTotES - outputPS1;
898  }
899  else if(outputPS1 == -1 && outputPS2 == 0 && esEEInterCalib_ != nullptr){
900  // ESP2 only working
901  double corrTotES = gammaprime*Alpha(etaEcal)*ePS2 * esEEInterCalib_->getGammaLow3() * esEEInterCalib_->getGammaLow2();
902  outputPS2 = gammaprime*Alpha(etaEcal)*ePS2 * esEEInterCalib_->getGammaLow3();
903  outputPS1 = corrTotES - outputPS2;
904  }
905  else{
906  // none working
907  outputPS1 = gammaprime*ePS1;
908  outputPS2 = gammaprime*Alpha(etaEcal)*ePS2;
909  }
910 
911  double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+outputPS1+outputPS2;
912 
913  //Correction of the residual energy dependency
914  constexpr double p0 = 1.00;
915  constexpr double p1 = 2.18;
916  constexpr double p2 =1.94;
917  constexpr double p3 =4.13;
918  constexpr double p4 =1.127;
919 
920  double corrfac=(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
921  outputPS1*=corrfac;
922  outputPS2*=corrfac;
923  double result = E*corrfac;
924 
925  return result;
926 }
927 
928 
929 // returns the corrected energy in the PS (1.65,2.6)
930 // only when (ePS1=0)&&(ePS2=0)
931 double
932 PFEnergyCalibration::EcorrPS_ePSNil(double eEcal,double eta) const {
933 
934  //Energy dependency
935  constexpr double p0= 1.02;
936  constexpr double p1= 0.165;
937  constexpr double p2= 6.5 ;
938  constexpr double p3= 2.1 ;
939 
940  //Eta dependency
941  constexpr double p4 = 1.02496e+00 ;
942  constexpr double p5 = -4.40176e-03 ;
943 
944  //so that <feta()> = 1
945  constexpr double norm = (p4+p5*(2.6+1.656)/2);
946 
947  double result = eEcal*(p0+p1*TMath::Exp(-TMath::Abs(eEcal-p3)/p2))*(p4+p5*eta)/norm;
948 
949  return result;
950 }
951 
952 
953 // returns the corrected energy in the area between the end of the PS and the end of the endcap (2.6,2.98)
954 double
956 
957  //Energy dependency
958  constexpr double p0 =1;
959  constexpr double p1 = 0.058;
960  constexpr double p2 =12.5;
961  constexpr double p3 =-1.05444e+00;
962  constexpr double p4 =-5.39557e+00;
963  constexpr double p5 =8.38444e+00;
964  constexpr double p6 = 6.10998e-01 ;
965 
966  //Eta dependency
967  constexpr double p7 =1.06161e+00;
968  constexpr double p8 = 0.41;
969  constexpr double p9 =2.918;
970  constexpr double p10 =0.0181;
971  constexpr double p11= 2.05;
972  constexpr double p12 =2.99;
973  constexpr double p13=0.0287;
974 
975  //so that <feta()> = 1
976  constexpr double norm=1.045;
977 
978  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;
979  return result;
980 }
981 
982 
983 
984 
985 // returns the corrected energy everywhere
986 // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
987 double
988 PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,
989  double eta,double phi,
990  bool crackCorrection ) const {
991 
992  constexpr double endBarrel=1.48;
993  constexpr double beginingPS=1.65;
994  constexpr double endPS=2.6;
995  constexpr double endEndCap=2.98;
996 
997  double result=0;
998 
999  eta=TMath::Abs(eta);
1000 
1001  if(eEcal>0){
1002  if(eta <= endBarrel) result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
1003  else if(eta <= beginingPS) result = EcorrZoneBeforePS(eEcal,eta);
1004  else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
1005  else if(eta < endPS) result = EcorrPS(eEcal,ePS1,ePS2,eta);
1006  else if(eta < endEndCap) result = EcorrZoneAfterPS(eEcal,eta);
1007  else result =eEcal;
1008  }
1009  else result = eEcal;// useful if eEcal=0 or eta>2.98
1010  //protection
1011  if(result<eEcal) result=eEcal;
1012  return result;
1013 }
1014 
1015 // returns the corrected energy everywhere
1016 // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
1017 double
1018 PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,double eta,double phi,double& ps1,double&ps2,bool crackCorrection) const {
1019 
1020  constexpr double endBarrel=1.48;
1021  constexpr double beginingPS=1.65;
1022  constexpr double endPS=2.6;
1023  constexpr double endEndCap=2.98;
1024 
1025  double result=0;
1026 
1027  eta=TMath::Abs(eta);
1028 
1029  if(eEcal>0){
1030  if(eta <= endBarrel) result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
1031  else if(eta <= beginingPS) result = EcorrZoneBeforePS(eEcal,eta);
1032  else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
1033  else if(eta < endPS) result = EcorrPS(eEcal,ePS1,ePS2,eta,ps1,ps2);
1034  else if(eta < endEndCap) result = EcorrZoneAfterPS(eEcal,eta);
1035  else result =eEcal;
1036  }
1037  else result = eEcal;// useful if eEcal=0 or eta>2.98
1038  // protection
1039  if(result<eEcal) result=eEcal;
1040  return result;
1041 }
std::unique_ptr< TF1 > faEtaBarrelH
const PerformancePayloadFromTFormula * pfCalibrations
double aEtaEndcapEH(double x) const
float getResult(PerformanceResult::ResultType, const BinningPointByMap &) const override
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
#define nullptr
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
#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