Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:24:01

0001 #include <map>
0002 #include <string>
0003 
0004 #include "TH1D.h"
0005 #include "TH2D.h"
0006 #include "TGraphErrors.h"
0007 
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0010 #include "FWCore/Utilities/interface/InputTag.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 
0016 #include "DataFormats/PatCandidates/interface/Jet.h"
0017 #include "PhysicsTools/PatUtils/interface/bJetSelector.h"
0018 #include "PhysicsTools/PatExamples/interface/BTagPerformance.h"
0019 #include "PhysicsTools/PatExamples/interface/PatBTagCommonHistos.h"
0020 
0021 class PatBTagAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0022 public:
0023   explicit PatBTagAnalyzer(const edm::ParameterSet &);
0024   ~PatBTagAnalyzer() override;
0025 
0026 private:
0027   void beginJob() override;
0028   void analyze(const edm::Event &, const edm::EventSetup &) override;
0029   void endJob() override;
0030 
0031   edm::EDGetTokenT<edm::View<pat::Jet> > jetToken_;
0032   edm::ParameterSet PatBjet_;
0033 
0034   std::string BTagpurity_;
0035   std::string BTagmethod_;
0036   std::string BTagdiscriminator_;
0037 
0038   bool BTagverbose;
0039   double BTagdisccut_;
0040   double BTagdiscmax_;
0041 
0042   std::string discname[10];
0043   std::string bname[10];
0044   std::string cname[10];
0045   BTagPerformance BTagPerf[10];
0046   std::map<int, std::string> udsgname;
0047 
0048   /// simple map to contain all histograms; histograms are booked in
0049   /// beginJob()
0050   std::map<std::string, TH1D *> histocontainer_;
0051   /// simple map to contain 2D  histograms; histograms are booked in
0052   /// beginJob()
0053   std::map<std::string, TH2D *> h2_;
0054   /// simple map to contain all graphs; graphs are booked in
0055   /// beginJob()
0056   std::map<std::string, TGraph *> graphcontainer_;
0057   /// simple map to contain all graphs; graphs are booked in
0058   /// beginJob()
0059   std::map<std::string, TGraphErrors *> grapherrorscontainer_;
0060 
0061   bJetSelector BTagger;
0062   PatBTagCommonHistos BTagHistograms;
0063 };
0064 
0065 PatBTagAnalyzer::PatBTagAnalyzer(const edm::ParameterSet &iConfig)
0066     : jetToken_(consumes<edm::View<pat::Jet> >(iConfig.getUntrackedParameter<edm::InputTag>("jetTag"))),
0067       PatBjet_(iConfig.getParameter<edm::ParameterSet>("BjetTag")),
0068       BTagpurity_(PatBjet_.getParameter<std::string>("purity")),
0069       BTagmethod_(PatBjet_.getUntrackedParameter<std::string>("tagger", "TC2")),
0070       BTagdiscriminator_(PatBjet_.getParameter<std::string>("discriminator")),
0071       BTagverbose(PatBjet_.getUntrackedParameter<bool>("verbose", false)),
0072       BTagdisccut_(PatBjet_.getUntrackedParameter<double>("mindiscriminatorcut", 5.0)),
0073       BTagdiscmax_(PatBjet_.getUntrackedParameter<double>("maxdiscriminatorcut", 15.0)),
0074       BTagger(iConfig.getParameter<edm::ParameterSet>("BJetOperatingPoints")),
0075       BTagHistograms(iConfig) {
0076   //now do what ever initialization is needed
0077   usesResource(TFileService::kSharedResource);
0078 }
0079 
0080 PatBTagAnalyzer::~PatBTagAnalyzer() {
0081   // do anything here that needs to be done at desctruction time
0082   // (e.g. close files, deallocate resources etc.)
0083 }
0084 
0085 // ------------ method called to for each event  ------------
0086 void PatBTagAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0087   using namespace edm;
0088 
0089   // first: get all objects from the event.
0090 
0091   edm::Handle<edm::View<pat::Jet> > jetHandle;
0092   iEvent.getByToken(jetToken_, jetHandle);
0093   edm::View<pat::Jet> jets = *jetHandle;  // get JETS
0094 
0095   // LOOP over all jets
0096 
0097   for (edm::View<pat::Jet>::const_iterator jet_iter = jets.begin(); jet_iter != jets.end(); ++jet_iter) {
0098     float bdiscriminator = jet_iter->bDiscriminator(BTagdiscriminator_);
0099     int flavor = jet_iter->partonFlavour();
0100     //
0101     // Fill in for performance standard pt(uncorrected) >10 and abs(eta)<2.4
0102     if (jet_iter->correctedJet("raw").pt() > 10 && fabs(jet_iter->eta()) < 2.4) {
0103       BTagPerf[0].Add(bdiscriminator, abs(flavor));
0104     }
0105 
0106     //Fill histograms
0107     BTagHistograms.Fill(jet_iter, "all");
0108     if (flavor == 0)
0109       BTagHistograms.Fill(jet_iter, "no_flavor");
0110     if (flavor == 5 || flavor == -5)
0111       BTagHistograms.Fill(jet_iter, "b");
0112     if (flavor == 4 || flavor == -4)
0113       BTagHistograms.Fill(jet_iter, "c");
0114     if ((-4 < flavor && flavor < 4 && flavor != 0) || (flavor == 21 || flavor == -21))
0115       BTagHistograms.Fill(jet_iter, "udsg");
0116 
0117   }  //end loop over jets
0118 }
0119 // ------------ method called once each job just before starting event loop  ------------
0120 void PatBTagAnalyzer::beginJob() {
0121   //
0122   // define some histograms using the framework tfileservice. Define the output file name in your .cfg.
0123   //
0124   edm::Service<TFileService> fs;
0125 
0126   TString suffix1 = "_test";
0127 
0128   //set performance variables collector
0129   for (int i = 0; i < 10; i++) {
0130     BTagPerf[i].Set(BTagmethod_);
0131     BTagPerf[i].SetMinDiscriminator(BTagdisccut_);
0132     BTagPerf[i].SetMaxDiscriminator(BTagdiscmax_);
0133   }
0134 
0135   histocontainer_["njets"] = fs->make<TH1D>("njets", "jet multiplicity for jets with p_{T} > 50 GeV/c", 10, 0, 10);
0136   // Std. 30 pt uncorr cut for performance
0137   discname[0] = "disc" + BTagmethod_ + "_udsg";
0138   bname[0] = "g" + BTagmethod_ + "_b";
0139   cname[0] = "g" + BTagmethod_ + "_c";
0140   udsgname[0] = "g" + BTagmethod_ + "_udsg";
0141 
0142   // 10 pt uncorr for performance + all,>0,>1,>2 tracks
0143   discname[1] = "Uncor10_disc" + BTagmethod_ + "_udsg";
0144   bname[1] = "Uncor10_g" + BTagmethod_ + "_b";
0145   cname[1] = "Uncor10_g" + BTagmethod_ + "_c";
0146   udsgname[1] = "Uncor10_g" + BTagmethod_ + "_udsg";
0147   discname[2] = "Uncor10t0_disc" + BTagmethod_ + "_udsg";
0148   bname[2] = "Uncor10t0_g" + BTagmethod_ + "_b";
0149   cname[2] = "Uncor10t0_g" + BTagmethod_ + "_c";
0150   udsgname[2] = "Uncor10t0_g" + BTagmethod_ + "_udsg";
0151   discname[3] = "Uncor10t1_disc" + BTagmethod_ + "_udsg";
0152   bname[3] = "Uncor10t1_g" + BTagmethod_ + "_b";
0153   cname[3] = "Uncor10t1_g" + BTagmethod_ + "_c";
0154   udsgname[3] = "Uncor10t1_g" + BTagmethod_ + "_udsg";
0155   discname[4] = "Uncor10t2_disc" + BTagmethod_ + "_udsg";
0156   bname[4] = "Uncor10t2_g" + BTagmethod_ + "_b";
0157   cname[4] = "Uncor10t2_g" + BTagmethod_ + "_c";
0158   udsgname[4] = "Uncor10t2_g" + BTagmethod_ + "_udsg";
0159 
0160   // 30 pt corr for performance + all,>0,>1,>2 tracks
0161   discname[5] = "Corr30_disc" + BTagmethod_ + "_udsg";
0162   bname[5] = "Corr30_g" + BTagmethod_ + "_b";
0163   cname[5] = "Corr30_g" + BTagmethod_ + "_c";
0164   udsgname[5] = "Corr30_g" + BTagmethod_ + "_udsg";
0165   discname[6] = "Corr30t0_disc" + BTagmethod_ + "_udsg";
0166   bname[6] = "Corr30t0_g" + BTagmethod_ + "_b";
0167   cname[6] = "Corr30t0_g" + BTagmethod_ + "_c";
0168   udsgname[6] = "Corr30t0_g" + BTagmethod_ + "_udsg";
0169   discname[7] = "Corr30t1_disc" + BTagmethod_ + "_udsg";
0170   bname[7] = "Corr30t1_g" + BTagmethod_ + "_b";
0171   cname[7] = "Corr30t1_g" + BTagmethod_ + "_c";
0172   udsgname[7] = "Corr30t1_g" + BTagmethod_ + "_udsg";
0173   discname[8] = "Corr30t2_disc" + BTagmethod_ + "_udsg";
0174   bname[8] = "Corr30t2_g" + BTagmethod_ + "_b";
0175   cname[8] = "Corr30t2_g" + BTagmethod_ + "_c";
0176   udsgname[8] = "Corr30t2_g" + BTagmethod_ + "_udsg";
0177 
0178   // check filter
0179   discname[9] = "check_disc" + BTagmethod_ + "_udsg";
0180   bname[9] = "check_g" + BTagmethod_ + "_b";
0181   cname[9] = "check_g" + BTagmethod_ + "_c";
0182   udsgname[9] = "check_g" + BTagmethod_ + "_udsg";
0183 
0184   for (int i = 1; i < 10; i++) {
0185     graphcontainer_[discname[i]] = fs->make<TGraph>(BTagPerf[i].GetN());
0186     graphcontainer_[discname[i]]->SetName(discname[i].c_str());
0187     grapherrorscontainer_[bname[i]] = fs->make<TGraphErrors>(BTagPerf[i].GetN());
0188     grapherrorscontainer_[bname[i]]->SetName(bname[i].c_str());
0189     grapherrorscontainer_[cname[i]] = fs->make<TGraphErrors>(BTagPerf[i].GetN());
0190     grapherrorscontainer_[cname[i]]->SetName(cname[i].c_str());
0191     grapherrorscontainer_[udsgname[i]] = fs->make<TGraphErrors>(BTagPerf[i].GetN());
0192     grapherrorscontainer_[udsgname[i]]->SetName(udsgname[i].c_str());
0193   }
0194   //Define histograms
0195   BTagHistograms.Set("all");
0196   BTagHistograms.Set("no_flavor");
0197   BTagHistograms.Set("b");
0198   BTagHistograms.Set("c");
0199   BTagHistograms.Set("udsg");
0200 
0201   // Set to save histogram errors
0202   BTagHistograms.Sumw2();
0203 }
0204 
0205 // ------------ method called once each job just after ending the event loop  ------------
0206 void PatBTagAnalyzer::endJob() {
0207   //ed++
0208   edm::Service<TFileService> fs;
0209 
0210   // Save performance plots as Tgraphs
0211 
0212   for (int i = 1; i < 10; i++) {
0213     BTagPerf[i].Eval();
0214     for (int n = 0; n < BTagPerf[i].GetN(); n++) {
0215       graphcontainer_[discname[i]]->SetPoint(
0216           n, BTagPerf[i].GetArray("udsg")[n], BTagPerf[i].GetArray("discriminator")[n]);
0217       grapherrorscontainer_[bname[i]]->SetPoint(n, BTagPerf[i].GetArray("b")[n], BTagPerf[i].GetArray("b")[n]);
0218       grapherrorscontainer_[cname[i]]->SetPoint(n, BTagPerf[i].GetArray("b")[n], BTagPerf[i].GetArray("c")[n]);
0219       grapherrorscontainer_[udsgname[i]]->SetPoint(n, BTagPerf[i].GetArray("b")[n], BTagPerf[i].GetArray("udsg")[n]);
0220       grapherrorscontainer_[bname[i]]->SetPointError(
0221           n, BTagPerf[i].GetArray("bErr")[n], BTagPerf[i].GetArray("bErr")[n]);
0222       grapherrorscontainer_[cname[i]]->SetPointError(
0223           n, BTagPerf[i].GetArray("bErr")[n], BTagPerf[i].GetArray("cErr")[n]);
0224       grapherrorscontainer_[udsgname[i]]->SetPointError(
0225           n, BTagPerf[i].GetArray("bErr")[n], BTagPerf[i].GetArray("udsgErr")[n]);
0226     }  //end for over BTagPerf[i] elements
0227     graphcontainer_[discname[i]]->SetTitle("Jet udsg-mistagging");
0228     grapherrorscontainer_[bname[i]]->SetTitle("Jet b-efficiency");
0229     grapherrorscontainer_[cname[i]]->SetTitle("Jet c-mistagging");
0230     grapherrorscontainer_[udsgname[i]]->SetTitle("discriminator vs udsg-mistagging");
0231   }  //end for over [i]
0232 
0233   // Save default cut performance plot
0234   BTagPerf[0].Eval();
0235 
0236   //  TFileDirectory TaggerDir = fs->mkdir(BTagmethod_);
0237   //  TGraphErrors *BTagger_b    = new TGraphErrors(BTagTool.GetN(),
0238   TGraphErrors *BTagger_b = fs->mkdir(BTagmethod_)
0239                                 .make<TGraphErrors>(BTagPerf[0].GetN(),
0240                                                     BTagPerf[0].GetArray("b").GetArray(),
0241                                                     BTagPerf[0].GetArray("b").GetArray(),
0242                                                     BTagPerf[0].GetArray("bErr").GetArray(),
0243                                                     BTagPerf[0].GetArray("bErr").GetArray());
0244 
0245   TGraphErrors *BTagger_c = new TGraphErrors(BTagPerf[0].GetN(),
0246                                              BTagPerf[0].GetArray("b").GetArray(),
0247                                              BTagPerf[0].GetArray("c").GetArray(),
0248                                              BTagPerf[0].GetArray("bErr").GetArray(),
0249                                              BTagPerf[0].GetArray("cErr").GetArray());
0250 
0251   TGraphErrors *BTagger_udsg = new TGraphErrors(BTagPerf[0].GetN(),
0252                                                 BTagPerf[0].GetArray("b").GetArray(),
0253                                                 BTagPerf[0].GetArray("udsg").GetArray(),
0254                                                 BTagPerf[0].GetArray("bErr").GetArray(),
0255                                                 BTagPerf[0].GetArray("udsgErr").GetArray());
0256   TGraph *discBTagger_udsg = new TGraph(
0257       BTagPerf[0].GetN(), BTagPerf[0].GetArray("udsg").GetArray(), BTagPerf[0].GetArray("discriminator").GetArray());
0258 
0259   BTagger_b->SetName(bname[0].c_str());
0260   BTagger_c->SetName(cname[0].c_str());
0261   BTagger_udsg->SetName(udsgname[0].c_str());
0262   discBTagger_udsg->SetName(discname[0].c_str());
0263 
0264   BTagger_b->SetTitle("Jet b-efficiency");
0265   BTagger_c->SetTitle("Jet c-mistagging");
0266   BTagger_udsg->SetTitle("Jet udsg-mistagging");
0267   discBTagger_udsg->SetTitle("discriminator vs udsg-mistagging");
0268 
0269   for (int i = 1; i < 10; i++) {
0270     graphcontainer_[discname[i]]->Write();
0271     grapherrorscontainer_[bname[i]]->Write();
0272     grapherrorscontainer_[cname[i]]->Write();
0273     grapherrorscontainer_[udsgname[i]]->Write();
0274   }
0275 
0276   BTagger_b->Write();
0277   BTagger_c->Write();
0278   BTagger_udsg->Write();
0279   discBTagger_udsg->Write();
0280 }
0281 
0282 //define this as a plug-in
0283 DEFINE_FWK_MODULE(PatBTagAnalyzer);