File indexing completed on 2024-04-06 12:27:55
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
0191 DEFINE_FWK_MODULE(rerunMVAIsolationOnMiniAOD_Phase2);
