File indexing completed on 2023-10-25 09:59:03
0001
0002
0003
0004
0005
0006
0007
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
0022 const EcalClusterEnergyCorrectionParameters *getParameters() const { return params_; }
0023
0024 void checkInit() const;
0025
0026
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
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
0043
0044 float EcalClusterEnergyCorrection::fEta(float energy, float eta, int algorithm) const {
0045
0046 if (algorithm != 0)
0047 return energy;
0048
0049 float ieta = fabs(eta) * (5 / 0.087);
0050 float p0 = (params_->params())[0];
0051 float p1 = (params_->params())[1];
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
0059 return correctedEnergy;
0060 }
0061
0062 float EcalClusterEnergyCorrection::fBrem(float e, float brem, int algorithm) const {
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075 int offset;
0076 if (algorithm == 0)
0077 offset = 0;
0078 else if (algorithm == 1)
0079 offset = 20;
0080 else {
0081
0082 return e;
0083 }
0084
0085
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
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
0118 return e / fCorr;
0119 }
0120
0121 float EcalClusterEnergyCorrection::fEtEta(float et, float eta, int algorithm) const {
0122
0123
0124
0125
0126
0127
0128
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
0142 return et;
0143 }
0144
0145
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) {
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
0167 if (fCorr < 0.5)
0168 fCorr = 0.5;
0169 if (fCorr > 1.5)
0170 fCorr = 1.5;
0171
0172
0173 return et / fCorr;
0174 }
0175
0176 float EcalClusterEnergyCorrection::getValue(const reco::SuperCluster &superCluster, const int mode) const {
0177
0178
0179
0180
0181
0182
0183
0184
0185
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
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
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
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
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");