BTau/src/BTagCalibrationReader.cc

0001 #include "CondTools/BTau/interface/BTagCalibrationReader.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003
0004 class BTagCalibrationReader::BTagCalibrationReaderImpl {
0005   friend class BTagCalibrationReader;
0006
0007 public:
0008   struct TmpEntry {
0009     float etaMin;
0010     float etaMax;
0011     float ptMin;
0012     float ptMax;
0013     float discrMin;
0014     float discrMax;
0015     TF1 func;
0016   };
0017
0018 private:
0019   BTagCalibrationReaderImpl(BTagEntry::OperatingPoint op,
0020                             const std::string &sysType,
0021                             const std::vector<std::string> &otherSysTypes = {});
0022
0023   void load(const BTagCalibration &c, BTagEntry::JetFlavor jf, std::string measurementType);
0024
0025   double eval(BTagEntry::JetFlavor jf, float eta, float pt, float discr) const;
0026
0027   double eval_auto_bounds(const std::string &sys, BTagEntry::JetFlavor jf, float eta, float pt, float discr) const;
0028
0029   std::pair<float, float> min_max_pt(BTagEntry::JetFlavor jf, float eta, float discr) const;
0030
0031   std::pair<float, float> min_max_eta(BTagEntry::JetFlavor jf, float discr) const;
0032
0033   BTagEntry::OperatingPoint op_;
0034   std::string sysType_;
0035   std::vector<std::vector<TmpEntry>> tmpData_;  // first index: jetFlavor
0036   std::vector<bool> useAbsEta_;                 // first index: jetFlavor
0037   std::map<std::string, std::shared_ptr<BTagCalibrationReaderImpl>> otherSysTypeReaders_;
0038 };
0039
0040 BTagCalibrationReader::BTagCalibrationReaderImpl::BTagCalibrationReaderImpl(
0041     BTagEntry::OperatingPoint op, const std::string &sysType, const std::vector<std::string> &otherSysTypes)
0042     : op_(op), sysType_(sysType), tmpData_(3), useAbsEta_(3, true) {
0043   for (const std::string &ost : otherSysTypes) {
0044     if (otherSysTypeReaders_.count(ost)) {
0045       throw cms::Exception("BTagCalibrationReader")
0046           << "Every otherSysType should only be given once. Duplicate: " << ost;
0047     }
0048     otherSysTypeReaders_[ost] = std::unique_ptr<BTagCalibrationReaderImpl>(new BTagCalibrationReaderImpl(op, ost));
0049   }
0050 }
0051
0052 void BTagCalibrationReader::BTagCalibrationReaderImpl::load(const BTagCalibration &c,
0053                                                             BTagEntry::JetFlavor jf,
0054                                                             std::string measurementType) {
0055   if (!tmpData_[jf].empty()) {
0056     throw cms::Exception("BTagCalibrationReader") << "Data for this jet-flavor is already loaded: " << jf;
0057   }
0058
0059   BTagEntry::Parameters params(op_, measurementType, sysType_);
0060   const std::vector<BTagEntry> &entries = c.getEntries(params);
0061
0062   for (const auto &be : entries) {
0063     if (be.params.jetFlavor != jf) {
0064       continue;
0065     }
0066
0067     TmpEntry te;
0068     te.etaMin = be.params.etaMin;
0069     te.etaMax = be.params.etaMax;
0070     te.ptMin = be.params.ptMin;
0071     te.ptMax = be.params.ptMax;
0072     te.discrMin = be.params.discrMin;
0073     te.discrMax = be.params.discrMax;
0074
0075     if (op_ == BTagEntry::OP_RESHAPING) {
0076       te.func = TF1("", be.formula.c_str(), be.params.discrMin, be.params.discrMax);
0077     } else {
0078       te.func = TF1("", be.formula.c_str(), be.params.ptMin, be.params.ptMax);
0079     }
0080
0081     tmpData_[be.params.jetFlavor].push_back(te);
0082     if (te.etaMin < 0) {
0083       useAbsEta_[be.params.jetFlavor] = false;
0084     }
0085   }
0086
0087   for (auto &p : otherSysTypeReaders_) {
0088     p.second->load(c, jf, measurementType);
0089   }
0090 }
0091
0092 double BTagCalibrationReader::BTagCalibrationReaderImpl::eval(BTagEntry::JetFlavor jf,
0093                                                               float eta,
0094                                                               float pt,
0095                                                               float discr) const {
0096   bool use_discr = (op_ == BTagEntry::OP_RESHAPING);
0097   if (useAbsEta_[jf] && eta < 0) {
0098     eta = -eta;
0099   }
0100
0101   // search linearly through eta, pt and discr ranges and eval
0102   // future: find some clever data structure based on intervals
0103   const auto &entries = tmpData_.at(jf);
0104   for (unsigned i = 0; i < entries.size(); ++i) {
0105     const auto &e = entries.at(i);
0106     if (e.etaMin <= eta && eta <= e.etaMax  // find eta
0107         && e.ptMin < pt && pt <= e.ptMax    // check pt
0108     ) {
0109       if (use_discr) {                                    // discr. reshaping?
0110         if (e.discrMin <= discr && discr < e.discrMax) {  // check discr
0111           return e.func.Eval(discr);
0112         }
0113       } else {
0114         return e.func.Eval(pt);
0115       }
0116     }
0117   }
0118
0119   return 0.;  // default value
0120 }
0121
0122 double BTagCalibrationReader::BTagCalibrationReaderImpl::eval_auto_bounds(
0123     const std::string &sys, BTagEntry::JetFlavor jf, float eta, float pt, float discr) const {
0124   auto sf_bounds_eta = min_max_eta(jf, discr);
0125   bool eta_is_out_of_bounds = false;
0126
0127   if (sf_bounds_eta.first < 0)
0128     sf_bounds_eta.first = -sf_bounds_eta.second;
0129
0130   if (useAbsEta_[jf] && eta < 0) {
0131     eta = -eta;
0132   }
0133
0134   if (eta <= sf_bounds_eta.first || eta > sf_bounds_eta.second) {
0135     eta_is_out_of_bounds = true;
0136   }
0137
0138   if (eta_is_out_of_bounds) {
0139     return 1.;
0140   }
0141
0142   auto sf_bounds = min_max_pt(jf, eta, discr);
0143   float pt_for_eval = pt;
0144   bool is_out_of_bounds = false;
0145
0146   if (pt <= sf_bounds.first) {
0147     pt_for_eval = sf_bounds.first + .0001;
0148     is_out_of_bounds = true;
0149   } else if (pt > sf_bounds.second) {
0150     pt_for_eval = sf_bounds.second - .0001;
0151     is_out_of_bounds = true;
0152   }
0153
0154   // get central SF (and maybe return)
0155   double sf = eval(jf, eta, pt_for_eval, discr);
0156   if (sys == sysType_) {
0157     return sf;
0158   }
0159
0160   // get sys SF (and maybe return)
0161   if (!otherSysTypeReaders_.count(sys)) {
0162     throw cms::Exception("BTagCalibrationReader") << "sysType not available (maybe not loaded?): " << sys;
0163   }
0164   double sf_err = otherSysTypeReaders_.at(sys)->eval(jf, eta, pt_for_eval, discr);
0165   if (!is_out_of_bounds) {
0166     return sf_err;
0167   }
0168
0169   // double uncertainty on out-of-bounds and return
0170   sf_err = sf + 2 * (sf_err - sf);
0171   return sf_err;
0172 }
0173
0174 std::pair<float, float> BTagCalibrationReader::BTagCalibrationReaderImpl::min_max_pt(BTagEntry::JetFlavor jf,
0175                                                                                      float eta,
0176                                                                                      float discr) const {
0177   bool use_discr = (op_ == BTagEntry::OP_RESHAPING);
0178   if (useAbsEta_[jf] && eta < 0) {
0179     eta = -eta;
0180   }
0181
0182   const auto &entries = tmpData_.at(jf);
0183   float min_pt = -1., max_pt = -1.;
0184   for (const auto &e : entries) {
0185     if (e.etaMin <= eta && eta <= e.etaMax  // find eta
0186     ) {
0187       if (min_pt < 0.) {  // init
0188         min_pt = e.ptMin;
0189         max_pt = e.ptMax;
0190         continue;
0191       }
0192
0193       if (use_discr) {                                    // discr. reshaping?
0194         if (e.discrMin <= discr && discr < e.discrMax) {  // check discr
0195           min_pt = min_pt < e.ptMin ? min_pt : e.ptMin;
0196           max_pt = max_pt > e.ptMax ? max_pt : e.ptMax;
0197         }
0198       } else {
0199         min_pt = min_pt < e.ptMin ? min_pt : e.ptMin;
0200         max_pt = max_pt > e.ptMax ? max_pt : e.ptMax;
0201       }
0202     }
0203   }
0204
0205   return std::make_pair(min_pt, max_pt);
0206 }
0207
0208 std::pair<float, float> BTagCalibrationReader::BTagCalibrationReaderImpl::min_max_eta(BTagEntry::JetFlavor jf,
0209                                                                                       float discr) const {
0210   bool use_discr = (op_ == BTagEntry::OP_RESHAPING);
0211
0212   const auto &entries = tmpData_.at(jf);
0213   float min_eta = 0., max_eta = 0.;
0214   for (const auto &e : entries) {
0215     if (use_discr) {                                    // discr. reshaping?
0216       if (e.discrMin <= discr && discr < e.discrMax) {  // check discr
0217         min_eta = min_eta < e.etaMin ? min_eta : e.etaMin;
0218         max_eta = max_eta > e.etaMax ? max_eta : e.etaMax;
0219       }
0220     } else {
0221       min_eta = min_eta < e.etaMin ? min_eta : e.etaMin;
0222       max_eta = max_eta > e.etaMax ? max_eta : e.etaMax;
0223     }
0224   }
0225
0226   return std::make_pair(min_eta, max_eta);
0227 }
0228
0229 BTagCalibrationReader::BTagCalibrationReader(BTagEntry::OperatingPoint op,
0230                                              const std::string &sysType,
0231                                              const std::vector<std::string> &otherSysTypes)
0232     : pimpl(new BTagCalibrationReaderImpl(op, sysType, otherSysTypes)) {}
0233
0234 void BTagCalibrationReader::load(const BTagCalibration &c,
0235                                  BTagEntry::JetFlavor jf,
0236                                  const std::string &measurementType) {
0237   pimpl->load(c, jf, measurementType);
0238 }
0239
0240 double BTagCalibrationReader::eval(BTagEntry::JetFlavor jf, float eta, float pt, float discr) const {
0241   return pimpl->eval(jf, eta, pt, discr);
0242 }
0243
0244 double BTagCalibrationReader::eval_auto_bounds(
0245     const std::string &sys, BTagEntry::JetFlavor jf, float eta, float pt, float discr) const {
0246   return pimpl->eval_auto_bounds(sys, jf, eta, pt, discr);
0247 }
0248
0249 std::pair<float, float> BTagCalibrationReader::min_max_pt(BTagEntry::JetFlavor jf, float eta, float discr) const {
0250   return pimpl->min_max_pt(jf, eta, discr);
0251 }