CMS 3D CMS Logo
EnergyUncertaintyElectronSpecific.cc
Go to the documentation of this file.
1 #include "RecoEgamma/EgammaElectronAlgos/interface/EnergyUncertaintyElectronSpecific.h"
2 
3 namespace {
4 
5  double computeEnergyUncertaintyGolden(double eta, double brem, double energy) {
6  double et = energy / cosh(eta);
7 
8  constexpr int nBinsEta = 5;
9  constexpr Double_t EtaBins[nBinsEta + 1] = {0.0, 0.4, 0.8, 1.5, 2.0, 2.5};
10 
11  constexpr int nBinsBrem = 6;
12  constexpr Double_t BremBins[nBinsBrem + 1] = {0.8, 1.0, 1.1, 1.2, 1.3, 1.5, 8.0};
13 
14  constexpr float par0[nBinsEta][nBinsBrem] = {
15  {0.00567891, 0.0065673, 0.00574742, 0.00542964, 0.00523293, 0.00547518},
16 
17  {0.00552517, 0.00611188, 0.0062729, 0.00574846, 0.00447373, 0.00595789},
18 
19  {0.00356679, 0.00503827, 0.00328016, 0.00592303, 0.00512479, 0.00484166},
20 
21  {0.0109195, 0.0102361, 0.0101576, 0.0120683, 0.0155326, 0.0225035},
22 
23  {0.0109632, 0.0103342, 0.0103486, 0.00862762, 0.0111448, 0.0146648}};
24 
25  static_assert(par0[0][0] == 0.00567891f);
26  static_assert(par0[0][1] == 0.0065673f);
27  static_assert(par0[1][3] == 0.00574846f);
28 
29  constexpr float par1[nBinsEta][nBinsBrem] = {{0.238685, 0.193642, 0.249171, 0.259997, 0.310505, 0.390506},
30 
31  {0.288736, 0.312303, 0.294717, 0.294491, 0.379178, 0.38164},
32 
33  {0.456456, 0.394912, 0.541713, 0.401744, 0.483151, 0.657995},
34 
35  {1.13803, 1.39866, 1.51353, 1.48587, 1.49732, 1.82363},
36 
37  {0.458212, 0.628761, 0.659144, 0.929563, 1.06724, 1.6427}};
38 
39  constexpr float par2[nBinsEta][nBinsBrem] = {{2.12035, 3.41493, 1.7401, 1.46234, 0.233226, -2.78168},
40 
41  {1.30552, 0.137905, 0.653793, 0.790746, -1.42584, -2.34653},
42 
43  {0.610716, 0.778879, -1.58577, 1.45098, -0.0985911, -3.47167},
44  {-3.48281, -6.4736, -8.03308, -7.55974, -7.98843, -10.1027},
45 
46  {0.995183, -2.42889, -2.14073, -6.27768, -7.68512, -13.3504}};
47 
48  Int_t iEtaSl = -1;
49  for (Int_t iEta = 0; iEta < nBinsEta; ++iEta) {
50  if (EtaBins[iEta] <= fabs(eta) && fabs(eta) < EtaBins[iEta + 1]) {
51  iEtaSl = iEta;
52  }
53  }
54 
55  Int_t iBremSl = -1;
56  for (Int_t iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
57  if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
58  iBremSl = iBrem;
59  }
60  }
61 
62  if (fabs(eta) > 2.5)
63  iEtaSl = nBinsEta - 1;
64  if (brem < BremBins[0])
65  iBremSl = 0;
66  if (brem > BremBins[nBinsBrem - 1])
67  iBremSl = nBinsBrem - 1;
68 
69  if (iEtaSl == -1) {
70  edm::LogError("BadRange") << "Bad eta value: " << eta << " in computeEnergyUncertaintyGolden";
71  return 0;
72  }
73 
74  if (iBremSl == -1) {
75  edm::LogError("BadRange") << "Bad brem value: " << brem << " in computeEnergyUncertaintyGolden";
76  return 0;
77  }
78 
79  float uncertainty = 0;
80  if (et < 5)
81  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]);
82  if (et > 100)
83  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (100 - par2[iEtaSl][iBremSl]);
84 
85  if (et > 5 && et < 100)
86  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]);
87 
88  return (uncertainty * energy);
89  }
90 
91  double computeEnergyUncertaintyBigbrem(double eta, double brem, double energy) {
92  const double et = energy / cosh(eta);
93 
94  constexpr int nBinsEta = 4;
95  constexpr Double_t EtaBins[nBinsEta + 1] = {0.0, 0.8, 1.5, 2.0, 2.5};
96 
97  constexpr int nBinsBrem = 1;
98  constexpr Double_t BremBins[nBinsBrem + 1] = {0.8, 8.0};
99 
100  constexpr float par0[nBinsEta][nBinsBrem] = {{0.00593389}, {0.00266954}, {0.00500623}, {0.00841038}};
101 
102  constexpr float par1[nBinsEta][nBinsBrem] = {{0.178275}, {0.811415}, {2.34018}, {1.06851}};
103 
104  constexpr float par2[nBinsEta][nBinsBrem] = {{-7.28273}, {-1.66063}, {-11.0129}, {-4.1259}};
105 
106  constexpr float par3[nBinsEta][nBinsBrem] = {{13.2632}, {1.03555}, {-0.200323}, {-0.0646195}};
107 
108  Int_t iEtaSl = -1;
109  for (Int_t iEta = 0; iEta < nBinsEta; ++iEta) {
110  if (EtaBins[iEta] <= fabs(eta) && fabs(eta) < EtaBins[iEta + 1]) {
111  iEtaSl = iEta;
112  }
113  }
114 
115  Int_t iBremSl = -1;
116  for (Int_t iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
117  if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
118  iBremSl = iBrem;
119  }
120  }
121 
122  if (fabs(eta) > 2.5)
123  iEtaSl = nBinsEta - 1;
124  if (brem < BremBins[0])
125  iBremSl = 0;
126  if (brem > BremBins[nBinsBrem - 1])
127  iBremSl = nBinsBrem - 1;
128 
129  if (iEtaSl == -1) {
130  edm::LogError("BadRange") << "Bad eta value: " << eta << " in computeEnergyUncertaintyBigbrem";
131  return 0;
132  }
133 
134  if (iBremSl == -1) {
135  edm::LogError("BadRange") << "Bad brem value: " << brem << " in computeEnergyUncertaintyBigbrem";
136  return 0;
137  }
138 
139  float uncertainty = 0;
140  if (et < 5)
141  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]) +
142  par3[iEtaSl][iBremSl] / ((5 - par2[iEtaSl][iBremSl]) * (5 - par2[iEtaSl][iBremSl]));
143  if (et > 100)
144  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (100 - par2[iEtaSl][iBremSl]) +
145  par3[iEtaSl][iBremSl] / ((100 - par2[iEtaSl][iBremSl]) * (100 - par2[iEtaSl][iBremSl]));
146 
147  if (et > 5 && et < 100)
148  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]) +
149  par3[iEtaSl][iBremSl] / ((et - par2[iEtaSl][iBremSl]) * (et - par2[iEtaSl][iBremSl]));
150 
151  return (uncertainty * energy);
152  }
153  double computeEnergyUncertaintyBadTrack(double eta, double brem, double energy) {
154  const double et = energy / cosh(eta);
155 
156  constexpr int nBinsEta = 4;
157  constexpr Double_t EtaBins[nBinsEta + 1] = {0.0, 0.7, 1.3, 1.8, 2.5};
158 
159  constexpr int nBinsBrem = 1;
160  constexpr Double_t BremBins[nBinsBrem + 1] = {0.8, 8.0};
161 
162  constexpr float par0[nBinsEta][nBinsBrem] = {{0.00601311}, {0.0059814}, {0.00953032}, {0.00728618}};
163 
164  constexpr float par1[nBinsEta][nBinsBrem] = {{0.390988}, {1.02668}, {2.27491}, {2.08268}};
165 
166  constexpr float par2[nBinsEta][nBinsBrem] = {{-4.11919}, {-2.87477}, {-7.61675}, {-8.66756}};
167 
168  constexpr float par3[nBinsEta][nBinsBrem] = {{4.61671}, {0.163447}, {-0.335786}, {-1.27831}};
169 
170  Int_t iEtaSl = -1;
171  for (Int_t iEta = 0; iEta < nBinsEta; ++iEta) {
172  if (EtaBins[iEta] <= fabs(eta) && fabs(eta) < EtaBins[iEta + 1]) {
173  iEtaSl = iEta;
174  }
175  }
176 
177  Int_t iBremSl = -1;
178  for (Int_t iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
179  if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
180  iBremSl = iBrem;
181  }
182  }
183 
184  if (fabs(eta) > 2.5)
185  iEtaSl = nBinsEta - 1;
186  if (brem < BremBins[0])
187  iBremSl = 0;
188  if (brem > BremBins[nBinsBrem - 1])
189  iBremSl = nBinsBrem - 1;
190 
191  if (iEtaSl == -1) {
192  edm::LogError("BadRange") << "Bad eta value: " << eta << " in computeEnergyUncertaintyBadTrack";
193  return 0;
194  }
195 
196  if (iBremSl == -1) {
197  edm::LogError("BadRange") << "Bad brem value: " << brem << " in computeEnergyUncertaintyBadTrack";
198  return 0;
199  }
200 
201  float uncertainty = 0;
202  if (et < 5)
203  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]) +
204  par3[iEtaSl][iBremSl] / ((5 - par2[iEtaSl][iBremSl]) * (5 - par2[iEtaSl][iBremSl]));
205  if (et > 100)
206  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (100 - par2[iEtaSl][iBremSl]) +
207  par3[iEtaSl][iBremSl] / ((100 - par2[iEtaSl][iBremSl]) * (100 - par2[iEtaSl][iBremSl]));
208 
209  if (et > 5 && et < 100)
210  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]) +
211  par3[iEtaSl][iBremSl] / ((et - par2[iEtaSl][iBremSl]) * (et - par2[iEtaSl][iBremSl]));
212 
213  return (uncertainty * energy);
214  }
215 
216  double computeEnergyUncertaintyShowering(double eta, double brem, double energy) {
217  const double et = energy / cosh(eta);
218 
219  constexpr int nBinsEta = 4;
220  constexpr Double_t EtaBins[nBinsEta + 1] = {0.0, 0.8, 1.2, 1.7, 2.5};
221 
222  constexpr int nBinsBrem = 5;
223  constexpr Double_t BremBins[nBinsBrem + 1] = {0.8, 1.8, 2.2, 3.0, 4.0, 8.0};
224 
225  constexpr float par0[nBinsEta][nBinsBrem] = {{0.0049351, 0.00566155, 0.0051397, 0.00468481, 0.00444475},
226 
227  {0.00201762, 0.00431475, 0.00501004, 0.00632666, 0.00636704},
228 
229  {-0.00729396, 0.00539783, 0.00608149, 0.00465335, 0.00642685},
230 
231  {0.0149449, 0.0216691, 0.0255957, 0.0206101, 0.0180508}};
232 
233  constexpr float par1[nBinsEta][nBinsBrem] = {{0.579925, 0.496137, 0.551947, 0.63011, 0.684261},
234 
235  {0.914762, 0.824483, 0.888521, 0.960241, 1.25728},
236 
237  {3.24295, 1.72935, 1.80606, 2.13562, 2.07592},
238 
239  {1.00448, 1.18393, 0.00775295, 2.59246, 3.1099}};
240 
241  constexpr float par2[nBinsEta][nBinsBrem] = {{-9.33987, -5.52543, -7.30079, -6.7722, -4.67614},
242 
243  {-4.48042, -5.02885, -4.77311, -3.36742, -5.53561},
244 
245  {-17.1458, -5.92807, -6.67563, -10.1105, -7.50257},
246 
247  {-2.09368, -4.56674, -44.2722, -13.1702, -13.6208}};
248 
249  constexpr float par3[nBinsEta][nBinsBrem] = {{1.62129, 1.19101, 1.89701, 1.81614, 1.64415},
250 
251  {-1.50473, -0.153502, -0.355145, -1.16499, -0.864123},
252 
253  {-4.69711, -2.18733, -0.922401, -0.230781, -2.91515},
254 
255  {0.455037, -0.601872, 241.516, -2.35024, -2.11069}};
256 
257  Int_t iEtaSl = -1;
258  for (Int_t iEta = 0; iEta < nBinsEta; ++iEta) {
259  if (EtaBins[iEta] <= fabs(eta) && fabs(eta) < EtaBins[iEta + 1]) {
260  iEtaSl = iEta;
261  }
262  }
263 
264  Int_t iBremSl = -1;
265  for (Int_t iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
266  if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
267  iBremSl = iBrem;
268  }
269  }
270 
271  if (fabs(eta) > 2.5)
272  iEtaSl = nBinsEta - 1;
273  if (brem < BremBins[0])
274  iBremSl = 0;
275  if (brem > BremBins[nBinsBrem - 1])
276  iBremSl = nBinsBrem - 1;
277 
278  if (iEtaSl == -1) {
279  edm::LogError("BadRange") << "Bad eta value: " << eta << " in computeEnergyUncertaintyShowering";
280  return 0;
281  }
282 
283  if (iBremSl == -1) {
284  edm::LogError("BadRange") << "Bad brem value: " << brem << " in computeEnergyUncertaintyShowering";
285  return 0;
286  }
287 
288  float uncertainty = 0;
289  if (et < 5)
290  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]) +
291  par3[iEtaSl][iBremSl] / ((5 - par2[iEtaSl][iBremSl]) * (5 - par2[iEtaSl][iBremSl]));
292  if (et > 100)
293  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (100 - par2[iEtaSl][iBremSl]) +
294  par3[iEtaSl][iBremSl] / ((100 - par2[iEtaSl][iBremSl]) * (100 - par2[iEtaSl][iBremSl]));
295 
296  if (et > 5 && et < 100)
297  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]) +
298  par3[iEtaSl][iBremSl] / ((et - par2[iEtaSl][iBremSl]) * (et - par2[iEtaSl][iBremSl]));
299 
300  return (uncertainty * energy);
301  }
302 
303  double computeEnergyUncertaintyCracks(double eta, double brem, double energy) {
304  const double et = energy / cosh(eta);
305 
306  constexpr int nBinsEta = 5;
307  constexpr Double_t EtaBins[nBinsEta + 1] = {0.0, 0.42, 0.78, 1.2, 1.52, 1.65};
308 
309  constexpr int nBinsBrem = 6;
310  constexpr Double_t BremBins[nBinsBrem + 1] = {0.8, 1.2, 1.5, 2.1, 3., 4, 8.0};
311 
312  constexpr float par0[nBinsEta][nBinsBrem] = {
313  {0.0139815, 0.00550839, 0.0108292, 0.00596201, -0.00498136, 0.000621696},
314 
315  {0.00467498, 0.00808463, 0.00546665, 0.00506318, 0.00608425, -4.45641e-06},
316 
317  {0.00971734, 0.00063951, -0.0121618, -0.00604365, 0.00492161, -0.00143907},
318 
319  {-0.0844907, -0.0592498, -0.0828631, -0.0740798, -0.0698045, -0.0699518},
320 
321  {-0.0999971, -0.0999996, -0.0989356, -0.0999965, -0.0833049, -0.020072}};
322 
323  constexpr float par1[nBinsEta][nBinsBrem] = {{0.569273, 0.674654, 0.523128, 1.02501, 1.75645, 0.955191},
324 
325  {0.697951, 0.580628, 0.814515, 0.819975, 0.829616, 1.18952},
326 
327  {3.79446, 2.47472, 5.12931, 3.42497, 1.84123, 2.3773},
328 
329  {19.9999, 10.4079, 16.6273, 15.9316, 15.4883, 14.7306},
330 
331  {15.9122, 18.5882, 19.9996, 19.9999, 18.2281, 8.1587}};
332 
333  constexpr float par2[nBinsEta][nBinsBrem] = {{-4.31243, -3.071, -2.56702, -7.74555, -21.3726, -6.2189},
334 
335  {-6.56009, -3.66067, -7.8275, -6.01641, -7.85456, -8.27071},
336 
337  {-49.9996, -25.0724, -49.985, -28.1932, -10.6485, -15.4014},
338 
339  {-39.9444, -25.1133, -49.9999, -50, -49.9998, -49.9998},
340 
341  {-30.1268, -42.6113, -46.6999, -47.074, -49.9995, -25.2897}};
342 
343  static_assert(par0[0][3] == 0.00596201f);
344  static_assert(par1[0][3] == 1.02501f);
345  static_assert(par2[0][3] == -7.74555f);
346 
347  static_assert(par0[2][4] == 0.00492161f);
348  static_assert(par1[2][4] == 1.84123f);
349  static_assert(par2[2][4] == -10.6485f);
350 
351  static_assert(par0[4][3] == -0.0999965f);
352  static_assert(par1[4][3] == 19.9999f);
353  static_assert(par2[4][3] == -47.074f);
354 
355  Int_t iEtaSl = -1;
356  for (Int_t iEta = 0; iEta < nBinsEta; ++iEta) {
357  if (EtaBins[iEta] <= fabs(eta) && fabs(eta) < EtaBins[iEta + 1]) {
358  iEtaSl = iEta;
359  }
360  }
361 
362  Int_t iBremSl = -1;
363  for (Int_t iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
364  if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
365  iBremSl = iBrem;
366  }
367  }
368 
369  if (fabs(eta) > 2.5)
370  iEtaSl = nBinsEta - 1;
371  if (brem < BremBins[0])
372  iBremSl = 0;
373  if (brem > BremBins[nBinsBrem - 1])
374  iBremSl = nBinsBrem - 1;
375 
376  if (iEtaSl == -1) {
377  edm::LogError("BadRange") << "Bad eta value: " << eta << " in computeEnergyUncertaintyCracks";
378  return 0;
379  }
380 
381  if (iBremSl == -1) {
382  edm::LogError("BadRange") << "Bad brem value: " << brem << " in computeEnergyUncertaintyCracks";
383  return 0;
384  }
385 
386  float uncertainty = 0;
387  if (et < 5)
388  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]);
389  if (et > 100)
390  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (100 - par2[iEtaSl][iBremSl]);
391 
392  if (et > 5 && et < 100)
393  uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]);
394 
395  return (uncertainty * energy);
396  }
397 
398 } // namespace
399 
400 using reco::GsfElectron;
401 
402 double egamma::electronEnergyUncertainty(GsfElectron::Classification c, double eta, double brem, double energy) {
403  if (c == GsfElectron::GOLDEN)
404  return computeEnergyUncertaintyGolden(eta, brem, energy);
405  if (c == GsfElectron::BIGBREM)
406  return computeEnergyUncertaintyBigbrem(eta, brem, energy);
407  if (c == GsfElectron::SHOWERING)
408  return computeEnergyUncertaintyShowering(eta, brem, energy);
409  if (c == GsfElectron::BADTRACK)
410  return computeEnergyUncertaintyBadTrack(eta, brem, energy);
411  if (c == GsfElectron::GAP)
412  return computeEnergyUncertaintyCracks(eta, brem, energy);
413  throw cms::Exception("GsfElectronAlgo|InternalError") << "unknown classification";
414 }
EgHLTOffHistBins_cfi.et
et
Definition: EgHLTOffHistBins_cfi.py:8
egamma::electronEnergyUncertainty
double electronEnergyUncertainty(reco::GsfElectron::Classification c, double eta, double brem, double energy)
Definition: EnergyUncertaintyElectronSpecific.cc:402
pfMETCorrectionType0_cfi.par3
par3
Definition: pfMETCorrectionType0_cfi.py:50
HltBtagPostValidation_cff.c
c
Definition: HltBtagPostValidation_cff.py:31
pfMETCorrectionType0_cfi.par1
par1
Definition: pfMETCorrectionType0_cfi.py:48
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
fitWZ.par0
par0
Definition: fitWZ.py:52
pfMETCorrectionType0_cfi.par2
par2
Definition: pfMETCorrectionType0_cfi.py:49
HCALHighEnergyHPDFilter_cfi.energy
energy
Definition: HCALHighEnergyHPDFilter_cfi.py:5
MuonErrorMatrixAnalyzer_cfi.EtaBins
EtaBins
Definition: MuonErrorMatrixAnalyzer_cfi.py:15
hgcalTowerMapProducer_cfi.nBinsEta
nBinsEta
Definition: hgcalTowerMapProducer_cfi.py:9
constexpr
#define constexpr
Definition: GCC11Compatibility.h:35
Generated for CMSSW Reference Manual by doxygen 1.8.11