Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-30 02:37:09

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 }