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