CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/RecoEgamma/EgammaPhotonAlgos/src/EnergyUncertaintyPhotonSpecific.cc

Go to the documentation of this file.
00001 #include "RecoEgamma/EgammaPhotonAlgos/interface/EnergyUncertaintyPhotonSpecific.h"
00002 #include "TMath.h"
00003 
00004 #include <iostream>
00005 
00006 
00007 EnergyUncertaintyPhotonSpecific::EnergyUncertaintyPhotonSpecific(  const edm::ParameterSet& config ) {
00008 
00009 
00010 
00011 }
00012 
00013 
00014 EnergyUncertaintyPhotonSpecific::~EnergyUncertaintyPhotonSpecific() {
00015 
00016 }
00017 
00018 
00019 
00020 void EnergyUncertaintyPhotonSpecific::init (  const edm::EventSetup& theEventSetup ) {
00021  
00022 }
00023 
00024 
00025 double EnergyUncertaintyPhotonSpecific::computePhotonEnergyUncertainty_lowR9(double eta, double brem, double energy){
00026 
00027  
00028 
00029   double et = energy/cosh(eta);
00030 
00031   const int nBinsEta=6;
00032   const double EtaBins[nBinsEta+1] = {0.0, 0.7, 1.15, 1.44, 1.56, 2.0, 2.5};
00033 
00034   const int nBinsBrem=6;
00035   const double BremBins  [nBinsBrem+1]= {0.8, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0};
00036 
00037 
00038   float par0[nBinsEta][nBinsBrem];
00039   float par1[nBinsEta][nBinsBrem];
00040   float par2[nBinsEta][nBinsBrem];
00041   float par3[nBinsEta][nBinsBrem];
00042 
00043   par0[0][0]=0.0232291;
00044   par1[0][0]=0;
00045   par2[0][0]=0;
00046   par3[0][0]=0;
00047 
00048   par0[0][1]=0.00703187;
00049   par1[0][1]=0.646644;
00050   par2[0][1]=-7.4698;
00051   par3[0][1]=5.53373e-08;
00052 
00053   par0[0][2]=0.00692465;
00054   par1[0][2]=0.292698;
00055   par2[0][2]=4.16907;
00056   par3[0][2]=5.61149e-06;
00057 
00058   par0[0][3]=0.00855993;
00059   par1[0][3]=0.280843;
00060   par2[0][3]=4.25527;
00061   par3[0][3]=9.6404e-07;
00062 
00063   par0[0][4]=0.00795058;
00064   par1[0][4]=0.370007;
00065   par2[0][4]=3.03429;
00066   par3[0][4]=4.43986e-07;
00067 
00068   par0[0][5]=0.0107494;
00069   par1[0][5]=0.276159;
00070   par2[0][5]=4.44532;
00071   par3[0][5]=2.58822e-06;
00072 
00073   par0[1][0]=0.0614866;
00074   par1[1][0]=0;
00075   par2[1][0]=0;
00076   par3[1][0]=0;
00077 
00078   par0[1][1]=0.00894211;
00079   par1[1][1]=0.466937;
00080   par2[1][1]=3.33434;
00081   par3[1][1]=0.000114835;
00082 
00083   par0[1][2]=0.0102959;
00084   par1[1][2]=0.313568;
00085   par2[1][2]=6.34301;
00086   par3[1][2]=2.86726e-07;
00087 
00088   par0[1][3]=0.0128934;
00089   par1[1][3]=0.302943;
00090   par2[1][3]=6.35598;
00091   par3[1][3]=0.00190694;
00092 
00093   par0[1][4]=0.0130199;
00094   par1[1][4]=0.505135;
00095   par2[1][4]=2.52964;
00096   par3[1][4]=0.120204;
00097 
00098   par0[1][5]=0.0180839;
00099   par1[1][5]=0.382134;
00100   par2[1][5]=5.3388;
00101   par3[1][5]=3.59921e-07;
00102 
00103   par0[2][0]=0.0291343;
00104   par1[2][0]=0;
00105   par2[2][0]=0;
00106   par3[2][0]=0;
00107 
00108   par0[2][1]=0.00876269;
00109   par1[2][1]=0.375159;
00110   par2[2][1]=7.11411;
00111   par3[2][1]=0.0438575;
00112 
00113   par0[2][2]=0.0120863;
00114   par1[2][2]=0.397635;
00115   par2[2][2]=5.97451;
00116   par3[2][2]=0.0469782;
00117 
00118   par0[2][3]=0.0112655;
00119   par1[2][3]=0.856565;
00120   par2[2][3]=-5.76122;
00121   par3[2][3]=4.99993;
00122 
00123   par0[2][4]=0.0168267;
00124   par1[2][4]=0.636468;
00125   par2[2][4]=-1.54548;
00126   par3[2][4]=4.99992;
00127 
00128   par0[2][5]=0.0168059;
00129   par1[2][5]=1.09268;
00130   par2[2][5]=-0.547554;
00131   par3[2][5]=0.0952985;
00132 
00133   par0[3][0]=0.158403;
00134   par1[3][0]=0;
00135   par2[3][0]=0;
00136   par3[3][0]=0;
00137 
00138   par0[3][1]=0.0717431;
00139   par1[3][1]=1.66981;
00140   par2[3][1]=6.86275;
00141   par3[3][1]=0.00543544;
00142 
00143   par0[3][2]=0.0385666;
00144   par1[3][2]=3.6319;
00145   par2[3][2]=-3.76633;
00146   par3[3][2]=6.56718e-05;
00147 
00148   par0[3][3]=0.0142631;
00149   par1[3][3]=8.85991;
00150   par2[3][3]=-32.6073;
00151   par3[3][3]=0.00119538;
00152 
00153   par0[3][4]=0.0421638;
00154   par1[3][4]=3.1289;
00155   par2[3][4]=-6.58653;
00156   par3[3][4]=1.10125e-05;
00157 
00158   par0[3][5]=0.046331;
00159   par1[3][5]=1.29951;
00160   par2[3][5]=1.76117;
00161   par3[3][5]=0.00204206;
00162 
00163   par0[4][0]=0.0483944;
00164   par1[4][0]=0;
00165   par2[4][0]=0;
00166   par3[4][0]=0;
00167 
00168   par0[4][1]=0.0168516;
00169   par1[4][1]=1.19617;
00170   par2[4][1]=-6.78666;
00171   par3[4][1]=4.98192;
00172 
00173   par0[4][2]=0.0243039;
00174   par1[4][2]=0.994626;
00175   par2[4][2]=-4.26073;
00176   par3[4][2]=4.99984;
00177 
00178   par0[4][3]=0.031795;
00179   par1[4][3]=0.875925;
00180   par2[4][3]=1.43183;
00181   par3[4][3]=0.0920944;
00182 
00183   par0[4][4]=0.0414953;
00184   par1[4][4]=0.654605;
00185   par2[4][4]=4.45367;
00186   par3[4][4]=0.030385;
00187 
00188   par0[4][5]=0.058031;
00189   par1[4][5]=0.292915;
00190   par2[4][5]=8.48307;
00191   par3[4][5]=0.0134321;
00192 
00193   par0[5][0]=0.107158;
00194   par1[5][0]=0;
00195   par2[5][0]=0;
00196   par3[5][0]=0;
00197 
00198   par0[5][1]=0.021685;
00199   par1[5][1]=0.574207;
00200   par2[5][1]=-0.566981;
00201   par3[5][1]=0.0120609;
00202 
00203   par0[5][2]=0.0196619;
00204   par1[5][2]=0.940217;
00205   par2[5][2]=-6.05845;
00206   par3[5][2]=0.000193818;
00207 
00208   par0[5][3]=0.0324734;
00209   par1[5][3]=0.574766;
00210   par2[5][3]=-5.23571;
00211   par3[5][3]=4.9419;
00212 
00213   par0[5][4]=0.0414953;
00214   par1[5][4]=0.654605;
00215   par2[5][4]=4.45367;
00216   par3[5][4]=0.030385;
00217 
00218   par0[5][5]=0.058031;
00219   par1[5][5]=0.292915;
00220   par2[5][5]=8.48307;
00221   par3[5][5]=0.0134321;
00222 
00223 
00224   int iEtaSl = -1;                                                                         
00225   for (int iEta = 0; iEta < nBinsEta; ++iEta){                                                                       
00226     if ( EtaBins[iEta] <= TMath::Abs(eta) && TMath::Abs(eta) <EtaBins[iEta+1] ){                         
00227       iEtaSl = iEta;                                                                                                                                     
00228     }                                                                                                                                                                    
00229   }
00230 
00231 
00232   int iBremSl = -1;                                                                        
00233   for (int iBrem = 0; iBrem < nBinsBrem; ++iBrem){                                                                       
00234     if ( BremBins[iBrem] <= brem && brem <BremBins[iBrem+1] ){                   
00235       iBremSl = iBrem;                                                                                                                                   
00236     }                                                                                                                                                                    
00237   }
00238 
00239   if (TMath::Abs(eta)>2.5) iEtaSl = nBinsEta-1;
00240   if (brem<BremBins[0]) iBremSl = 0;
00241   if (brem>BremBins[nBinsBrem-1]) iBremSl = nBinsBrem-1;
00242 
00243   float uncertainty = 0;
00244   if (et<5) uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl]/(5-par2[iEtaSl][iBremSl]) + par3[iEtaSl][iBremSl]/((5-par2[iEtaSl][iBremSl])*(5-par2[iEtaSl][iBremSl]));
00245   if (et>200) uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl]/(200-par2[iEtaSl][iBremSl]) + par3[iEtaSl][iBremSl]/((200-par2[iEtaSl][iBremSl])*(200-par2[iEtaSl][iBremSl]));
00246 
00247   if (et>5 && et<200) uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl]/(et-par2[iEtaSl][iBremSl]) + par3[iEtaSl][iBremSl]/((et-par2[iEtaSl][iBremSl])*(et-par2[iEtaSl][iBremSl]));
00248 
00249   return (uncertainty*energy);
00250 
00251 }
00252 
00253 double EnergyUncertaintyPhotonSpecific::computePhotonEnergyUncertainty_highR9(double eta, double brem, double energy){
00254 
00255 
00256 
00257   double et = energy/cosh(eta);
00258 
00259   const int nBinsEta=6;
00260   const double EtaBins[nBinsEta+1] = {0.0, 0.7, 1.15, 1.44, 1.56, 2.0, 2.5};
00261 
00262   const int nBinsBrem=6;
00263   const double BremBins  [nBinsBrem+1]= {0.8, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0};
00264 
00265   float par0[nBinsEta][nBinsBrem];
00266   float par1[nBinsEta][nBinsBrem];
00267   float par2[nBinsEta][nBinsBrem];
00268   float par3[nBinsEta][nBinsBrem];
00269 
00270   par0[0][0]=0.00806753;
00271   par1[0][0]=0.143754;
00272   par2[0][0]=-0.00368104;
00273   par3[0][0]=0.219829;
00274 
00275   par0[0][1]=0.00899298;
00276   par1[0][1]=0.10159;
00277   par2[0][1]=4.70884;
00278   par3[0][1]=9.07419e-08;
00279 
00280   par0[0][2]=0;
00281   par1[0][2]=0;
00282   par2[0][2]=0;
00283   par3[0][2]=0;
00284 
00285   par0[0][3]=0;
00286   par1[0][3]=0;
00287   par2[0][3]=0;
00288   par3[0][3]=0;
00289 
00290   par0[0][4]=0;
00291   par1[0][4]=0;
00292   par2[0][4]=0;
00293   par3[0][4]=0;
00294 
00295   par0[0][5]=0;
00296   par1[0][5]=0;
00297   par2[0][5]=0;
00298   par3[0][5]=0;
00299 
00300   par0[1][0]=0.00880649;
00301   par1[1][0]=0.0716169;
00302   par2[1][0]=5.23856;
00303   par3[1][0]=0.00632907;
00304 
00305   par0[1][1]=0.00972275;
00306   par1[1][1]=0.0752675;
00307   par2[1][1]=3.35623;
00308   par3[1][1]=2.49397e-07;
00309 
00310   par0[1][2]=0;
00311   par1[1][2]=0;
00312   par2[1][2]=0;
00313   par3[1][2]=0;
00314 
00315   par0[1][3]=0;
00316   par1[1][3]=0;
00317   par2[1][3]=0;
00318   par3[1][3]=0;
00319 
00320   par0[1][4]=0;
00321   par1[1][4]=0;
00322   par2[1][4]=0;
00323   par3[1][4]=0;
00324 
00325   par0[1][5]=0;
00326   par1[1][5]=0;
00327   par2[1][5]=0;
00328   par3[1][5]=0;
00329 
00330   par0[2][0]=0.0101474;
00331   par1[2][0]=-0.332171;
00332   par2[2][0]=-31.8456;
00333   par3[2][0]=22.543;
00334 
00335   par0[2][1]=0.0109109;
00336   par1[2][1]=0.0425903;
00337   par2[2][1]=6.52561;
00338   par3[2][1]=2.18593e-08;
00339 
00340   par0[2][2]=0;
00341   par1[2][2]=0;
00342   par2[2][2]=0;
00343   par3[2][2]=0;
00344 
00345   par0[2][3]=0;
00346   par1[2][3]=0;
00347   par2[2][3]=0;
00348   par3[2][3]=0;
00349 
00350   par0[2][4]=0;
00351   par1[2][4]=0;
00352   par2[2][4]=0;
00353   par3[2][4]=0;
00354 
00355   par0[2][5]=0;
00356   par1[2][5]=0;
00357   par2[2][5]=0;
00358   par3[2][5]=0;
00359 
00360   par0[3][0]=0.00343003;
00361   par1[3][0]=11.5791;
00362   par2[3][0]=-112.084;
00363   par3[3][0]=-863.968;
00364 
00365   par0[3][1]=0.0372159;
00366   par1[3][1]=1.44028;
00367   par2[3][1]=-40;
00368   par3[3][1]=0.00102639;
00369 
00370   par0[3][2]=0;
00371   par1[3][2]=0;
00372   par2[3][2]=0;
00373   par3[3][2]=0;
00374 
00375   par0[3][3]=0;
00376   par1[3][3]=0;
00377   par2[3][3]=0;
00378   par3[3][3]=0;
00379 
00380   par0[3][4]=0;
00381   par1[3][4]=0;
00382   par2[3][4]=0;
00383   par3[3][4]=0;
00384 
00385   par0[3][5]=0;
00386   par1[3][5]=0;
00387   par2[3][5]=0;
00388   par3[3][5]=0;
00389 
00390   par0[4][0]=0.0192411;
00391   par1[4][0]=0.0511006;
00392   par2[4][0]=7.56304;
00393   par3[4][0]=0.00331583;
00394 
00395   par0[4][1]=0.0195124;
00396   par1[4][1]=0.104321;
00397   par2[4][1]=5.71476;
00398   par3[4][1]=6.12472e-06;
00399 
00400   par0[4][2]=0;
00401   par1[4][2]=0;
00402   par2[4][2]=0;
00403   par3[4][2]=0;
00404 
00405   par0[4][3]=0;
00406   par1[4][3]=0;
00407   par2[4][3]=0;
00408   par3[4][3]=0;
00409 
00410   par0[4][4]=0;
00411   par1[4][4]=0;
00412   par2[4][4]=0;
00413   par3[4][4]=0;
00414 
00415   par0[4][5]=0;
00416   par1[4][5]=0;
00417   par2[4][5]=0;
00418   par3[4][5]=0;
00419 
00420   par0[5][0]=0.0203644;
00421   par1[5][0]=-0.050789;
00422   par2[5][0]=-7.96854;
00423   par3[5][0]=4.71223;
00424 
00425   par0[5][1]=0.0198718;
00426   par1[5][1]=0.106859;
00427   par2[5][1]=3.54235;
00428   par3[5][1]=6.89631e-06;
00429 
00430   par0[5][2]=0;
00431   par1[5][2]=0;
00432   par2[5][2]=0;
00433   par3[5][2]=0;
00434 
00435   par0[5][3]=0;
00436   par1[5][3]=0;
00437   par2[5][3]=0;
00438   par3[5][3]=0;
00439 
00440   par0[5][4]=0;
00441   par1[5][4]=0;
00442   par2[5][4]=0;
00443   par3[5][4]=0;
00444 
00445   par0[5][5]=0;
00446   par1[5][5]=0;
00447   par2[5][5]=0;
00448   par3[5][5]=0;
00449 
00450 
00451   int iEtaSl = -1;                                                                         
00452   for (int iEta = 0; iEta < nBinsEta; ++iEta){                                                                       
00453     if ( EtaBins[iEta] <= TMath::Abs(eta) && TMath::Abs(eta) <EtaBins[iEta+1] ){                         
00454       iEtaSl = iEta;                                                                                                                                     
00455     }                                                                                                                                                                    
00456   }
00457 
00458 
00459   int iBremSl = -1;                                                                        
00460   for (int iBrem = 0; iBrem < nBinsBrem; ++iBrem){                                                                       
00461     if ( BremBins[iBrem] <= brem && brem <BremBins[iBrem+1] ){                   
00462       iBremSl = iBrem;                                                                                                                                   
00463     }                                                                                                                                                                    
00464   }
00465 
00466   if (TMath::Abs(eta)>2.5) iEtaSl = nBinsEta-1;
00467   if (brem<BremBins[0]) iBremSl = 0;
00468   if (brem>BremBins[nBinsBrem-1]) iBremSl = nBinsBrem-1;
00469   if (brem>2) iBremSl = 1;
00470 
00471   float uncertainty = 0;
00472   if (et<5) uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl]/(5-par2[iEtaSl][iBremSl]) + par3[iEtaSl][iBremSl]/((5-par2[iEtaSl][iBremSl])*(5-par2[iEtaSl][iBremSl]));
00473   if (et>200) uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl]/(200-par2[iEtaSl][iBremSl]) + par3[iEtaSl][iBremSl]/((200-par2[iEtaSl][iBremSl])*(200-par2[iEtaSl][iBremSl]));
00474 
00475   if (et>5 && et<200) uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl]/(et-par2[iEtaSl][iBremSl]) + par3[iEtaSl][iBremSl]/((et-par2[iEtaSl][iBremSl])*(et-par2[iEtaSl][iBremSl]));
00476 
00477   return (uncertainty*energy);
00478 
00479 }