Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:27:14

0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "DataFormats/JetReco/interface/GenJet.h"
0005 #include "DataFormats/PatCandidates/interface/Tau.h"
0006 #include "FWCore/ServiceRegistry/interface/Service.h"
0007 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0008 
0009 #include "TH1D.h"
0010 #include "TH2D.h"
0011 #include "TGraphAsymmErrors.h"
0012 
0013 class rerunMVAIsolationOnMiniAOD_Phase2 : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0014 public:
0015   explicit rerunMVAIsolationOnMiniAOD_Phase2(const edm::ParameterSet &);
0016 
0017 private:
0018   virtual void analyze(const edm::Event &, const edm::EventSetup &) override;
0019   virtual void endJob() override;
0020 
0021   edm::EDGetTokenT<std::vector<pat::Tau>> tauToken_;
0022 
0023   TH1D *h_raw_bkg, *h_raw_sig;
0024 
0025   TH1D *h_pt_denom_bkg, *h_pt_denom_sig;
0026   TH1D *h_pt_num_bkg[7], *h_pt_num_sig[7];
0027   TGraphAsymmErrors *eff[7], *fak[7];
0028 
0029   TH1D *h_inc_denom_bkg, *h_inc_num_bkg;
0030   TH1D *h_inc_denom_sig, *h_inc_num_sig;
0031   TGraphAsymmErrors *inceff, *incfak;
0032 
0033   TH2D *h2_denom_bkg, *h2_denom_sig;
0034   TH2D *h2_num_bkg[7], *h2_num_sig[7];
0035   TH2D *eff2[7], *fak2[7];
0036 
0037   bool haveGenJets;
0038   edm::EDGetTokenT<std::vector<reco::GenJet>> genJetToken_;
0039 };
0040 
0041 rerunMVAIsolationOnMiniAOD_Phase2::rerunMVAIsolationOnMiniAOD_Phase2(const edm::ParameterSet &iConfig) {
0042   tauToken_ = consumes<pat::TauCollection>(iConfig.getParameter<edm::InputTag>("tauCollection"));
0043 
0044   haveGenJets = false;
0045   if (iConfig.existsAs<edm::InputTag>("genJetCollection")) {
0046     haveGenJets = true;
0047     genJetToken_ = consumes<std::vector<reco::GenJet>>(iConfig.getParameter<edm::InputTag>("genJetCollection"));
0048   }
0049 
0050   edm::Service<TFileService> fileService;
0051   h_raw_sig = fileService->make<TH1D>("h_raw_sig", ";BDT output;taus / 0.02", 110, -1.1, 1.1);
0052   h_raw_bkg = fileService->make<TH1D>("h_raw_bkg", ";BDT output;taus / 0.02", 110, -1.1, 1.1);
0053 
0054   h_inc_num_sig = fileService->make<TH1D>("h_inc_num_sig", ";working point;taus / bin", 7, 0.5, 7.5);
0055   h_inc_denom_sig = fileService->make<TH1D>("h_inc_denom_sig", ";working point;taus / bin", 7, 0.5, 7.5);
0056   h_inc_num_bkg = fileService->make<TH1D>("h_inc_num_bkg", ";working point;taus / bin", 7, 0.5, 7.5);
0057   h_inc_denom_bkg = fileService->make<TH1D>("h_inc_denom_bkg", ";working point;taus / bin", 7, 0.5, 7.5);
0058   inceff = fileService->make<TGraphAsymmErrors>(7);
0059   inceff->SetName("inceff");
0060   incfak = fileService->make<TGraphAsymmErrors>(7);
0061   incfak->SetName("incfak");
0062 
0063   const int n = 10;
0064   const double x[n + 1] = {20., 25.4196, 32.3079, 41.0627, 52.19, 66.3325, 84.3074, 107.153, 136.19, 173.095, 220.};
0065   h_pt_denom_sig = fileService->make<TH1D>("h_pt_denom_sig", ";p_{T} [GeV];taus / bin", n, x);
0066   h_pt_denom_bkg = fileService->make<TH1D>("h_pt_denom_bkg", ";p_{T} [GeV];taus / bin", n, x);
0067   for (int i = 0; i < 7; ++i) {
0068     const TString tag = TString::Itoa(i, 10);
0069     h_pt_num_sig[i] = fileService->make<TH1D>("h_pt_num_sig_" + tag, ";p_{T} [GeV];taus / bin", n, x);
0070     h_pt_num_bkg[i] = fileService->make<TH1D>("h_pt_num_bkg_" + tag, ";p_{T} [GeV];taus / bin", n, x);
0071     eff[i] = fileService->make<TGraphAsymmErrors>(10);
0072     eff[i]->SetName("eff_" + TString::Itoa(i, 10));
0073     fak[i] = fileService->make<TGraphAsymmErrors>(10);
0074     fak[i]->SetName("fak_" + TString::Itoa(i, 10));
0075   }
0076   const double xeta[4] = {0., 1.5, 2.3, 3.};
0077   const double xpt[3] = {20., 50., 220.};
0078   h2_denom_sig = fileService->make<TH2D>("h2_denom_sig", ";|#eta|;p_{T} [GeV]", 3, xeta, 2, xpt);
0079   h2_denom_bkg = fileService->make<TH2D>("h2_denom_bkg", ";|#eta|;p_{T} [GeV]", 3, xeta, 2, xpt);
0080   for (int i = 0; i < 7; ++i) {
0081     const TString tag = TString::Itoa(i, 10);
0082     h2_num_sig[i] = fileService->make<TH2D>("h2_num_sig_" + tag, ";|#eta|;p_{T} [GeV]", 3, xeta, 2, xpt);
0083     h2_num_sig[i]->Sumw2();
0084     h2_num_bkg[i] = fileService->make<TH2D>("h2_num_bkg_" + tag, ";|#eta|;p_{T} [GeV]", 3, xeta, 2, xpt);
0085     h2_num_bkg[i]->Sumw2();
0086     eff2[i] = fileService->make<TH2D>("eff2_" + tag, ";|#eta|;p_{T} [GeV]", 3, xeta, 2, xpt);
0087     eff2[i]->Sumw2();
0088     fak2[i] = fileService->make<TH2D>("fak2_" + tag, ";|#eta|;p_{T} [GeV]", 3, xeta, 2, xpt);
0089     fak2[i]->Sumw2();
0090   }
0091 }
0092 
0093 void rerunMVAIsolationOnMiniAOD_Phase2::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0094   edm::Handle<pat::TauCollection> taus;
0095   iEvent.getByToken(tauToken_, taus);
0096 
0097   const bool isRealData = iEvent.isRealData();
0098 
0099   edm::Handle<std::vector<reco::GenJet>> genJets;
0100   if (!isRealData && haveGenJets) {
0101     iEvent.getByToken(genJetToken_, genJets);
0102   } else {
0103     haveGenJets = false;
0104   }
0105 
0106   for (auto i = taus->begin(); i != taus->end(); ++i) {
0107     const double pt = i->pt();
0108     const double eta = std::abs(i->eta());
0109     if (pt < 20. || pt >= 220. || eta >= 3.)
0110       continue;
0111     if (!i->tauID("decayModeFindingNewDMs"))
0112       continue;
0113 
0114     bool isPU = false;
0115     if (haveGenJets) {
0116       double mindr_jet = 9.;
0117       for (auto j = genJets->begin(); j != genJets->end(); ++j) {
0118         if (reco::deltaR(*i, *j) < mindr_jet) {
0119           mindr_jet = reco::deltaR(*i, *j);
0120         }
0121       }
0122       if (mindr_jet >= 0.4)
0123         isPU = true;
0124     }
0125     if (isPU)
0126       continue;
0127 
0128     bool isSig = false;
0129     if (!isRealData) {
0130       if (i->genJet()) {
0131         isSig = true;
0132       }
0133     }
0134 
0135     const double byIsolationMVAPhase2raw = i->tauID("byIsolationMVADBnewDMwLTPhase2raw");
0136     double wp[7];
0137     wp[0] = i->tauID("byVVLooseIsolationMVADBnewDMwLTPhase2");
0138     wp[1] = i->tauID("byVLooseIsolationMVADBnewDMwLTPhase2");
0139     wp[2] = i->tauID("byLooseIsolationMVADBnewDMwLTPhase2");
0140     wp[3] = i->tauID("byMediumIsolationMVADBnewDMwLTPhase2");
0141     wp[4] = i->tauID("byTightIsolationMVADBnewDMwLTPhase2");
0142     wp[5] = i->tauID("byVTightIsolationMVADBnewDMwLTPhase2");
0143     wp[6] = i->tauID("byVVTightIsolationMVADBnewDMwLTPhase2");
0144 
0145     if (isSig) {
0146       h_raw_sig->Fill(byIsolationMVAPhase2raw);
0147       h_pt_denom_sig->Fill(pt);
0148       h2_denom_sig->Fill(eta, pt);
0149       for (int j = 0; j < 7; ++j) {
0150         h_inc_denom_sig->Fill(j + 1);
0151         if (wp[j]) {
0152           h_inc_num_sig->Fill(j + 1);
0153           h_pt_num_sig[j]->Fill(pt);
0154           h2_num_sig[j]->Fill(eta, pt);
0155         }
0156       }
0157     } else {
0158       h_raw_bkg->Fill(byIsolationMVAPhase2raw);
0159       h_pt_denom_bkg->Fill(pt);
0160       h2_denom_bkg->Fill(eta, pt);
0161       for (int j = 0; j < 7; ++j) {
0162         h_inc_denom_bkg->Fill(j + 1);
0163         if (wp[j]) {
0164           h_inc_num_bkg->Fill(j + 1);
0165           h_pt_num_bkg[j]->Fill(pt);
0166           h2_num_bkg[j]->Fill(eta, pt);
0167         }
0168       }
0169     }
0170   }
0171 }
0172 
0173 void rerunMVAIsolationOnMiniAOD_Phase2::endJob() {
0174   inceff->Divide(h_inc_num_sig, h_inc_denom_sig, "e0");
0175   inceff->SetTitle(";working point;tagging efficiency");
0176   incfak->Divide(h_inc_num_bkg, h_inc_denom_bkg, "e0");
0177   incfak->SetTitle(";working point;tagging efficiency");
0178   for (int i = 0; i < 7; ++i) {
0179     eff[i]->Divide(h_pt_num_sig[i], h_pt_denom_sig, "e0");
0180     eff[i]->SetTitle(";p_{T} [GeV];tagging efficiency");
0181     fak[i]->Divide(h_pt_num_bkg[i], h_pt_denom_bkg, "e0");
0182     fak[i]->SetTitle(";p_{T} [GeV];tagging efficiency");
0183     eff2[i]->Divide(h2_num_sig[i], h2_denom_sig, 1., 1., "B");
0184     eff2[i]->SetTitle("tagging efficiency;|eta|;p_{T} [GeV]");
0185     fak2[i]->Divide(h2_num_bkg[i], h2_denom_bkg, 1., 1., "B");
0186     fak2[i]->SetTitle("tagging efficiency;|eta|;p_{T} [GeV]");
0187   }
0188 }
0189 
0190 //define this as a plug-in
0191 DEFINE_FWK_MODULE(rerunMVAIsolationOnMiniAOD_Phase2);