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