NanoAOD/plugins/BTagSFProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Original Author:  Spandan Mondal
0004 //         Created:  Tue, 13 Mar 2018 09:26:52 GMT
0005 //
0006 //
0007 // NOTE: this file is deprecated and only kept for compatibility when running the
0008 // frozen nanoAOD V10 config
0009
0010 // system include files
0011 #include <memory>
0012
0013 // user include files
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/stream/EDProducer.h"
0016
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/MakerMacros.h"
0019
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Utilities/interface/StreamID.h"
0022 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0024
0025 #include "DataFormats/PatCandidates/interface/Jet.h"
0026 #include "DataFormats/NanoAOD/interface/FlatTable.h"
0027
0028 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0029 #include "DataFormats/Common/interface/ValueMap.h"
0030
0031 #include "CondFormats/BTauObjects/interface/BTagCalibration.h"
0032 #include "CondTools/BTau/interface/BTagCalibrationReader.h"
0033
0034 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0035
0036 #include <vector>
0037 #include <string>
0038
0039 class BTagSFProducer : public edm::stream::EDProducer<> {
0040 public:
0041   BTagSFProducer(const edm::ParameterSet& iConfig)
0042       : src_(consumes<std::vector<pat::Jet>>(iConfig.getParameter<edm::InputTag>("src"))),
0043         cut_(iConfig.getParameter<std::string>("cut")),
0044         discNames_(iConfig.getParameter<std::vector<std::string>>("discNames")),
0045         discShortNames_(iConfig.getParameter<std::vector<std::string>>("discShortNames")),
0046         weightFiles_(iConfig.getParameter<std::vector<std::string>>("weightFiles")),
0047         operatingPoints_(iConfig.getParameter<std::vector<std::string>>("operatingPoints")),
0048         measurementTypesB_(iConfig.getParameter<std::vector<std::string>>("measurementTypesB")),
0049         measurementTypesC_(iConfig.getParameter<std::vector<std::string>>("measurementTypesC")),
0050         measurementTypesUDSG_(iConfig.getParameter<std::vector<std::string>>("measurementTypesUDSG")),
0051         sysTypes_(iConfig.getParameter<std::vector<std::string>>("sysTypes")) {
0052     produces<nanoaod::FlatTable>();
0053
0054     nDiscs = discNames_.size();
0055     assert(discShortNames_.size() == nDiscs && weightFiles_.size() == nDiscs && operatingPoints_.size() == nDiscs &&
0056            measurementTypesB_.size() == nDiscs && measurementTypesC_.size() == nDiscs &&
0057            measurementTypesUDSG_.size() == nDiscs && sysTypes_.size() == nDiscs);
0058
0059     bool validate = iConfig.getUntrackedParameter<bool>("validate");
0060     for (unsigned int iDisc = 0; iDisc < nDiscs; ++iDisc) {
0061       if (weightFiles_[iDisc] != "unavailable") {
0062         // setup calibration
0063         BTagCalibration calib;
0064         edm::FileInPath fip(weightFiles_[iDisc]);
0065         calib = BTagCalibration(discShortNames_[iDisc], fip.fullPath(), validate);
0066
0067         // determine op
0068         std::string opname;
0069         if (operatingPoints_[iDisc] == "0" || operatingPoints_[iDisc] == "loose") {
0070           op = BTagEntry::OP_LOOSE;
0071           opname = "loose";
0072         } else if (operatingPoints_[iDisc] == "1" || operatingPoints_[iDisc] == "medium") {
0073           op = BTagEntry::OP_MEDIUM;
0074           opname = "medium";
0075         } else if (operatingPoints_[iDisc] == "2" || operatingPoints_[iDisc] == "tight") {
0076           op = BTagEntry::OP_TIGHT;
0077           opname = "tight";
0078         } else if (operatingPoints_[iDisc] == "3" || operatingPoints_[iDisc] == "reshaping") {
0079           op = BTagEntry::OP_RESHAPING;
0080           opname = "discriminator reshaping";
0081         }
0082
0083         // setup reader
0084         BTagCalibrationReader reader;
0085         reader = BTagCalibrationReader(op, sysTypes_[iDisc]);
0086         reader.load(calib, BTagEntry::FLAV_B, measurementTypesB_[iDisc]);
0087         reader.load(calib, BTagEntry::FLAV_C, measurementTypesC_[iDisc]);
0088         reader.load(calib, BTagEntry::FLAV_UDSG, measurementTypesUDSG_[iDisc]);
0089
0090         //calibs.push_back(calib);
0091         readers.push_back(reader);
0092
0093         // report
0094         edm::LogInfo("BTagSFProducer") << "Loaded " + discShortNames_[iDisc] + " SFs from weight file " +
0095                                               weightFiles_[iDisc] + " with\noperating point: " + opname +
0096                                               ",\nmeasurement type: B=" + measurementTypesB_[iDisc] +
0097                                               ", C=" + measurementTypesC_[iDisc] +
0098                                               ", UDSG=" + measurementTypesUDSG_[iDisc] +
0099                                               ",\nsystematic type: " + sysTypes_[iDisc] + ".\n"
0100                                        << std::endl;
0101
0102         // find if multiple MiniAOD branches need to be summed up (e.g., DeepCSV b+bb) and separate them using '+' delimiter from config
0103         std::stringstream dName(discNames_[iDisc]);
0104         std::string branch;
0105         std::vector<std::string> branches;
0106         while (std::getline(dName, branch, '+')) {
0107           branches.push_back(branch);
0108         }
0109         inBranchNames.push_back(branches);
0110       } else {
0111         //BTagCalibration calib;
0112         BTagCalibrationReader reader;
0113         //calibs.push_back(calib);            //dummy, so that index of vectors still match
0114         readers.push_back(reader);  //dummy, so that index of vectors still match
0115         std::vector<std::string> branches;
0116         branches.push_back("");
0117         inBranchNames.push_back(branches);
0118
0119         // report
0120         edm::LogWarning("BTagSFProducer")
0121             << "Skipped loading BTagCalibration for " + discShortNames_[iDisc] +
0122                    " as it was marked as unavailable in the configuration file. Event weights will not be stored.\n"
0123             << std::endl;
0124       }
0125     }
0126   }
0127
0128   ~BTagSFProducer() override{};
0129
0130   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0131     edm::ParameterSetDescription desc;
0132
0133     desc.add<edm::InputTag>("src")->setComment("input AK4 jet collection");
0134     desc.add<std::string>("cut")->setComment("minimum pT and maximum eta cuts for jets");
0135     desc.add<std::vector<std::string>>("discNames")->setComment("name of b-tag discriminator branch in MiniAOD");
0136     desc.add<std::vector<std::string>>("discShortNames")->setComment("common name of discriminator");
0137     desc.add<std::vector<std::string>>("weightFiles")->setComment("path to the .csv file containing the SFs");
0138     desc.add<std::vector<std::string>>("operatingPoints")
0139         ->setComment("loose = 0, medium = 1, tight = 2, disriminator reshaping = 3");
0140     desc.add<std::vector<std::string>>("measurementTypesB")
0141         ->setComment("e.g. \"ttbar\", \"comb\", \"incl\", \"iterativefit\" for b jets");
0142     desc.add<std::vector<std::string>>("measurementTypesC")
0143         ->setComment("e.g. \"ttbar\", \"comb\", \"incl\", \"iterativefit\" for c jets");
0144     desc.add<std::vector<std::string>>("measurementTypesUDSG")
0145         ->setComment("e.g. \"ttbar\", \"comb\", \"incl\", \"iterativefit\" for light jets");
0146     desc.add<std::vector<std::string>>("sysTypes")
0147         ->setComment(
0148             "\"up\", \"central\", \"down\", but arbitrary strings possible, like \"up_generator\" or \"up_jec\"");
0149     desc.addUntracked<bool>("validate", false)->setComment("validate the function expressions in the weightFiles");
0150     descriptions.add("BTagWeightTable", desc);
0151   }
0152
0153 private:
0154   void produce(edm::Event&, const edm::EventSetup&) override;
0155
0156   edm::EDGetTokenT<std::vector<pat::Jet>> src_;
0157   const StringCutObjectSelector<pat::Jet> cut_;
0158
0159   std::vector<std::string> discNames_;
0160   std::vector<std::string> discShortNames_;
0161   std::vector<std::string> weightFiles_;
0162   std::vector<std::string> operatingPoints_;
0163   std::vector<std::string> measurementTypesB_;
0164   std::vector<std::string> measurementTypesC_;
0165   std::vector<std::string> measurementTypesUDSG_;
0166   std::vector<std::string> sysTypes_;
0167
0168   BTagEntry::OperatingPoint op;
0169   std::vector<std::vector<std::string>> inBranchNames;
0170   //std::vector<BTagCalibration> calibs;
0171   std::vector<BTagCalibrationReader> readers;
0172   unsigned int nDiscs;
0173 };
0174
0175 // ------------ method called to produce the data  ------------
0176 void BTagSFProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0177   using namespace edm;
0178   using namespace std;
0179
0180   Handle<pat::JetCollection> jets;
0181   iEvent.getByToken(src_, jets);
0182
0183   double pt;
0184   double eta;
0185   int flavour;
0186   double bdisc;
0187   double SF;
0188
0189   double EventWt;
0190
0191   auto out = std::make_unique<nanoaod::FlatTable>(1, "btagWeight", true);
0192   out->setDoc("b-tagging event weights");
0193
0194   for (unsigned int iDisc = 0; iDisc < nDiscs; ++iDisc) {  // loop over b-tagging algorithms
0195
0196     if (weightFiles_[iDisc] != "unavailable") {
0197       EventWt = 1.;
0198       for (const pat::Jet& jet : *jets) {  // loop over jets and accumulate product of SF for each jet
0199         pt = jet.pt();
0200         eta = jet.eta();
0201         bdisc = 0.;
0202
0203         if (op == BTagEntry::OP_RESHAPING) {
0204           for (const string& inBranch :
0205                inBranchNames[iDisc]) {  //sum up the discriminator values if multiple, e.g. DeepCSV b+bb
0206             bdisc += jet.bDiscriminator(inBranch);
0207           }
0208         }
0209
0210         flavour = jet.hadronFlavour();
0211
0212         if (cut_(jet)) {             //multiply SF of only the jets that pass the cut
0213           if (fabs(flavour) == 5) {  // b jets
0214             SF = readers[iDisc].eval_auto_bounds(sysTypes_[iDisc], BTagEntry::FLAV_B, eta, pt, bdisc);
0215           } else if (fabs(flavour) == 4) {  // c jets
0216             SF = readers[iDisc].eval_auto_bounds(sysTypes_[iDisc], BTagEntry::FLAV_C, eta, pt, bdisc);
0217           } else {  // others
0218             SF = readers[iDisc].eval_auto_bounds(sysTypes_[iDisc], BTagEntry::FLAV_UDSG, eta, pt, bdisc);
0219           }
0220         } else {
0221           SF = 1.;
0222         }
0223
0224         if (SF == 0.) {  // default value of SF is set to 1 in case BTagCalibration returns 0
0225           //no need to log this as could be pretty common, leaving the cout commented in case this is needed by the author for simple debugging
0226           //cout << discShortNames_[iDisc]+" SF not found for jet with pT="+to_string(pt)+", eta="+to_string(eta)+", discValue="+to_string(bdisc)+", flavour="+to_string(flavour) +". Setting SF to 1." << endl;
0227           SF = 1.;
0228         }
0229
0230         EventWt *= SF;
0231       }
0232
0233       out->addColumnValue<float>(discShortNames_[iDisc], EventWt, "b-tag event weight for " + discShortNames_[iDisc]);
0234     }
0235   }
0236
0237   iEvent.put(move(out));
0238 }
0239
0240 //define this as a plug-in
0241 DEFINE_FWK_MODULE(BTagSFProducer);