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