Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:26:27

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   //calibChrisClean.C calibration parameters bhumika Nov, 2018
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   //added by Bhumika on 2 august 2018
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         //getStatusCode() == 0 => active channel
0170         // apply correction if all recHits are dead
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   // Use calorimetric energy as true energy for neutral particles
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   // Barrel calibration
0208   if (absEta < 1.48) {
0209     // The energy correction
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     // Protection against negative calibration
0215     if (a < -0.25 || b < -0.25) {
0216       a = 1.;
0217       b = 1.;
0218       thresh = 0.;
0219     }
0220 
0221     // The new estimate of the true energy
0222     t = min(999.9, max(tt, thresh + a * e + b * h));
0223 
0224     // The angular correction
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     // Endcap calibration
0239   } else {
0240     // The energy correction
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     // Protection against negative calibration
0246     if (a < -0.25 || b < -0.25) {
0247       a = 1.;
0248       b = 1.;
0249       thresh = 0.;
0250     }
0251 
0252     // The new estimate of the true energy
0253     t = min(999.9, max(tt, thresh + a * e + b * h));
0254 
0255     // The angular correction
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     //t = min(999.9,max(tt, thresh + etaCorrE*a*e + etaCorrH*b*h));
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   // Protection
0286   if (e < 0. || h < 0.) {
0287     // Some protection against crazy calibration
0288     if (e < 0.)
0289       e = ee;
0290     if (h < 0.)
0291       h = hh;
0292   }
0293 
0294   // And that's it !
0295 }
0296 
0297 // The calibration functions
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 //added by Bhumika Kansal on 3 august 2018
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 ////             CORRECTION OF PHOTONS' ENERGY             ////
0599 ////                                                       ////
0600 ////              Material effect: No tracker              ////
0601 ////       Tuned on CMSSW_2_1_0_pre4, Full Sim events      ////
0602 ////                                                       ////
0603 ///////////////////////////////////////////////////////////////
0604 ////                                                       ////
0605 ////            Jonathan Biteau - June 2008                ////
0606 ////                                                       ////
0607 ///////////////////////////////////////////////////////////////
0608 
0609 ///////////////////////////////////////////////////////////////
0610 ////                                                       ////
0611 ////  USEFUL FUNCTIONS FOR THE CORRECTION IN THE BARREL    ////
0612 ////                                                       ////
0613 ///////////////////////////////////////////////////////////////
0614 
0615 //useful to compute the signed distance to the closest crack in the barrel
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;  // 3.14159265358979323846;
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   //Location of the 18 phi-cracks
0636   const std::vector<double> cPhi = fillcPhi();
0637 }  // namespace
0638 
0639 //compute the unsigned distance to the closest phi-crack in the barrel
0640 double PFEnergyCalibration::dCrackPhi(double phi, double eta) const {
0641   //Shift of this location if eta<0
0642   constexpr double delta_cPhi = 0.00638;
0643 
0644   double m;  //the result
0645 
0646   //the location is shifted
0647   if (eta < 0)
0648     phi += delta_cPhi;
0649 
0650   if (phi >= -pi && phi <= pi) {
0651     //the problem of the extrema
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     //between these extrema...
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.;  //if there is a problem, we assum that we are in a crack
0672     std::cout << "Problem in dminphi" << std::endl;
0673   }
0674   if (eta < 0)
0675     m = -m;  //because of the disymetry
0676   return m;
0677 }
0678 
0679 // corrects the effect of phi-cracks
0680 double PFEnergyCalibration::CorrPhi(double phi, double eta) const {
0681   // we use 3 gaussians to correct the phi-cracks effect
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 // corrects the effect of  |eta|-cracks
0703 double PFEnergyCalibration::CorrEta(double eta) const {
0704   // we use a gaussian with a screwness for each of the 5 |eta|-cracks
0705   constexpr double a[] = {6.13349e-01, 5.08146e-01, 4.44480e-01, 3.3487e-01, 7.65627e-01};     // amplitude
0706   constexpr double m[] = {-1.79514e-02, 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00};   // mean
0707   constexpr double s[] = {7.92382e-03, 3.06028e-03, 3.36139e-03, 3.94521e-03, 8.63950e-04};    // sigma
0708   constexpr double sa[] = {1.27228e+01, 3.81517e-02, 1.63507e-01, -6.56480e-02, 1.87160e-01};  // screwness amplitude
0709   constexpr double ss[] = {5.48753e-02, -1.00223e-02, 2.22866e-03, 4.26288e-04, 2.67937e-03};  // screwness sigma
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 //corrects the global behaviour in the barrel
0720 double PFEnergyCalibration::CorrBarrel(double E, double eta) const {
0721   //Energy dependency
0722   /*
0723   //YM Parameters 52XX:
0724   constexpr double p0=1.00000e+00;
0725   constexpr double p1=3.27753e+01;
0726   constexpr double p2=2.28552e-02;
0727   constexpr double p3=3.06139e+00;
0728   constexpr double p4=2.25135e-01;
0729   constexpr double p5=1.47824e+00;
0730   constexpr double p6=1.09e-02;
0731   constexpr double p7=4.19343e+01;
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   //Eta dependency
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 ////  USEFUL FUNCTIONS FOR THE CORRECTION IN THE ENDCAPS   ////
0754 ////  Parameters tuned for:                                ////
0755 ////          dR(ClustersPS1,ClusterEcal) < 0.08           ////
0756 ////          dR(ClustersPS2,ClusterEcal) < 0.13           ////
0757 ////                                                       ////
0758 ///////////////////////////////////////////////////////////////
0759 
0760 //Alpha, Beta, Gamma give the weight of each sub-detector (PS layer1, PS layer2 and Ecal) in the areas of the endcaps where there is a PS
0761 // Etot = Beta*eEcal + Gamma*(ePS1 + Alpha*ePS2)
0762 
0763 double PFEnergyCalibration::Alpha(double eta) const {
0764   //Energy dependency
0765   constexpr double p0 = 5.97621e-01;
0766 
0767   //Eta dependency
0768   constexpr double p1 = -1.86407e-01;
0769   constexpr double p2 = 3.85197e-01;
0770 
0771   //so that <feta()> = 1
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   //Energy dependency
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   //Eta dependency
0787   constexpr double p4 = 1.02496e+00;
0788   constexpr double p5 = -4.40176e-03;
0789 
0790   //so that <feta()> = 1
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   //Energy dependency
0799   constexpr double p0 = 2.49752e-02;
0800 
0801   //Eta dependency
0802   constexpr double p1 = 6.48816e-02;
0803   constexpr double p2 = -1.59517e-02;
0804 
0805   //so that <feta()> = 1
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 ////   THE CORRECTIONS IN THE BARREL AND IN THE ENDCAPS    ////
0816 ////                                                       ////
0817 ///////////////////////////////////////////////////////////////
0818 
0819 // returns the corrected energy in the barrel (0,1.48)
0820 // Global Behaviour, phi and eta cracks are taken into account
0821 double PFEnergyCalibration::EcorrBarrel(double E, double eta, double phi, bool crackCorrection) const {
0822   // double result = E*CorrBarrel(E,eta)*CorrEta(eta)*CorrPhi(phi,eta);
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 // returns the corrected energy in the area between the barrel and the PS (1.48,1.65)
0830 double PFEnergyCalibration::EcorrZoneBeforePS(double E, double eta) const {
0831   //Energy dependency
0832   constexpr double p0 = 1;
0833   constexpr double p1 = 0.18;
0834   constexpr double p2 = 8.;
0835 
0836   //Eta dependency
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   //so that <feta()> = 1
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 // returns the corrected energy in the PS (1.65,2.6)
0852 // only when (ePS1>0)||(ePS2>0)
0853 double PFEnergyCalibration::EcorrPS(double eEcal, double ePS1, double ePS2, double etaEcal) const {
0854   // gives the good weights to each subdetector
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   //Correction of the residual energy dependency
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 // returns the corrected energy in the PS (1.65,2.6)
0871 // only when (ePS1>0)||(ePS2>0)
0872 double PFEnergyCalibration::EcorrPS(
0873     double eEcal, double ePS1, double ePS2, double etaEcal, double& outputPS1, double& outputPS2) const {
0874   // gives the good weights to each subdetector
0875   double gammaprime = Gamma(etaEcal) / 9e-5;
0876 
0877   if (outputPS1 == 0 && outputPS2 == 0 && esEEInterCalib_ != nullptr) {
0878     // both ES planes working
0879     // scaling factor accounting for data-mc
0880     outputPS1 = gammaprime * ePS1 * esEEInterCalib_->getGammaLow0();
0881     outputPS2 = gammaprime * Alpha(etaEcal) * ePS2 * esEEInterCalib_->getGammaLow3();
0882   } else if (outputPS1 == 0 && outputPS2 == -1 && esEEInterCalib_ != nullptr) {
0883     // ESP1 only working
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     // ESP2 only working
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     // none working
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   //Correction of the residual energy dependency
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 // returns the corrected energy in the PS (1.65,2.6)
0917 // only when (ePS1=0)&&(ePS2=0)
0918 double PFEnergyCalibration::EcorrPS_ePSNil(double eEcal, double eta) const {
0919   //Energy dependency
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   //Eta dependency
0926   constexpr double p4 = 1.02496e+00;
0927   constexpr double p5 = -4.40176e-03;
0928 
0929   //so that <feta()> = 1
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 // returns the corrected energy in the area between the end of the PS and the end of the endcap (2.6,2.98)
0938 double PFEnergyCalibration::EcorrZoneAfterPS(double E, double eta) const {
0939   //Energy dependency
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   //Eta dependency
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   //so that <feta()> = 1
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 // returns the corrected energy everywhere
0966 // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
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;  // useful if eEcal=0 or eta>2.98
0993   //protection
0994   if (result < eEcal)
0995     result = eEcal;
0996   return result;
0997 }
0998 
0999 // returns the corrected energy everywhere
1000 // this work should be improved between 1.479 and 1.52 (junction barrel-endcap)
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;  // useful if eEcal=0 or eta>2.98
1033   // protection
1034   if (result < eEcal)
1035     result = eEcal;
1036   return result;
1037 }