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 
524 //compute the unsigned distance to the closest phi-crack in the barrel
525 double
526 PFEnergyCalibration::dCrackPhi(double phi, double eta){
527 
528  constexpr double pi= M_PI;// 3.14159265358979323846;
529 
530  //Location of the 18 phi-cracks
531  static std::vector<double> cPhi;
532  if(cPhi.size()==0)
533  {
534  cPhi.resize(18,0);
535  cPhi[0]=2.97025;
536  for(unsigned i=1;i<=17;++i) cPhi[i]=cPhi[0]-2*i*pi/18;
537  }
538 
539  //Shift of this location if eta<0
540  constexpr double delta_cPhi=0.00638;
541 
542  double m; //the result
543 
544  //the location is shifted
545  if(eta<0) phi +=delta_cPhi;
546 
547  if (phi>=-pi && phi<=pi){
548 
549  //the problem of the extrema
550  if (phi<cPhi[17] || phi>=cPhi[0]){
551  if (phi<0) phi+= 2*pi;
552  m = minimum(phi -cPhi[0],phi-cPhi[17]-2*pi);
553  }
554 
555  //between these extrema...
556  else{
557  bool OK = false;
558  unsigned i=16;
559  while(!OK){
560  if (phi<cPhi[i]){
561  m=minimum(phi-cPhi[i+1],phi-cPhi[i]);
562  OK=true;
563  }
564  else i-=1;
565  }
566  }
567  }
568  else{
569  m=0.; //if there is a problem, we assum that we are in a crack
570  std::cout<<"Problem in dminphi"<<std::endl;
571  }
572  if(eta<0) m=-m; //because of the disymetry
573  return m;
574 }
575 
576 // corrects the effect of phi-cracks
577 double
578 PFEnergyCalibration::CorrPhi(double phi, double eta) {
579 
580  // we use 3 gaussians to correct the phi-cracks effect
581  constexpr double p1= 5.59379e-01;
582  constexpr double p2= -1.26607e-03;
583  constexpr double p3= 9.61133e-04;
584 
585  constexpr double p4= 1.81691e-01;
586  constexpr double p5= -4.97535e-03;
587  constexpr double p6= 1.31006e-03;
588 
589  constexpr double p7= 1.38498e-01;
590  constexpr double p8= 1.18599e-04;
591  constexpr double p9= 2.01858e-03;
592 
593 
594  double dminphi = dCrackPhi(phi,eta);
595 
596  double result = (1+p1*TMath::Gaus(dminphi,p2,p3)+p4*TMath::Gaus(dminphi,p5,p6)+p7*TMath::Gaus(dminphi,p8,p9));
597 
598  return result;
599 }
600 
601 
602 // corrects the effect of |eta|-cracks
603 double
604 PFEnergyCalibration::CorrEta(double eta){
605 
606  // we use a gaussian with a screwness for each of the 5 |eta|-cracks
607  constexpr double a[] = {6.13349e-01, 5.08146e-01, 4.44480e-01, 3.3487e-01, 7.65627e-01}; // amplitude
608  constexpr double m[] = {-1.79514e-02, 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00}; // mean
609  constexpr double s[] = {7.92382e-03, 3.06028e-03, 3.36139e-03, 3.94521e-03, 8.63950e-04}; // sigma
610  constexpr double sa[] = {1.27228e+01, 3.81517e-02, 1.63507e-01, -6.56480e-02, 1.87160e-01}; // screwness amplitude
611  constexpr double ss[] = {5.48753e-02, -1.00223e-02, 2.22866e-03, 4.26288e-04, 2.67937e-03}; // screwness sigma
612  double result = 1;
613 
614  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]));
615 
616  return result;
617 }
618 
619 
620 //corrects the global behaviour in the barrel
621 double
622 PFEnergyCalibration::CorrBarrel(double E, double eta) {
623 
624  //Energy dependency
625  /*
626  //YM Parameters 52XX:
627  constexpr double p0=1.00000e+00;
628  constexpr double p1=3.27753e+01;
629  constexpr double p2=2.28552e-02;
630  constexpr double p3=3.06139e+00;
631  constexpr double p4=2.25135e-01;
632  constexpr double p5=1.47824e+00;
633  constexpr double p6=1.09e-02;
634  constexpr double p7=4.19343e+01;
635  */
636  constexpr double p0 = 0.9944;
637  constexpr double p1 = 9.827;
638  constexpr double p2 = 1.503;
639  constexpr double p3 = 1.196;
640  constexpr double p4 = 0.3349;
641  constexpr double p5 = 0.89;
642  constexpr double p6 = 0.004361;
643  constexpr double p7 = 51.51;
644  //Eta dependency
645  constexpr double p8=2.705593e-03;
646 
647  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);
648 
649  return result;
650 }
651 
652 
653 
662 
663 
664 //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
665 // Etot = Beta*eEcal + Gamma*(ePS1 + Alpha*ePS2)
666 
667 double
668 PFEnergyCalibration::Alpha(double eta) {
669 
670  //Energy dependency
671  constexpr double p0 = 5.97621e-01;
672 
673  //Eta dependency
674  constexpr double p1 =-1.86407e-01;
675  constexpr double p2 = 3.85197e-01;
676 
677  //so that <feta()> = 1
678  constexpr double norm = (p1+p2*(2.6+1.656)/2);
679 
680  double result = p0*(p1+p2*eta)/norm;
681 
682  return result;
683 }
684 
685 double
686 PFEnergyCalibration::Beta(double E, double eta) {
687 
688  //Energy dependency
689  constexpr double p0 = 0.032;
690  constexpr double p1 = 9.70394e-02;
691  constexpr double p2 = 2.23072e+01;
692  constexpr double p3 = 100;
693 
694  //Eta dependency
695  constexpr double p4 = 1.02496e+00 ;
696  constexpr double p5 = -4.40176e-03 ;
697 
698  //so that <feta()> = 1
699  constexpr double norm = (p4+p5*(2.6+1.656)/2);
700 
701  double result = (1.0012+p0*TMath::Exp(-E/p3)+p1*TMath::Exp(-E/p2))*(p4+p5*eta)/norm;
702  return result;
703 }
704 
705 
706 double
707 PFEnergyCalibration::Gamma(double etaEcal) {
708 
709  //Energy dependency
710  constexpr double p0 = 2.49752e-02;
711 
712  //Eta dependency
713  constexpr double p1 = 6.48816e-02;
714  constexpr double p2 = -1.59517e-02;
715 
716  //so that <feta()> = 1
717  constexpr double norm = (p1+p2*(2.6+1.656)/2);
718 
719  double result = p0*(p1+p2*etaEcal)/norm;
720 
721  return result;
722 }
723 
724 
725 
731 
732 
733 // returns the corrected energy in the barrel (0,1.48)
734 // Global Behaviour, phi and eta cracks are taken into account
735 double
736 PFEnergyCalibration::EcorrBarrel(double E, double eta, double phi,
737  bool crackCorrection ){
738 
739  // double result = E*CorrBarrel(E,eta)*CorrEta(eta)*CorrPhi(phi,eta);
740  double correction = crackCorrection ? std::max(CorrEta(eta),CorrPhi(phi,eta)) : 1.;
741  double result = E * CorrBarrel(E,eta) * correction;
742 
743  return result;
744 }
745 
746 
747 // returns the corrected energy in the area between the barrel and the PS (1.48,1.65)
748 double
749 PFEnergyCalibration::EcorrZoneBeforePS(double E, double eta){
750 
751  //Energy dependency
752  constexpr double p0 =1;
753  constexpr double p1 =0.18;
754  constexpr double p2 =8.;
755 
756  //Eta dependency
757  constexpr double p3 =0.3;
758  constexpr double p4 =1.11;
759  constexpr double p5 =0.025;
760  constexpr double p6 =1.49;
761  constexpr double p7 =0.6;
762 
763  //so that <feta()> = 1
764  constexpr double norm = 1.21;
765 
766  double result = E*(p0+p1*TMath::Exp(-E/p2))*(p3+p4*TMath::Gaus(eta,p6,p5)+p7*eta)/norm;
767 
768  return result;
769 }
770 
771 
772 // returns the corrected energy in the PS (1.65,2.6)
773 // only when (ePS1>0)||(ePS2>0)
774 double
775 PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal) {
776 
777  // gives the good weights to each subdetector
778  double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+Gamma(etaEcal)*(ePS1+Alpha(etaEcal)*ePS2)/9e-5 ;
779 
780  //Correction of the residual energy dependency
781  constexpr double p0 = 1.00;
782  constexpr double p1 = 2.18;
783  constexpr double p2 =1.94;
784  constexpr double p3 =4.13;
785  constexpr double p4 =1.127;
786 
787  double result = E*(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
788 
789  return result;
790 }
791 
792 // returns the corrected energy in the PS (1.65,2.6)
793 // only when (ePS1>0)||(ePS2>0)
794 double
795 PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal,double & outputPS1, double & outputPS2) {
796 
797  // gives the good weights to each subdetector
798  double gammaprime=Gamma(etaEcal)/9e-5;
799  outputPS1=gammaprime*ePS1;
800  outputPS2=gammaprime*Alpha(etaEcal)*ePS2;
801  double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+outputPS1+outputPS2;
802 
803  //Correction of the residual energy dependency
804  constexpr double p0 = 1.00;
805  constexpr double p1 = 2.18;
806  constexpr double p2 =1.94;
807  constexpr double p3 =4.13;
808  constexpr double p4 =1.127;
809 
810  double corrfac=(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
811  outputPS1*=corrfac;
812  outputPS2*=corrfac;
813  double result = E*corrfac;
814 
815  return result;
816 }
817 
818 
819 // returns the corrected energy in the PS (1.65,2.6)
820 // only when (ePS1=0)&&(ePS2=0)
821 double
822 PFEnergyCalibration::EcorrPS_ePSNil(double eEcal,double eta){
823 
824  //Energy dependency
825  constexpr double p0= 1.02;
826  constexpr double p1= 0.165;
827  constexpr double p2= 6.5 ;
828  constexpr double p3= 2.1 ;
829 
830  //Eta dependency
831  constexpr double p4 = 1.02496e+00 ;
832  constexpr double p5 = -4.40176e-03 ;
833 
834  //so that <feta()> = 1
835  constexpr double norm = (p4+p5*(2.6+1.656)/2);
836 
837  double result = eEcal*(p0+p1*TMath::Exp(-TMath::Abs(eEcal-p3)/p2))*(p4+p5*eta)/norm;
838 
839  return result;
840 }
841 
842 
843 // returns the corrected energy in the area between the end of the PS and the end of the endcap (2.6,2.98)
844 double
845 PFEnergyCalibration::EcorrZoneAfterPS(double E, double eta){
846 
847  //Energy dependency
848  constexpr double p0 =1;
849  constexpr double p1 = 0.058;
850  constexpr double p2 =12.5;
851  constexpr double p3 =-1.05444e+00;
852  constexpr double p4 =-5.39557e+00;
853  constexpr double p5 =8.38444e+00;
854  constexpr double p6 = 6.10998e-01 ;
855 
856  //Eta dependency
857  constexpr double p7 =1.06161e+00;
858  constexpr double p8 = 0.41;
859  constexpr double p9 =2.918;
860  constexpr double p10 =0.0181;
861  constexpr double p11= 2.05;
862  constexpr double p12 =2.99;
863  constexpr double p13=0.0287;
864 
865  //so that <feta()> = 1
866  constexpr double norm=1.045;
867 
868  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;
869  return result;
870 }
871 
872 
873 
874 
875 // returns the corrected energy everywhere
876 // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
877 double
878 PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,
879  double eta,double phi,
880  bool crackCorrection ) {
881 
882  constexpr double endBarrel=1.48;
883  constexpr double beginingPS=1.65;
884  constexpr double endPS=2.6;
885  constexpr double endEndCap=2.98;
886 
887  double result=0;
888 
889  eta=TMath::Abs(eta);
890 
891  if(eEcal>0){
892  if(eta <= endBarrel) result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
893  else if(eta <= beginingPS) result = EcorrZoneBeforePS(eEcal,eta);
894  else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
895  else if(eta < endPS) result = EcorrPS(eEcal,ePS1,ePS2,eta);
896  else if(eta < endEndCap) result = EcorrZoneAfterPS(eEcal,eta);
897  else result =eEcal;
898  }
899  else result = eEcal;// useful if eEcal=0 or eta>2.98
900  //protection
901  if(result<eEcal) result=eEcal;
902  return result;
903 }
904 
905 // returns the corrected energy everywhere
906 // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
907 double
908 PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,double eta,double phi,double& ps1,double&ps2,bool crackCorrection) {
909 
910  constexpr double endBarrel=1.48;
911  constexpr double beginingPS=1.65;
912  constexpr double endPS=2.6;
913  constexpr double endEndCap=2.98;
914 
915  double result=0;
916 
917  eta=TMath::Abs(eta);
918 
919  if(eEcal>0){
920  if(eta <= endBarrel) result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
921  else if(eta <= beginingPS) result = EcorrZoneBeforePS(eEcal,eta);
922  else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
923  else if(eta < endPS) result = EcorrPS(eEcal,ePS1,ePS2,eta,ps1,ps2);
924  else if(eta < endEndCap) result = EcorrZoneAfterPS(eEcal,eta);
925  else result =eEcal;
926  }
927  else result = eEcal;// useful if eEcal=0 or eta>2.98
928  // protection
929  if(result<eEcal) result=eEcal;
930  return result;
931 }
PerformancePayloadFromTFormula::getResult
float getResult(PerformanceResult::ResultType, BinningPointByMap) const
Definition: PerformancePayloadFromTFormula.cc:8
i
int i
Definition: DBlmapReader.cc:9
PFEnergyCalibration::Beta
double Beta(double E, double eta)
Definition: PFEnergyCalibration.cc:686
PFEnergyCalibration::EcorrZoneAfterPS
double EcorrZoneAfterPS(double E, double eta)
Definition: PFEnergyCalibration.cc:845
PFEnergyCalibration::Alpha
double Alpha(double eta)
Definition: PFEnergyCalibration.cc:668
PFEnergyCalibration::pfCalibrations
const PerformancePayloadFromTFormula * pfCalibrations
Definition: PFEnergyCalibration.h:83
reco::OK
pair< int, edm::FunctionWithDict > OK
Definition: findMethod.cc:70
PFEnergyCalibration::bBarrel
double bBarrel(double x) const
Definition: PFEnergyCalibration.cc:241
reco::PFCluster
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:42
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:736
PFEnergyCalibration::EcorrZoneBeforePS
double EcorrZoneBeforePS(double E, double eta)
Definition: PFEnergyCalibration.cc:749
PFEnergyCalibration::EcorrPS_ePSNil
double EcorrPS_ePSNil(double eEcal, double eta)
Definition: PFEnergyCalibration.cc:822
PFEnergyCalibration::CorrEta
double CorrEta(double eta)
Definition: PFEnergyCalibration.cc:604
min
#define min(a, b)
Definition: mlp_lapack.h:161
eta
T eta() const
Definition: Basic3DVectorLD.h:179
operator<<
std::ostream & operator<<(std::ostream &out, const ALILine &li)
Definition: ALILine.cc:187
PFEnergyCalibration::energyEmHad
void energyEmHad(double t, double &e, double &h, double eta, double phi) const
Definition: PFEnergyCalibration.cc:108
PFEnergyCalibration::Gamma
double Gamma(double etaEcal)
Definition: PFEnergyCalibration.cc:707
calib
MVATrainerComputer * calib
Definition: MVATrainer.cc:64
PFEnergyCalibration::CorrPhi
double CorrPhi(double phi, double eta)
Definition: PFEnergyCalibration.cc:578
PFEnergyCalibration::CorrBarrel
double CorrBarrel(double E, double eta)
Definition: PFEnergyCalibration.cc:622
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
reco::PFCluster::calculatePositionREP
void calculatePositionREP()
computes posrep_ once and for all
Definition: PFCluster.h:79
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:76
p2
double p2[4]
Definition: TauolaWrapper.h:90
PFEnergyCalibration::dCrackPhi
double dCrackPhi(double phi, double eta)
Definition: PFEnergyCalibration.cc:526
reco::PFCluster::energy
double energy() const
cluster energy
Definition: PFCluster.h:73
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
M_PI
#define M_PI
Definition: BFit3D.cc:3
BinningPointByMap::insert
bool insert(BinningVariables::BinningVariablesType, float)
Definition: BinningPointByMap.cc:3
h
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
b
double b
Definition: hdecay.h:120
PFEnergyCalibration::EcorrPS
double EcorrPS(double eEcal, double ePS1, double ePS2, double etaEcal)
Definition: PFEnergyCalibration.cc:775
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:878
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
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
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
constexpr
#define constexpr
Definition: GCC11Compatibility.h:36
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