PatExamples/plugins/PatBJetTagAnalyzer.cc

0001 #include <iostream>
0002 #include <cmath>
0003 #include <vector>
0004 #include <string>
0005
0006 #include <TH1.h>
0007
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015
0016 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0017
0018 #include "DataFormats/PatCandidates/interface/Jet.h"
0019
0020 class PatBJetTagAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0021 public:
0022   /// constructor and destructor
0023   PatBJetTagAnalyzer(const edm::ParameterSet &params);
0024   ~PatBJetTagAnalyzer() override;
0025
0026   // virtual methods called from base class EDAnalyzer
0027   void beginJob() override;
0028   void analyze(const edm::Event &event, const edm::EventSetup &es) override;
0029
0030 private:
0031   // configuration parameters
0032   edm::EDGetTokenT<pat::JetCollection> jetsToken_;
0033
0034   double jetPtCut_;   // minimum (uncorrected) jet energy
0035   double jetEtaCut_;  // maximum |eta| for jet
0036
0037   // jet flavour constants
0038
0039   enum Flavour { ALL_JETS = 0, UDSG_JETS, C_JETS, B_JETS, NONID_JETS, N_JET_TYPES };
0040
0041   TH1 *flavours_;
0042
0043   // one group of plots per jet flavour;
0044   struct Plots {
0045     TH1 *discrTC, *discrSSV, *discrCSV;
0046   } plots_[N_JET_TYPES];
0047 };
0048
0049 PatBJetTagAnalyzer::PatBJetTagAnalyzer(const edm::ParameterSet &params)
0050     : jetsToken_(consumes<pat::JetCollection>(params.getParameter<edm::InputTag>("jets"))),
0051       jetPtCut_(params.getParameter<double>("jetPtCut")),
0052       jetEtaCut_(params.getParameter<double>("jetEtaCut")) {
0053   usesResource(TFileService::kSharedResource);
0054 }
0055
0056 PatBJetTagAnalyzer::~PatBJetTagAnalyzer() {}
0057
0058 void PatBJetTagAnalyzer::beginJob() {
0059   // retrieve handle to auxiliary service
0060   //  used for storing histograms into ROOT file
0061   edm::Service<TFileService> fs;
0062
0063   flavours_ = fs->make<TH1F>("flavours", "jet flavours", 5, 0, 5);
0064
0065   // book histograms for all jet flavours
0066   for (unsigned int i = 0; i < N_JET_TYPES; i++) {
0067     Plots &plots = plots_[i];
0068     const char *flavour, *name;
0069
0070     switch ((Flavour)i) {
0071       case ALL_JETS:
0072         flavour = "all jets";
0073         name = "all";
0074         break;
0075       case UDSG_JETS:
0076         flavour = "light flavour jets";
0077         name = "udsg";
0078         break;
0079       case C_JETS:
0080         flavour = "charm jets";
0081         name = "c";
0082         break;
0083       case B_JETS:
0084         flavour = "bottom jets";
0085         name = "b";
0086         break;
0087       default:
0088         flavour = "unidentified jets";
0089         name = "ni";
0090         break;
0091     }
0092
0093     plots.discrTC = fs->make<TH1F>(
0094         Form("discrTC_%s", name), Form("track counting (\"high efficiency\") in %s", flavour), 100, 0, 20);
0095     plots.discrSSV =
0096         fs->make<TH1F>(Form("discrSSV_%s", name), Form("simple secondary vertex in %s", flavour), 100, 0, 10);
0097     plots.discrCSV =
0098         fs->make<TH1F>(Form("discrCSV_%s", name), Form("combined secondary vertex in %s", flavour), 100, 0, 1);
0099   }
0100 }
0101
0102 void PatBJetTagAnalyzer::analyze(const edm::Event &event, const edm::EventSetup &es) {
0103   // handle to the jets collection
0104   edm::Handle<pat::JetCollection> jetsHandle;
0105   event.getByToken(jetsToken_, jetsHandle);
0106
0107   // now go through all jets
0108   for (pat::JetCollection::const_iterator jet = jetsHandle->begin(); jet != jetsHandle->end(); ++jet) {
0109     // only look at jets that pass the pt and eta cut
0110     if (jet->pt() < jetPtCut_ || std::abs(jet->eta()) > jetEtaCut_)
0111       continue;
0112
0113     Flavour flavour;
0114     // find out the jet flavour (differs between quark and anti-quark)
0115     switch (std::abs(jet->partonFlavour())) {
0116       case 1:
0117       case 2:
0118       case 3:
0119       case 21:
0120         flavour = UDSG_JETS;
0121         break;
0122       case 4:
0123         flavour = C_JETS;
0124         break;
0125       case 5:
0126         flavour = B_JETS;
0127         break;
0128       default:
0129         flavour = NONID_JETS;
0130     }
0131
0132     // simply count the number of accepted jets
0133     flavours_->Fill(ALL_JETS);
0134     flavours_->Fill(flavour);
0135
0136     double discrTC = jet->bDiscriminator("trackCountingHighEffBJetTags");
0137     double discrSSV = jet->bDiscriminator("simpleSecondaryVertexBJetTags");
0138     double discrCSV = jet->bDiscriminator("combinedSecondaryVertexBJetTags");
0139
0140     plots_[ALL_JETS].discrTC->Fill(discrTC);
0141     plots_[flavour].discrTC->Fill(discrTC);
0142
0143     plots_[ALL_JETS].discrSSV->Fill(discrSSV);
0144     plots_[flavour].discrSSV->Fill(discrSSV);
0145
0146     plots_[ALL_JETS].discrCSV->Fill(discrCSV);
0147     plots_[flavour].discrCSV->Fill(discrCSV);
0148   }
0149 }
0150
0151 #include "FWCore/Framework/interface/MakerMacros.h"
0152
0153 DEFINE_FWK_MODULE(PatBJetTagAnalyzer);