EgammaElectronAlgos/src/ecalClusterEnergyUncertaintyElectronSpecific.cc

0001 #include "RecoEgamma/EgammaElectronAlgos/interface/ecalClusterEnergyUncertaintyElectronSpecific.h"
0002 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0003
0004 float egamma::ecalClusterEnergyUncertaintyElectronSpecific(reco::SuperCluster const& superCluster) {
0005   float en = superCluster.energy();
0006   float eta = fabs(superCluster.eta());
0007   float et = en / cosh(eta);
0008   float brem = superCluster.etaWidth() != 0 ? superCluster.phiWidth() / superCluster.etaWidth() : 0;
0009
0010   const int nBinsEta = 6;
0011   const float EtaBins[nBinsEta + 1] = {0.0, 0.7, 1.15, 1.44, 1.56, 2.0, 2.5};
0012
0013   const int nBinsBrem = 6;
0014   const float BremBins[nBinsBrem + 1] = {0.8, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0};
0015
0016   float par0[nBinsEta][nBinsBrem];
0017   float par1[nBinsEta][nBinsBrem];
0018   float par2[nBinsEta][nBinsBrem];
0019   float par3[nBinsEta][nBinsBrem];
0020
0021   par0[0][0] = 0.00640519;
0022   par1[0][0] = 0.257578;
0023   par2[0][0] = 1.72437;
0024   par3[0][0] = 4.04686e-06;
0025
0026   par0[0][1] = 0.00709569;
0027   par1[0][1] = 0.279844;
0028   par2[0][1] = 1.13789;
0029   par3[0][1] = 1.16239e-05;
0030
0031   par0[0][2] = 0.0075544;
0032   par1[0][2] = 0.341346;
0033   par2[0][2] = 0.513396;
0034   par3[0][2] = 2.90054e-06;
0035
0036   par0[0][3] = 0.00659365;
0037   par1[0][3] = 0.517649;
0038   par2[0][3] = -3.1847;
0039   par3[0][3] = 7.37152e-07;
0040
0041   par0[0][4] = 0.00771696;
0042   par1[0][4] = 0.492897;
0043   par2[0][4] = -1.42222;
0044   par3[0][4] = 0.000358677;
0045
0046   par0[0][5] = 0.00561532;
0047   par1[0][5] = 0.655138;
0048   par2[0][5] = -3.29839;
0049   par3[0][5] = 6.25898e-07;
0050
0051   par0[1][0] = 0.00273646;
0052   par1[1][0] = 0.714568;
0053   par2[1][0] = -4.82956;
0054   par3[1][0] = 4.45878e-07;
0055
0056   par0[1][1] = 0.00679797;
0057   par1[1][1] = 0.472856;
0058   par2[1][1] = -0.281699;
0059   par3[1][1] = 5.46479e-05;
0060
0061   par0[1][2] = 0.00845532;
0062   par1[1][2] = 0.611624;
0063   par2[1][2] = -1.10104;
0064   par3[1][2] = 1.16803e-05;
0065
0066   par0[1][3] = 0.00831068;
0067   par1[1][3] = 0.853653;
0068   par2[1][3] = -4.23761;
0069   par3[1][3] = 2.61247e-05;
0070
0071   par0[1][4] = 0.00845457;
0072   par1[1][4] = 0.984985;
0073   par2[1][4] = -5.19548;
0074   par3[1][4] = 2.05044e-07;
0075
0076   par0[1][5] = 0.0110227;
0077   par1[1][5] = 1.00356;
0078   par2[1][5] = -4.31936;
0079   par3[1][5] = 0.14384;
0080
0081   par0[2][0] = -0.00192618;
0082   par1[2][0] = 1.69986;
0083   par2[2][0] = -16.4355;
0084   par3[2][0] = 1.94946e-06;
0085
0086   par0[2][1] = 0.0067622;
0087   par1[2][1] = 0.792209;
0088   par2[2][1] = -1.18521;
0089   par3[2][1] = 0.066577;
0090
0091   par0[2][2] = 0.00761595;
0092   par1[2][2] = 1.03058;
0093   par2[2][2] = -4.17237;
0094   par3[2][2] = 0.168543;
0095
0096   par0[2][3] = 0.0119179;
0097   par1[2][3] = 0.910145;
0098   par2[2][3] = -2.14122;
0099   par3[2][3] = 0.00342264;
0100
0101   par0[2][4] = 0.0139921;
0102   par1[2][4] = 1.01488;
0103   par2[2][4] = -2.46637;
0104   par3[2][4] = 0.0458434;
0105
0106   par0[2][5] = 0.013724;
0107   par1[2][5] = 1.49078;
0108   par2[2][5] = -6.60661;
0109   par3[2][5] = 0.297821;
0110
0111   par0[3][0] = -0.00197909;
0112   par1[3][0] = 4.40696;
0113   par2[3][0] = -4.88737;
0114   par3[3][0] = 4.99999;
0115
0116   par0[3][1] = 0.0340196;
0117   par1[3][1] = 3.86278;
0118   par2[3][1] = -10.899;
0119   par3[3][1] = 0.130098;
0120
0121   par0[3][2] = 0.0102397;
0122   par1[3][2] = 8.99643;
0123   par2[3][2] = -31.5122;
0124   par3[3][2] = 0.00118335;
0125
0126   par0[3][3] = 0.0110891;
0127   par1[3][3] = 8.01794;
0128   par2[3][3] = -21.9038;
0129   par3[3][3] = 0.000245975;
0130
0131   par0[3][4] = 0.0328931;
0132   par1[3][4] = 4.73441;
0133   par2[3][4] = -12.1148;
0134   par3[3][4] = 3.01721e-05;
0135
0136   par0[3][5] = 0.0395614;
0137   par1[3][5] = 3.54327;
0138   par2[3][5] = -12.6514;
0139   par3[3][5] = 0.119761;
0140
0141   par0[4][0] = 0.0121809;
0142   par1[4][0] = 0.965608;
0143   par2[4][0] = -4.19667;
0144   par3[4][0] = 0.129896;
0145
0146   par0[4][1] = 0.0168951;
0147   par1[4][1] = 1.0218;
0148   par2[4][1] = -4.03078;
0149   par3[4][1] = 0.374291;
0150
0151   par0[4][2] = 0.0213549;
0152   par1[4][2] = 1.29613;
0153   par2[4][2] = -4.89024;
0154   par3[4][2] = 0.0297165;
0155
0156   par0[4][3] = 0.0262602;
0157   par1[4][3] = 1.41674;
0158   par2[4][3] = -5.94928;
0159   par3[4][3] = 0.19298;
0160
0161   par0[4][4] = 0.0334892;
0162   par1[4][4] = 1.48572;
0163   par2[4][4] = -5.3175;
0164   par3[4][4] = 0.0157013;
0165
0166   par0[4][5] = 0.0347093;
0167   par1[4][5] = 1.63127;
0168   par2[4][5] = -7.27426;
0169   par3[4][5] = 0.201164;
0170
0171   par0[5][0] = 0.0185321;
0172   par1[5][0] = 0.255205;
0173   par2[5][0] = 1.56798;
0174   par3[5][0] = 5.07655e-11;
0175
0176   par0[5][1] = 0.0182718;
0177   par1[5][1] = 0.459086;
0178   par2[5][1] = -0.48198;
0179   par3[5][1] = 0.00114946;
0180
0181   par0[5][2] = 0.0175505;
0182   par1[5][2] = 0.92848;
0183   par2[5][2] = -4.52737;
0184   par3[5][2] = 0.154827;
0185
0186   par0[5][3] = 0.0233833;
0187   par1[5][3] = 0.804105;
0188   par2[5][3] = -3.75131;
0189   par3[5][3] = 2.84172;
0190
0191   par0[5][4] = 0.0334892;
0192   par1[5][4] = 1.48572;
0193   par2[5][4] = -5.3175;
0194   par3[5][4] = 0.0157013;
0195
0196   par0[5][5] = 0.0347093;
0197   par1[5][5] = 1.63127;
0198   par2[5][5] = -7.27426;
0199   par3[5][5] = 0.201164;
0200
0201   int iEtaSl = -1;
0202   for (int iEta = 0; iEta < nBinsEta; ++iEta) {
0203     if (EtaBins[iEta] <= eta && eta < EtaBins[iEta + 1]) {
0204       iEtaSl = iEta;
0205     }
0206   }
0207
0208   int iBremSl = -1;
0209   for (int iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
0210     if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
0211       iBremSl = iBrem;
0212     }
0213   }
0214
0215   //this code is confusing as it has no concept of under and overflow bins
0216   //we will use Et as an example but also applies to eta
0217   //underflow is 1st bin (naively currently labeled as 0 to 0.7, its really <0.7)
0218   //overflow is the final bin (naviely currently labeled as 5 to 10, its really >=5)
0219   //logic: if brem is 0<=brem <0.7 it will be already set to the 1st bin, this checks if its <0
0220   //logic: if brem is 5<=brem<10 it will be set to the last bin so this then checks if its >5 at which point
0221   //it also assigns it to the last bin. The value of 5 will have already been assigned in the for
0222   //loop above to the last bin so its okay that its a >5 test
0223   if (eta > EtaBins[nBinsEta - 1])
0224     iEtaSl = nBinsEta - 1;
0225   if (brem < BremBins[0])
0226     iBremSl = 0;
0227   if (brem > BremBins[nBinsBrem - 1])
0228     iBremSl = nBinsBrem - 1;
0229
0230   float uncertainty = 0;
0231   if (et <= 5)
0232     uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]) +
0233                   par3[iEtaSl][iBremSl] / ((5 - par2[iEtaSl][iBremSl]) * (5 - par2[iEtaSl][iBremSl]));
0234   if (et >= 200)
0235     uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (200 - par2[iEtaSl][iBremSl]) +
0236                   par3[iEtaSl][iBremSl] / ((200 - par2[iEtaSl][iBremSl]) * (200 - par2[iEtaSl][iBremSl]));
0237
0238   if (et > 5 && et < 200)
0239     uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]) +
0240                   par3[iEtaSl][iBremSl] / ((et - par2[iEtaSl][iBremSl]) * (et - par2[iEtaSl][iBremSl]));
0241
0242   return (uncertainty * en);
0243 }