Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-07 22:33:46

0001 #include "RecoEgamma/EgammaElectronAlgos/interface/EnergyUncertaintyElectronSpecific.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 
0004 namespace {
0005 
0006   double computeEnergyUncertaintyGolden(double eta, double brem, double energy) {
0007     double et = energy / cosh(eta);
0008 
0009     constexpr int nBinsEta = 5;
0010     constexpr double EtaBins[nBinsEta + 1] = {0.0, 0.4, 0.8, 1.5, 2.0, 2.5};
0011 
0012     constexpr int nBinsBrem = 6;
0013     constexpr double BremBins[nBinsBrem + 1] = {0.8, 1.0, 1.1, 1.2, 1.3, 1.5, 8.0};
0014 
0015     constexpr float par0[nBinsEta][nBinsBrem] = {
0016         {0.00567891, 0.0065673, 0.00574742, 0.00542964, 0.00523293, 0.00547518},
0017 
0018         {0.00552517, 0.00611188, 0.0062729, 0.00574846, 0.00447373, 0.00595789},
0019 
0020         {0.00356679, 0.00503827, 0.00328016, 0.00592303, 0.00512479, 0.00484166},
0021 
0022         {0.0109195, 0.0102361, 0.0101576, 0.0120683, 0.0155326, 0.0225035},
0023 
0024         {0.0109632, 0.0103342, 0.0103486, 0.00862762, 0.0111448, 0.0146648}};
0025 
0026     static_assert(par0[0][0] == 0.00567891f);
0027     static_assert(par0[0][1] == 0.0065673f);
0028     static_assert(par0[1][3] == 0.00574846f);
0029 
0030     constexpr float par1[nBinsEta][nBinsBrem] = {{0.238685, 0.193642, 0.249171, 0.259997, 0.310505, 0.390506},
0031 
0032                                                  {0.288736, 0.312303, 0.294717, 0.294491, 0.379178, 0.38164},
0033 
0034                                                  {0.456456, 0.394912, 0.541713, 0.401744, 0.483151, 0.657995},
0035 
0036                                                  {1.13803, 1.39866, 1.51353, 1.48587, 1.49732, 1.82363},
0037 
0038                                                  {0.458212, 0.628761, 0.659144, 0.929563, 1.06724, 1.6427}};
0039 
0040     constexpr float par2[nBinsEta][nBinsBrem] = {{2.12035, 3.41493, 1.7401, 1.46234, 0.233226, -2.78168},
0041 
0042                                                  {1.30552, 0.137905, 0.653793, 0.790746, -1.42584, -2.34653},
0043 
0044                                                  {0.610716, 0.778879, -1.58577, 1.45098, -0.0985911, -3.47167},
0045                                                  {-3.48281, -6.4736, -8.03308, -7.55974, -7.98843, -10.1027},
0046 
0047                                                  {0.995183, -2.42889, -2.14073, -6.27768, -7.68512, -13.3504}};
0048 
0049     Int_t iEtaSl = -1;
0050     for (Int_t iEta = 0; iEta < nBinsEta; ++iEta) {
0051       if (EtaBins[iEta] <= fabs(eta) && fabs(eta) < EtaBins[iEta + 1]) {
0052         iEtaSl = iEta;
0053       }
0054     }
0055 
0056     Int_t iBremSl = -1;
0057     for (Int_t iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
0058       if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
0059         iBremSl = iBrem;
0060       }
0061     }
0062 
0063     if (fabs(eta) > 2.5)
0064       iEtaSl = nBinsEta - 1;
0065     if (brem < BremBins[0])
0066       iBremSl = 0;
0067     if (brem > BremBins[nBinsBrem - 1])
0068       iBremSl = nBinsBrem - 1;
0069 
0070     if (iEtaSl == -1) {
0071       edm::LogError("BadRange") << "Bad eta value: " << eta << " in computeEnergyUncertaintyGolden";
0072       return 0;
0073     }
0074 
0075     if (iBremSl == -1) {
0076       edm::LogError("BadRange") << "Bad brem value: " << brem << " in computeEnergyUncertaintyGolden";
0077       return 0;
0078     }
0079 
0080     float uncertainty = 0;
0081     if (et < 5)
0082       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]);
0083     if (et > 100)
0084       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (100 - par2[iEtaSl][iBremSl]);
0085 
0086     if (et > 5 && et < 100)
0087       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]);
0088 
0089     return (uncertainty * energy);
0090   }
0091 
0092   double computeEnergyUncertaintyBigbrem(double eta, double brem, double energy) {
0093     const double et = energy / cosh(eta);
0094 
0095     constexpr int nBinsEta = 4;
0096     constexpr double EtaBins[nBinsEta + 1] = {0.0, 0.8, 1.5, 2.0, 2.5};
0097 
0098     constexpr int nBinsBrem = 1;
0099     constexpr double BremBins[nBinsBrem + 1] = {0.8, 8.0};
0100 
0101     constexpr float par0[nBinsEta][nBinsBrem] = {{0.00593389}, {0.00266954}, {0.00500623}, {0.00841038}};
0102 
0103     constexpr float par1[nBinsEta][nBinsBrem] = {{0.178275}, {0.811415}, {2.34018}, {1.06851}};
0104 
0105     constexpr float par2[nBinsEta][nBinsBrem] = {{-7.28273}, {-1.66063}, {-11.0129}, {-4.1259}};
0106 
0107     constexpr float par3[nBinsEta][nBinsBrem] = {{13.2632}, {1.03555}, {-0.200323}, {-0.0646195}};
0108 
0109     Int_t iEtaSl = -1;
0110     for (Int_t iEta = 0; iEta < nBinsEta; ++iEta) {
0111       if (EtaBins[iEta] <= fabs(eta) && fabs(eta) < EtaBins[iEta + 1]) {
0112         iEtaSl = iEta;
0113       }
0114     }
0115 
0116     Int_t iBremSl = -1;
0117     for (Int_t iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
0118       if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
0119         iBremSl = iBrem;
0120       }
0121     }
0122 
0123     if (fabs(eta) > 2.5)
0124       iEtaSl = nBinsEta - 1;
0125     if (brem < BremBins[0])
0126       iBremSl = 0;
0127     if (brem > BremBins[nBinsBrem - 1])
0128       iBremSl = nBinsBrem - 1;
0129 
0130     if (iEtaSl == -1) {
0131       edm::LogError("BadRange") << "Bad eta value: " << eta << " in computeEnergyUncertaintyBigbrem";
0132       return 0;
0133     }
0134 
0135     if (iBremSl == -1) {
0136       edm::LogError("BadRange") << "Bad brem value: " << brem << " in computeEnergyUncertaintyBigbrem";
0137       return 0;
0138     }
0139 
0140     float uncertainty = 0;
0141     if (et < 5)
0142       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]) +
0143                     par3[iEtaSl][iBremSl] / ((5 - par2[iEtaSl][iBremSl]) * (5 - par2[iEtaSl][iBremSl]));
0144     if (et > 100)
0145       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (100 - par2[iEtaSl][iBremSl]) +
0146                     par3[iEtaSl][iBremSl] / ((100 - par2[iEtaSl][iBremSl]) * (100 - par2[iEtaSl][iBremSl]));
0147 
0148     if (et > 5 && et < 100)
0149       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]) +
0150                     par3[iEtaSl][iBremSl] / ((et - par2[iEtaSl][iBremSl]) * (et - par2[iEtaSl][iBremSl]));
0151 
0152     return (uncertainty * energy);
0153   }
0154   double computeEnergyUncertaintyBadTrack(double eta, double brem, double energy) {
0155     const double et = energy / cosh(eta);
0156 
0157     constexpr int nBinsEta = 4;
0158     constexpr double EtaBins[nBinsEta + 1] = {0.0, 0.7, 1.3, 1.8, 2.5};
0159 
0160     constexpr int nBinsBrem = 1;
0161     constexpr double BremBins[nBinsBrem + 1] = {0.8, 8.0};
0162 
0163     constexpr float par0[nBinsEta][nBinsBrem] = {{0.00601311}, {0.0059814}, {0.00953032}, {0.00728618}};
0164 
0165     constexpr float par1[nBinsEta][nBinsBrem] = {{0.390988}, {1.02668}, {2.27491}, {2.08268}};
0166 
0167     constexpr float par2[nBinsEta][nBinsBrem] = {{-4.11919}, {-2.87477}, {-7.61675}, {-8.66756}};
0168 
0169     constexpr float par3[nBinsEta][nBinsBrem] = {{4.61671}, {0.163447}, {-0.335786}, {-1.27831}};
0170 
0171     Int_t iEtaSl = -1;
0172     for (Int_t iEta = 0; iEta < nBinsEta; ++iEta) {
0173       if (EtaBins[iEta] <= fabs(eta) && fabs(eta) < EtaBins[iEta + 1]) {
0174         iEtaSl = iEta;
0175       }
0176     }
0177 
0178     Int_t iBremSl = -1;
0179     for (Int_t iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
0180       if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
0181         iBremSl = iBrem;
0182       }
0183     }
0184 
0185     if (fabs(eta) > 2.5)
0186       iEtaSl = nBinsEta - 1;
0187     if (brem < BremBins[0])
0188       iBremSl = 0;
0189     if (brem > BremBins[nBinsBrem - 1])
0190       iBremSl = nBinsBrem - 1;
0191 
0192     if (iEtaSl == -1) {
0193       edm::LogError("BadRange") << "Bad eta value: " << eta << " in computeEnergyUncertaintyBadTrack";
0194       return 0;
0195     }
0196 
0197     if (iBremSl == -1) {
0198       edm::LogError("BadRange") << "Bad brem value: " << brem << " in computeEnergyUncertaintyBadTrack";
0199       return 0;
0200     }
0201 
0202     float uncertainty = 0;
0203     if (et < 5)
0204       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]) +
0205                     par3[iEtaSl][iBremSl] / ((5 - par2[iEtaSl][iBremSl]) * (5 - par2[iEtaSl][iBremSl]));
0206     if (et > 100)
0207       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (100 - par2[iEtaSl][iBremSl]) +
0208                     par3[iEtaSl][iBremSl] / ((100 - par2[iEtaSl][iBremSl]) * (100 - par2[iEtaSl][iBremSl]));
0209 
0210     if (et > 5 && et < 100)
0211       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]) +
0212                     par3[iEtaSl][iBremSl] / ((et - par2[iEtaSl][iBremSl]) * (et - par2[iEtaSl][iBremSl]));
0213 
0214     return (uncertainty * energy);
0215   }
0216 
0217   double computeEnergyUncertaintyShowering(double eta, double brem, double energy) {
0218     const double et = energy / cosh(eta);
0219 
0220     constexpr int nBinsEta = 4;
0221     constexpr double EtaBins[nBinsEta + 1] = {0.0, 0.8, 1.2, 1.7, 2.5};
0222 
0223     constexpr int nBinsBrem = 5;
0224     constexpr double BremBins[nBinsBrem + 1] = {0.8, 1.8, 2.2, 3.0, 4.0, 8.0};
0225 
0226     constexpr float par0[nBinsEta][nBinsBrem] = {{0.0049351, 0.00566155, 0.0051397, 0.00468481, 0.00444475},
0227 
0228                                                  {0.00201762, 0.00431475, 0.00501004, 0.00632666, 0.00636704},
0229 
0230                                                  {-0.00729396, 0.00539783, 0.00608149, 0.00465335, 0.00642685},
0231 
0232                                                  {0.0149449, 0.0216691, 0.0255957, 0.0206101, 0.0180508}};
0233 
0234     constexpr float par1[nBinsEta][nBinsBrem] = {{0.579925, 0.496137, 0.551947, 0.63011, 0.684261},
0235 
0236                                                  {0.914762, 0.824483, 0.888521, 0.960241, 1.25728},
0237 
0238                                                  {3.24295, 1.72935, 1.80606, 2.13562, 2.07592},
0239 
0240                                                  {1.00448, 1.18393, 0.00775295, 2.59246, 3.1099}};
0241 
0242     constexpr float par2[nBinsEta][nBinsBrem] = {{-9.33987, -5.52543, -7.30079, -6.7722, -4.67614},
0243 
0244                                                  {-4.48042, -5.02885, -4.77311, -3.36742, -5.53561},
0245 
0246                                                  {-17.1458, -5.92807, -6.67563, -10.1105, -7.50257},
0247 
0248                                                  {-2.09368, -4.56674, -44.2722, -13.1702, -13.6208}};
0249 
0250     constexpr float par3[nBinsEta][nBinsBrem] = {{1.62129, 1.19101, 1.89701, 1.81614, 1.64415},
0251 
0252                                                  {-1.50473, -0.153502, -0.355145, -1.16499, -0.864123},
0253 
0254                                                  {-4.69711, -2.18733, -0.922401, -0.230781, -2.91515},
0255 
0256                                                  {0.455037, -0.601872, 241.516, -2.35024, -2.11069}};
0257 
0258     Int_t iEtaSl = -1;
0259     for (Int_t iEta = 0; iEta < nBinsEta; ++iEta) {
0260       if (EtaBins[iEta] <= fabs(eta) && fabs(eta) < EtaBins[iEta + 1]) {
0261         iEtaSl = iEta;
0262       }
0263     }
0264 
0265     Int_t iBremSl = -1;
0266     for (Int_t iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
0267       if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
0268         iBremSl = iBrem;
0269       }
0270     }
0271 
0272     if (fabs(eta) > 2.5)
0273       iEtaSl = nBinsEta - 1;
0274     if (brem < BremBins[0])
0275       iBremSl = 0;
0276     if (brem > BremBins[nBinsBrem - 1])
0277       iBremSl = nBinsBrem - 1;
0278 
0279     if (iEtaSl == -1) {
0280       edm::LogError("BadRange") << "Bad eta value: " << eta << " in computeEnergyUncertaintyShowering";
0281       return 0;
0282     }
0283 
0284     if (iBremSl == -1) {
0285       edm::LogError("BadRange") << "Bad brem value: " << brem << " in computeEnergyUncertaintyShowering";
0286       return 0;
0287     }
0288 
0289     float uncertainty = 0;
0290     if (et < 5)
0291       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]) +
0292                     par3[iEtaSl][iBremSl] / ((5 - par2[iEtaSl][iBremSl]) * (5 - par2[iEtaSl][iBremSl]));
0293     if (et > 100)
0294       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (100 - par2[iEtaSl][iBremSl]) +
0295                     par3[iEtaSl][iBremSl] / ((100 - par2[iEtaSl][iBremSl]) * (100 - par2[iEtaSl][iBremSl]));
0296 
0297     if (et > 5 && et < 100)
0298       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]) +
0299                     par3[iEtaSl][iBremSl] / ((et - par2[iEtaSl][iBremSl]) * (et - par2[iEtaSl][iBremSl]));
0300 
0301     return (uncertainty * energy);
0302   }
0303 
0304   double computeEnergyUncertaintyCracks(double eta, double brem, double energy) {
0305     const double et = energy / cosh(eta);
0306 
0307     constexpr int nBinsEta = 5;
0308     constexpr double EtaBins[nBinsEta + 1] = {0.0, 0.42, 0.78, 1.2, 1.52, 1.65};
0309 
0310     constexpr int nBinsBrem = 6;
0311     constexpr double BremBins[nBinsBrem + 1] = {0.8, 1.2, 1.5, 2.1, 3., 4, 8.0};
0312 
0313     constexpr float par0[nBinsEta][nBinsBrem] = {
0314         {0.0139815, 0.00550839, 0.0108292, 0.00596201, -0.00498136, 0.000621696},
0315 
0316         {0.00467498, 0.00808463, 0.00546665, 0.00506318, 0.00608425, -4.45641e-06},
0317 
0318         {0.00971734, 0.00063951, -0.0121618, -0.00604365, 0.00492161, -0.00143907},
0319 
0320         {-0.0844907, -0.0592498, -0.0828631, -0.0740798, -0.0698045, -0.0699518},
0321 
0322         {-0.0999971, -0.0999996, -0.0989356, -0.0999965, -0.0833049, -0.020072}};
0323 
0324     constexpr float par1[nBinsEta][nBinsBrem] = {{0.569273, 0.674654, 0.523128, 1.02501, 1.75645, 0.955191},
0325 
0326                                                  {0.697951, 0.580628, 0.814515, 0.819975, 0.829616, 1.18952},
0327 
0328                                                  {3.79446, 2.47472, 5.12931, 3.42497, 1.84123, 2.3773},
0329 
0330                                                  {19.9999, 10.4079, 16.6273, 15.9316, 15.4883, 14.7306},
0331 
0332                                                  {15.9122, 18.5882, 19.9996, 19.9999, 18.2281, 8.1587}};
0333 
0334     constexpr float par2[nBinsEta][nBinsBrem] = {{-4.31243, -3.071, -2.56702, -7.74555, -21.3726, -6.2189},
0335 
0336                                                  {-6.56009, -3.66067, -7.8275, -6.01641, -7.85456, -8.27071},
0337 
0338                                                  {-49.9996, -25.0724, -49.985, -28.1932, -10.6485, -15.4014},
0339 
0340                                                  {-39.9444, -25.1133, -49.9999, -50, -49.9998, -49.9998},
0341 
0342                                                  {-30.1268, -42.6113, -46.6999, -47.074, -49.9995, -25.2897}};
0343 
0344     static_assert(par0[0][3] == 0.00596201f);
0345     static_assert(par1[0][3] == 1.02501f);
0346     static_assert(par2[0][3] == -7.74555f);
0347 
0348     static_assert(par0[2][4] == 0.00492161f);
0349     static_assert(par1[2][4] == 1.84123f);
0350     static_assert(par2[2][4] == -10.6485f);
0351 
0352     static_assert(par0[4][3] == -0.0999965f);
0353     static_assert(par1[4][3] == 19.9999f);
0354     static_assert(par2[4][3] == -47.074f);
0355 
0356     Int_t iEtaSl = -1;
0357     for (Int_t iEta = 0; iEta < nBinsEta; ++iEta) {
0358       if (EtaBins[iEta] <= fabs(eta) && fabs(eta) < EtaBins[iEta + 1]) {
0359         iEtaSl = iEta;
0360       }
0361     }
0362 
0363     Int_t iBremSl = -1;
0364     for (Int_t iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
0365       if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
0366         iBremSl = iBrem;
0367       }
0368     }
0369 
0370     if (fabs(eta) > 2.5)
0371       iEtaSl = nBinsEta - 1;
0372     if (brem < BremBins[0])
0373       iBremSl = 0;
0374     if (brem > BremBins[nBinsBrem - 1])
0375       iBremSl = nBinsBrem - 1;
0376 
0377     if (iEtaSl == -1) {
0378       edm::LogError("BadRange") << "Bad eta value: " << eta << " in computeEnergyUncertaintyCracks";
0379       return 0;
0380     }
0381 
0382     if (iBremSl == -1) {
0383       edm::LogError("BadRange") << "Bad brem value: " << brem << " in computeEnergyUncertaintyCracks";
0384       return 0;
0385     }
0386 
0387     float uncertainty = 0;
0388     if (et < 5)
0389       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]);
0390     if (et > 100)
0391       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (100 - par2[iEtaSl][iBremSl]);
0392 
0393     if (et > 5 && et < 100)
0394       uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]);
0395 
0396     return (uncertainty * energy);
0397   }
0398 
0399 }  // namespace
0400 
0401 using reco::GsfElectron;
0402 
0403 double egamma::electronEnergyUncertainty(GsfElectron::Classification c, double eta, double brem, double energy) {
0404   if (c == GsfElectron::GOLDEN)
0405     return computeEnergyUncertaintyGolden(eta, brem, energy);
0406   if (c == GsfElectron::BIGBREM)
0407     return computeEnergyUncertaintyBigbrem(eta, brem, energy);
0408   if (c == GsfElectron::SHOWERING)
0409     return computeEnergyUncertaintyShowering(eta, brem, energy);
0410   if (c == GsfElectron::BADTRACK)
0411     return computeEnergyUncertaintyBadTrack(eta, brem, energy);
0412   if (c == GsfElectron::GAP)
0413     return computeEnergyUncertaintyCracks(eta, brem, energy);
0414   throw cms::Exception("GsfElectronAlgo|InternalError") << "unknown classification";
0415 }