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