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