File indexing completed on 2024-04-06 12:27:26
0001
0002 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
0003 #include "CondFormats/PhysicsToolsObjects/interface/PerformancePayloadFromTFormula.h"
0004
0005 #include "CondFormats/ESObjects/interface/ESEEIntercalibConstants.h"
0006
0007 #include <TMath.h>
0008 #include <cmath>
0009 #include <vector>
0010 #include <map>
0011 #include <algorithm>
0012 #include <numeric>
0013
0014 using namespace std;
0015
0016 PFEnergyCalibration::PFEnergyCalibration() {
0017
0018 faBarrel = std::make_unique<TF1>(
0019 "faBarrel", "[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))", 0., 1000.);
0020 faBarrel->SetParameter(0, -30.7141);
0021 faBarrel->SetParameter(1, 31.7583);
0022 faBarrel->SetParameter(2, 4.40594);
0023 faBarrel->SetParameter(3, 1.70914);
0024 faBarrel->SetParameter(4, 0.0613696);
0025 faBarrel->SetParameter(5, 0.000104857);
0026 faBarrel->SetParameter(6, -1.38927);
0027 faBarrel->SetParameter(7, -0.743082);
0028 fbBarrel = std::make_unique<TF1>(
0029 "fbBarrel", "[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))", 0., 1000.);
0030 fbBarrel->SetParameter(0, 2.25366);
0031 fbBarrel->SetParameter(1, 0.537715);
0032 fbBarrel->SetParameter(2, -4.81375);
0033 fbBarrel->SetParameter(3, 12.109);
0034 fbBarrel->SetParameter(4, 1.80577);
0035 fbBarrel->SetParameter(5, 0.187919);
0036 fbBarrel->SetParameter(6, -6.26234);
0037 fbBarrel->SetParameter(7, -0.607392);
0038 fcBarrel = std::make_unique<TF1>(
0039 "fcBarrel", "[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))", 0., 1000.);
0040 fcBarrel->SetParameter(0, 1.5125962);
0041 fcBarrel->SetParameter(1, 0.855057);
0042 fcBarrel->SetParameter(2, -6.04199);
0043 fcBarrel->SetParameter(3, 2.08229);
0044 fcBarrel->SetParameter(4, 0.592266);
0045 fcBarrel->SetParameter(5, 0.0291232);
0046 fcBarrel->SetParameter(6, 0.364802);
0047 fcBarrel->SetParameter(7, -1.50142);
0048 faEtaBarrelEH = std::make_unique<TF1>("faEtaBarrelEH", "[0]+[1]*exp(-x/[2])", 0., 1000.);
0049 faEtaBarrelEH->SetParameter(0, 0.0185555);
0050 faEtaBarrelEH->SetParameter(1, -0.0470674);
0051 faEtaBarrelEH->SetParameter(2, 396.959);
0052 fbEtaBarrelEH = std::make_unique<TF1>("fbEtaBarrelEH", "[0]+[1]*exp(-x/[2])", 0., 1000.);
0053 fbEtaBarrelEH->SetParameter(0, 0.0396458);
0054 fbEtaBarrelEH->SetParameter(1, 0.114128);
0055 fbEtaBarrelEH->SetParameter(2, 251.405);
0056 faEtaBarrelH = std::make_unique<TF1>("faEtaBarrelH", "[0]+[1]*x", 0., 1000.);
0057 faEtaBarrelH->SetParameter(0, 0.00434994);
0058 faEtaBarrelH->SetParameter(1, -5.16564e-06);
0059 fbEtaBarrelH = std::make_unique<TF1>("fbEtaBarrelH", "[0]+[1]*exp(-x/[2])", 0., 1000.);
0060 fbEtaBarrelH->SetParameter(0, -0.0232604);
0061 fbEtaBarrelH->SetParameter(1, 0.0937525);
0062 fbEtaBarrelH->SetParameter(2, 34.9935);
0063
0064 faEndcap = std::make_unique<TF1>(
0065 "faEndcap", "[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))", 0., 1000.);
0066 faEndcap->SetParameter(0, 1.17227);
0067 faEndcap->SetParameter(1, 13.1489);
0068 faEndcap->SetParameter(2, -29.1672);
0069 faEndcap->SetParameter(3, 0.604223);
0070 faEndcap->SetParameter(4, 0.0426363);
0071 faEndcap->SetParameter(5, 3.30898e-15);
0072 faEndcap->SetParameter(6, 0.165293);
0073 faEndcap->SetParameter(7, -7.56786);
0074 fbEndcap = std::make_unique<TF1>(
0075 "fbEndcap", "[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))", 0., 1000.);
0076 fbEndcap->SetParameter(0, -0.974251);
0077 fbEndcap->SetParameter(1, 1.61733);
0078 fbEndcap->SetParameter(2, 0.0629183);
0079 fbEndcap->SetParameter(3, 7.78495);
0080 fbEndcap->SetParameter(4, -0.774289);
0081 fbEndcap->SetParameter(5, 7.81399e-05);
0082 fbEndcap->SetParameter(6, 0.139116);
0083 fbEndcap->SetParameter(7, -4.25551);
0084 fcEndcap = std::make_unique<TF1>(
0085 "fcEndcap", "[0]+((([1]+([2]/sqrt(x)))*exp(-(x^[6]/[3])))-([4]*exp(-(x^[7]/[5]))))", 0., 1000.);
0086 fcEndcap->SetParameter(0, 1.01863);
0087 fcEndcap->SetParameter(1, 1.29787);
0088 fcEndcap->SetParameter(2, -3.97293);
0089 fcEndcap->SetParameter(3, 21.7805);
0090 fcEndcap->SetParameter(4, 0.810195);
0091 fcEndcap->SetParameter(5, 0.234134);
0092 fcEndcap->SetParameter(6, 1.42226);
0093 fcEndcap->SetParameter(7, -0.0997326);
0094 faEtaEndcapEH = std::make_unique<TF1>("faEtaEndcapEH", "[0]+[1]*exp(-x/[2])", 0., 1000.);
0095 faEtaEndcapEH->SetParameter(0, 0.0112692);
0096 faEtaEndcapEH->SetParameter(1, -2.68063);
0097 faEtaEndcapEH->SetParameter(2, 2.90973);
0098 fbEtaEndcapEH = std::make_unique<TF1>("fbEtaEndcapEH", "[0]+[1]*exp(-x/[2])", 0., 1000.);
0099 fbEtaEndcapEH->SetParameter(0, -0.0192991);
0100 fbEtaEndcapEH->SetParameter(1, -0.265);
0101 fbEtaEndcapEH->SetParameter(2, 80.5502);
0102 faEtaEndcapH = std::make_unique<TF1>("faEtaEndcapH", "[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))", 0., 1000.);
0103 faEtaEndcapH->SetParameter(0, -0.0106029);
0104 faEtaEndcapH->SetParameter(1, -0.692207);
0105 faEtaEndcapH->SetParameter(2, 0.0542991);
0106 faEtaEndcapH->SetParameter(3, -0.171435);
0107 faEtaEndcapH->SetParameter(4, -61.2277);
0108 fbEtaEndcapH = std::make_unique<TF1>("fbEtaEndcapH", "[0]+[1]*exp(-x/[2])+[3]*[3]*exp(-x*x/([4]*[4]))", 0., 1000.);
0109 fbEtaEndcapH->SetParameter(0, 0.0214894);
0110 fbEtaEndcapH->SetParameter(1, -0.266704);
0111 fbEtaEndcapH->SetParameter(2, 5.2112);
0112 fbEtaEndcapH->SetParameter(3, 0.303578);
0113 fbEtaEndcapH->SetParameter(4, -104.367);
0114
0115
0116
0117 fcEtaBarrelH = std::make_unique<TF1>("fcEtaBarrelH", "[3]*((x-[0])^[1])+[2]", 0., 1000.);
0118 fcEtaBarrelH->SetParameter(0, 0);
0119 fcEtaBarrelH->SetParameter(1, 2);
0120 fcEtaBarrelH->SetParameter(2, 0);
0121 fcEtaBarrelH->SetParameter(3, 1);
0122
0123 fcEtaEndcapH = std::make_unique<TF1>("fcEtaEndcapH", "[3]*((x-[0])^[1])+[2]", 0., 1000.);
0124 fcEtaEndcapH->SetParameter(0, 0);
0125 fcEtaEndcapH->SetParameter(1, 0);
0126 fcEtaEndcapH->SetParameter(2, 0.05);
0127 fcEtaEndcapH->SetParameter(3, 0);
0128
0129 fdEtaEndcapH = std::make_unique<TF1>("fdEtaEndcapH", "[3]*((x-[0])^[1])+[2]", 0., 1000.);
0130 fdEtaEndcapH->SetParameter(0, 1.5);
0131 fdEtaEndcapH->SetParameter(1, 4);
0132 fdEtaEndcapH->SetParameter(2, -1.1);
0133 fdEtaEndcapH->SetParameter(3, 1.0);
0134
0135 fcEtaBarrelEH = std::make_unique<TF1>("fcEtaBarrelEH", "[3]*((x-[0])^[1])+[2]", 0., 1000.);
0136 fcEtaBarrelEH->SetParameter(0, 0);
0137 fcEtaBarrelEH->SetParameter(1, 2);
0138 fcEtaBarrelEH->SetParameter(2, 0);
0139 fcEtaBarrelEH->SetParameter(3, 1);
0140
0141 fcEtaEndcapEH = std::make_unique<TF1>("fcEtaEndcapEH", "[3]*((x-[0])^[1])+[2]", 0., 1000.);
0142 fcEtaEndcapEH->SetParameter(0, 0);
0143 fcEtaEndcapEH->SetParameter(1, 0);
0144 fcEtaEndcapEH->SetParameter(2, 0);
0145 fcEtaEndcapEH->SetParameter(3, 0);
0146
0147 fdEtaEndcapEH = std::make_unique<TF1>("fdEtaEndcapEH", "[3]*((x-[0])^[1])+[2]", 0., 1000.);
0148 fdEtaEndcapEH->SetParameter(0, 1.5);
0149 fdEtaEndcapEH->SetParameter(1, 2.0);
0150 fdEtaEndcapEH->SetParameter(2, 0.6);
0151 fdEtaEndcapEH->SetParameter(3, 1.0);
0152 }
0153
0154 PFEnergyCalibration::CalibratedEndcapPFClusterEnergies PFEnergyCalibration::calibrateEndcapClusterEnergies(
0155 reco::PFCluster const& eeCluster,
0156 std::vector<reco::PFCluster const*> const& psClusterPointers,
0157 ESChannelStatus const& channelStatus,
0158 bool applyCrackCorrections) const {
0159 double ps1_energy_sum = 0.;
0160 double ps2_energy_sum = 0.;
0161 bool condP1 = true;
0162 bool condP2 = true;
0163
0164 for (auto const& psclus : psClusterPointers) {
0165 bool cond = true;
0166 for (auto const& recH : psclus->recHitFractions()) {
0167 auto strip = recH.recHitRef()->detId();
0168 if (strip != ESDetId(0)) {
0169
0170
0171 if (channelStatus.getMap().find(strip)->getStatusCode() == 0) {
0172 cond = false;
0173 break;
0174 }
0175 }
0176 }
0177
0178 if (psclus->layer() == PFLayer::PS1) {
0179 ps1_energy_sum += psclus->energy();
0180 condP1 &= cond;
0181 } else if (psclus->layer() == PFLayer::PS2) {
0182 ps2_energy_sum += psclus->energy();
0183 condP2 &= cond;
0184 }
0185 }
0186
0187 double ePS1 = condP1 ? -1. : 0.;
0188 double ePS2 = condP2 ? -1. : 0.;
0189
0190 double cluscalibe = energyEm(eeCluster, ps1_energy_sum, ps2_energy_sum, ePS1, ePS2, applyCrackCorrections);
0191
0192 return {cluscalibe, ePS1, ePS2};
0193 }
0194
0195 void PFEnergyCalibration::energyEmHad(double t, double& e, double& h, double eta, double phi) const {
0196
0197 const double tt = t;
0198 const double ee = e;
0199 const double hh = h;
0200 double etaCorrE = 1.;
0201 double etaCorrH = 1.;
0202 auto absEta = std::abs(eta);
0203 t = min(999.9, max(tt, e + h));
0204 if (t < 1.)
0205 return;
0206
0207
0208 if (absEta < 1.48) {
0209
0210 double a = e > 0. ? aBarrel(t) : 1.;
0211 double b = e > 0. ? bBarrel(t) : cBarrel(t);
0212 double thresh = e > 0. ? threshE : threshH;
0213
0214
0215 if (a < -0.25 || b < -0.25) {
0216 a = 1.;
0217 b = 1.;
0218 thresh = 0.;
0219 }
0220
0221
0222 t = min(999.9, max(tt, thresh + a * e + b * h));
0223
0224
0225 if (e > 0. && thresh > 0.) {
0226 etaCorrE = 1.0 + aEtaBarrelEH(t) + 1.3 * bEtaBarrelEH(t) * cEtaBarrelEH(absEta);
0227 etaCorrH = 1.0;
0228 } else {
0229 etaCorrE = 1.0 + aEtaBarrelH(t) + 1.3 * bEtaBarrelH(t) * cEtaBarrelH(absEta);
0230 etaCorrH = 1.0 + aEtaBarrelH(t) + bEtaBarrelH(t) * cEtaBarrelH(absEta);
0231 }
0232 if (e > 0. && thresh > 0.)
0233 e = h > 0. ? threshE - threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
0234 if (h > 0. && thresh > 0.) {
0235 h = threshH + etaCorrH * b * h;
0236 }
0237
0238
0239 } else {
0240
0241 double a = e > 0. ? aEndcap(t) : 1.;
0242 double b = e > 0. ? bEndcap(t) : cEndcap(t);
0243 double thresh = e > 0. ? threshE : threshH;
0244
0245
0246 if (a < -0.25 || b < -0.25) {
0247 a = 1.;
0248 b = 1.;
0249 thresh = 0.;
0250 }
0251
0252
0253 t = min(999.9, max(tt, thresh + a * e + b * h));
0254
0255
0256 const double dEta = std::abs(absEta - 1.5);
0257 const double etaPow = dEta * dEta * dEta * dEta;
0258
0259 if (e > 0. && thresh > 0.) {
0260 if (absEta < 2.5) {
0261 etaCorrE = 1. + aEtaEndcapEH(t) + bEtaEndcapEH(t) * cEtaEndcapEH(absEta);
0262 } else {
0263 etaCorrE = 1. + aEtaEndcapEH(t) + 1.3 * bEtaEndcapEH(t) * dEtaEndcapEH(absEta);
0264 }
0265
0266 etaCorrH = 1. + aEtaEndcapEH(t) + bEtaEndcapEH(t) * (0.04 + etaPow);
0267 } else {
0268 etaCorrE = 1.;
0269 if (absEta < 2.5) {
0270 etaCorrH = 1. + aEtaEndcapH(t) + bEtaEndcapH(t) * cEtaEndcapH(absEta);
0271 } else {
0272 etaCorrH = 1. + aEtaEndcapH(t) + bEtaEndcapH(t) * dEtaEndcapH(absEta);
0273 }
0274 }
0275
0276
0277
0278 if (e > 0. && thresh > 0.)
0279 e = h > 0. ? threshE - threshH + etaCorrE * a * e : threshE + etaCorrE * a * e;
0280 if (h > 0. && thresh > 0.) {
0281 h = threshH + etaCorrH * b * h;
0282 }
0283 }
0284
0285
0286 if (e < 0. || h < 0.) {
0287
0288 if (e < 0.)
0289 e = ee;
0290 if (h < 0.)
0291 h = hh;
0292 }
0293
0294
0295 }
0296
0297
0298 double PFEnergyCalibration::aBarrel(double x) const {
0299 if (pfCalibrations) {
0300 BinningPointByMap point;
0301 point.insert(BinningVariables::JetEt, x);
0302 return pfCalibrations->getResult(PerformanceResult::PFfa_BARREL, point);
0303
0304 } else {
0305 return faBarrel->Eval(x);
0306 }
0307 }
0308
0309 double PFEnergyCalibration::bBarrel(double x) const {
0310 if (pfCalibrations) {
0311 BinningPointByMap point;
0312 point.insert(BinningVariables::JetEt, x);
0313 return pfCalibrations->getResult(PerformanceResult::PFfb_BARREL, point);
0314
0315 } else {
0316 return fbBarrel->Eval(x);
0317 }
0318 }
0319
0320 double PFEnergyCalibration::cBarrel(double x) const {
0321 if (pfCalibrations) {
0322 BinningPointByMap point;
0323 point.insert(BinningVariables::JetEt, x);
0324 return pfCalibrations->getResult(PerformanceResult::PFfc_BARREL, point);
0325
0326 } else {
0327 return fcBarrel->Eval(x);
0328 }
0329 }
0330
0331 double PFEnergyCalibration::aEtaBarrelEH(double x) const {
0332 if (pfCalibrations) {
0333 BinningPointByMap point;
0334 point.insert(BinningVariables::JetEt, x);
0335 return pfCalibrations->getResult(PerformanceResult::PFfaEta_BARRELEH, point);
0336
0337 } else {
0338 return faEtaBarrelEH->Eval(x);
0339 }
0340 }
0341
0342 double PFEnergyCalibration::bEtaBarrelEH(double x) const {
0343 if (pfCalibrations) {
0344 BinningPointByMap point;
0345 point.insert(BinningVariables::JetEt, x);
0346 return pfCalibrations->getResult(PerformanceResult::PFfbEta_BARRELEH, point);
0347
0348 } else {
0349 return fbEtaBarrelEH->Eval(x);
0350 }
0351 }
0352
0353 double PFEnergyCalibration::aEtaBarrelH(double x) const {
0354 if (pfCalibrations) {
0355 BinningPointByMap point;
0356 point.insert(BinningVariables::JetEt, x);
0357 return pfCalibrations->getResult(PerformanceResult::PFfaEta_BARRELH, point);
0358
0359 } else {
0360 return faEtaBarrelH->Eval(x);
0361 }
0362 }
0363
0364 double PFEnergyCalibration::bEtaBarrelH(double x) const {
0365 if (pfCalibrations) {
0366 BinningPointByMap point;
0367 point.insert(BinningVariables::JetEt, x);
0368 return pfCalibrations->getResult(PerformanceResult::PFfbEta_BARRELH, point);
0369
0370 } else {
0371 return fbEtaBarrelH->Eval(x);
0372 }
0373 }
0374
0375 double PFEnergyCalibration::aEndcap(double x) const {
0376 if (pfCalibrations) {
0377 BinningPointByMap point;
0378 point.insert(BinningVariables::JetEt, x);
0379 return pfCalibrations->getResult(PerformanceResult::PFfa_ENDCAP, point);
0380
0381 } else {
0382 return faEndcap->Eval(x);
0383 }
0384 }
0385
0386 double PFEnergyCalibration::bEndcap(double x) const {
0387 if (pfCalibrations) {
0388 BinningPointByMap point;
0389 point.insert(BinningVariables::JetEt, x);
0390 return pfCalibrations->getResult(PerformanceResult::PFfb_ENDCAP, point);
0391
0392 } else {
0393 return fbEndcap->Eval(x);
0394 }
0395 }
0396
0397 double PFEnergyCalibration::cEndcap(double x) const {
0398 if (pfCalibrations) {
0399 BinningPointByMap point;
0400 point.insert(BinningVariables::JetEt, x);
0401 return pfCalibrations->getResult(PerformanceResult::PFfc_ENDCAP, point);
0402
0403 } else {
0404 return fcEndcap->Eval(x);
0405 }
0406 }
0407
0408 double PFEnergyCalibration::aEtaEndcapEH(double x) const {
0409 if (pfCalibrations) {
0410 BinningPointByMap point;
0411 point.insert(BinningVariables::JetEt, x);
0412 return pfCalibrations->getResult(PerformanceResult::PFfaEta_ENDCAPEH, point);
0413
0414 } else {
0415 return faEtaEndcapEH->Eval(x);
0416 }
0417 }
0418
0419 double PFEnergyCalibration::bEtaEndcapEH(double x) const {
0420 if (pfCalibrations) {
0421 BinningPointByMap point;
0422 point.insert(BinningVariables::JetEt, x);
0423 return pfCalibrations->getResult(PerformanceResult::PFfbEta_ENDCAPEH, point);
0424
0425 } else {
0426 return fbEtaEndcapEH->Eval(x);
0427 }
0428 }
0429
0430 double PFEnergyCalibration::aEtaEndcapH(double x) const {
0431 if (pfCalibrations) {
0432 BinningPointByMap point;
0433 point.insert(BinningVariables::JetEt, x);
0434 return pfCalibrations->getResult(PerformanceResult::PFfaEta_ENDCAPH, point);
0435
0436 } else {
0437 return faEtaEndcapH->Eval(x);
0438 }
0439 }
0440
0441 double PFEnergyCalibration::bEtaEndcapH(double x) const {
0442 if (pfCalibrations) {
0443 BinningPointByMap point;
0444 point.insert(BinningVariables::JetEt, x);
0445 return pfCalibrations->getResult(PerformanceResult::PFfbEta_ENDCAPH, point);
0446 } else {
0447 return fbEtaEndcapH->Eval(x);
0448 }
0449 }
0450
0451
0452
0453 double PFEnergyCalibration::cEtaBarrelH(double x) const {
0454 if (pfCalibrations) {
0455 BinningPointByMap point;
0456 point.insert(BinningVariables::JetEt, x);
0457 return pfCalibrations->getResult(PerformanceResult::PFfcEta_BARRELH, point);
0458
0459 } else {
0460 return fcEtaBarrelH->Eval(x);
0461 }
0462 }
0463 double PFEnergyCalibration::cEtaEndcapH(double x) const {
0464 if (pfCalibrations) {
0465 BinningPointByMap point;
0466 point.insert(BinningVariables::JetEt, x);
0467 return pfCalibrations->getResult(PerformanceResult::PFfcEta_ENDCAPH, point);
0468
0469 } else {
0470 return fcEtaEndcapH->Eval(x);
0471 }
0472 }
0473
0474 double PFEnergyCalibration::dEtaEndcapH(double x) const {
0475 if (pfCalibrations) {
0476 BinningPointByMap point;
0477 point.insert(BinningVariables::JetEt, x);
0478 return pfCalibrations->getResult(PerformanceResult::PFfdEta_ENDCAPH, point);
0479
0480 } else {
0481 return fdEtaEndcapH->Eval(x);
0482 }
0483 }
0484
0485 double PFEnergyCalibration::cEtaBarrelEH(double x) const {
0486 if (pfCalibrations) {
0487 BinningPointByMap point;
0488 point.insert(BinningVariables::JetEt, x);
0489 return pfCalibrations->getResult(PerformanceResult::PFfcEta_BARRELEH, point);
0490
0491 } else {
0492 return fcEtaBarrelEH->Eval(x);
0493 }
0494 }
0495
0496 double PFEnergyCalibration::cEtaEndcapEH(double x) const {
0497 if (pfCalibrations) {
0498 BinningPointByMap point;
0499 point.insert(BinningVariables::JetEt, x);
0500 return pfCalibrations->getResult(PerformanceResult::PFfcEta_ENDCAPEH, point);
0501
0502 } else {
0503 return fcEtaEndcapEH->Eval(x);
0504 }
0505 }
0506
0507 double PFEnergyCalibration::dEtaEndcapEH(double x) const {
0508 if (pfCalibrations) {
0509 BinningPointByMap point;
0510 point.insert(BinningVariables::JetEt, x);
0511 return pfCalibrations->getResult(PerformanceResult::PFfdEta_ENDCAPEH, point);
0512
0513 } else {
0514 return fdEtaEndcapEH->Eval(x);
0515 }
0516 }
0517
0518 double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
0519 double ePS1,
0520 double ePS2,
0521 bool crackCorrection) const {
0522 return Ecorr(clusterEcal.energy(), ePS1, ePS2, clusterEcal.eta(), clusterEcal.phi(), crackCorrection);
0523 }
0524
0525 double PFEnergyCalibration::energyEm(const reco::PFCluster& clusterEcal,
0526 double ePS1,
0527 double ePS2,
0528 double& ps1,
0529 double& ps2,
0530 bool crackCorrection) const {
0531 return Ecorr(clusterEcal.energy(), ePS1, ePS2, clusterEcal.eta(), clusterEcal.phi(), ps1, ps2, crackCorrection);
0532 }
0533
0534 std::ostream& operator<<(std::ostream& out, const PFEnergyCalibration& calib) {
0535 if (!out)
0536 return out;
0537
0538 out << "PFEnergyCalibration -- " << endl;
0539
0540 if (calib.pfCalibrations) {
0541 static const std::map<std::string, PerformanceResult::ResultType> functType = {
0542 {"PFfa_BARREL", PerformanceResult::PFfa_BARREL},
0543 {"PFfa_ENDCAP", PerformanceResult::PFfa_ENDCAP},
0544 {"PFfb_BARREL", PerformanceResult::PFfb_BARREL},
0545 {"PFfb_ENDCAP", PerformanceResult::PFfb_ENDCAP},
0546 {"PFfc_BARREL", PerformanceResult::PFfc_BARREL},
0547 {"PFfc_ENDCAP", PerformanceResult::PFfc_ENDCAP},
0548 {"PFfaEta_BARRELH", PerformanceResult::PFfaEta_BARRELH},
0549 {"PFfaEta_ENDCAPH", PerformanceResult::PFfaEta_ENDCAPH},
0550 {"PFfbEta_BARRELH", PerformanceResult::PFfbEta_BARRELH},
0551 {"PFfbEta_ENDCAPH", PerformanceResult::PFfbEta_ENDCAPH},
0552 {"PFfaEta_BARRELEH", PerformanceResult::PFfaEta_BARRELEH},
0553 {"PFfaEta_ENDCAPEH", PerformanceResult::PFfaEta_ENDCAPEH},
0554 {"PFfbEta_BARRELEH", PerformanceResult::PFfbEta_BARRELEH},
0555 {"PFfbEta_ENDCAPEH", PerformanceResult::PFfbEta_ENDCAPEH},
0556 {"PFfcEta_BARRELH", PerformanceResult::PFfcEta_BARRELH},
0557 {"PFfcEta_ENDCAPH", PerformanceResult::PFfcEta_ENDCAPH},
0558 {"PFfdEta_ENDCAPH", PerformanceResult::PFfdEta_ENDCAPH},
0559 {"PFfcEta_BARRELEH", PerformanceResult::PFfcEta_BARRELEH},
0560 {"PFfcEta_ENDCAPEH", PerformanceResult::PFfcEta_ENDCAPEH},
0561 {"PFfdEta_ENDCAPEH", PerformanceResult::PFfdEta_ENDCAPEH}
0562
0563 };
0564
0565 for (std::map<std::string, PerformanceResult::ResultType>::const_iterator func = functType.begin();
0566 func != functType.end();
0567 ++func) {
0568 cout << "Function: " << func->first << endl;
0569 PerformanceResult::ResultType fType = func->second;
0570 calib.pfCalibrations->printFormula(fType);
0571 }
0572
0573 } else {
0574 std::cout << "Default calibration functions : " << std::endl;
0575
0576 calib.faBarrel->Print();
0577 calib.fbBarrel->Print();
0578 calib.fcBarrel->Print();
0579 calib.faEtaBarrelEH->Print();
0580 calib.fbEtaBarrelEH->Print();
0581 calib.faEtaBarrelH->Print();
0582 calib.fbEtaBarrelH->Print();
0583 calib.faEndcap->Print();
0584 calib.fbEndcap->Print();
0585 calib.fcEndcap->Print();
0586 calib.faEtaEndcapEH->Print();
0587 calib.fbEtaEndcapEH->Print();
0588 calib.faEtaEndcapH->Print();
0589 calib.fbEtaEndcapH->Print();
0590
0591 }
0592
0593 return out;
0594 }
0595
0596
0597
0598
0599
0600
0601
0602
0603
0604
0605
0606
0607
0608
0609
0610
0611
0612
0613
0614
0615
0616 double PFEnergyCalibration::minimum(double a, double b) const {
0617 if (std::abs(b) < std::abs(a))
0618 a = b;
0619 return a;
0620 }
0621
0622 namespace {
0623 constexpr double pi = M_PI;
0624
0625 std::vector<double> fillcPhi() {
0626 std::vector<double> retValue;
0627 retValue.resize(18, 0);
0628 retValue[0] = 2.97025;
0629 for (unsigned i = 1; i <= 17; ++i)
0630 retValue[i] = retValue[0] - 2 * i * pi / 18;
0631
0632 return retValue;
0633 }
0634
0635
0636 const std::vector<double> cPhi = fillcPhi();
0637 }
0638
0639
0640 double PFEnergyCalibration::dCrackPhi(double phi, double eta) const {
0641
0642 constexpr double delta_cPhi = 0.00638;
0643
0644 double m;
0645
0646
0647 if (eta < 0)
0648 phi += delta_cPhi;
0649
0650 if (phi >= -pi && phi <= pi) {
0651
0652 if (phi < cPhi[17] || phi >= cPhi[0]) {
0653 if (phi < 0)
0654 phi += 2 * pi;
0655 m = minimum(phi - cPhi[0], phi - cPhi[17] - 2 * pi);
0656 }
0657
0658
0659 else {
0660 bool OK = false;
0661 unsigned i = 16;
0662 while (!OK) {
0663 if (phi < cPhi[i]) {
0664 m = minimum(phi - cPhi[i + 1], phi - cPhi[i]);
0665 OK = true;
0666 } else
0667 i -= 1;
0668 }
0669 }
0670 } else {
0671 m = 0.;
0672 std::cout << "Problem in dminphi" << std::endl;
0673 }
0674 if (eta < 0)
0675 m = -m;
0676 return m;
0677 }
0678
0679
0680 double PFEnergyCalibration::CorrPhi(double phi, double eta) const {
0681
0682 constexpr double p1 = 5.59379e-01;
0683 constexpr double p2 = -1.26607e-03;
0684 constexpr double p3 = 9.61133e-04;
0685
0686 constexpr double p4 = 1.81691e-01;
0687 constexpr double p5 = -4.97535e-03;
0688 constexpr double p6 = 1.31006e-03;
0689
0690 constexpr double p7 = 1.38498e-01;
0691 constexpr double p8 = 1.18599e-04;
0692 constexpr double p9 = 2.01858e-03;
0693
0694 double dminphi = dCrackPhi(phi, eta);
0695
0696 double result =
0697 (1 + p1 * TMath::Gaus(dminphi, p2, p3) + p4 * TMath::Gaus(dminphi, p5, p6) + p7 * TMath::Gaus(dminphi, p8, p9));
0698
0699 return result;
0700 }
0701
0702
0703 double PFEnergyCalibration::CorrEta(double eta) const {
0704
0705 constexpr double a[] = {6.13349e-01, 5.08146e-01, 4.44480e-01, 3.3487e-01, 7.65627e-01};
0706 constexpr double m[] = {-1.79514e-02, 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00};
0707 constexpr double s[] = {7.92382e-03, 3.06028e-03, 3.36139e-03, 3.94521e-03, 8.63950e-04};
0708 constexpr double sa[] = {1.27228e+01, 3.81517e-02, 1.63507e-01, -6.56480e-02, 1.87160e-01};
0709 constexpr double ss[] = {5.48753e-02, -1.00223e-02, 2.22866e-03, 4.26288e-04, 2.67937e-03};
0710 double result = 1;
0711
0712 for (unsigned i = 0; i <= 4; i++)
0713 result += a[i] * TMath::Gaus(eta, m[i], s[i]) *
0714 (1 + sa[i] * TMath::Sign(1., eta - m[i]) * TMath::Exp(-std::abs(eta - m[i]) / ss[i]));
0715
0716 return result;
0717 }
0718
0719
0720 double PFEnergyCalibration::CorrBarrel(double E, double eta) const {
0721
0722
0723
0724
0725
0726
0727
0728
0729
0730
0731
0732
0733 constexpr double p0 = 0.9944;
0734 constexpr double p1 = 9.827;
0735 constexpr double p2 = 1.503;
0736 constexpr double p3 = 1.196;
0737 constexpr double p4 = 0.3349;
0738 constexpr double p5 = 0.89;
0739 constexpr double p6 = 0.004361;
0740 constexpr double p7 = 51.51;
0741
0742 constexpr double p8 = 2.705593e-03;
0743
0744 double result =
0745 (p0 + 1 / (p1 + p2 * TMath::Power(E, p3)) + p4 * TMath::Exp(-E / p5) + p6 * TMath::Exp(-E * E / (p7 * p7))) *
0746 (1 + p8 * eta * eta);
0747
0748 return result;
0749 }
0750
0751
0752
0753
0754
0755
0756
0757
0758
0759
0760
0761
0762
0763 double PFEnergyCalibration::Alpha(double eta) const {
0764
0765 constexpr double p0 = 5.97621e-01;
0766
0767
0768 constexpr double p1 = -1.86407e-01;
0769 constexpr double p2 = 3.85197e-01;
0770
0771
0772 constexpr double norm = (p1 + p2 * (2.6 + 1.656) / 2);
0773
0774 double result = p0 * (p1 + p2 * eta) / norm;
0775
0776 return result;
0777 }
0778
0779 double PFEnergyCalibration::Beta(double E, double eta) const {
0780
0781 constexpr double p0 = 0.032;
0782 constexpr double p1 = 9.70394e-02;
0783 constexpr double p2 = 2.23072e+01;
0784 constexpr double p3 = 100;
0785
0786
0787 constexpr double p4 = 1.02496e+00;
0788 constexpr double p5 = -4.40176e-03;
0789
0790
0791 constexpr double norm = (p4 + p5 * (2.6 + 1.656) / 2);
0792
0793 double result = (1.0012 + p0 * TMath::Exp(-E / p3) + p1 * TMath::Exp(-E / p2)) * (p4 + p5 * eta) / norm;
0794 return result;
0795 }
0796
0797 double PFEnergyCalibration::Gamma(double etaEcal) const {
0798
0799 constexpr double p0 = 2.49752e-02;
0800
0801
0802 constexpr double p1 = 6.48816e-02;
0803 constexpr double p2 = -1.59517e-02;
0804
0805
0806 constexpr double norm = (p1 + p2 * (2.6 + 1.656) / 2);
0807
0808 double result = p0 * (p1 + p2 * etaEcal) / norm;
0809
0810 return result;
0811 }
0812
0813
0814
0815
0816
0817
0818
0819
0820
0821 double PFEnergyCalibration::EcorrBarrel(double E, double eta, double phi, bool crackCorrection) const {
0822
0823 double correction = crackCorrection ? std::max(CorrEta(eta), CorrPhi(phi, eta)) : 1.;
0824 double result = E * CorrBarrel(E, eta) * correction;
0825
0826 return result;
0827 }
0828
0829
0830 double PFEnergyCalibration::EcorrZoneBeforePS(double E, double eta) const {
0831
0832 constexpr double p0 = 1;
0833 constexpr double p1 = 0.18;
0834 constexpr double p2 = 8.;
0835
0836
0837 constexpr double p3 = 0.3;
0838 constexpr double p4 = 1.11;
0839 constexpr double p5 = 0.025;
0840 constexpr double p6 = 1.49;
0841 constexpr double p7 = 0.6;
0842
0843
0844 constexpr double norm = 1.21;
0845
0846 double result = E * (p0 + p1 * TMath::Exp(-E / p2)) * (p3 + p4 * TMath::Gaus(eta, p6, p5) + p7 * eta) / norm;
0847
0848 return result;
0849 }
0850
0851
0852
0853 double PFEnergyCalibration::EcorrPS(double eEcal, double ePS1, double ePS2, double etaEcal) const {
0854
0855 double E = Beta(1.0155 * eEcal + 0.025 * (ePS1 + 0.5976 * ePS2) / 9e-5, etaEcal) * eEcal +
0856 Gamma(etaEcal) * (ePS1 + Alpha(etaEcal) * ePS2) / 9e-5;
0857
0858
0859 constexpr double p0 = 1.00;
0860 constexpr double p1 = 2.18;
0861 constexpr double p2 = 1.94;
0862 constexpr double p3 = 4.13;
0863 constexpr double p4 = 1.127;
0864
0865 double result = E * (p0 + p1 * TMath::Exp(-E / p2) - p3 * TMath::Exp(-E / p4));
0866
0867 return result;
0868 }
0869
0870
0871
0872 double PFEnergyCalibration::EcorrPS(
0873 double eEcal, double ePS1, double ePS2, double etaEcal, double& outputPS1, double& outputPS2) const {
0874
0875 double gammaprime = Gamma(etaEcal) / 9e-5;
0876
0877 if (outputPS1 == 0 && outputPS2 == 0 && esEEInterCalib_ != nullptr) {
0878
0879
0880 outputPS1 = gammaprime * ePS1 * esEEInterCalib_->getGammaLow0();
0881 outputPS2 = gammaprime * Alpha(etaEcal) * ePS2 * esEEInterCalib_->getGammaLow3();
0882 } else if (outputPS1 == 0 && outputPS2 == -1 && esEEInterCalib_ != nullptr) {
0883
0884 double corrTotES = gammaprime * ePS1 * esEEInterCalib_->getGammaLow0() * esEEInterCalib_->getGammaLow1();
0885 outputPS1 = gammaprime * ePS1 * esEEInterCalib_->getGammaLow0();
0886 outputPS2 = corrTotES - outputPS1;
0887 } else if (outputPS1 == -1 && outputPS2 == 0 && esEEInterCalib_ != nullptr) {
0888
0889 double corrTotES =
0890 gammaprime * Alpha(etaEcal) * ePS2 * esEEInterCalib_->getGammaLow3() * esEEInterCalib_->getGammaLow2();
0891 outputPS2 = gammaprime * Alpha(etaEcal) * ePS2 * esEEInterCalib_->getGammaLow3();
0892 outputPS1 = corrTotES - outputPS2;
0893 } else {
0894
0895 outputPS1 = gammaprime * ePS1;
0896 outputPS2 = gammaprime * Alpha(etaEcal) * ePS2;
0897 }
0898
0899 double E = Beta(1.0155 * eEcal + 0.025 * (ePS1 + 0.5976 * ePS2) / 9e-5, etaEcal) * eEcal + outputPS1 + outputPS2;
0900
0901
0902 constexpr double p0 = 1.00;
0903 constexpr double p1 = 2.18;
0904 constexpr double p2 = 1.94;
0905 constexpr double p3 = 4.13;
0906 constexpr double p4 = 1.127;
0907
0908 double corrfac = (p0 + p1 * TMath::Exp(-E / p2) - p3 * TMath::Exp(-E / p4));
0909 outputPS1 *= corrfac;
0910 outputPS2 *= corrfac;
0911 double result = E * corrfac;
0912
0913 return result;
0914 }
0915
0916
0917
0918 double PFEnergyCalibration::EcorrPS_ePSNil(double eEcal, double eta) const {
0919
0920 constexpr double p0 = 1.02;
0921 constexpr double p1 = 0.165;
0922 constexpr double p2 = 6.5;
0923 constexpr double p3 = 2.1;
0924
0925
0926 constexpr double p4 = 1.02496e+00;
0927 constexpr double p5 = -4.40176e-03;
0928
0929
0930 constexpr double norm = (p4 + p5 * (2.6 + 1.656) / 2);
0931
0932 double result = eEcal * (p0 + p1 * TMath::Exp(-std::abs(eEcal - p3) / p2)) * (p4 + p5 * eta) / norm;
0933
0934 return result;
0935 }
0936
0937
0938 double PFEnergyCalibration::EcorrZoneAfterPS(double E, double eta) const {
0939
0940 constexpr double p0 = 1;
0941 constexpr double p1 = 0.058;
0942 constexpr double p2 = 12.5;
0943 constexpr double p3 = -1.05444e+00;
0944 constexpr double p4 = -5.39557e+00;
0945 constexpr double p5 = 8.38444e+00;
0946 constexpr double p6 = 6.10998e-01;
0947
0948
0949 constexpr double p7 = 1.06161e+00;
0950 constexpr double p8 = 0.41;
0951 constexpr double p9 = 2.918;
0952 constexpr double p10 = 0.0181;
0953 constexpr double p11 = 2.05;
0954 constexpr double p12 = 2.99;
0955 constexpr double p13 = 0.0287;
0956
0957
0958 constexpr double norm = 1.045;
0959
0960 double result = E * (p0 + p1 * TMath::Exp(-(E - p3) / p2) + 1 / (p4 + p5 * TMath::Power(E, p6))) *
0961 (p7 + p8 * TMath::Gaus(eta, p9, p10) + p11 * TMath::Gaus(eta, p12, p13)) / norm;
0962 return result;
0963 }
0964
0965
0966
0967 double PFEnergyCalibration::Ecorr(
0968 double eEcal, double ePS1, double ePS2, double eta, double phi, bool crackCorrection) const {
0969 constexpr double endBarrel = 1.48;
0970 constexpr double beginingPS = 1.65;
0971 constexpr double endPS = 2.6;
0972 constexpr double endEndCap = 2.98;
0973
0974 double result = 0;
0975
0976 eta = std::abs(eta);
0977
0978 if (eEcal > 0) {
0979 if (eta <= endBarrel)
0980 result = EcorrBarrel(eEcal, eta, phi, crackCorrection);
0981 else if (eta <= beginingPS)
0982 result = EcorrZoneBeforePS(eEcal, eta);
0983 else if ((eta < endPS) && ePS1 == 0 && ePS2 == 0)
0984 result = EcorrPS_ePSNil(eEcal, eta);
0985 else if (eta < endPS)
0986 result = EcorrPS(eEcal, ePS1, ePS2, eta);
0987 else if (eta < endEndCap)
0988 result = EcorrZoneAfterPS(eEcal, eta);
0989 else
0990 result = eEcal;
0991 } else
0992 result = eEcal;
0993
0994 if (result < eEcal)
0995 result = eEcal;
0996 return result;
0997 }
0998
0999
1000
1001 double PFEnergyCalibration::Ecorr(double eEcal,
1002 double ePS1,
1003 double ePS2,
1004 double eta,
1005 double phi,
1006 double& ps1,
1007 double& ps2,
1008 bool crackCorrection) const {
1009 constexpr double endBarrel = 1.48;
1010 constexpr double beginingPS = 1.65;
1011 constexpr double endPS = 2.6;
1012 constexpr double endEndCap = 2.98;
1013
1014 double result = 0;
1015
1016 eta = std::abs(eta);
1017
1018 if (eEcal > 0) {
1019 if (eta <= endBarrel)
1020 result = EcorrBarrel(eEcal, eta, phi, crackCorrection);
1021 else if (eta <= beginingPS)
1022 result = EcorrZoneBeforePS(eEcal, eta);
1023 else if ((eta < endPS) && ePS1 == 0 && ePS2 == 0)
1024 result = EcorrPS_ePSNil(eEcal, eta);
1025 else if (eta < endPS)
1026 result = EcorrPS(eEcal, ePS1, ePS2, eta, ps1, ps2);
1027 else if (eta < endEndCap)
1028 result = EcorrZoneAfterPS(eEcal, eta);
1029 else
1030 result = eEcal;
1031 } else
1032 result = eEcal;
1033
1034 if (result < eEcal)
1035 result = eEcal;
1036 return result;
1037 }