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