Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:47 

0001 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronEnergyCorrector.h"
0002 #include "RecoEgamma/EgammaElectronAlgos/interface/EnergyUncertaintyElectronSpecific.h"
0003 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h"
0004 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionBaseClass.h"
0005 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 
0008 /****************************************************************************
0009  *
0010  * Classification based eta corrections for the ecal cluster energy
0011  *
0012  * \author Federico Ferri - INFN Milano, Bicocca university
0013  * \author Ivica Puljak - FESB, Split
0014  * \author Stephanie Baffioni - Laboratoire Leprince-Ringuet - École polytechnique, CNRS/IN2P3
0015  *
0016  ****************************************************************************/
0017 
0018 namespace {
0019 
0020   inline float energyError(float E, float* par) {
0021     return sqrt(pow(par[0] / sqrt(E), 2) + pow(par[1] / E, 2) + pow(par[2], 2));
0022   }
0023 
0024   // main correction function
0025   // new corrections: taken from EcalClusterCorrectionObjectSpecific.cc (N. Chanon et al.)
0026   // this is to prepare for class based corrections, for the time being the parameters are the same as for the SC corrections
0027   // code fully duplicated here, to be improved; electron means algorithm==0 and mode==0
0028 
0029   float fEta(float energy, float eta, int algorithm) {
0030     // corrections for electrons
0031     if (algorithm != 0) {
0032       edm::LogWarning("ElectronEnergyCorrector::fEta") << "algorithm should be 0 for electrons !";
0033       return energy;
0034     }
0035     //std::cout << "fEta function" << std::endl;
0036 
0037     // this correction is setup only for EB
0038     if (algorithm != 0)
0039       return energy;
0040 
0041     float ieta = fabs(eta) * (5 / 0.087);
0042     // bypass the DB reading for the time being
0043     //float p0 = (params_->params())[0];  // should be 40.2198
0044     //float p1 = (params_->params())[1];  // should be -3.03103e-6
0045     float p0 = 40.2198;
0046     float p1 = -3.03103e-6;
0047 
0048     //std::cout << "ieta=" << ieta << std::endl;
0049 
0050     float correctedEnergy = energy;
0051     if (ieta < p0)
0052       correctedEnergy = energy;
0053     else
0054       correctedEnergy = energy / (1.0 + p1 * (ieta - p0) * (ieta - p0));
0055     //std::cout << "ECEC fEta = " << correctedEnergy << std::endl;
0056     return correctedEnergy;
0057   }
0058 
0059   float fBremEta(float sigmaPhiSigmaEta, float eta, int algorithm, reco::GsfElectron::Classification cl) {
0060     // corrections for electrons
0061     if (algorithm != 0) {
0062       edm::LogWarning("ElectronEnergyCorrector::fBremEta") << "algorithm should be 0 for electrons !";
0063       return 1.;
0064     }
0065 
0066     const float etaCrackMin = 1.44;
0067     const float etaCrackMax = 1.56;
0068 
0069     //STD
0070     const int nBinsEta = 14;
0071     float leftEta[nBinsEta] = {0.02, 0.25, 0.46, 0.81, 0.91, 1.01, 1.16, etaCrackMax, 1.653, 1.8, 2.0, 2.2, 2.3, 2.4};
0072     float rightEta[nBinsEta] = {0.25, 0.42, 0.77, 0.91, 1.01, 1.13, etaCrackMin, 1.653, 1.8, 2.0, 2.2, 2.3, 2.4, 2.5};
0073 
0074     // eta = 0
0075     if (std::abs(eta) < leftEta[0]) {
0076       eta = 0.02;
0077     }
0078 
0079     // outside acceptance
0080     if (std::abs(eta) >= rightEta[nBinsEta - 1]) {
0081       eta = 2.49;
0082     }  //if (DBG) std::cout << " WARNING [applyScCorrections]: std::abs(eta)  >=  rightEta[nBinsEta-1] " << std::endl;}
0083 
0084     int tmpEta = -1;
0085     for (int iEta = 0; iEta < nBinsEta; ++iEta) {
0086       if (leftEta[iEta] <= std::abs(eta) && std::abs(eta) < rightEta[iEta]) {
0087         tmpEta = iEta;
0088       }
0089     }
0090 
0091     float xcorr[nBinsEta][reco::GsfElectron::GAP + 1] = {{1.00227, 1.00227, 1.00227, 1.00227, 1.00227},
0092                                                          {1.00252, 1.00252, 1.00252, 1.00252, 1.00252},
0093                                                          {1.00225, 1.00225, 1.00225, 1.00225, 1.00225},
0094                                                          {1.00159, 1.00159, 1.00159, 1.00159, 1.00159},
0095                                                          {0.999475, 0.999475, 0.999475, 0.999475, 0.999475},
0096                                                          {0.997203, 0.997203, 0.997203, 0.997203, 0.997203},
0097                                                          {0.993886, 0.993886, 0.993886, 0.993886, 0.993886},
0098                                                          {0.971262, 0.971262, 0.971262, 0.971262, 0.971262},
0099                                                          {0.975922, 0.975922, 0.975922, 0.975922, 0.975922},
0100                                                          {0.979087, 0.979087, 0.979087, 0.979087, 0.979087},
0101                                                          {0.98495, 0.98495, 0.98495, 0.98495, 0.98495},
0102                                                          {0.98781, 0.98781, 0.98781, 0.98781, 0.98781},
0103                                                          {0.989546, 0.989546, 0.989546, 0.989546, 0.989546},
0104                                                          {0.989638, 0.989638, 0.989638, 0.989638, 0.989638}};
0105 
0106     float par0[nBinsEta][reco::GsfElectron::GAP + 1] = {{0.994949, 1.00718, 1.00718, 1.00556, 1.00718},
0107                                                         {1.009, 1.00713, 1.00713, 1.00248, 1.00713},
0108                                                         {0.999395, 1.00641, 1.00641, 1.00293, 1.00641},
0109                                                         {0.988662, 1.00761, 1.001, 0.99972, 1.00761},
0110                                                         {0.998443, 1.00682, 1.001, 1.00282, 1.00682},
0111                                                         {1.00285, 1.0073, 1.001, 1.00396, 1.0073},
0112                                                         {0.993053, 1.00462, 1.01341, 1.00184, 1.00462},
0113                                                         {1.10561, 0.972798, 1.02835, 0.995218, 0.972798},
0114                                                         {0.893741, 0.981672, 0.98982, 1.01712, 0.981672},
0115                                                         {0.911123, 0.98251, 1.03466, 1.00824, 0.98251},
0116                                                         {0.981931, 0.986123, 0.954295, 1.0202, 0.986123},
0117                                                         {0.905634, 0.990124, 0.928934, 0.998492, 0.990124},
0118                                                         {0.919343, 0.990187, 0.967526, 0.963923, 0.990187},
0119                                                         {0.844783, 0.99372, 0.923808, 0.953001, 0.99372}};
0120 
0121     float par1[nBinsEta][reco::GsfElectron::GAP + 1] = {
0122         {0.0111034, -0.00187886, -0.00187886, -0.00289304, -0.00187886},
0123         {-0.00969012, -0.00227574, -0.00227574, -0.00182187, -0.00227574},
0124         {0.00389454, -0.00259935, -0.00259935, -0.00211059, -0.00259935},
0125         {0.017095, -0.00433692, -0.00302335, -0.00241385, -0.00433692},
0126         {-0.00049009, -0.00551324, -0.00302335, -0.00532352, -0.00551324},
0127         {-0.00252723, -0.00799669, -0.00302335, -0.00823109, -0.00799669},
0128         {0.00332567, -0.00870057, -0.0170581, -0.011482, -0.00870057},
0129         {-0.213285, -0.000771577, -0.036007, -0.0187526, -0.000771577},
0130         {0.1741, -0.00202028, -0.0233995, -0.0302066, -0.00202028},
0131         {0.152794, 0.00441308, -0.0468563, -0.0158817, 0.00441308},
0132         {0.0351465, 0.00832913, 0.0358028, -0.0233262, 0.00832913},
0133         {0.185781, 0.00742879, 0.08858, 0.00568078, 0.00742879},
0134         {0.153088, 0.0094608, 0.0489979, 0.0491897, 0.0094608},
0135         {0.296681, 0.00560406, 0.106492, 0.0652007, 0.00560406}};
0136 
0137     float par2[nBinsEta][reco::GsfElectron::GAP + 1] = {{-0.00330844, 0, 0, 5.62441e-05, 0},
0138                                                         {0.00329373, 0, 0, -0.000113883, 0},
0139                                                         {-0.00104661, 0, 0, -0.000152794, 0},
0140                                                         {-0.0060409, 0, -0.000257724, -0.000202099, 0},
0141                                                         {-0.000742866, 0, -0.000257724, -2.06003e-05, 0},
0142                                                         {-0.00205425, 0, -0.000257724, 3.84179e-05, 0},
0143                                                         {-0.00350757, 0, 0.000590483, 0.000323723, 0},
0144                                                         {0.0794596, -0.00276696, 0.00205854, 0.000356716, -0.00276696},
0145                                                         {-0.092436, -0.00471028, 0.00062096, 0.00088347, -0.00471028},
0146                                                         {-0.0855029, -0.00809139, 0.00284102, -0.00366903, -0.00809139},
0147                                                         {-0.0306209, -0.00944584, -0.0145892, -0.00176969, -0.00944584},
0148                                                         {-0.0996414, -0.00960462, -0.0328264, -0.00983844, -0.00960462},
0149                                                         {-0.0784107, -0.010172, -0.0256722, -0.0215133, -0.010172},
0150                                                         {-0.145815, -0.00943169, -0.0414525, -0.027087, -0.00943169}};
0151 
0152     float sigmaPhiSigmaEtaMin[reco::GsfElectron::GAP + 1] = {0.8, 0.8, 0.8, 0.8, 0.8};
0153     float sigmaPhiSigmaEtaMax[reco::GsfElectron::GAP + 1] = {5., 5., 5., 5., 5.};
0154     float sigmaPhiSigmaEtaFit[reco::GsfElectron::GAP + 1] = {1.2, 1.2, 1.2, 1.2, 1.2};
0155 
0156     // extra protections
0157     // fix sigmaPhiSigmaEta boundaries
0158     if (sigmaPhiSigmaEta < sigmaPhiSigmaEtaMin[cl]) {
0159       sigmaPhiSigmaEta = sigmaPhiSigmaEtaMin[cl];
0160     }
0161     if (sigmaPhiSigmaEta > sigmaPhiSigmaEtaMax[cl]) {
0162       sigmaPhiSigmaEta = sigmaPhiSigmaEtaMax[cl];
0163     }
0164 
0165     // In eta cracks/gaps
0166     if (tmpEta == -1)  // need to interpolate
0167     {
0168       float tmpInter = 1;
0169       for (int iEta = 0; iEta < nBinsEta - 1; ++iEta) {
0170         if (rightEta[iEta] <= std::abs(eta) && std::abs(eta) < leftEta[iEta + 1]) {
0171           if (sigmaPhiSigmaEta >= sigmaPhiSigmaEtaFit[cl]) {
0172             tmpInter =
0173                 (par0[iEta][cl] + sigmaPhiSigmaEta * par1[iEta][cl] +
0174                  sigmaPhiSigmaEta * sigmaPhiSigmaEta * par2[iEta][cl] + par0[iEta + 1][cl] +
0175                  sigmaPhiSigmaEta * par1[iEta + 1][cl] + sigmaPhiSigmaEta * sigmaPhiSigmaEta * par2[iEta + 1][cl]) /
0176                 2.;
0177           } else
0178             tmpInter = (xcorr[iEta][cl] + xcorr[iEta + 1][cl]) / 2.;
0179         }
0180       }
0181       return tmpInter;
0182     }
0183 
0184     if (sigmaPhiSigmaEta >= sigmaPhiSigmaEtaFit[cl]) {
0185       return par0[tmpEta][cl] + sigmaPhiSigmaEta * par1[tmpEta][cl] +
0186              sigmaPhiSigmaEta * sigmaPhiSigmaEta * par2[tmpEta][cl];
0187     } else {
0188       return xcorr[tmpEta][cl];
0189     }
0190 
0191     return 1.;
0192   }
0193 
0194   float fEt(float ET, int algorithm, reco::GsfElectron::Classification cl) {
0195     if (algorithm == 0)  //Electrons EB
0196     {
0197       const float parClassIndep[5] = {0.97213, 0.999528, 5.61192e-06, 0.0143269, -17.1776};
0198       const float par[reco::GsfElectron::GAP + 1][5] = {
0199           {0.974327, 0.996127, 5.99401e-05, 0.159813, -3.80392},
0200           {0.97213, 0.999528, 5.61192e-06, 0.0143269, -17.1776},
0201           {0.940666, 0.988894, 0.00017474, 0.25603, -4.58153},
0202           {0.969526, 0.98572, 0.000193842, 4.21548, -1.37159},
0203           {parClassIndep[0], parClassIndep[1], parClassIndep[2], parClassIndep[3], parClassIndep[4]}};
0204       if (ET > 200) {
0205         ET = 200;
0206       }
0207       if (ET > 100) {
0208         return (parClassIndep[1] + ET * parClassIndep[2]) * (1 - parClassIndep[3] * exp(ET / parClassIndep[4]));
0209       }
0210       if (ET >= 10) {
0211         return (par[cl][1] + ET * par[cl][2]) * (1 - par[cl][3] * exp(ET / par[cl][4]));
0212       }
0213       if (ET >= 5) {
0214         return par[cl][0];
0215       }
0216       return 1.;
0217     } else if (algorithm == 1)  //Electrons EE
0218     {
0219       float par[reco::GsfElectron::GAP + 1][5] = {{0.930081, 0.996683, 3.54079e-05, 0.0460187, -23.2461},
0220                                                   {0.930081, 0.996683, 3.54079e-05, 0.0460187, -23.2461},
0221                                                   {0.930081, 0.996683, 3.54079e-05, 0.0460187, -23.2461},
0222                                                   {0.930081, 0.996683, 3.54079e-05, 0.0460187, -23.2461},
0223                                                   {0.930081, 0.996683, 3.54079e-05, 0.0460187, -23.2461}};
0224       if (ET > 200) {
0225         ET = 200;
0226       }
0227       if (ET < 5) {
0228         return 1.;
0229       }
0230       if (5 <= ET && ET < 10) {
0231         return par[cl][0];
0232       }
0233       if (10 <= ET && ET <= 200) {
0234         return (par[cl][1] + ET * par[cl][2]) * (1 - par[cl][3] * exp(ET / par[cl][4]));
0235       }
0236     } else {
0237       edm::LogWarning("ElectronEnergyCorrector::fEt") << "algorithm should be 0 or 1 for electrons !";
0238     }
0239     return 1.;
0240   }
0241 
0242   float fEnergy(float E, int algorithm, reco::GsfElectron::Classification cl) {
0243     if (algorithm == 0)  // Electrons EB
0244     {
0245       return 1.;
0246     } else if (algorithm == 1)  // Electrons EE
0247     {
0248       float par0[reco::GsfElectron::GAP + 1] = {400, 400, 400, 400, 400};
0249       float par1[reco::GsfElectron::GAP + 1] = {0.999545, 0.982475, 0.986217, 0.996763, 0.982475};
0250       float par2[reco::GsfElectron::GAP + 1] = {1.26568e-05, 4.95413e-05, 5.02161e-05, 2.8485e-06, 4.95413e-05};
0251       float par3[reco::GsfElectron::GAP + 1] = {0.0696757, 0.16886, 0.115317, 0.12098, 0.16886};
0252       float par4[reco::GsfElectron::GAP + 1] = {-54.3468, -30.1517, -26.3805, -62.0538, -30.1517};
0253 
0254       if (E > par0[cl]) {
0255         E = par0[cl];
0256       }
0257       if (E < 0) {
0258         return 1.;
0259       }
0260       if (0 <= E && E <= par0[cl]) {
0261         return (par1[cl] + E * par2[cl]) * (1 - par3[cl] * exp(E / par4[cl]));
0262       }
0263     } else {
0264       edm::LogWarning("ElectronEnergyCorrector::fEnergy") << "algorithm should be 0 or 1 for electrons !";
0265     }
0266     return 1.;
0267   }
0268 
0269 }  // namespace
0270 
0271 double egamma::classBasedElectronEnergyUncertainty(reco::GsfElectron const& electron) {
0272   double ecalEnergy = electron.correctedEcalEnergy();
0273   double eleEta = electron.superCluster()->eta();
0274   double brem = electron.superCluster()->etaWidth() / electron.superCluster()->phiWidth();
0275   return egamma::electronEnergyUncertainty(electron.classification(), eleEta, brem, ecalEnergy);
0276 }
0277 
0278 double egamma::simpleElectronEnergyUncertainty(reco::GsfElectron const& electron) {
0279   double error = 999.;
0280   double ecalEnergy = electron.correctedEcalEnergy();
0281 
0282   if (electron.isEB()) {
0283     float parEB[3] = {5.24e-02, 2.01e-01, 1.00e-02};
0284     error = ecalEnergy * energyError(ecalEnergy, parEB);
0285   } else if (electron.isEE()) {
0286     float parEE[3] = {1.46e-01, 9.21e-01, 1.94e-03};
0287     error = ecalEnergy * energyError(ecalEnergy, parEE);
0288   } else {
0289     edm::LogWarning("ElectronEnergyCorrector::simpleParameterizationUncertainty")
0290         << "nor barrel neither endcap electron !";
0291   }
0292 
0293   return error;
0294 }
0295 
0296 float egamma::classBasedElectronEnergy(reco::GsfElectron const& electron,
0297                                        reco::BeamSpot const& bs,
0298                                        EcalClusterFunctionBaseClass const& crackCorrectionFunction) {
0299   auto elClass = electron.classification();
0300 
0301   // new corrections from N. Chanon et al., taken from EcalClusterCorrectionObjectSpecific.cc
0302   float corr = 1.;
0303   float corr2 = 1.;
0304   float energy = electron.superCluster()->energy();
0305 
0306   //int subdet = electron.superCluster()->seed()->hitsAndFractions()[0].first.subdetId();
0307 
0308   if (electron.isEB()) {
0309     float cetacorr = fEta(electron.superCluster()->rawEnergy(), electron.superCluster()->eta(), 0) /
0310                      electron.superCluster()->rawEnergy();
0311     energy = electron.superCluster()->rawEnergy() * cetacorr;  //previously in CMSSW
0312     //energy = superCluster.rawEnergy()*fEta(e5x5, superCluster.seed()->eta(), 0)/e5x5;
0313   } else if (electron.isEE()) {
0314     energy = electron.superCluster()->rawEnergy() + electron.superCluster()->preshowerEnergy();
0315   } else {
0316     edm::LogWarning("ElectronEnergyCorrector::classBasedParameterizationEnergy")
0317         << "nor barrel neither endcap electron !";
0318   }
0319 
0320   corr = fBremEta(electron.superCluster()->phiWidth() / electron.superCluster()->etaWidth(),
0321                   electron.superCluster()->eta(),
0322                   0,
0323                   elClass);
0324 
0325   float et = energy * std::sin(2 * std::atan(std::exp(-electron.superCluster()->eta()))) / corr;
0326 
0327   if (electron.isEB()) {
0328     corr2 = corr * fEt(et, 0, elClass);
0329   } else if (electron.isEE()) {
0330     corr2 = corr * fEnergy(energy / corr, 1, elClass);
0331   } else {
0332     edm::LogWarning("ElectronEnergyCorrector::classBasedParameterizationEnergy")
0333         << "nor barrel neither endcap electron !";
0334   }
0335 
0336   float newEnergy = energy / corr2;
0337 
0338   // cracks
0339   double crackcor = 1.;
0340   for (reco::CaloCluster_iterator cIt = electron.superCluster()->clustersBegin();
0341        cIt != electron.superCluster()->clustersEnd();
0342        ++cIt) {
0343     const reco::CaloClusterPtr cc = *cIt;
0344     crackcor *= (electron.superCluster()->rawEnergy() + cc->energy() * (crackCorrectionFunction.getValue(*cc) - 1.)) /
0345                 electron.superCluster()->rawEnergy();
0346   }
0347   newEnergy *= crackcor;
0348 
0349   return newEnergy;
0350 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!