EgammaCoreTools/plugins/EcalClusterEnergyCorrection.cc

0001 /** \class EcalClusterEnergyCorrection
0002   *  Function that provides supercluster energy correction due to Bremsstrahlung loss
0003   *
0004   *  $Id: EcalClusterEnergyCorrection.h
0005   *  $Date:
0006   *  $Revision:
0007   *  \author Yurii Maravin, KSU, March 2009
0008   */
0009
0010 #include "CondFormats/DataRecord/interface/EcalClusterEnergyCorrectionParametersRcd.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/ConsumesCollector.h"
0014 #include "CondFormats/EcalObjects/interface/EcalClusterEnergyCorrectionParameters.h"
0015 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionBaseClass.h"
0016
0017 class EcalClusterEnergyCorrection : public EcalClusterFunctionBaseClass {
0018 public:
0019   EcalClusterEnergyCorrection(const edm::ParameterSet &, edm::ConsumesCollector iC) : paramsToken_(iC.esConsumes()) {}
0020
0021   // get/set explicit methods for parameters
0022   const EcalClusterEnergyCorrectionParameters *getParameters() const { return params_; }
0023   // check initialization
0024   void checkInit() const;
0025
0026   // compute the correction
0027   float getValue(const reco::SuperCluster &, const int mode) const override;
0028   float getValue(const reco::BasicCluster &, const EcalRecHitCollection &) const override { return 0.; };
0029
0030   // set parameters
0031   void init(const edm::EventSetup &es) override;
0032
0033 private:
0034   float fEta(float e, float eta, int algorithm) const;
0035   float fBrem(float e, float eta, int algorithm) const;
0036   float fEtEta(float et, float eta, int algorithm) const;
0037
0038   const edm::ESGetToken<EcalClusterEnergyCorrectionParameters, EcalClusterEnergyCorrectionParametersRcd> paramsToken_;
0039   const EcalClusterEnergyCorrectionParameters *params_ = nullptr;
0040 };
0041
0042 // Shower leakage corrections developed by Jungzhie et al. using TB data
0043 // Developed for EB only!
0044 float EcalClusterEnergyCorrection::fEta(float energy, float eta, int algorithm) const {
0045   // this correction is setup only for EB
0046   if (algorithm != 0)
0047     return energy;
0048
0049   float ieta = fabs(eta) * (5 / 0.087);
0050   float p0 = (params_->params())[0];  // should be 40.2198
0051   float p1 = (params_->params())[1];  // should be -3.03103e-6
0052
0053   float correctedEnergy = energy;
0054   if (ieta < p0)
0055     correctedEnergy = energy;
0056   else
0057     correctedEnergy = energy / (1.0 + p1 * (ieta - p0) * (ieta - p0));
0058   //std::cout << "ECEC fEta = " << correctedEnergy << std::endl;
0059   return correctedEnergy;
0060 }
0061
0062 float EcalClusterEnergyCorrection::fBrem(float e, float brem, int algorithm) const {
0063   // brem == phiWidth/etaWidth of the SuperCluster
0064   // e  == energy of the SuperCluster
0065   // first parabola (for br > threshold)
0066   // p0 + p1*x + p2*x^2
0067   // second parabola (for br <= threshold)
0068   // ax^2 + bx + c, make y and y' the same in threshold
0069   // y = p0 + p1*threshold + p2*threshold^2
0070   // yprime = p1 + 2*p2*threshold
0071   // a = p3
0072   // b = yprime - 2*a*threshold
0073   // c = y - a*threshold^2 - b*threshold
0074
0075   int offset;
0076   if (algorithm == 0)
0077     offset = 0;
0078   else if (algorithm == 1)
0079     offset = 20;
0080   else {
0081     // not supported, produce no correction
0082     return e;
0083   }
0084
0085   //Make No Corrections if brem is invalid!
0086   if (brem == 0)
0087     return e;
0088
0089   float bremLowThr = (params_->params())[2 + offset];
0090   float bremHighThr = (params_->params())[3 + offset];
0091   if (brem < bremLowThr)
0092     brem = bremLowThr;
0093   if (brem > bremHighThr)
0094     brem = bremHighThr;
0095
0096   // Parameters provided in cfg file
0097   float p0 = (params_->params())[4 + offset];
0098   float p1 = (params_->params())[5 + offset];
0099   float p2 = (params_->params())[6 + offset];
0100   float p3 = (params_->params())[7 + offset];
0101   float p4 = (params_->params())[8 + offset];
0102   //
0103   float threshold = p4;
0104
0105   float y = p0 * threshold * threshold + p1 * threshold + p2;
0106   float yprime = 2 * p0 * threshold + p1;
0107   float a = p3;
0108   float b = yprime - 2 * a * threshold;
0109   float c = y - a * threshold * threshold - b * threshold;
0110
0111   float fCorr = 1;
0112   if (brem < threshold)
0113     fCorr = p0 * brem * brem + p1 * brem + p2;
0114   else
0115     fCorr = a * brem * brem + b * brem + c;
0116
0117   //std::cout << "ECEC fBrem " << e/fCorr << std::endl;
0118   return e / fCorr;
0119 }
0120
0121 float EcalClusterEnergyCorrection::fEtEta(float et, float eta, int algorithm) const {
0122   // et -- Et of the SuperCluster (with respect to (0,0,0))
0123   // eta -- eta of the SuperCluster
0124
0125   //std::cout << "fEtEta, mode = " << algorithm << std::endl;
0126   //std::cout << "ECEC: p0    " << (params_->params())[9]  << " " << (params_->params())[10] << " " << (params_->params())[11] << " " << (params_->params())[12] << std::endl;
0127   //std::cout << "ECEC: p1    " << (params_->params())[13] << " " << (params_->params())[14] << " " << (params_->params())[15] << " " << (params_->params())[16] << std::endl;
0128   //std::cout << "ECEC: fcorr " << (params_->params())[17] << " " << (params_->params())[18] << " " << (params_->params())[19] << std::endl;
0129
0130   float fCorr = 0.;
0131   int offset;
0132   if (algorithm == 0)
0133     offset = 0;
0134   else if (algorithm == 1)
0135     offset = 20;
0136   else if (algorithm == 10)
0137     offset = 28;
0138   else if (algorithm == 11)
0139     offset = 39;
0140   else {
0141     // not supported, produce no correction
0142     return et;
0143   }
0144
0145   // Barrel
0146   if (algorithm == 0 || algorithm == 10) {
0147     float p0 = (params_->params())[9 + offset] +
0148                (params_->params())[10 + offset] / (et + (params_->params())[11 + offset]) +
0149                (params_->params())[12 + offset] / (et * et);
0150     float p1 = (params_->params())[13 + offset] +
0151                (params_->params())[14 + offset] / (et + (params_->params())[15 + offset]) +
0152                (params_->params())[16 + offset] / (et * et);
0153
0154     fCorr = p0 + p1 * atan((params_->params())[17 + offset] * ((params_->params())[18 + offset] - fabs(eta))) +
0155             (params_->params())[19 + offset] * fabs(eta);
0156
0157   } else if (algorithm == 1 || algorithm == 11) {  // Endcap
0158     float p0 = (params_->params())[9 + offset] + (params_->params())[10 + offset] / sqrt(et);
0159     float p1 = (params_->params())[11 + offset] + (params_->params())[12 + offset] / sqrt(et);
0160     float p2 = (params_->params())[13 + offset] + (params_->params())[14 + offset] / sqrt(et);
0161     float p3 = (params_->params())[15 + offset] + (params_->params())[16 + offset] / sqrt(et);
0162
0163     fCorr = p0 + p1 * fabs(eta) + p2 * eta * eta + p3 / fabs(eta);
0164   }
0165
0166   // cap the correction at 50%
0167   if (fCorr < 0.5)
0168     fCorr = 0.5;
0169   if (fCorr > 1.5)
0170     fCorr = 1.5;
0171
0172   //std::cout << "ECEC fEtEta " << et/fCorr << std::endl;
0173   return et / fCorr;
0174 }
0175
0176 float EcalClusterEnergyCorrection::getValue(const reco::SuperCluster &superCluster, const int mode) const {
0177   // mode flags:
0178   // hybrid modes: 1 - return f(eta) correction in GeV
0179   //               2 - return f(eta) + f(brem) correction
0180   //               3 - return f(eta) + f(brem) + f(et, eta) correction
0181   // multi5x5:     4 - return f(brem) correction
0182   //               5 - return f(brem) + f(et, eta) correction
0183
0184   // special cases: mode = 10 -- return f(et, eta) correction with respect to already corrected SC in barrel
0185   //                mode = 11 -- return f(et, eta) correction with respect to already corrected SC in endcap
0186
0187   checkInit();
0188
0189   float eta = fabs(superCluster.eta());
0190   float brem = superCluster.phiWidth() / superCluster.etaWidth();
0191   int algorithm = -1;
0192
0193   if (mode <= 3 || mode == 10) {
0194     // algorithm: hybrid
0195     algorithm = 0;
0196
0197     float energy = superCluster.rawEnergy();
0198     if (mode == 10) {
0199       algorithm = 10;
0200       energy = superCluster.energy();
0201     }
0202     float correctedEnergy = fEta(energy, eta, algorithm);
0203
0204     if (mode == 1) {
0205       return correctedEnergy - energy;
0206
0207     } else {
0208       // now apply F(brem)
0209       correctedEnergy = fBrem(correctedEnergy, brem, algorithm);
0210       if (mode == 2) {
0211         return correctedEnergy - energy;
0212       }
0213
0214       float correctedEt = correctedEnergy / cosh(eta);
0215       correctedEt = fEtEta(correctedEt, eta, algorithm);
0216       correctedEnergy = correctedEt * cosh(eta);
0217       return correctedEnergy - energy;
0218     }
0219   } else if (mode == 4 || mode == 5 || mode == 11) {
0220     algorithm = 1;
0221
0222     float energy = superCluster.rawEnergy() + superCluster.preshowerEnergy();
0223     if (mode == 11) {
0224       algorithm = 11;
0225       energy = superCluster.energy();
0226     }
0227
0228     float correctedEnergy = fBrem(energy, brem, algorithm);
0229     if (mode == 4) {
0230       return correctedEnergy - energy;
0231     }
0232
0233     float correctedEt = correctedEnergy / cosh(eta);
0234     correctedEt = fEtEta(correctedEt, eta, algorithm);
0235     correctedEnergy = correctedEt * cosh(eta);
0236     return correctedEnergy - energy;
0237
0238   } else {
0239     // perform no correction
0240     return 0;
0241   }
0242 }
0243
0244 void EcalClusterEnergyCorrection::init(const edm::EventSetup &es) { params_ = &es.getData(paramsToken_); }
0245
0246 void EcalClusterEnergyCorrection::checkInit() const {
0247   if (!params_) {
0248     // non-initialized function parameters: throw exception
0249     throw cms::Exception("EcalClusterEnergyCorrection::checkInit()")
0250         << "Trying to access an uninitialized crack correction function.\n"
0251            "Please call `init( edm::EventSetup &)' before any use of the function.\n";
0252   }
0253 }
0254
0255 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionFactory.h"
0256 DEFINE_EDM_PLUGIN(EcalClusterFunctionFactory, EcalClusterEnergyCorrection, "EcalClusterEnergyCorrection");