File indexing completed on 2024-04-06 12:24:47
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 }
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 }
