File indexing completed on 2024-04-06 12:24:58
0001 #include "RecoEgamma/EgammaPhotonAlgos/interface/EnergyUncertaintyPhotonSpecific.h"
0002
0003 #include <iostream>
0004
0005 EnergyUncertaintyPhotonSpecific::EnergyUncertaintyPhotonSpecific(const edm::ParameterSet& config) {}
0006
0007 EnergyUncertaintyPhotonSpecific::~EnergyUncertaintyPhotonSpecific() {}
0008
0009 void EnergyUncertaintyPhotonSpecific::init(const edm::EventSetup& theEventSetup) {}
0010
0011 double EnergyUncertaintyPhotonSpecific::computePhotonEnergyUncertainty_lowR9(double eta, double brem, double energy) {
0012 double et = energy / cosh(eta);
0013
0014 constexpr int nBinsEta = 6;
0015 const double EtaBins[nBinsEta + 1] = {0.0, 0.7, 1.15, 1.44, 1.56, 2.0, 2.5};
0016
0017 constexpr int nBinsBrem = 6;
0018 const double BremBins[nBinsBrem + 1] = {0.8, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0};
0019
0020 constexpr std::array<std::array<float, nBinsBrem>, nBinsEta> par0 = {
0021 {{{0.0232291, 0.00703187, 0.00692465, 0.00855993, 0.00795058, 0.0107494}},
0022 {{0.0614866, 0.00894211, 0.0102959, 0.0128934, 0.0130199, 0.0180839}},
0023 {{0.0291343, 0.00876269, 0.0120863, 0.0112655, 0.0168267, 0.0168059}},
0024 {{0.158403, 0.0717431, 0.0385666, 0.0142631, 0.0421638, 0.046331}},
0025 {{0.0483944, 0.0168516, 0.0243039, 0.031795, 0.0414953, 0.058031}},
0026 {{0.107158, 0.021685, 0.0196619, 0.0324734, 0.0414953, 0.058031}}}};
0027
0028 constexpr std::array<std::array<float, nBinsBrem>, nBinsEta> par1 = {
0029 {{{0., 0.646644, 0.292698, 0.280843, 0.370007, 0.276159}},
0030 {{0., 0.466937, 0.313568, 0.302943, 0.505135, 0.382134}},
0031 {{0., 0.375159, 0.397635, 0.856565, 0.636468, 1.09268}},
0032 {{0., 1.66981, 3.6319, 8.85991, 3.1289, 1.29951}},
0033 {{0., 1.19617, 0.994626, 0.875925, 0.654605, 0.292915}},
0034 {{0., 0.574207, 0.940217, 0.574766, 0.654605, 0.292915}}}};
0035
0036 constexpr std::array<std::array<float, nBinsBrem>, nBinsEta> par2 = {
0037 {{{0., -7.4698, 4.16907, 4.25527, 3.03429, 4.44532}},
0038 {{0., 3.33434, 6.34301, 6.35598, 2.52964, 5.3388}},
0039 {{0., 7.11411, 5.97451, -5.76122, -1.54548, -0.547554}},
0040 {{0., 6.86275, -3.76633, -32.6073, -6.58653, 1.76117}},
0041 {{0., -6.78666, -4.26073, 1.43183, 4.45367, 8.48307}},
0042 {{0., -0.566981, -6.05845, -5.23571, 4.45367, 8.48307}}}};
0043
0044 constexpr std::array<std::array<float, nBinsBrem>, nBinsEta> par3 = {
0045 {{{0., 5.53373e-08, 5.61149e-06, 9.6404e-07, 4.43986e-07, 2.58822e-06}},
0046 {{0., 0.000114835, 2.86726e-07, 0.00190694, 0.120204, 3.59921e-07}},
0047 {{0., 0.0438575, 0.0469782, 4.99993, 4.99992, 0.0952985}},
0048 {{0., 0.00543544, 6.56718e-05, 0.00119538, 1.10125e-05, 0.00204206}},
0049 {{0., 4.98192, 4.99984, 0.0920944, 0.030385, 0.0134321}},
0050 {{0., 0.0120609, 0.000193818, 4.9419, 0.030385, 0.0134321}}}};
0051
0052 int iEtaSl = -1;
0053 for (int iEta = 0; iEta < nBinsEta; ++iEta) {
0054 if (EtaBins[iEta] <= std::abs(eta) && std::abs(eta) < EtaBins[iEta + 1]) {
0055 iEtaSl = iEta;
0056 }
0057 }
0058
0059 int iBremSl = -1;
0060 for (int iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
0061 if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
0062 iBremSl = iBrem;
0063 }
0064 }
0065
0066 if (std::abs(eta) > 2.5)
0067 iEtaSl = nBinsEta - 1;
0068 if (brem < BremBins[0])
0069 iBremSl = 0;
0070 if (brem > BremBins[nBinsBrem - 1])
0071 iBremSl = nBinsBrem - 1;
0072
0073 float uncertainty = 0;
0074 if (iBremSl >= 0 && iBremSl < nBinsBrem && iEtaSl >= 0 && iEtaSl < nBinsEta) {
0075 if (et < 5)
0076 uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]) +
0077 par3[iEtaSl][iBremSl] / ((5 - par2[iEtaSl][iBremSl]) * (5 - par2[iEtaSl][iBremSl]));
0078 if (et > 200)
0079 uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (200 - par2[iEtaSl][iBremSl]) +
0080 par3[iEtaSl][iBremSl] / ((200 - par2[iEtaSl][iBremSl]) * (200 - par2[iEtaSl][iBremSl]));
0081
0082 if (et > 5 && et < 200)
0083 uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]) +
0084 par3[iEtaSl][iBremSl] / ((et - par2[iEtaSl][iBremSl]) * (et - par2[iEtaSl][iBremSl]));
0085 }
0086
0087 return (uncertainty * energy);
0088 }
0089
0090 double EnergyUncertaintyPhotonSpecific::computePhotonEnergyUncertainty_highR9(double eta, double brem, double energy) {
0091 double et = energy / cosh(eta);
0092
0093 constexpr int nBinsEta = 6;
0094 const double EtaBins[nBinsEta + 1] = {0.0, 0.7, 1.15, 1.44, 1.56, 2.0, 2.5};
0095
0096 constexpr int nBinsBrem = 2;
0097 const double BremBins[nBinsBrem + 1] = {0.8, 1.0, 2.0};
0098
0099 constexpr std::array<std::array<float, nBinsBrem>, nBinsEta> par0 = {{{{0.00806753, 0.00899298}},
0100 {{0.00880649, 0.00972275}},
0101 {{0.0101474, 0.0109109}},
0102 {{0.00343003, 0.0372159}},
0103 {{0.0192411, 0.0195124}},
0104 {{0.0203644, 0.0198718}}}};
0105
0106 constexpr std::array<std::array<float, nBinsBrem>, nBinsEta> par1 = {{{{0.143754, 0.10159}},
0107 {{0.0716169, 0.0752675}},
0108 {{-0.332171, 0.0425903}},
0109 {{11.5791, 1.44028}},
0110 {{0.0511006, 0.104321}},
0111 {{-0.050789, 0.106859}}}};
0112
0113 constexpr std::array<std::array<float, nBinsBrem>, nBinsEta> par2 = {{{{-0.00368104, 4.70884}},
0114 {{5.23856, 3.35623}},
0115 {{-31.8456, 6.52561}},
0116 {{-112.084, -40.}},
0117 {{7.56304, 5.71476}},
0118 {{-7.96854, 3.54235}}}};
0119
0120 constexpr std::array<std::array<float, nBinsBrem>, nBinsEta> par3 = {{{{0.219829, 9.07419e-08}},
0121 {{0.00632907, 2.49397e-07}},
0122 {{22.543, 2.18593e-08}},
0123 {{-863.968, 0.00102639}},
0124 {{0.00331583, 6.12472e-06}},
0125 {{4.71223, 6.89631e-06}}}};
0126
0127 int iEtaSl = -1;
0128 for (int iEta = 0; iEta < nBinsEta; ++iEta) {
0129 if (EtaBins[iEta] <= std::abs(eta) && std::abs(eta) < EtaBins[iEta + 1]) {
0130 iEtaSl = iEta;
0131 }
0132 }
0133
0134 int iBremSl = -1;
0135 for (int iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
0136 if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
0137 iBremSl = iBrem;
0138 }
0139 }
0140
0141 if (std::abs(eta) > 2.5)
0142 iEtaSl = nBinsEta - 1;
0143 if (brem < BremBins[0])
0144 iBremSl = 0;
0145 if (brem > BremBins[nBinsBrem - 1])
0146 iBremSl = nBinsBrem - 1;
0147
0148 float uncertainty = 0;
0149 if (iBremSl >= 0 && iBremSl < nBinsBrem && iEtaSl >= 0 && iEtaSl < nBinsEta) {
0150 if (et < 5)
0151 uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]) +
0152 par3[iEtaSl][iBremSl] / ((5 - par2[iEtaSl][iBremSl]) * (5 - par2[iEtaSl][iBremSl]));
0153 else if (et > 200)
0154 uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (200 - par2[iEtaSl][iBremSl]) +
0155 par3[iEtaSl][iBremSl] / ((200 - par2[iEtaSl][iBremSl]) * (200 - par2[iEtaSl][iBremSl]));
0156 else if (et >= 5 && et <= 200)
0157 uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]) +
0158 par3[iEtaSl][iBremSl] / ((et - par2[iEtaSl][iBremSl]) * (et - par2[iEtaSl][iBremSl]));
0159 }
0160
0161 return (uncertainty * energy);
0162 }
