00001 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
00002 #include "CondFormats/PhysicsToolsObjects/interface/PerformancePayloadFromTFormula.h"
00003 #include <TMath.h>
00004 #include <math.h>
00005 #include <vector>
00006 #include <TF1.h>
00007 #include <map>
00008
00009 using namespace std;
00010
00011 PFEnergyCalibration::PFEnergyCalibration() : pfCalibrations(0)
00012 {
00013 initializeCalibrationFunctions();
00014 }
00015
00016 PFEnergyCalibration::~PFEnergyCalibration()
00017 {
00018
00019 delete faBarrel;
00020 delete fbBarrel;
00021 delete fcBarrel;
00022 delete faEtaBarrel;
00023 delete fbEtaBarrel;
00024 delete faEndcap;
00025 delete fbEndcap;
00026 delete fcEndcap;
00027 delete faEtaEndcap;
00028 delete fbEtaEndcap;
00029
00030 }
00031
00032 void
00033 PFEnergyCalibration::initializeCalibrationFunctions() {
00034
00035
00036
00037 threshE = 3.5;
00038 threshH = 2.5;
00039
00040
00041 faBarrel = new TF1("faBarrel","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
00042 fbBarrel = new TF1("fbBarrel","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
00043 fcBarrel = new TF1("fcBarrel","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
00044 faEtaBarrel = new TF1("faEtaBarrel","[0]+[1]*exp(-x/[2])",1.,1000.);
00045 fbEtaBarrel = new TF1("fbEtaBarrel","[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",1.,1000.);
00046 faBarrel->SetParameter(0,1.15665);
00047 fbBarrel->SetParameter(0,0.994603);
00048 fcBarrel->SetParameter(0,0.956544);
00049 faEtaBarrel->SetParameter(0,0.014664);
00050 fbEtaBarrel->SetParameter(0,0.00975451);
00051 faBarrel->SetParameter(1,0.165627);
00052 fbBarrel->SetParameter(1,0.13632);
00053 fcBarrel->SetParameter(1,0.0857207);
00054 faEtaBarrel->SetParameter(1,-0.0426776);
00055 fbEtaBarrel->SetParameter(1,0.102247);
00056 faBarrel->SetParameter(2,0.827718);
00057 fbBarrel->SetParameter(2,-0.758013);
00058 fcBarrel->SetParameter(2,-0.44347);
00059 faEtaBarrel->SetParameter(2,431.054);
00060 fbEtaBarrel->SetParameter(2,436.21);
00061 faBarrel->SetParameter(3,231.339);
00062 fbBarrel->SetParameter(3,183.627);
00063 fcBarrel->SetParameter(3,63.3479);
00064 faBarrel->SetParameter(4,2.45332);
00065 fbBarrel->SetParameter(4,1);
00066 fcBarrel->SetParameter(4,1.24174);
00067 faBarrel->SetParameter(5,29.6603);
00068 fbBarrel->SetParameter(5,39.6784);
00069 fcBarrel->SetParameter(5,12.322);
00070
00071
00072 faEndcap = new TF1("faEndcap","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
00073 fbEndcap = new TF1("fbEndcap","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
00074 fcEndcap = new TF1("fcEndcap","[0]+([1]+[2]/sqrt(x))*exp(-x/[3])-[4]*exp(-x*x/[5])",1.,1000.);
00075 faEtaEndcap = new TF1("faEtaEndcap","[0]+[1]*exp(-x/[2])",1.,1000.);
00076 fbEtaEndcap = new TF1("fbEtaEndcap","[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))",1.,1000.);
00077 faEndcap->SetParameter(0,1.1272);
00078 fbEndcap->SetParameter(0,0.982824);
00079 fcEndcap->SetParameter(0,0.950244);
00080 faEtaEndcap->SetParameter(0,-0.000582903);
00081 fbEtaEndcap->SetParameter(0,0.0267319);
00082 faEndcap->SetParameter(1,0.258536);
00083 fbEndcap->SetParameter(1,0.0977533);
00084 fcEndcap->SetParameter(1,0.00564779);
00085 faEtaEndcap->SetParameter(1,-0.000482148);
00086 fbEtaEndcap->SetParameter(1,-0.554552);
00087 faEndcap->SetParameter(2,0.808071);
00088 fbEndcap->SetParameter(2,0.155416);
00089 fcEndcap->SetParameter(2,0.227162);
00090 faEtaEndcap->SetParameter(2,209.466);
00091 fbEtaEndcap->SetParameter(2,1.71188);
00092 faEndcap->SetParameter(3,214.039);
00093 fbEndcap->SetParameter(3,240.379);
00094 fcEndcap->SetParameter(3,207.786);
00095 fbEtaEndcap->SetParameter(3,0.235834);
00096 faEndcap->SetParameter(4,2);
00097 fbEndcap->SetParameter(4,1.2);
00098 fcEndcap->SetParameter(4,1.32824);
00099 fbEtaEndcap->SetParameter(4,-135.431);
00100 faEndcap->SetParameter(5,47.2602);
00101 fbEndcap->SetParameter(5,78.3083);
00102 fcEndcap->SetParameter(5,22.1825);
00103
00104 }
00105
00106 void
00107 PFEnergyCalibration::energyEmHad(double t, double& e, double&h, double eta, double phi) const {
00108
00109
00110
00111 double tt = t;
00112 double ee = e;
00113 double hh = h;
00114 double a = 1.;
00115 double b = 1.;
00116 double etaCorrE = 1.;
00117 double etaCorrH = 1.;
00118 t = min(999.9,max(tt,e+h));
00119 if ( t < 1. ) return;
00120
00121
00122 if ( fabs(eta) < 1.48 ) {
00123
00124
00125 a = e>0. ? aBarrel(t) : 1.;
00126 b = e>0. ? bBarrel(t) : cBarrel(t);
00127 double thresh = e > 0. ? threshE : threshH;
00128
00129
00130 if ( a < -0.25 || b < -0.25 ) {
00131 a = 1.;
00132 b = 1.;
00133 thresh = 0.;
00134 }
00135
00136
00137 t = min(999.9,max(tt, thresh+a*e+b*h));
00138
00139
00140 etaCorrE = 1. + aEtaBarrel(t) + bEtaBarrel(t)*fabs(eta)*fabs(eta);
00141 etaCorrH = 1.;
00142
00143 t = max(tt, thresh+etaCorrE*a*e+etaCorrH*b*h);
00144
00145 if ( e > 0. && thresh > 0. )
00146 e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
00147 if ( h > 0. && thresh > 0. )
00148 h = threshH + etaCorrH * b * h;
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165 } else {
00166
00167
00168 a = e>0. ? aEndcap(t) : 1.;
00169 b = e>0. ? bEndcap(t) : cEndcap(t);
00170 double thresh = e > 0. ? threshE : threshH;
00171
00172 if ( a < -0.25 || b < -0.25 ) {
00173 a = 1.;
00174 b = 1.;
00175 thresh = 0.;
00176 }
00177
00178
00179 t = min(999.9,max(tt, thresh+a*e+b*h));
00180
00181
00182 double dEta = fabs ( fabs(eta) - 1.5 );
00183 double etaPow = dEta * dEta * dEta * dEta;
00184
00185 etaCorrE = 1. + aEtaEndcap(t) + bEtaEndcap(t)*etaPow;
00186 etaCorrH = 1. + aEtaEndcap(t) + bEtaEndcap(t)*etaPow;
00187
00188
00189
00190
00191
00192
00193
00194
00195 t = min(999.9,max(tt, thresh + etaCorrE*a*e + etaCorrH*b*h));
00196
00197 if ( e > 0. && thresh > 0. )
00198 e = h > 0. ? threshE-threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
00199 if ( h > 0. && thresh > 0. )
00200 h = threshH + b * etaCorrH * h;
00201
00202
00203 }
00204
00205
00206 if ( e < 0. || h < 0. ) {
00207
00208
00209
00210
00211
00212
00213 if ( e < 0. ) e = ee;
00214 if ( h < 0. ) h = hh;
00215 }
00216
00217
00218
00219
00220 }
00221
00222
00223 double
00224 PFEnergyCalibration::aBarrel(double x) const {
00225
00226 if ( pfCalibrations ) {
00227
00228 BinningPointByMap point;
00229 point.insert(BinningVariables::JetEt, x);
00230 return pfCalibrations->getResult(PerformanceResult::PFfa_BARREL,point);
00231
00232 } else {
00233
00234 return faBarrel->Eval(x);
00235
00236 }
00237 }
00238
00239 double
00240 PFEnergyCalibration::bBarrel(double x) const {
00241
00242 if ( pfCalibrations ) {
00243
00244 BinningPointByMap point;
00245 point.insert(BinningVariables::JetEt, x);
00246 return pfCalibrations->getResult(PerformanceResult::PFfb_BARREL,point);
00247
00248 } else {
00249
00250 return fbBarrel->Eval(x);
00251
00252 }
00253 }
00254
00255 double
00256 PFEnergyCalibration::cBarrel(double x) const {
00257
00258 if ( pfCalibrations ) {
00259
00260 BinningPointByMap point;
00261 point.insert(BinningVariables::JetEt, x);
00262 return pfCalibrations->getResult(PerformanceResult::PFfc_BARREL,point);
00263
00264 } else {
00265
00266 return fcBarrel->Eval(x);
00267
00268 }
00269 }
00270
00271 double
00272 PFEnergyCalibration::aEtaBarrel(double x) const {
00273
00274 if ( pfCalibrations ) {
00275
00276 BinningPointByMap point;
00277 point.insert(BinningVariables::JetEt, x);
00278 return pfCalibrations->getResult(PerformanceResult::PFfaEta_BARREL,point);
00279
00280 } else {
00281
00282 return faEtaBarrel->Eval(x);
00283
00284 }
00285 }
00286
00287 double
00288 PFEnergyCalibration::bEtaBarrel(double x) const {
00289
00290 if ( pfCalibrations ) {
00291
00292 BinningPointByMap point;
00293 point.insert(BinningVariables::JetEt, x);
00294 return pfCalibrations->getResult(PerformanceResult::PFfbEta_BARREL,point);
00295
00296 } else {
00297
00298 return fbEtaBarrel->Eval(x);
00299
00300 }
00301 }
00302
00303 double
00304 PFEnergyCalibration::aEndcap(double x) const {
00305
00306 if ( pfCalibrations ) {
00307
00308 BinningPointByMap point;
00309 point.insert(BinningVariables::JetEt, x);
00310 return pfCalibrations->getResult(PerformanceResult::PFfa_ENDCAP,point);
00311
00312 } else {
00313
00314 return faEndcap->Eval(x);
00315
00316 }
00317 }
00318
00319 double
00320 PFEnergyCalibration::bEndcap(double x) const {
00321
00322 if ( pfCalibrations ) {
00323
00324 BinningPointByMap point;
00325 point.insert(BinningVariables::JetEt, x);
00326 return pfCalibrations->getResult(PerformanceResult::PFfb_ENDCAP,point);
00327
00328 } else {
00329
00330 return fbEndcap->Eval(x);
00331
00332 }
00333 }
00334
00335 double
00336 PFEnergyCalibration::cEndcap(double x) const {
00337
00338 if ( pfCalibrations ) {
00339
00340 BinningPointByMap point;
00341 point.insert(BinningVariables::JetEt, x);
00342 return pfCalibrations->getResult(PerformanceResult::PFfc_ENDCAP,point);
00343
00344 } else {
00345
00346 return fcEndcap->Eval(x);
00347
00348 }
00349 }
00350
00351 double
00352 PFEnergyCalibration::aEtaEndcap(double x) const {
00353
00354 if ( pfCalibrations ) {
00355
00356 BinningPointByMap point;
00357 point.insert(BinningVariables::JetEt, x);
00358 return pfCalibrations->getResult(PerformanceResult::PFfaEta_ENDCAP,point);
00359
00360 } else {
00361
00362 return faEtaEndcap->Eval(x);
00363
00364 }
00365 }
00366
00367 double
00368 PFEnergyCalibration::bEtaEndcap(double x) const {
00369
00370 if ( pfCalibrations ) {
00371
00372 BinningPointByMap point;
00373 point.insert(BinningVariables::JetEt, x);
00374 return pfCalibrations->getResult(PerformanceResult::PFfbEta_ENDCAP,point);
00375
00376 } else {
00377
00378 return fbEtaEndcap->Eval(x);
00379
00380 }
00381 }
00382
00383
00384 double
00385 PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
00386 std::vector<double> &EclustersPS1,
00387 std::vector<double> &EclustersPS2,
00388 bool crackCorrection ){
00389 double eEcal = clusterEcal.energy();
00390
00391 reco::PFCluster myPFCluster=clusterEcal;
00392 myPFCluster.calculatePositionREP();
00393 double eta = myPFCluster.positionREP().eta();
00394 double phi = myPFCluster.positionREP().phi();
00395
00396 double ePS1 = 0;
00397 double ePS2 = 0;
00398
00399 for(unsigned i=0;i<EclustersPS1.size();i++) ePS1 += EclustersPS1[i];
00400 for(unsigned i=0;i<EclustersPS2.size();i++) ePS2 += EclustersPS2[i];
00401
00402 double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi, crackCorrection);
00403 if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
00404 return calibrated;
00405 }
00406
00407 double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
00408 std::vector<double> &EclustersPS1,
00409 std::vector<double> &EclustersPS2,
00410 double& ps1,double& ps2,
00411 bool crackCorrection){
00412 double eEcal = clusterEcal.energy();
00413
00414 reco::PFCluster myPFCluster=clusterEcal;
00415 myPFCluster.calculatePositionREP();
00416 double eta = myPFCluster.positionREP().eta();
00417 double phi = myPFCluster.positionREP().phi();
00418
00419 double ePS1 = 0;
00420 double ePS2 = 0;
00421
00422 for(unsigned i=0;i<EclustersPS1.size();i++) ePS1 += EclustersPS1[i];
00423 for(unsigned i=0;i<EclustersPS2.size();i++) ePS2 += EclustersPS2[i];
00424
00425 double calibrated = Ecorr(eEcal,ePS1,ePS2,eta,phi,ps1,ps2,crackCorrection);
00426 if(eEcal!=0 && calibrated==0) std::cout<<"Eecal = "<<eEcal<<" eta = "<<eta<<" phi = "<<phi<<std::endl;
00427 return calibrated;
00428 }
00429
00430
00431 std::ostream& operator<<(std::ostream& out,
00432 const PFEnergyCalibration& calib) {
00433
00434 if(!out ) return out;
00435
00436 out<<"PFEnergyCalibration -- "<<endl;
00437
00438 if ( calib.pfCalibrations ) {
00439
00440 std::cout << "Functions taken from the global tags : " << std::endl;
00441
00442 static std::map<std::string, PerformanceResult::ResultType> functType;
00443
00444 functType["PFfa_BARREL"] = PerformanceResult::PFfa_BARREL;
00445 functType["PFfa_ENDCAP"] = PerformanceResult::PFfa_ENDCAP;
00446 functType["PFfb_BARREL"] = PerformanceResult::PFfb_BARREL;
00447 functType["PFfb_ENDCAP"] = PerformanceResult::PFfb_ENDCAP;
00448 functType["PFfc_BARREL"] = PerformanceResult::PFfc_BARREL;
00449 functType["PFfc_ENDCAP"] = PerformanceResult::PFfc_ENDCAP;
00450 functType["PFfaEta_BARREL"] = PerformanceResult::PFfaEta_BARREL;
00451 functType["PFfaEta_ENDCAP"] = PerformanceResult::PFfaEta_ENDCAP;
00452 functType["PFfbEta_BARREL"] = PerformanceResult::PFfbEta_BARREL;
00453 functType["PFfbEta_ENDCAP"] = PerformanceResult::PFfbEta_ENDCAP;
00454
00455 for(std::map<std::string,PerformanceResult::ResultType>::const_iterator
00456 func = functType.begin();
00457 func != functType.end();
00458 ++func) {
00459
00460 cout << "Function: " << func->first << endl;
00461 PerformanceResult::ResultType fType = func->second;
00462 calib.pfCalibrations->printFormula(fType);
00463 }
00464
00465 } else {
00466
00467 std::cout << "Default calibration functions : " << std::endl;
00468
00469 calib.faBarrel->Print();
00470 calib.fbBarrel->Print();
00471 calib.fcBarrel->Print();
00472 calib.faEtaBarrel->Print();
00473 calib.fbEtaBarrel->Print();
00474 calib.faEndcap->Print();
00475 calib.fbEndcap->Print();
00476 calib.fcEndcap->Print();
00477 calib.faEtaEndcap->Print();
00478 calib.fbEtaEndcap->Print();
00479 }
00480
00481 return out;
00482 }
00483
00484
00485
00486
00499
00500
00501
00507
00508
00509
00510 double
00511 PFEnergyCalibration::minimum(double a,double b){
00512 if(TMath::Abs(b)<TMath::Abs(a)) a=b;
00513 return a;
00514 }
00515
00516
00517
00518 double
00519 PFEnergyCalibration::dCrackPhi(double phi, double eta){
00520
00521 static double pi= M_PI;
00522
00523
00524 static std::vector<double> cPhi;
00525 if(cPhi.size()==0)
00526 {
00527 cPhi.resize(18,0);
00528 cPhi[0]=2.97025;
00529 for(unsigned i=1;i<=17;++i) cPhi[i]=cPhi[0]-2*i*pi/18;
00530 }
00531
00532
00533 static double delta_cPhi=0.00638;
00534
00535 double m;
00536
00537
00538 if(eta<0) phi +=delta_cPhi;
00539
00540 if (phi>=-pi && phi<=pi){
00541
00542
00543 if (phi<cPhi[17] || phi>=cPhi[0]){
00544 if (phi<0) phi+= 2*pi;
00545 m = minimum(phi -cPhi[0],phi-cPhi[17]-2*pi);
00546 }
00547
00548
00549 else{
00550 bool OK = false;
00551 unsigned i=16;
00552 while(!OK){
00553 if (phi<cPhi[i]){
00554 m=minimum(phi-cPhi[i+1],phi-cPhi[i]);
00555 OK=true;
00556 }
00557 else i-=1;
00558 }
00559 }
00560 }
00561 else{
00562 m=0.;
00563 std::cout<<"Problem in dminphi"<<std::endl;
00564 }
00565 if(eta<0) m=-m;
00566 return m;
00567 }
00568
00569
00570 double
00571 PFEnergyCalibration::CorrPhi(double phi, double eta) {
00572
00573
00574 static double p1= 5.59379e-01;
00575 static double p2= -1.26607e-03;
00576 static double p3= 9.61133e-04;
00577
00578 static double p4= 1.81691e-01;
00579 static double p5= -4.97535e-03;
00580 static double p6= 1.31006e-03;
00581
00582 static double p7= 1.38498e-01;
00583 static double p8= 1.18599e-04;
00584 static double p9= 2.01858e-03;
00585
00586
00587 double dminphi = dCrackPhi(phi,eta);
00588
00589 double result = (1+p1*TMath::Gaus(dminphi,p2,p3)+p4*TMath::Gaus(dminphi,p5,p6)+p7*TMath::Gaus(dminphi,p8,p9));
00590
00591 return result;
00592 }
00593
00594
00595
00596 double
00597 PFEnergyCalibration::CorrEta(double eta){
00598
00599
00600 static std::vector<double> a;
00601 static std::vector<double> m;
00602 static std::vector<double> s;
00603 static std::vector<double> sa;
00604 static std::vector<double> ss;
00605
00606 if(a.size()==0)
00607 {
00608 a.push_back(6.13349e-01) ;a.push_back(5.08146e-01) ;a.push_back(4.44480e-01) ;a.push_back(3.3487e-01) ;a.push_back(7.65627e-01) ;
00609 m.push_back(-1.79514e-02);m.push_back(4.44747e-01) ;m.push_back(7.92824e-01) ;m.push_back(1.14090e+00) ;m.push_back(1.47464e+00) ;
00610 s.push_back(7.92382e-03) ;s.push_back(3.06028e-03) ;s.push_back(3.36139e-03) ;s.push_back(3.94521e-03) ;s.push_back(8.63950e-04) ;
00611 sa.push_back(1.27228e+01);sa.push_back(3.81517e-02) ;sa.push_back(1.63507e-01);sa.push_back(-6.56480e-02);sa.push_back(1.87160e-01);
00612 ss.push_back(5.48753e-02);ss.push_back(-1.00223e-02);ss.push_back(2.22866e-03);ss.push_back(4.26288e-04) ;ss.push_back(2.67937e-03);
00613 }
00614 double result = 1;
00615
00616 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]));
00617
00618 return result;
00619 }
00620
00621
00622
00623 double
00624 PFEnergyCalibration::CorrBarrel(double E, double eta) {
00625
00626
00627 static double p0=1.00000e+00;
00628 static double p1=3.27753e+01;
00629 static double p2=2.28552e-02;
00630 static double p3=3.06139e+00;
00631 static double p4=2.25135e-01;
00632 static double p5=1.47824e+00;
00633 static double p6=1.09e-02;
00634 static double p7=4.19343e+01;
00635
00636
00637 static double p8=2.705593e-03;
00638
00639 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);
00640
00641 return result;
00642 }
00643
00644
00645
00654
00655
00656
00657
00658
00659 double
00660 PFEnergyCalibration::Alpha(double eta) {
00661
00662
00663 static double p0 = 5.97621e-01;
00664
00665
00666 static double p1 =-1.86407e-01;
00667 static double p2 = 3.85197e-01;
00668
00669
00670 static double norm = (p1+p2*(2.6+1.656)/2);
00671
00672 double result = p0*(p1+p2*eta)/norm;
00673
00674 return result;
00675 }
00676
00677 double
00678 PFEnergyCalibration::Beta(double E, double eta) {
00679
00680
00681 static double p0 = 0.032;
00682 static double p1 = 9.70394e-02;
00683 static double p2 = 2.23072e+01;
00684 static double p3 = 100;
00685
00686
00687 static double p4 = 1.02496e+00 ;
00688 static double p5 = -4.40176e-03 ;
00689
00690
00691 static double norm = (p4+p5*(2.6+1.656)/2);
00692
00693 double result = (1.0012+p0*TMath::Exp(-E/p3)+p1*TMath::Exp(-E/p2))*(p4+p5*eta)/norm;
00694 return result;
00695 }
00696
00697
00698 double
00699 PFEnergyCalibration::Gamma(double etaEcal) {
00700
00701
00702 static double p0 = 2.49752e-02;
00703
00704
00705 static double p1 = 6.48816e-02;
00706 static double p2 = -1.59517e-02;
00707
00708
00709 static double norm = (p1+p2*(2.6+1.656)/2);
00710
00711 double result = p0*(p1+p2*etaEcal)/norm;
00712
00713 return result;
00714 }
00715
00716
00717
00723
00724
00725
00726
00727 double
00728 PFEnergyCalibration::EcorrBarrel(double E, double eta, double phi,
00729 bool crackCorrection ){
00730
00731
00732 double correction = crackCorrection ? std::max(CorrEta(eta),CorrPhi(phi,eta)) : 1.;
00733 double result = E * CorrBarrel(E,eta) * correction;
00734
00735 return result;
00736 }
00737
00738
00739
00740 double
00741 PFEnergyCalibration::EcorrZoneBeforePS(double E, double eta){
00742
00743
00744 static double p0 =1;
00745 static double p1 =0.18;
00746 static double p2 =8.;
00747
00748
00749 static double p3 =0.3;
00750 static double p4 =1.11;
00751 static double p5 =0.025;
00752 static double p6 =1.49;
00753 static double p7 =0.6;
00754
00755
00756 static double norm = 1.21;
00757
00758 double result = E*(p0+p1*TMath::Exp(-E/p2))*(p3+p4*TMath::Gaus(eta,p6,p5)+p7*eta)/norm;
00759
00760 return result;
00761 }
00762
00763
00764
00765
00766 double
00767 PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal) {
00768
00769
00770 double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+Gamma(etaEcal)*(ePS1+Alpha(etaEcal)*ePS2)/9e-5 ;
00771
00772
00773 static double p0 = 1.00;
00774 static double p1 = 2.18;
00775 static double p2 =1.94;
00776 static double p3 =4.13;
00777 static double p4 =1.127;
00778
00779 double result = E*(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
00780
00781 return result;
00782 }
00783
00784
00785
00786 double
00787 PFEnergyCalibration::EcorrPS(double eEcal,double ePS1,double ePS2,double etaEcal,double & outputPS1, double & outputPS2) {
00788
00789
00790 double gammaprime=Gamma(etaEcal)/9e-5;
00791 outputPS1=gammaprime*ePS1;
00792 outputPS2=gammaprime*Alpha(etaEcal)*ePS2;
00793 double E = Beta(1.0155*eEcal+0.025*(ePS1+0.5976*ePS2)/9e-5,etaEcal)*eEcal+outputPS1+outputPS2;
00794
00795
00796 static double p0 = 1.00;
00797 static double p1 = 2.18;
00798 static double p2 =1.94;
00799 static double p3 =4.13;
00800 static double p4 =1.127;
00801
00802 double corrfac=(p0+p1*TMath::Exp(-E/p2)-p3*TMath::Exp(-E/p4));
00803 outputPS1*=corrfac;
00804 outputPS2*=corrfac;
00805 double result = E*corrfac;
00806
00807 return result;
00808 }
00809
00810
00811
00812
00813 double
00814 PFEnergyCalibration::EcorrPS_ePSNil(double eEcal,double eta){
00815
00816
00817 static double p0= 1.02;
00818 static double p1= 0.165;
00819 static double p2= 6.5 ;
00820 static double p3= 2.1 ;
00821
00822
00823 static double p4 = 1.02496e+00 ;
00824 static double p5 = -4.40176e-03 ;
00825
00826
00827 static double norm = (p4+p5*(2.6+1.656)/2);
00828
00829 double result = eEcal*(p0+p1*TMath::Exp(-TMath::Abs(eEcal-p3)/p2))*(p4+p5*eta)/norm;
00830
00831 return result;
00832 }
00833
00834
00835
00836 double
00837 PFEnergyCalibration::EcorrZoneAfterPS(double E, double eta){
00838
00839
00840 static double p0 =1;
00841 static double p1 = 0.058;
00842 static double p2 =12.5;
00843 static double p3 =-1.05444e+00;
00844 static double p4 =-5.39557e+00;
00845 static double p5 =8.38444e+00;
00846 static double p6 = 6.10998e-01 ;
00847
00848
00849 static double p7 =1.06161e+00;
00850 static double p8 = 0.41;
00851 static double p9 =2.918;
00852 static double p10 =0.0181;
00853 static double p11= 2.05;
00854 static double p12 =2.99;
00855 static double p13=0.0287;
00856
00857
00858 static double norm=1.045;
00859
00860 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;
00861 return result;
00862 }
00863
00864
00865
00866
00867
00868
00869 double
00870 PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,
00871 double eta,double phi,
00872 bool crackCorrection ) {
00873
00874 static double endBarrel=1.48;
00875 static double beginingPS=1.65;
00876 static double endPS=2.6;
00877 static double endEndCap=2.98;
00878
00879 double result=0;
00880
00881 eta=TMath::Abs(eta);
00882
00883 if(eEcal>0){
00884 if(eta <= endBarrel) result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
00885 else if(eta <= beginingPS) result = EcorrZoneBeforePS(eEcal,eta);
00886 else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
00887 else if(eta < endPS) result = EcorrPS(eEcal,ePS1,ePS2,eta);
00888 else if(eta < endEndCap) result = EcorrZoneAfterPS(eEcal,eta);
00889 else result =eEcal;
00890 }
00891 else result = eEcal;
00892
00893 if(result<eEcal) result=eEcal;
00894 return result;
00895 }
00896
00897
00898
00899 double
00900 PFEnergyCalibration::Ecorr(double eEcal,double ePS1,double ePS2,double eta,double phi,double& ps1,double&ps2,bool crackCorrection) {
00901
00902 static double endBarrel=1.48;
00903 static double beginingPS=1.65;
00904 static double endPS=2.6;
00905 static double endEndCap=2.98;
00906
00907 double result=0;
00908
00909 eta=TMath::Abs(eta);
00910
00911 if(eEcal>0){
00912 if(eta <= endBarrel) result = EcorrBarrel(eEcal,eta,phi,crackCorrection);
00913 else if(eta <= beginingPS) result = EcorrZoneBeforePS(eEcal,eta);
00914 else if((eta < endPS) && ePS1==0 && ePS2==0) result = EcorrPS_ePSNil(eEcal,eta);
00915 else if(eta < endPS) result = EcorrPS(eEcal,ePS1,ePS2,eta,ps1,ps2);
00916 else if(eta < endEndCap) result = EcorrZoneAfterPS(eEcal,eta);
00917 else result =eEcal;
00918 }
00919 else result = eEcal;
00920
00921 if(result<eEcal) result=eEcal;
00922 return result;
00923 }