CMS 3D CMS Logo

PFEnergyCalibration.cc
Go to the documentation of this file.
1 
4 
6 
7 #include <TMath.h>
8 #include <cmath>
9 #include <vector>
10 #include <TF1.h>
11 #include <map>
12 #include <algorithm>
13 #include <numeric>
14 
15 using namespace std;
16 
17 PFEnergyCalibration::PFEnergyCalibration() : pfCalibrations(nullptr), esEEInterCalib_(nullptr)
18 {
20 }
21 
23 {
24 
25 }
26 
27 void
29 
30  threshE = 3.5;
31  threshH = 2.5;
32 
33 
34  //calibChrisClean.C calibration parameters bhumika Nov, 2018
35  faBarrel = std::make_unique<TF1>("faBarrel","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
36  faBarrel->SetParameter(0,-30.7141);
37  faBarrel->SetParameter(1,31.7583);
38  faBarrel->SetParameter(2,4.40594);
39  faBarrel->SetParameter(3,1.70914);
40  faBarrel->SetParameter(4,0.0613696);
41  faBarrel->SetParameter(5,0.000104857);
42  faBarrel->SetParameter(6,-1.38927);
43  faBarrel->SetParameter(7,-0.743082);
44  fbBarrel = std::make_unique<TF1>("fbBarrel","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
45  fbBarrel->SetParameter(0,2.25366);
46  fbBarrel->SetParameter(1,0.537715);
47  fbBarrel->SetParameter(2,-4.81375);
48  fbBarrel->SetParameter(3,12.109);
49  fbBarrel->SetParameter(4,1.80577);
50  fbBarrel->SetParameter(5,0.187919);
51  fbBarrel->SetParameter(6,-6.26234);
52  fbBarrel->SetParameter(7,-0.607392);
53  fcBarrel = std::make_unique<TF1>("fcBarrel","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
54  fcBarrel->SetParameter(0,1.5125962);
55  fcBarrel->SetParameter(1,0.855057);
56  fcBarrel->SetParameter(2,-6.04199);
57  fcBarrel->SetParameter(3,2.08229);
58  fcBarrel->SetParameter(4,0.592266);
59  fcBarrel->SetParameter(5,0.0291232);
60  fcBarrel->SetParameter(6,0.364802);
61  fcBarrel->SetParameter(7,-1.50142);
62  faEtaBarrelEH = std::make_unique<TF1>("faEtaBarrelEH","[0]+[1]*exp(-x/[2])",0.,1000.);
63  faEtaBarrelEH->SetParameter(0,0.0185555);
64  faEtaBarrelEH->SetParameter(1,-0.0470674);
65  faEtaBarrelEH->SetParameter(2,396.959);
66  fbEtaBarrelEH = std::make_unique<TF1>("fbEtaBarrelEH","[0]+[1]*exp(-x/[2])",0.,1000.);
67  fbEtaBarrelEH->SetParameter(0,0.0396458);
68  fbEtaBarrelEH->SetParameter(1,0.114128);
69  fbEtaBarrelEH->SetParameter(2,251.405);
70  faEtaBarrelH = std::make_unique<TF1>("faEtaBarrelH","[0]+[1]*x",0.,1000.);
71  faEtaBarrelH->SetParameter(0,0.00434994);
72  faEtaBarrelH->SetParameter(1,-5.16564e-06);
73  fbEtaBarrelH = std::make_unique<TF1>("fbEtaBarrelH","[0]+[1]*exp(-x/[2])",0.,1000.);
74  fbEtaBarrelH->SetParameter(0,-0.0232604);
75  fbEtaBarrelH->SetParameter(1,0.0937525);
76  fbEtaBarrelH->SetParameter(2,34.9935);
77 
78  faEndcap = std::make_unique<TF1>("faEndcap","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
79  faEndcap->SetParameter(0,1.17227);
80  faEndcap->SetParameter(1,13.1489);
81  faEndcap->SetParameter(2,-29.1672);
82  faEndcap->SetParameter(3,0.604223);
83  faEndcap->SetParameter(4,0.0426363);
84  faEndcap->SetParameter(5,3.30898e-15);
85  faEndcap->SetParameter(6,0.165293);
86  faEndcap->SetParameter(7,-7.56786);
87  fbEndcap = std::make_unique<TF1>("fbEndcap","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
88  fbEndcap->SetParameter(0,-0.974251);
89  fbEndcap->SetParameter(1,1.61733);
90  fbEndcap->SetParameter(2,0.0629183);
91  fbEndcap->SetParameter(3,7.78495);
92  fbEndcap->SetParameter(4,-0.774289);
93  fbEndcap->SetParameter(5,7.81399e-05);
94  fbEndcap->SetParameter(6,0.139116);
95  fbEndcap->SetParameter(7,-4.25551);
96  fcEndcap = std::make_unique<TF1>("fcEndcap","[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))",0.,1000.);
97  fcEndcap->SetParameter(0,1.01863);
98  fcEndcap->SetParameter(1,1.29787);
99  fcEndcap->SetParameter(2,-3.97293);
100  fcEndcap->SetParameter(3,21.7805);
101  fcEndcap->SetParameter(4,0.810195);
102  fcEndcap->SetParameter(5,0.234134);
103  fcEndcap->SetParameter(6,1.42226);
104  fcEndcap->SetParameter(7,-0.0997326);
105  faEtaEndcapEH = std::make_unique<TF1>("faEtaEndcapEH","[0]+[1]*exp(-x/[2])",0.,1000.);
106  faEtaEndcapEH->SetParameter(0,0.0112692);
107  faEtaEndcapEH->SetParameter(1,-2.68063);
108  faEtaEndcapEH->SetParameter(2,2.90973);
109  fbEtaEndcapEH = std::make_unique<TF1>("fbEtaEndcapEH","[0]+[1]*exp(-x/[2])",0.,1000.);
110  fbEtaEndcapEH->SetParameter(0,-0.0192991);
111  fbEtaEndcapEH->SetParameter(1,-0.265);
112  fbEtaEndcapEH->SetParameter(2,80.5502);
113  faEtaEndcapH = std::make_unique<TF1>("faEtaEndcapH","[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",0.,1000.);
114  faEtaEndcapH->SetParameter(0,-0.0106029);
115  faEtaEndcapH->SetParameter(1,-0.692207);
116  faEtaEndcapH->SetParameter(2,0.0542991);
117  faEtaEndcapH->SetParameter(3,-0.171435);
118  faEtaEndcapH->SetParameter(4,-61.2277);
119  fbEtaEndcapH = std::make_unique<TF1>("fbEtaEndcapH","[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",0.,1000.);
120  fbEtaEndcapH->SetParameter(0,0.0214894);
121  fbEtaEndcapH->SetParameter(1,-0.266704);
122  fbEtaEndcapH->SetParameter(2,5.2112);
123  fbEtaEndcapH->SetParameter(3,0.303578);
124  fbEtaEndcapH->SetParameter(4,-104.367);
125 
126  //added by Bhumika on 2 august 2018
127 
128  fcEtaBarrelH = std::make_unique<TF1>("fcEtaBarrelH","[3]*((x-[0])^[1])+[2]",0.,1000.);
129  fcEtaBarrelH->SetParameter(0,0);
130  fcEtaBarrelH->SetParameter(1,2);
131  fcEtaBarrelH->SetParameter(2,0);
132  fcEtaBarrelH->SetParameter(3,1);
133 
134  fcEtaEndcapH = std::make_unique<TF1>("fcEtaEndcapH","[3]*((x-[0])^[1])+[2]",0.,1000.);
135  fcEtaEndcapH->SetParameter(0,0);
136  fcEtaEndcapH->SetParameter(1,0);
137  fcEtaEndcapH->SetParameter(2,0.05);
138  fcEtaEndcapH->SetParameter(3,0);
139 
140  fdEtaEndcapH = std::make_unique<TF1>("fdEtaEndcapH","[3]*((x-[0])^[1])+[2]",0.,1000.);
141  fdEtaEndcapH->SetParameter(0,1.5);
142  fdEtaEndcapH->SetParameter(1,4);
143  fdEtaEndcapH->SetParameter(2,-1.1);
144  fdEtaEndcapH->SetParameter(3,1.0);
145 
146  fcEtaBarrelEH = std::make_unique<TF1>("fcEtaBarrelEH","[3]*((x-[0])^[1])+[2]",0.,1000.);
147  fcEtaBarrelEH->SetParameter(0,0);
148  fcEtaBarrelEH->SetParameter(1,2);
149  fcEtaBarrelEH->SetParameter(2,0);
150  fcEtaBarrelEH->SetParameter(3,1);
151 
152  fcEtaEndcapEH = std::make_unique<TF1>("fcEtaEndcapEH","[3]*((x-[0])^[1])+[2]",0.,1000.);
153  fcEtaEndcapEH->SetParameter(0,0);
154  fcEtaEndcapEH->SetParameter(1,0);
155  fcEtaEndcapEH->SetParameter(2,0);
156  fcEtaEndcapEH->SetParameter(3,0);
157 
158  fdEtaEndcapEH = std::make_unique<TF1>("fdEtaEndcapEH","[3]*((x-[0])^[1])+[2]",0.,1000.);
159  fdEtaEndcapEH->SetParameter(0,1.5);
160  fdEtaEndcapEH->SetParameter(1,2.0);
161  fdEtaEndcapEH->SetParameter(2,0.6);
162  fdEtaEndcapEH->SetParameter(3,1.0);
163 
164 
165 
166 }
167 
168 void
169 PFEnergyCalibration::energyEmHad(double t, double& e, double&h, double eta, double phi) const {
170 
171  // Use calorimetric energy as true energy for neutral particles
172  double tt = t;
173  double ee = e;
174  double hh = h;
175  double a = 1.;
176  double b = 1.;
177  double etaCorrE = 1.;
178  double etaCorrH = 1.;
179  auto absEta = std::abs(eta);
180  t = min(999.9,max(tt,e+h));
181  if ( t < 1. ) return;
182 
183  // Barrel calibration
184  if ( absEta < 1.48 ) {
185  // The energy correction
186  a = e>0. ? aBarrel(t) : 1.;
187  b = e>0. ? bBarrel(t) : cBarrel(t);
188  double thresh = e > 0. ? threshE : threshH;
189 
190  // Protection against negative calibration
191  if ( a < -0.25 || b < -0.25 ) {
192  a = 1.;
193  b = 1.;
194  thresh = 0.;
195  }
196 
197  // The new estimate of the true energy
198  t = min(999.9,max(tt, thresh+a*e+b*h));
199 
200  // The angular correction
201  if ( e > 0. && thresh > 0. ) {
202  etaCorrE = 1.0 + aEtaBarrelEH(t) + 1.3*bEtaBarrelEH(t)*cEtaBarrelEH(absEta);
203  etaCorrH = 1.0;
204  } else {
205  etaCorrE = 1.0 + aEtaBarrelH(t) + 1.3*bEtaBarrelH(t)*cEtaBarrelH(absEta);
206  etaCorrH = 1.0 + aEtaBarrelH(t) + bEtaBarrelH(t)*cEtaBarrelH(absEta);
207 
208 
209  }
210  if ( e > 0. && thresh > 0. )
211  e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
212  if ( h > 0. && thresh > 0. ) {
213  h = threshH + etaCorrH * b * h;
214  }
215 
216  // Endcap calibration
217  } else {
218  // The energy correction
219  a = e>0. ? aEndcap(t) : 1.;
220  b = e>0. ? bEndcap(t) : cEndcap(t);
221  double thresh = e > 0. ? threshE : threshH;
222 
223  // Protection against negative calibration
224  if ( a < -0.25 || b < -0.25 ) {
225  a = 1.;
226  b = 1.;
227  thresh = 0.;
228  }
229 
230  // The new estimate of the true energy
231  t = min(999.9,max(tt, thresh+a*e+b*h));
232 
233  // The angular correction
234  double dEta = std::abs( absEta - 1.5 );
235  double etaPow = dEta * dEta * dEta * dEta;
236 
237 
238  if ( e > 0. && thresh > 0. ) {
239  if(absEta<2.5) {
240  etaCorrE = 1. + aEtaEndcapEH(t) + bEtaEndcapEH(t)*cEtaEndcapEH(absEta);
241  }
242  else {
243 
244  etaCorrE = 1. + aEtaEndcapEH(t) + 1.3*bEtaEndcapEH(t)*dEtaEndcapEH(absEta);
245  }
246 
247  etaCorrH = 1. + aEtaEndcapEH(t) + bEtaEndcapEH(t)*(0.04 + etaPow);
248  } else {
249  etaCorrE = 1.;
250  if(absEta<2.5) {
251  etaCorrH = 1. + aEtaEndcapH(t) + bEtaEndcapH(t)*cEtaEndcapH(absEta);
252  }
253  else {
254  etaCorrH = 1. + aEtaEndcapH(t) + bEtaEndcapH(t)*dEtaEndcapH(absEta);
255  }
256  }
257 
258  //t = min(999.9,max(tt, thresh + etaCorrE*a*e + etaCorrH*b*h));
259 
260  if ( e > 0. && thresh > 0. )
261  e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
262  if ( h > 0. && thresh > 0. ) {
263  h = threshH + etaCorrH * b * h;
264  }
265  }
266 
267  // Protection
268  if ( e < 0. || h < 0. ) {
269 
270  // Some protection against crazy calibration
271  if ( e < 0. ) e = ee;
272  if ( h < 0. ) h = hh;
273  }
274 
275  // And that's it !
276 
277 
278 }
279 
280 // The calibration functions
281 double
283 
284  if ( pfCalibrations ) {
288 
289  } else {
290 
291  return faBarrel->Eval(x);
292 
293  }
294 }
295 
296 double
298 
299  if ( pfCalibrations ) {
300 
304 
305  } else {
306 
307  return fbBarrel->Eval(x);
308 
309  }
310 }
311 
312 double
314 
315  if ( pfCalibrations ) {
316 
320 
321  } else {
322 
323  return fcBarrel->Eval(x);
324 
325  }
326 }
327 
328 double
330 
331  if ( pfCalibrations ) {
332 
336 
337  } else {
338 
339  return faEtaBarrelEH->Eval(x);
340 
341  }
342 }
343 
344 double
346 
347  if ( pfCalibrations ) {
348 
352 
353  } else {
354 
355  return fbEtaBarrelEH->Eval(x);
356 
357  }
358 }
359 
360 double
362 
363  if ( pfCalibrations ) {
364 
368 
369  } else {
370 
371  return faEtaBarrelH->Eval(x);
372 
373  }
374 }
375 
376 double
378 
379  if ( pfCalibrations ) {
380 
384 
385  } else {
386 
387  return fbEtaBarrelH->Eval(x);
388 
389  }
390 }
391 
392 double
394 
395  if ( pfCalibrations ) {
396 
400 
401  } else {
402 
403  return faEndcap->Eval(x);
404 
405  }
406 }
407 
408 double
410 
411  if ( pfCalibrations ) {
412 
416 
417  } else {
418 
419  return fbEndcap->Eval(x);
420 
421  }
422 }
423 
424 double
426 
427  if ( pfCalibrations ) {
428 
432 
433  } else {
434 
435  return fcEndcap->Eval(x);
436 
437  }
438 }
439 
440 double
442 
443  if ( pfCalibrations ) {
444 
448 
449  } else {
450 
451  return faEtaEndcapEH->Eval(x);
452 
453  }
454 }
455 
456 double
458 
459  if ( pfCalibrations ) {
460 
464 
465  } else {
466 
467  return fbEtaEndcapEH->Eval(x);
468 
469  }
470 }
471 
472 double
474 
475  if ( pfCalibrations ) {
476 
480 
481  } else {
482 
483  return faEtaEndcapH->Eval(x);
484 
485  }
486 }
487 
488 double
490 
491  if ( pfCalibrations ) {
492 
496  } else {
497  return fbEtaEndcapH->Eval(x);
498 
499  }
500 }
501 
502 //added by Bhumika Kansal on 3 august 2018
503 
504 double
506  if ( pfCalibrations ) {
510 
511  } else {
512  return fcEtaBarrelH->Eval(x);
513  }
514 }
515 double
517  if ( pfCalibrations ) {
521 
522  } else{
523 
524  return fcEtaEndcapH->Eval(x);
525  }
526 }
527 
528 double
530  if ( pfCalibrations ) {
534 
535  } else{
536 
537  return fdEtaEndcapH->Eval(x);
538  }
539 }
540 
541 double
543  if ( pfCalibrations ) {
547 
548  } else{
549 
550  return fcEtaBarrelEH->Eval(x);
551  }
552 }
553 
554 double
556  if ( pfCalibrations ) {
560 
561  } else{
562 
563  return fcEtaEndcapEH->Eval(x);
564  }
565 }
566 
567 double
569  if ( pfCalibrations ) {
573 
574  } else{
575 
576  return fdEtaEndcapEH->Eval(x);
577  }
578 }
579 //
580 double
582  std::vector<double> &EclustersPS1,
583  std::vector<double> &EclustersPS2,
584  bool crackCorrection ) const {
585  double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
586  double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
587  return energyEm(clusterEcal, ePS1, ePS2, crackCorrection);
588 }
589 
590 double
592  double ePS1,
593  double ePS2,
594  bool crackCorrection ) const {
595  double eEcal = clusterEcal.energy();
596  //temporaty ugly fix
597  reco::PFCluster myPFCluster=clusterEcal;
598  myPFCluster.calculatePositionREP();
599  double eta = myPFCluster.positionREP().eta();
600  double phi = myPFCluster.positionREP().phi();
601 
602  double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi, crackCorrection);
603  // if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
604  return calibrated;
605 }
606 
608  std::vector<double> &EclustersPS1,
609  std::vector<double> &EclustersPS2,
610  double& ps1,double& ps2,
611  bool crackCorrection) const {
612  double ePS1(std::accumulate(EclustersPS1.begin(), EclustersPS1.end(), 0.0));
613  double ePS2(std::accumulate(EclustersPS2.begin(), EclustersPS2.end(), 0.0));
614  return energyEm(clusterEcal, ePS1, ePS2, ps1, ps2, crackCorrection);
615 }
617  double ePS1, double ePS2,
618  double& ps1,double& ps2,
619  bool crackCorrection) const {
620  double eEcal = clusterEcal.energy();
621  //temporaty ugly fix
622  reco::PFCluster myPFCluster=clusterEcal;
623  myPFCluster.calculatePositionREP();
624  double eta = myPFCluster.positionREP().eta();
625  double phi = myPFCluster.positionREP().phi();
626 
627  double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi,ps1,ps2,crackCorrection);
628  // if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
629  return calibrated;
630 }
631 
632 
633 std::ostream& operator<<(std::ostream& out,
634  const PFEnergyCalibration& calib) {
635 
636  if(!out ) return out;
637 
638  out<<"PFEnergyCalibration -- "<<endl;
639 
640  if ( calib.pfCalibrations ) {
641 
642  static const std::map<std::string, PerformanceResult::ResultType> functType = {
643  {"PFfa_BARREL", PerformanceResult::PFfa_BARREL},
644  {"PFfa_ENDCAP", PerformanceResult::PFfa_ENDCAP},
645  {"PFfb_BARREL", PerformanceResult::PFfb_BARREL},
646  {"PFfb_ENDCAP", PerformanceResult::PFfb_ENDCAP},
647  {"PFfc_BARREL", PerformanceResult::PFfc_BARREL},
648  {"PFfc_ENDCAP", PerformanceResult::PFfc_ENDCAP},
649  {"PFfaEta_BARRELH", PerformanceResult::PFfaEta_BARRELH},
650  {"PFfaEta_ENDCAPH", PerformanceResult::PFfaEta_ENDCAPH},
651  {"PFfbEta_BARRELH", PerformanceResult::PFfbEta_BARRELH},
652  {"PFfbEta_ENDCAPH", PerformanceResult::PFfbEta_ENDCAPH},
653  {"PFfaEta_BARRELEH", PerformanceResult::PFfaEta_BARRELEH},
654  {"PFfaEta_ENDCAPEH", PerformanceResult::PFfaEta_ENDCAPEH},
655  {"PFfbEta_BARRELEH", PerformanceResult::PFfbEta_BARRELEH},
656  {"PFfbEta_ENDCAPEH", PerformanceResult::PFfbEta_ENDCAPEH},
657  {"PFfcEta_BARRELH", PerformanceResult::PFfcEta_BARRELH},
658  {"PFfcEta_ENDCAPH", PerformanceResult::PFfcEta_ENDCAPH},
659  {"PFfdEta_ENDCAPH", PerformanceResult::PFfdEta_ENDCAPH},
660  {"PFfcEta_BARRELEH", PerformanceResult::PFfcEta_BARRELEH},
661  {"PFfcEta_ENDCAPEH", PerformanceResult::PFfcEta_ENDCAPEH},
662  {"PFfdEta_ENDCAPEH", PerformanceResult::PFfdEta_ENDCAPEH}
663 
664  };
665 
666 
667  for(std::map<std::string,PerformanceResult::ResultType>::const_iterator
668  func = functType.begin();
669  func != functType.end();
670  ++func) {
671 
672  cout << "Function: " << func->first << endl;
673  PerformanceResult::ResultType fType = func->second;
674  calib.pfCalibrations->printFormula(fType);
675  }
676 
677  } else {
678 
679  std::cout << "Default calibration functions : " << std::endl;
680 
681  calib.faBarrel->Print();
682  calib.fbBarrel->Print();
683  calib.fcBarrel->Print();
684  calib.faEtaBarrelEH->Print();
685  calib.fbEtaBarrelEH->Print();
686  calib.faEtaBarrelH->Print();
687  calib.fbEtaBarrelH->Print();
688  calib.faEndcap->Print();
689  calib.fbEndcap->Print();
690  calib.fcEndcap->Print();
691  calib.faEtaEndcapEH->Print();
692  calib.fbEtaEndcapEH->Print();
693  calib.faEtaEndcapH->Print();
694  calib.fbEtaEndcapH->Print();
695  //
696 
697  }
698 
699  return out;
700 }
701 
702 
703 
704 
717 
718 
719 
725 
726 
727 //useful to compute the signed distance to the closest crack in the barrel
728 double
729 PFEnergyCalibration::minimum(double a,double b) const {
730  if(TMath::Abs(b)<TMath::Abs(a)) a=b;
731  return a;
732 }
733 
734 namespace {
735  constexpr double pi= M_PI;// 3.14159265358979323846;
736 
737 
738  std::vector<double> fillcPhi() {
739  std::vector<double> retValue;
740  retValue.resize(18,0);
741  retValue[0]=2.97025;
742  for(unsigned i=1;i<=17;++i) retValue[i]=retValue[0]-2*i*pi/18;
743 
744  return retValue;
745  }
746 
747  //Location of the 18 phi-cracks
748  const std::vector<double> cPhi = fillcPhi();
749 }
750 
751 //compute the unsigned distance to the closest phi-crack in the barrel
752 double
753 PFEnergyCalibration::dCrackPhi(double phi, double eta) const {
754 
755 
756  //Shift of this location if eta<0
757  constexpr double delta_cPhi=0.00638;
758 
759  double m; //the result
760 
761  //the location is shifted
762  if(eta<0) phi +=delta_cPhi;
763 
764  if (phi>=-pi && phi<=pi){
765 
766  //the problem of the extrema
767  if (phi<cPhi[17] || phi>=cPhi[0]){
768  if (phi<0) phi+= 2*pi;
769  m = minimum(phi -cPhi[0],phi-cPhi[17]-2*pi);
770  }
771 
772  //between these extrema...
773  else{
774  bool OK = false;
775  unsigned i=16;
776  while(!OK){
777  if (phi<cPhi[i]){
778  m=minimum(phi-cPhi[i+1],phi-cPhi[i]);
779  OK=true;
780  }
781  else i-=1;
782  }
783  }
784  }
785  else{
786  m=0.; //if there is a problem, we assum that we are in a crack
787  std::cout<<"Problem in dminphi"<<std::endl;
788  }
789  if(eta<0) m=-m; //because of the disymetry
790  return m;
791 }
792 
793 // corrects the effect of phi-cracks
794 double
795 PFEnergyCalibration::CorrPhi(double phi, double eta) const {
796 
797  // we use 3 gaussians to correct the phi-cracks effect
798  constexpr double p1= 5.59379e-01;
799  constexpr double p2= -1.26607e-03;
800  constexpr double p3= 9.61133e-04;
801 
802  constexpr double p4= 1.81691e-01;
803  constexpr double p5= -4.97535e-03;
804  constexpr double p6= 1.31006e-03;
805 
806  constexpr double p7= 1.38498e-01;
807  constexpr double p8= 1.18599e-04;
808  constexpr double p9= 2.01858e-03;
809 
810 
811  double dminphi = dCrackPhi(phi,eta);
812 
813  double result = (1+p1*TMath::Gaus(dminphi,p2,p3)+p4*TMath::Gaus(dminphi,p5,p6)+p7*TMath::Gaus(dminphi,p8,p9));
814 
815  return result;
816 }
817 
818 
819 // corrects the effect of |eta|-cracks
820 double
822 
823  // we use a gaussian with a screwness for each of the 5 |eta|-cracks
824  constexpr double a[] = {6.13349e-01, 5.08146e-01, 4.44480e-01, 3.3487e-01, 7.65627e-01}; // amplitude
825  constexpr double m[] = {-1.79514e-02, 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00}; // mean
826  constexpr double s[] = {7.92382e-03, 3.06028e-03, 3.36139e-03, 3.94521e-03, 8.63950e-04}; // sigma
827  constexpr double sa[] = {1.27228e+01, 3.81517e-02, 1.63507e-01, -6.56480e-02, 1.87160e-01}; // screwness amplitude
828  constexpr double ss[] = {5.48753e-02, -1.00223e-02, 2.22866e-03, 4.26288e-04, 2.67937e-03}; // screwness sigma
829  double result = 1;
830 
831  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]));
832 
833  return result;
834 }
835 
836 
837 //corrects the global behaviour in the barrel
838 double
839 PFEnergyCalibration::CorrBarrel(double E, double eta) const {
840 
841  //Energy dependency
842  /*
843  //YM Parameters 52XX:
844  constexpr double p0=1.00000e+00;
845  constexpr double p1=3.27753e+01;
846  constexpr double p2=2.28552e-02;
847  constexpr double p3=3.06139e+00;
848  constexpr double p4=2.25135e-01;
849  constexpr double p5=1.47824e+00;
850  constexpr double p6=1.09e-02;
851  constexpr double p7=4.19343e+01;
852  */
853  constexpr double p0 = 0.9944;
854  constexpr double p1 = 9.827;
855  constexpr double p2 = 1.503;
856  constexpr double p3 = 1.196;
857  constexpr double p4 = 0.3349;
858  constexpr double p5 = 0.89;
859  constexpr double p6 = 0.004361;
860  constexpr double p7 = 51.51;
861  //Eta dependency
862  constexpr double p8=2.705593e-03;
863 
864  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);
865 
866  return result;
867 }
868 
869 
870 
879 
880 
881 //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
882 // Etot = Beta*eEcal + Gamma*(ePS1 + Alpha*ePS2)
883 
884 double
886 
887  //Energy dependency
888  constexpr double p0 = 5.97621e-01;
889 
890  //Eta dependency
891  constexpr double p1 =-1.86407e-01;
892  constexpr double p2 = 3.85197e-01;
893 
894  //so that <feta()> = 1
895  constexpr double norm = (p1+p2*(2.6+1.656)/2);
896 
897  double result = p0*(p1+p2*eta)/norm;
898 
899  return result;
900 }
901 
902 double
903 PFEnergyCalibration::Beta(double E, double eta) const {
904 
905  //Energy dependency
906  constexpr double p0 = 0.032;
907  constexpr double p1 = 9.70394e-02;
908  constexpr double p2 = 2.23072e+01;
909  constexpr double p3 = 100;
910 
911  //Eta dependency
912  constexpr double p4 = 1.02496e+00 ;
913  constexpr double p5 = -4.40176e-03 ;
914 
915  //so that <feta()> = 1
916  constexpr double norm = (p4+p5*(2.6+1.656)/2);
917 
918  double result = (1.0012+p0*TMath::Exp(-E/p3)+p1*TMath::Exp(-E/p2))*(p4+p5*eta)/norm;
919  return result;
920 }
921 
922 
923 double
924 PFEnergyCalibration::Gamma(double etaEcal) const {
925 
926  //Energy dependency
927  constexpr double p0 = 2.49752e-02;
928 
929  //Eta dependency
930  constexpr double p1 = 6.48816e-02;
931  constexpr double p2 = -1.59517e-02;
932 
933  //so that <feta()> = 1
934  constexpr double norm = (p1+p2*(2.6+1.656)/2);
935 
936  double result = p0*(p1+p2*etaEcal)/norm;
937 
938  return result;
939 }
940 
941 
942 
948 
949 
950 // returns the corrected energy in the barrel (0,1.48)
951 // Global Behaviour, phi and eta cracks are taken into account
952 double
953 PFEnergyCalibration::EcorrBarrel(double E, double eta, double phi,
954  bool crackCorrection ) const {
955 
956  // double result = E*CorrBarrel(E,eta)*CorrEta(eta)*CorrPhi(phi,eta);
957  double correction = crackCorrection ? std::max(CorrEta(eta),CorrPhi(phi,eta)) : 1.;
958  double result = E * CorrBarrel(E,eta) * correction;
959 
960  return result;
961 }
962 
963 
964 // returns the corrected energy in the area between the barrel and the PS (1.48,1.65)
965 double
967 
968  //Energy dependency
969  constexpr double p0 =1;
970  constexpr double p1 =0.18;
971  constexpr double p2 =8.;
972 
973  //Eta dependency
974  constexpr double p3 =0.3;
975  constexpr double p4 =1.11;
976  constexpr double p5 =0.025;
977  constexpr double p6 =1.49;
978  constexpr double p7 =0.6;
979 
980  //so that <feta()> = 1
981  constexpr double norm = 1.21;
982 
983  double result = E*(p0+p1*TMath::Exp(-E/p2))*(p3+p4*TMath::Gaus(eta,p6,p5)+p7*eta)/norm;
984 
985  return result;
986 }
987 
988 
989 // returns the corrected energy in the PS (1.65,2.6)
990 // only when (ePS1>0)||(ePS2>0)
991 double
992 PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal) const {
993 
994  // gives the good weights to each subdetector
995  double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+Gamma(etaEcal)*(ePS1+Alpha(etaEcal)*ePS2)/9e-5 ;
996 
997  //Correction of the residual energy dependency
998  constexpr double p0 = 1.00;
999  constexpr double p1 = 2.18;
1000  constexpr double p2 =1.94;
1001  constexpr double p3 =4.13;
1002  constexpr double p4 =1.127;
1003 
1004  double result = E*(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
1005 
1006  return result;
1007 }
1008 
1009 // returns the corrected energy in the PS (1.65,2.6)
1010 // only when (ePS1>0)||(ePS2>0)
1011 double
1012 PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal,double & outputPS1, double & outputPS2) const {
1013 
1014  // gives the good weights to each subdetector
1015  double gammaprime=Gamma(etaEcal)/9e-5;
1016 
1017  if(outputPS1 == 0 && outputPS2 == 0 && esEEInterCalib_ != nullptr){
1018  // both ES planes working
1019  // scaling factor accounting for data-mc
1020  outputPS1=gammaprime*ePS1 * esEEInterCalib_->getGammaLow0();
1021  outputPS2=gammaprime*Alpha(etaEcal)*ePS2 * esEEInterCalib_->getGammaLow3();
1022  }
1023  else if(outputPS1 == 0 && outputPS2 == -1 && esEEInterCalib_ != nullptr){
1024  // ESP1 only working
1025  double corrTotES = gammaprime*ePS1 * esEEInterCalib_->getGammaLow0() * esEEInterCalib_->getGammaLow1();
1026  outputPS1 = gammaprime*ePS1 * esEEInterCalib_->getGammaLow0();
1027  outputPS2 = corrTotES - outputPS1;
1028  }
1029  else if(outputPS1 == -1 && outputPS2 == 0 && esEEInterCalib_ != nullptr){
1030  // ESP2 only working
1031  double corrTotES = gammaprime*Alpha(etaEcal)*ePS2 * esEEInterCalib_->getGammaLow3() * esEEInterCalib_->getGammaLow2();
1032  outputPS2 = gammaprime*Alpha(etaEcal)*ePS2 * esEEInterCalib_->getGammaLow3();
1033  outputPS1 = corrTotES - outputPS2;
1034  }
1035  else{
1036  // none working
1037  outputPS1 = gammaprime*ePS1;
1038  outputPS2 = gammaprime*Alpha(etaEcal)*ePS2;
1039  }
1040 
1041  double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+outputPS1+outputPS2;
1042 
1043  //Correction of the residual energy dependency
1044  constexpr double p0 = 1.00;
1045  constexpr double p1 = 2.18;
1046  constexpr double p2 =1.94;
1047  constexpr double p3 =4.13;
1048  constexpr double p4 =1.127;
1049 
1050  double corrfac=(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
1051  outputPS1*=corrfac;
1052  outputPS2*=corrfac;
1053  double result = E*corrfac;
1054 
1055  return result;
1056 }
1057 
1058 
1059 // returns the corrected energy in the PS (1.65,2.6)
1060 // only when (ePS1=0)&&(ePS2=0)
1061 double
1062 PFEnergyCalibration::EcorrPS_ePSNil(double eEcal,double eta) const {
1063 
1064  //Energy dependency
1065  constexpr double p0= 1.02;
1066  constexpr double p1= 0.165;
1067  constexpr double p2= 6.5 ;
1068  constexpr double p3= 2.1 ;
1069 
1070  //Eta dependency
1071  constexpr double p4 = 1.02496e+00 ;
1072  constexpr double p5 = -4.40176e-03 ;
1073 
1074  //so that <feta()> = 1
1075  constexpr double norm = (p4+p5*(2.6+1.656)/2);
1076 
1077  double result = eEcal*(p0+p1*TMath::Exp(-TMath::Abs(eEcal-p3)/p2))*(p4+p5*eta)/norm;
1078 
1079  return result;
1080 }
1081 
1082 
1083 // returns the corrected energy in the area between the end of the PS and the end of the endcap (2.6,2.98)
1084 double
1086 
1087  //Energy dependency
1088  constexpr double p0 =1;
1089  constexpr double p1 = 0.058;
1090  constexpr double p2 =12.5;
1091  constexpr double p3 =-1.05444e+00;
1092  constexpr double p4 =-5.39557e+00;
1093  constexpr double p5 =8.38444e+00;
1094  constexpr double p6 = 6.10998e-01 ;
1095 
1096  //Eta dependency
1097  constexpr double p7 =1.06161e+00;
1098  constexpr double p8 = 0.41;
1099  constexpr double p9 =2.918;
1100  constexpr double p10 =0.0181;
1101  constexpr double p11= 2.05;
1102  constexpr double p12 =2.99;
1103  constexpr double p13=0.0287;
1104 
1105  //so that <feta()> = 1
1106  constexpr double norm=1.045;
1107 
1108  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;
1109  return result;
1110 }
1111 
1112 
1113 
1114 
1115 // returns the corrected energy everywhere
1116 // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
1117 double
1118 PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,
1119  double eta,double phi,
1120  bool crackCorrection ) const {
1121 
1122  constexpr double endBarrel=1.48;
1123  constexpr double beginingPS=1.65;
1124  constexpr double endPS=2.6;
1125  constexpr double endEndCap=2.98;
1126 
1127  double result=0;
1128 
1129  eta=TMath::Abs(eta);
1130 
1131  if(eEcal>0){
1132  if(eta <= endBarrel) result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
1133  else if(eta <= beginingPS) result = EcorrZoneBeforePS(eEcal,eta);
1134  else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
1135  else if(eta < endPS) result = EcorrPS(eEcal,ePS1,ePS2,eta);
1136  else if(eta < endEndCap) result = EcorrZoneAfterPS(eEcal,eta);
1137  else result =eEcal;
1138  }
1139  else result = eEcal;// useful if eEcal=0 or eta>2.98
1140  //protection
1141  if(result<eEcal) result=eEcal;
1142  return result;
1143 }
1144 
1145 // returns the corrected energy everywhere
1146 // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
1147 double
1148 PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,double eta,double phi,double& ps1,double&ps2,bool crackCorrection) const {
1149 
1150  constexpr double endBarrel=1.48;
1151  constexpr double beginingPS=1.65;
1152  constexpr double endPS=2.6;
1153  constexpr double endEndCap=2.98;
1154 
1155  double result=0;
1156 
1157  eta=TMath::Abs(eta);
1158 
1159  if(eEcal>0){
1160  if(eta <= endBarrel) result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
1161  else if(eta <= beginingPS) result = EcorrZoneBeforePS(eEcal,eta);
1162  else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
1163  else if(eta < endPS) result = EcorrPS(eEcal,ePS1,ePS2,eta,ps1,ps2);
1164  else if(eta < endEndCap) result = EcorrZoneAfterPS(eEcal,eta);
1165  else result =eEcal;
1166  }
1167  else result = eEcal;// useful if eEcal=0 or eta>2.98
1168  // protection
1169  if(result<eEcal) result=eEcal;
1170  return result;
1171 }
std::unique_ptr< TF1 > faEtaBarrelH
std::unique_ptr< TF1 > fcEtaEndcapEH
const PerformancePayloadFromTFormula * pfCalibrations
double aEtaEndcapEH(double x) const
float getResult(PerformanceResult::ResultType, const BinningPointByMap &) const override
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
std::unique_ptr< TF1 > faEtaBarrelEH
double bBarrel(double x) const
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:47
double cEtaEndcapEH(double x) const
T Sign(T A, T B)
Definition: MathUtil.h:54
double aEtaBarrelEH(double x) const
double aEndcap(double x) const
double bEndcap(double x) const
const ESEEIntercalibConstants * esEEInterCalib_
double bEtaEndcapH(double x) const
#define nullptr
double Ecorr(double eEcal, double ePS1, double ePS2, double eta, double phi, bool crackCorrection=true) const
std::unique_ptr< TF1 > fbEtaBarrelEH
double aEtaBarrelH(double x) const
double EcorrZoneBeforePS(double E, double eta) const
double cEtaBarrelH(double x) const
double Alpha(double eta) const
void energyEmHad(double t, double &e, double &h, double eta, double phi) const
const Double_t pi
double cEtaBarrelEH(double x) const
std::unique_ptr< TF1 > fcEtaBarrelEH
void calculatePositionREP()
computes posrep_ once and for all
Definition: PFCluster.h:100
double bEtaBarrelEH(double x) const
double p4[4]
Definition: TauolaWrapper.h:92
std::unique_ptr< TF1 > fbBarrel
std::unique_ptr< TF1 > faEtaEndcapEH
friend std::ostream & operator<<(std::ostream &out, const PFEnergyCalibration &calib)
double cEndcap(double x) const
const REPPoint & positionREP() const
cluster position: rho, eta, phi
Definition: PFCluster.h:97
T Abs(T a)
Definition: MathUtil.h:49
double Beta(double E, double eta) const
std::unique_ptr< TF1 > faBarrel
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::unique_ptr< TF1 > faEndcap
T min(T a, T b)
Definition: MathUtil.h:58
double dEtaEndcapEH(double x) const
std::unique_ptr< TF1 > fcEtaBarrelH
double p2[4]
Definition: TauolaWrapper.h:90
double energy() const
cluster energy
Definition: PFCluster.h:82
#define M_PI
double bEtaEndcapEH(double x) const
double EcorrZoneAfterPS(double E, double eta) const
std::pair< int, edm::FunctionWithDict > OK
Definition: findMethod.cc:136
bool insert(BinningVariables::BinningVariablesType, float)
double bEtaBarrelH(double x) const
double dEtaEndcapH(double x) const
std::unique_ptr< TF1 > fbEtaBarrelH
std::unique_ptr< TF1 > fdEtaEndcapEH
double b
Definition: hdecay.h:120
std::unique_ptr< TF1 > faEtaEndcapH
double CorrBarrel(double E, double eta) const
double Gamma(double etaEcal) const
double energyEm(const reco::PFCluster &clusterEcal, std::vector< double > &EclustersPS1, std::vector< double > &EclustersPS2, bool crackCorrection=true) const
double dCrackPhi(double phi, double eta) const
double p1[4]
Definition: TauolaWrapper.h:89
double a
Definition: hdecay.h:121
std::unique_ptr< TF1 > fcEtaEndcapH
std::unique_ptr< TF1 > fbEndcap
std::unique_ptr< TF1 > fcBarrel
void printFormula(PerformanceResult::ResultType res) const
std::unique_ptr< TF1 > fdEtaEndcapH
double CorrEta(double eta) const
double aEtaEndcapH(double x) const
double aBarrel(double x) const
double CorrPhi(double phi, double eta) const
double EcorrPS_ePSNil(double eEcal, double eta) const
std::unique_ptr< TF1 > fbEtaEndcapEH
double cBarrel(double x) const
double cEtaEndcapH(double x) const
double EcorrBarrel(double E, double eta, double phi, bool crackCorrection=true) const
std::unique_ptr< TF1 > fcEndcap
#define constexpr
*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
std::unique_ptr< TF1 > fbEtaEndcapH
double EcorrPS(double eEcal, double ePS1, double ePS2, double etaEcal) const
double p3[4]
Definition: TauolaWrapper.h:91
double minimum(double a, double b) const