EgammaTools/src/EnergyScaleCorrection.cc

0001 #include "RecoEgamma/EgammaTools/interface/EnergyScaleCorrection.h"
0002
0003 #include "FWCore/ParameterSet/interface/FileInPath.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include <cstdlib>
0007 #include <cmath>
0008 #include <iomanip>
0009 #include <algorithm>
0010 #include <sstream>
0011 #include <iterator>
0012
0013 EnergyScaleCorrection::EnergyScaleCorrection(const std::string& correctionFileName, unsigned int genSeed)
0014     : smearingType_(ECALELF) {
0015   if (!correctionFileName.empty()) {
0016     std::string filename = correctionFileName + "_scales.dat";
0017     readScalesFromFile(filename);
0018     if (scales_.empty()) {
0019       throw cms::Exception("EnergyScaleCorrection") << "scale correction map empty";
0020     }
0021   }
0022
0023   if (!correctionFileName.empty()) {
0024     std::string filename = correctionFileName + "_smearings.dat";
0025     readSmearingsFromFile(filename);
0026     if (smearings_.empty()) {
0027       throw cms::Exception("EnergyScaleCorrection") << "smearing correction map empty";
0028     }
0029   }
0030 }
0031
0032 float EnergyScaleCorrection::scaleCorr(unsigned int runNumber,
0033                                        double et,
0034                                        double eta,
0035                                        double r9,
0036                                        unsigned int gainSeed,
0037                                        std::bitset<kErrNrBits> uncBitMask) const {
0038   const ScaleCorrection* scaleCorr = getScaleCorr(runNumber, et, eta, r9, gainSeed);
0039   if (scaleCorr != nullptr)
0040     return scaleCorr->scale();
0041   else
0042     return kDefaultScaleVal_;
0043 }
0044
0045 float EnergyScaleCorrection::scaleCorrUncert(unsigned int runNumber,
0046                                              double et,
0047                                              double eta,
0048                                              double r9,
0049                                              unsigned int gainSeed,
0050                                              std::bitset<kErrNrBits> uncBitMask) const {
0051   const ScaleCorrection* scaleCorr = getScaleCorr(runNumber, et, eta, r9, gainSeed);
0052   if (scaleCorr != nullptr)
0053     return scaleCorr->scaleErr(uncBitMask);
0054   else
0055     return 0.;
0056 }
0057
0058 float EnergyScaleCorrection::smearingSigma(
0059     int runnr, double et, double eta, double r9, unsigned int gainSeed, ParamSmear par, float nSigma) const {
0060   if (par == kRho)
0061     return smearingSigma(runnr, et, eta, r9, gainSeed, nSigma, 0.);
0062   if (par == kPhi)
0063     return smearingSigma(runnr, et, eta, r9, gainSeed, 0., nSigma);
0064   return smearingSigma(runnr, et, eta, r9, gainSeed, 0., 0.);
0065 }
0066
0067 float EnergyScaleCorrection::smearingSigma(
0068     int runnr, double et, double eta, double r9, unsigned int gainSeed, float nrSigmaRho, float nrSigmaPhi) const {
0069   const SmearCorrection* smearCorr = getSmearCorr(runnr, et, eta, r9, gainSeed);
0070
0071   if (smearCorr != nullptr)
0072     return smearCorr->sigma(nrSigmaRho, nrSigmaPhi);
0073   else
0074     return kDefaultSmearVal_;
0075 }
0076
0077 const EnergyScaleCorrection::ScaleCorrection* EnergyScaleCorrection::getScaleCorr(
0078     unsigned int runnr, double et, double eta, double r9, unsigned int gainSeed) const {
0079   // buld the category based on the values of the object
0080   CorrectionCategory category(runnr, et, eta, r9, gainSeed);
0081   auto result = scales_.find(category);
0082
0083   if (result == scales_.end()) {
0084     edm::LogInfo("EnergyScaleCorrection")
0085         << "Scale category not found: " << category << " Returning uncorrected value.";
0086     return nullptr;
0087   }
0088
0089   //validate the result, just to be sure
0090   if (!result->first.inCategory(runnr, et, eta, r9, gainSeed)) {
0091     throw cms::Exception("LogicError") << " error found scale category " << result->first
0092                                        << " that does not contain run " << runnr << " et " << et << " eta " << eta
0093                                        << " r9 " << r9 << " gain seed " << gainSeed;
0094   }
0095   return &result->second;
0096 }
0097
0098 const EnergyScaleCorrection::SmearCorrection* EnergyScaleCorrection::getSmearCorr(
0099     unsigned int runnr, double et, double eta, double r9, unsigned int gainSeed) const {
0100   // buld the category based on the values of the object
0101   CorrectionCategory category(runnr, et, eta, r9, gainSeed);
0102   auto result = smearings_.find(category);
0103
0104   if (result == smearings_.end()) {
0105     edm::LogInfo("EnergyScaleCorrection")
0106         << "Smear category not found: " << category << " Returning uncorrected value.";
0107     return nullptr;
0108   }
0109
0110   //validate the result, just to be sure
0111   if (!result->first.inCategory(runnr, et, eta, r9, gainSeed)) {
0112     throw cms::Exception("LogicError") << " error found smear category " << result->first
0113                                        << " that does not contain run " << runnr << " et " << et << " eta " << eta
0114                                        << " r9 " << r9 << " gain seed " << gainSeed;
0115   }
0116   return &result->second;
0117 }
0118
0119 void EnergyScaleCorrection::addScale(const std::string& category,
0120                                      int runMin,
0121                                      int runMax,
0122                                      double energyScale,
0123                                      double energyScaleErrStat,
0124                                      double energyScaleErrSyst,
0125                                      double energyScaleErrGain) {
0126   CorrectionCategory cat(category, runMin, runMax);  // build the category from the string
0127   auto result = scales_.equal_range(cat);
0128   if (result.first != result.second) {
0129     throw cms::Exception("ConfigError") << "Category already defined! " << cat;
0130   }
0131
0132   ScaleCorrection corr(energyScale, energyScaleErrStat, energyScaleErrSyst, energyScaleErrGain);
0133   scales_.insert(result.first, {cat, corr});  //use a hint where to insert
0134 }
0135
0136 void EnergyScaleCorrection::addScale(int runMin,
0137                                      int runMax,
0138                                      double etaMin,
0139                                      double etaMax,
0140                                      double r9Min,
0141                                      double r9Max,
0142                                      double etMin,
0143                                      double etMax,
0144                                      unsigned int gain,
0145                                      double energyScale,
0146                                      double energyScaleErrStat,
0147                                      double energyScaleErrSyst,
0148                                      double energyScaleErrGain) {
0149   CorrectionCategory cat(runMin, runMax, etaMin, etaMax, r9Min, r9Max, etMin, etMax, gain);
0150
0151   auto result = scales_.equal_range(cat);
0152   if (result.first != result.second) {
0153     throw cms::Exception("ConfigError") << "Category already defined! " << cat;
0154   }
0155
0156   ScaleCorrection corr(energyScale, energyScaleErrStat, energyScaleErrSyst, energyScaleErrGain);
0157   scales_.insert(result.first, {cat, corr});
0158 }
0159
0160 void EnergyScaleCorrection::addSmearing(const std::string& category,
0161                                         int runMin,
0162                                         int runMax,
0163                                         double rho,
0164                                         double errRho,
0165                                         double phi,
0166                                         double errPhi,
0167                                         double eMean,
0168                                         double errEMean) {
0169   CorrectionCategory cat(category);
0170   auto res = smearings_.equal_range(cat);
0171
0172   if (res.first != res.second) {
0173     throw cms::Exception("EnergyScaleCorrection") << "Smearing category already defined " << cat;
0174   }
0175
0176   SmearCorrection corr(rho, errRho, phi, errPhi, eMean, errEMean);
0177   smearings_.insert(res.first, {cat, corr});  //use a hint from res
0178 }
0179
0180 void EnergyScaleCorrection::setSmearingType(FileFormat value) {
0181   if (value <= 1) {
0182     smearingType_ = value;
0183   } else {
0184     smearingType_ = UNKNOWN;
0185   }
0186 }
0187
0188 void EnergyScaleCorrection::readScalesFromFile(const std::string& filename) {
0189   std::ifstream file(edm::FileInPath(filename).fullPath().c_str());
0190
0191   if (!file.good()) {
0192     throw cms::Exception("EnergyScaleCorrection") << "file " << filename << " not readable.";
0193   }
0194
0195   int runMin, runMax;
0196   std::string category, region2;
0197   double energyScale, energyScaleErr, energyScaleErrStat, energyScaleErrSyst, energyScaleErrGain;
0198
0199   double etaMin;      ///< Min eta value for the bin
0200   double etaMax;      ///< Max eta value for the bin
0201   double r9Min;       ///< Min R9 vaule for the bin
0202   double r9Max;       ///< Max R9 value for the bin
0203   double etMin;       ///< Min Et value for the bin
0204   double etMax;       ///< Max Et value for the bin
0205   unsigned int gain;  ///< 12, 6, 1, 61 (double gain switch)
0206
0207   // TO count the #columns in that txt file and decide based on that the version to read
0208   std::string line;
0209   std::stringstream stream;
0210   getline(file, line);
0211   stream.clear();
0212   stream << line;
0213
0214   int ncols = std::distance(std::istream_iterator<std::string>(stream), std::istream_iterator<std::string>());
0215
0216   file.seekg(0, std::ios::beg);
0217
0218   if (ncols == 9) {
0219     for (file >> category; file.good(); file >> category) {
0220       file >> region2 >> runMin >> runMax >> energyScale >> energyScaleErr >> energyScaleErrStat >>
0221           energyScaleErrSyst >> energyScaleErrGain;
0222       addScale(category, runMin, runMax, energyScale, energyScaleErrStat, energyScaleErrSyst, energyScaleErrGain);
0223     }
0224   } else {
0225     if (file.peek() == 'r')
0226       file.ignore(1000, 10);
0227
0228     for (file >> runMin; file.good(); file >> runMin) {
0229       file >> runMax >> etaMin >> etaMax >> r9Min >> r9Max >> etMin >> etMax >> gain >> energyScale >> energyScaleErr;
0230       file.ignore(1000, 10);
0231       energyScaleErrStat = energyScaleErr;
0232       energyScaleErrSyst = 0;
0233       energyScaleErrGain = 0;
0234
0235       addScale(runMin,
0236                runMax,
0237                etaMin,
0238                etaMax,
0239                r9Min,
0240                r9Max,
0241                etMin,
0242                etMax,
0243                gain,
0244                energyScale,
0245                energyScaleErrStat,
0246                energyScaleErrSyst,
0247                energyScaleErrGain);
0248     }
0249   }
0250
0251   file.close();
0252   return;
0253 }
0254
0255 //also more or less untouched function from the orginal package
0256 void EnergyScaleCorrection::readSmearingsFromFile(const std::string& filename) {
0257   std::ifstream file(edm::FileInPath(filename).fullPath().c_str());
0258   if (!file.good()) {
0259     throw cms::Exception("EnergyScaleCorrection") << "file " << filename << " not readable";
0260   }
0261
0262   int runMin = 0;
0263   int runMax = 900000;
0264   int unused = 0;
0265   std::string category, region2;
0266   double rho, phi, eMean, errRho, errPhi, errEMean;
0267   double etaMin, etaMax, r9Min, r9Max;
0268   std::string phiString, errPhiString;
0269
0270   while (file.peek() != EOF && file.good()) {
0271     if (file.peek() == 10) {  // 10 = \n
0272       file.get();
0273       continue;
0274     }
0275
0276     if (file.peek() == 35) {  // 35 = #
0277       file.ignore(1000, 10);  // ignore the rest of the line until \n
0278       continue;
0279     }
0280
0281     if (smearingType_ == UNKNOWN) {  // trying to guess: not recommended
0282       throw cms::Exception("ConfigError") << "unknown smearing type";
0283
0284     } else if (smearingType_ == GLOBE) {
0285       file >> category >> unused >> etaMin >> etaMax >> r9Min >> r9Max >> runMin >> runMax >> eMean >> errEMean >>
0286           rho >> errRho >> phi >> errPhi;
0287
0288       addSmearing(category, runMin, runMax, rho, errRho, phi, errPhi, eMean, errEMean);
0289
0290     } else if (smearingType_ == ECALELF) {
0291       file >> category >> eMean >> errEMean >> rho >> errRho >> phiString >> errPhiString;
0292
0293       if (phiString == "M_PI_2")
0294         phi = M_PI_2;
0295       else
0296         phi = std::stod(phiString);
0297
0298       if (errPhiString == "M_PI_2")
0299         errPhi = M_PI_2;
0300       else
0301         errPhi = std::stod(errPhiString);
0302
0303       addSmearing(category, runMin, runMax, rho, errRho, phi, errPhi, eMean, errEMean);
0304
0305     } else {
0306       file >> category >> rho >> phi;
0307       errRho = errPhi = eMean = errEMean = 0;
0308       addSmearing(category, runMin, runMax, rho, errRho, phi, errPhi, eMean, errEMean);
0309     }
0310   }
0311
0312   file.close();
0313   return;
0314 }
0315
0316 std::ostream& EnergyScaleCorrection::ScaleCorrection::print(std::ostream& os) const {
0317   os << "( " << scale_ << " +/- " << scaleErrStat_ << " +/- " << scaleErrSyst_ << " +/- " << scaleErrGain_ << ")";
0318   return os;
0319 }
0320
0321 float EnergyScaleCorrection::ScaleCorrection::scaleErr(const std::bitset<kErrNrBits>& uncBitMask) const {
0322   double totErr(0);
0323   auto pow2 = [](const double& x) { return x * x; };
0324
0325   if (uncBitMask.test(kErrStatBitNr))
0326     totErr += pow2(scaleErrStat_);
0327   if (uncBitMask.test(kErrSystBitNr))
0328     totErr += pow2(scaleErrSyst_);
0329   if (uncBitMask.test(kErrGainBitNr))
0330     totErr += pow2(scaleErrGain_);
0331
0332   return std::sqrt(totErr);
0333 }
0334
0335 std::ostream& EnergyScaleCorrection::SmearCorrection::print(std::ostream& os) const {
0336   os << rho_ << " +/- " << rhoErr_ << "\t" << phi_ << " +/- " << phiErr_ << "\t" << eMean_ << " +/- " << eMeanErr_;
0337   return os;
0338 }
0339
0340 //here be dragons
0341 //this function is nasty and needs to be replaced
0342 EnergyScaleCorrection::CorrectionCategory::CorrectionCategory(const std::string& category, int runnrMin, int runnrMax)
0343     : runMin_(runnrMin),
0344       runMax_(runnrMax),
0345       etaMin_(0),
0346       etaMax_(3),
0347       r9Min_(-1),
0348       r9Max_(999),
0349       etMin_(0),
0350       etMax_(9999999),
0351       gain_(0) {
0352   size_t p1, p2;  // boundary
0353
0354   // eta region
0355   p1 = category.find("absEta_");
0356   if (category.find("absEta_0_1") != std::string::npos) {
0357     etaMin_ = 0;
0358     etaMax_ = 1;
0359   } else if (category.find("absEta_1_1.4442") != std::string::npos) {
0360     etaMin_ = 1;
0361     etaMax_ = 1.479;
0362   } else if (category.find("absEta_1.566_2") != std::string::npos) {
0363     etaMin_ = 1.479;
0364     etaMax_ = 2;
0365   } else if (category.find("absEta_2_2.5") != std::string::npos) {
0366     etaMin_ = 2;
0367     etaMax_ = 3;
0368   } else {
0369     if (p1 != std::string::npos) {
0370       p1 = category.find('_', p1);
0371       p2 = category.find('_', p1 + 1);
0372       etaMin_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1));
0373       p1 = p2;
0374       p2 = category.find('-', p1);
0375       etaMax_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1));
0376     }
0377   }
0378
0379   if (category.find("EBlowEta") != std::string::npos) {
0380     etaMin_ = 0;
0381     etaMax_ = 1;
0382   };
0383   if (category.find("EBhighEta") != std::string::npos) {
0384     etaMin_ = 1;
0385     etaMax_ = 1.479;
0386   };
0387   if (category.find("EElowEta") != std::string::npos) {
0388     etaMin_ = 1.479;
0389     etaMax_ = 2;
0390   };
0391   if (category.find("EEhighEta") != std::string::npos) {
0392     etaMin_ = 2;
0393     etaMax_ = 7;
0394   };
0395
0396   // Et region
0397   p1 = category.find("-Et_");
0398
0399   if (p1 != std::string::npos) {
0400     p1 = category.find('_', p1);
0401     p2 = category.find('_', p1 + 1);
0402     etMin_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1));
0403     p1 = p2;
0404     p2 = category.find('-', p1);
0405     etMax_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1));
0406   }
0407
0408   if (category.find("gold") != std::string::npos || category.find("Gold") != std::string::npos ||
0409       category.find("highR9") != std::string::npos) {
0410     r9Min_ = 0.94;
0411     r9Max_ = std::numeric_limits<float>::max();
0412   } else if (category.find("bad") != std::string::npos || category.find("Bad") != std::string::npos ||
0413              category.find("lowR9") != std::string::npos) {
0414     r9Min_ = -1;
0415     r9Max_ = 0.94;
0416   };
0417   // R9 region
0418   p1 = category.find("-R9");
0419   if (p1 != std::string::npos) {
0420     p1 = category.find('_', p1);
0421     p2 = category.find('_', p1 + 1);
0422     r9Min_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1));
0423     // If there is one value, just set lower bound
0424     if (p2 != std::string::npos) {
0425       p1 = p2;
0426       p2 = category.find('-', p1);
0427       r9Max_ = std::stof(category.substr(p1 + 1, p2 - p1 - 1));
0428       if (r9Max_ >= 1.0)
0429         r9Max_ = std::numeric_limits<float>::max();
0430     }
0431   }
0432   //------------------------------
0433   p1 = category.find("gainEle_");  // Position of first character
0434   if (p1 != std::string::npos) {
0435     p1 += 8;                      // Position of character after _
0436     p2 = category.find('-', p1);  // Position of - or end of string
0437     gain_ = std::stoul(category.substr(p1, p2 - p1), nullptr);
0438   }
0439   //so turns out the code does an inclusive X<=Y<=Z search for bins
0440   //which is what we want for run numbers
0441   //however then the problem is when we get a value exactly at the bin boundary
0442   //for the et/eta/r9 which then gives multiple bins
0443   //so we just decrement the maxValues ever so slightly to ensure that they are different
0444   //from the next bins min value
0445   etMax_ = std::nextafterf(etMax_, std::numeric_limits<float>::min());
0446   etaMax_ = std::nextafterf(etaMax_, std::numeric_limits<float>::min());
0447   r9Max_ = std::nextafterf(r9Max_, std::numeric_limits<float>::min());
0448 }
0449
0450 ///for the new file format
0451 EnergyScaleCorrection::CorrectionCategory::CorrectionCategory(unsigned int runMin,
0452                                                               unsigned int runMax,
0453                                                               float etaMin,
0454                                                               float etaMax,
0455                                                               float r9Min,
0456                                                               float r9Max,
0457                                                               float etMin,
0458                                                               float etMax,
0459                                                               unsigned int gainSeed)
0460     : runMin_(runMin),
0461       runMax_(runMax),
0462       etaMin_(etaMin),
0463       etaMax_(etaMax),
0464       r9Min_(r9Min),
0465       r9Max_(r9Max),
0466       etMin_(etMin),
0467       etMax_(etMax),
0468       gain_(gainSeed) {
0469   ///Same logic as the above constructor to avoid problems at the bin
0470   ///boundary of et/eta/R9 - just decrement the maxValues
0471   ///ever so slightly to ensure that they are different
0472   ///from the next bins min value
0473   etMax_ = std::nextafterf(etMax_, std::numeric_limits<float>::min());
0474   etaMax_ = std::nextafterf(etaMax_, std::numeric_limits<float>::min());
0475   r9Max_ = std::nextafterf(r9Max_, std::numeric_limits<float>::min());
0476 };
0477
0478 bool EnergyScaleCorrection::CorrectionCategory::inCategory(
0479     const unsigned int runnr, const float et, const float eta, const float r9, const unsigned int gainSeed) const {
0480   return runnr >= runMin_ && runnr <= runMax_ && et >= etMin_ && et <= etMax_ && eta >= etaMin_ && eta <= etaMax_ &&
0481          r9 >= r9Min_ && r9 <= r9Max_ && (gain_ == 0 || gainSeed == gain_);
0482 }
0483
0484 bool EnergyScaleCorrection::CorrectionCategory::operator<(const EnergyScaleCorrection::CorrectionCategory& b) const {
0485   if (runMin_ < b.runMin_ && runMax_ < b.runMax_)
0486     return true;
0487   if (runMax_ > b.runMax_ && runMin_ > b.runMin_)
0488     return false;
0489
0490   if (etaMin_ < b.etaMin_ && etaMax_ < b.etaMax_)
0491     return true;
0492   if (etaMax_ > b.etaMax_ && etaMin_ > b.etaMin_)
0493     return false;
0494
0495   if (r9Min_ < b.r9Min_ && r9Max_ < b.r9Max_)
0496     return true;
0497   if (r9Max_ > b.r9Max_ && r9Min_ > b.r9Min_)
0498     return false;
0499
0500   if (etMin_ < b.etMin_ && etMax_ < b.etMax_)
0501     return true;
0502   if (etMax_ > b.etMax_ && etMin_ > b.etMin_)
0503     return false;
0504
0505   if (gain_ == 0 || b.gain_ == 0)
0506     return false;  // if corrections are not categorized in gain then default gain value should always return false in order to have a match with the category
0507   if (gain_ < b.gain_)
0508     return true;
0509   else
0510     return false;
0511   return false;
0512 }
0513
0514 std::ostream& EnergyScaleCorrection::CorrectionCategory::print(std::ostream& os) const {
0515   os << runMin_ << " " << runMax_ << "\t" << etaMin_ << " " << etaMax_ << "\t" << r9Min_ << " " << r9Max_ << "\t"
0516      << etMin_ << " " << etMax_ << "\t" << gain_;
0517   return os;
0518 }