diff --git a/DeepNtuplizer/data/EfficienciesAndSF_ID_GH.root b/DeepNtuplizer/data/EfficienciesAndSF_ID_GH.root new file mode 100644 index 00000000000..7c59c95fb43 Binary files /dev/null and b/DeepNtuplizer/data/EfficienciesAndSF_ID_GH.root differ diff --git a/DeepNtuplizer/data/EfficienciesAndSF_ISO_GH.root b/DeepNtuplizer/data/EfficienciesAndSF_ISO_GH.root new file mode 100644 index 00000000000..7227b8c7e11 Binary files /dev/null and b/DeepNtuplizer/data/EfficienciesAndSF_ISO_GH.root differ diff --git a/DeepNtuplizer/data/EfficienciesAndSF_Period4.root b/DeepNtuplizer/data/EfficienciesAndSF_Period4.root new file mode 100644 index 00000000000..6806d8ad161 Binary files /dev/null and b/DeepNtuplizer/data/EfficienciesAndSF_Period4.root differ diff --git a/DeepNtuplizer/data/Tracking_EfficienciesAndSF_BCDEFGH.root b/DeepNtuplizer/data/Tracking_EfficienciesAndSF_BCDEFGH.root new file mode 100644 index 00000000000..c29c506ba39 Binary files /dev/null and b/DeepNtuplizer/data/Tracking_EfficienciesAndSF_BCDEFGH.root differ diff --git a/DeepNtuplizer/data/egammaEffi.txt_EGM2D.root b/DeepNtuplizer/data/egammaEffi.txt_EGM2D.root new file mode 100644 index 00000000000..265448065aa Binary files /dev/null and b/DeepNtuplizer/data/egammaEffi.txt_EGM2D.root differ diff --git a/DeepNtuplizer/data/pileup_MC_2016.root b/DeepNtuplizer/data/pileup_MC_2016.root new file mode 100644 index 00000000000..4a2a782096c Binary files /dev/null and b/DeepNtuplizer/data/pileup_MC_2016.root differ diff --git a/DeepNtuplizer/data/pileup_data_2016.root b/DeepNtuplizer/data/pileup_data_2016.root new file mode 100644 index 00000000000..db1564ce9a6 Binary files /dev/null and b/DeepNtuplizer/data/pileup_data_2016.root differ diff --git a/DeepNtuplizer/interface/ntuple_JetInfo.h b/DeepNtuplizer/interface/ntuple_JetInfo.h index c99a0b7117d..efc9ed40027 100644 --- a/DeepNtuplizer/interface/ntuple_JetInfo.h +++ b/DeepNtuplizer/interface/ntuple_JetInfo.h @@ -21,8 +21,9 @@ class ntuple_JetInfo: public ntuple_content{ ntuple_JetInfo():ntuple_content(), gluonReduction_(0), useherwcompat_matching_(false), - isherwig_(false) -{} + isherwig_(false), + isData_(false) + {} void getInput(const edm::ParameterSet& iConfig); void initBranches(TTree* ); @@ -130,6 +131,8 @@ class ntuple_JetInfo: public ntuple_content{ bool useherwcompat_matching_; bool isherwig_; + bool isData_; + bool removeUndefined_; /////////branches @@ -167,6 +170,7 @@ class ntuple_JetInfo: public ntuple_content{ int isPhysLeptonicB_; int isPhysLeptonicB_C_; int isPhysTau_; + int isRealData_; // global variables float npv_; diff --git a/DeepNtuplizer/interface/ntuple_eventInfo.h b/DeepNtuplizer/interface/ntuple_eventInfo.h new file mode 100644 index 00000000000..0a25217eca4 --- /dev/null +++ b/DeepNtuplizer/interface/ntuple_eventInfo.h @@ -0,0 +1,133 @@ +/* + * ntuple_eventInfo.h + * + * Created on: 20 Feb 2018 + * Author: dwalter + */ + +#ifndef DEEPNTUPLES_DEEPNTUPLIZER_INTERFACE_NTUPLE_EVENTINFO_H_ +#define DEEPNTUPLES_DEEPNTUPLIZER_INTERFACE_NTUPLE_EVENTINFO_H_ + +#include "ntuple_content.h" +#include "triggerInfo.h" + +#include +#include +#include +#include +#include +#include "TGraphAsymmErrors.h" +#include + + + +#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" + +void readHistoFromGraph(TGraphAsymmErrors* graph, TH1D** h, TString name); +void initializeScalefactor(std::vector dirs, std::vector names, std::vector* sfHists, std::vector periods); +void initializeScalefactor(std::vector dirs, std::vector names, std::vector* sfHists, std::vector periods); +double getScalefactor(double x, double y, std::vector hists, unsigned int period, bool takeOverflowX = false); +double getScalefactor(double x, std::vector hists, unsigned int period); + + +/* + * For MC weights such as pileup, lhe, ... later: lepton scalefactors + */ +class ntuple_eventInfo: public ntuple_content{ +public: + ntuple_eventInfo():ntuple_content(), + isData_(false) +{} + + void getInput(const edm::ParameterSet& iConfig); + void initBranches(TTree* ); + void readEvent(const edm::Event& iEvent); + + //use either of these functions + + bool fillBranches(const pat::Jet &, const size_t& jetidx, const edm::View * coll=0); + + void setLHEToken(edm::EDGetTokenT lheToken) { + lheToken_ = lheToken; + } + void setMuonsToken(edm::EDGetTokenT > muonToken){ + muonToken_ = muonToken; + } + void setElectronsToken(edm::EDGetTokenT > electronToken){ + electronToken_ = electronToken; + } + void setTriggerToken(edm::EDGetTokenT triggerToken){ + triggerToken_ = triggerToken; + } + + + +private: + + edm::EDGetTokenT lheToken_; + edm::EDGetTokenT > muonToken_; + edm::EDGetTokenT > electronToken_; + edm::EDGetTokenT triggerToken_; + + edm::Handle lheInfo; + edm::Handle > muons; + edm::Handle > electrons; + + + bool isData_; + bool useLHEWeights_; + + std::vector lumis; + std::vector periods; + + double sigma; + unsigned int nEvents; + + std::string pupDataDir_; + std::string pupMCDir_; + + //std::string sfTrigger_mu_BCDEF_Dir_; + //std::string sfTrigger_mu_BCDEF_Name_; + //std::string sfTrigger_mu_GH_Dir_; + //std::string sfTrigger_mu_GH_Name_; + + std::vector sfTrigger_mu_Dir_; + std::vector sfTrigger_mu_Name_; + std::vector sfTrigger_e_Dir_; + std::vector sfTrigger_e_Name_; + std::vector sfTrigger_emu_Dir_; + std::vector sfTrigger_emu_Name_; + std::vector sfMuonId_Dir_; + std::vector sfMuonId_Name_; + std::vector sfMuonIso_Dir_; + std::vector sfMuonIso_Name_; + std::vector sfElIdAndIso_Dir_; + std::vector sfElIdAndIso_Name_; + std::vector sfMuonTracking_Dir_; + std::vector sfMuonTracking_Name_; + + std::vector sfTrigger_mu_Hist; + std::vector sfTrigger_e_Hist; + std::vector sfTrigger_emu_Hist; + std::vector sfMuonId_Hist; + std::vector sfMuonIso_Hist; + std::vector sfElIdAndIso_Hist; + + std::vector sfMuonTracking_Hist; + + std::vector triggers; + + // global variables + + float ntrueInt_; + std::vector pupWeights; + + /////////branches + float event_weight_; + + + +}; + + +#endif /* DEEPNTUPLES_DEEPNTUPLIZER_INTERFACE_NTUPLE_EVENTINFO_H_ */ diff --git a/DeepNtuplizer/interface/sorting_modules.h b/DeepNtuplizer/interface/sorting_modules.h index dbd16082401..cb0fb3fc719 100644 --- a/DeepNtuplizer/interface/sorting_modules.h +++ b/DeepNtuplizer/interface/sorting_modules.h @@ -100,10 +100,23 @@ class sortingClass{ } static bool compareByABCInv(const sortingClass& a, const sortingClass& b){ - return !compareByABC(a,b); + + compareResult tmpres=compare(a,b,0); + if(tmpres==cmp_smaller) return false; + if(tmpres==cmp_greater) return true; + + tmpres=compare(a,b,1); + if(tmpres==cmp_smaller) return false; + if(tmpres==cmp_greater) return true; + + tmpres=compare(a,b,2); + if(tmpres==cmp_smaller) return false; + if(tmpres==cmp_greater) return true; + + return false; } - //private: + //private: float sortValA,sortValB,sortValC; diff --git a/DeepNtuplizer/interface/triggerInfo.h b/DeepNtuplizer/interface/triggerInfo.h new file mode 100644 index 00000000000..74ba7725f53 --- /dev/null +++ b/DeepNtuplizer/interface/triggerInfo.h @@ -0,0 +1,36 @@ +/* + * treeReader.h + * + * Created on: 17 May 2018 + * Author: dwalter + */ + +#ifndef DEEPNTUPLES_DEEPNTUPLIZER_INTERFACE_TRIGGERINFO_H_ +#define DEEPNTUPLES_DEEPNTUPLIZER_INTERFACE_TRIGGERINFO_H_ + + +#include +#include +#include "FWCore/Common/interface/TriggerNames.h" +#include "DataFormats/Common/interface/TriggerResults.h" +#include "DataFormats/PatCandidates/interface/TriggerObjectStandAlone.h" +#include "DataFormats/PatCandidates/interface/PackedTriggerPrescales.h" +#include "FWCore/Framework/interface/Event.h" + +class TriggerInfo { + +public: + TriggerInfo(const edm::Event& iEvent, + const edm::EDGetTokenT& triggerBitsToken); + bool IsTriggered(std::string triggername) const ; + bool Exists(std::string triggername) const ; + bool IsAnyTriggered(std::vector< std::string > triggers) const ; + bool IsAnyTriggered(std::string triggerstring) const; + void ListTriggers() ; + std::map GetTriggers() const; +private: + std::map triggers; + std::vector triggernames; +}; + +#endif diff --git a/DeepNtuplizer/plots/plot_event_ttdilep.py b/DeepNtuplizer/plots/plot_event_ttdilep.py new file mode 100644 index 00000000000..e1214f3ae93 --- /dev/null +++ b/DeepNtuplizer/plots/plot_event_ttdilep.py @@ -0,0 +1,651 @@ + +import pdb +from ROOT import TH1F, TFile, TCanvas, gROOT, gStyle, THStack, TLegend, TGaxis, gPad, TPad, TGraph, TLine +import numpy as np +import os +from array import array +import datetime +from DataFormats.FWLite import Events, Handle + +gROOT.SetBatch() # don't pop up canvases +gROOT.SetStyle('Plain') # white background +gStyle.SetFillStyle(0) +TGaxis.SetMaxDigits(4) # Force scientific notation for numbers with more than 4 digits +#gStyle.SetBorderSize(0) +#gStyle.SetTextFont(43) +#gStyle.SetTextSize(14) + + + +class pileupweight: + """ DOC """ + + def __init__(self, file_pileup_MC, file_pileup_Data): + hist_mc = file_pileup_MC.Get("pileup") + hist_data = file_pileup_Data.Get("pileup") + self.pupWeights = [0] + hist_mc.Scale(1. / hist_mc.Integral()) + hist_data.Scale(1. / hist_data.Integral()) + for bin in range(1, hist_data.GetNbinsX() + 2): + if hist_mc.GetBinContent(bin) != 0: + w = hist_data.GetBinContent(bin) / hist_mc.GetBinContent(bin) + self.pupWeights.append(w) + else: + self.pupWeights.append(1) + + def getPileupWeight(self,nInteractions): + if nInteractions >= len(self.pupWeights): + return 0 + return self.pupWeights[nInteractions] + +class scaleFactor(): + """ + class for isolation and id scale factors + hist is the 2d histogram with the scale factor with x and y + """ + def __init__(self,hist): + self.hist = hist + self.xaxis = self.hist.GetXaxis() + self.yaxis = self.hist.GetYaxis() + self.binx = 0 + self.biny = 0 + self.scalefactor = 1. + self.maxBinX = self.xaxis.GetLast() + self.maxBinY = self.yaxis.GetLast() + self.minX = self.xaxis.GetXmin() + self.maxX = self.xaxis.GetXmax() + self.minY = self.yaxis.GetXmin() + self.maxY = self.yaxis.GetXmax() + + def getScalefactor(self,x,y): + + self.binx = self.xaxis.FindBin(x) + self.biny = self.yaxis.FindBin(y) + if x >= self.maxX: self.binx = self.maxBinX + if y >= self.maxY: self.biny = self.maxBinY + self.scalefactor = self.hist.GetBinContent(self.binx, self.biny) + return self.scalefactor + + +class scaleFactor_tracking(): + """ + class for tracking scale factors + the tracking scale factors are saved in a TGraphAsymmErrors object which has to be written in a histogram first + """ + def __init__(self,graph, name): + self.graph = graph + self.npoints = self.graph.GetN() + self.x_centers = self.graph.GetX() + self.y_centers = self.graph.GetY() + self.x_lows = np.zeros(self.npoints) + self.x_highs = np.zeros(self.npoints) + self.y_lows = np.zeros(self.npoints) + self.y_highs = np.zeros(self.npoints) + self.x_edges = np.zeros(self.npoints +1) + + for i in range(0,self.npoints): + self.x_lows[i] = self.graph.GetErrorXlow(i) + self.x_highs[i] = self.graph.GetErrorXhigh(i) + self.y_lows[i] = self.graph.GetErrorYlow(i) + self.y_lows[i] = self.graph.GetErrorYhigh(i) + + self.x_edges[i] = self.x_centers[i] - self.x_lows[i] + self.x_edges[self.npoints] = self.x_centers[self.npoints - 1] + self.x_highs[self.npoints-1] + self.hist = TH1F(name, name,self.npoints,self.x_edges) + + for i in range(0,self.npoints): + self.hist.SetBinContent(i+1, self.y_centers[i]) + self.hist.SetBinError(i+1, max(self.y_lows[i], self.y_highs[i])) + + self.bin = 0 + self.axis = self.hist.GetXaxis() + self.scalefactor = 1. + + def getScalefactor(self,x): + self.bin = self.axis.FindBin(abs(x)) + self.scalefactor = self.hist.GetBinContent(self.bin) + return self.scalefactor + + +class process: + """ DOC """ + datadir = "/afs/desy.de/user/d/dwalter/CMSSW_8_0_29/src/DeepNTuples/DeepNtuplizer/data/" + pupw = pileupweight( + TFile(datadir + "MyMCPileupHist.root"), + TFile(datadir + "MyDataPileupHist.root") + ) + # emu trigger scalefactors from + # https://gitlab.cern.ch/ttH/reference/blob/955cff0b09f2c95bc480ae0bf3145aab9ce08fcc/definitions/Moriond17.md#64-trigger-scale-factors + # https://gitlab.cern.ch/ttH/reference/blob/955cff0b09f2c95bc480ae0bf3145aab9ce08fcc/definitions/Moriond17.md#41-triggers + file_sf_emu_trigger = TFile(datadir + "triggerSummary_emu_ReReco2016_ttH.root") + h_sf_emu_trigger = file_sf_emu_trigger.Get("scalefactor_eta2d_with_syst") + sf_emu_trigger = scaleFactor(h_sf_emu_trigger) + + # muon scalefactor hists from https://twiki.cern.ch/twiki/bin/viewauth/CMS/MuonWorkInProgressAndPagResults + file_sf_muon_id_BCDEF = TFile(datadir + "EfficienciesAndSF_ID_BCDEF.root") + file_sf_muon_iso_BCDEF = TFile(datadir + "EfficienciesAndSF_ISO_BCDEF.root") + file_sf_muon_id_GH = TFile(datadir + "EfficienciesAndSF_ID_GH.root") + file_sf_muon_iso_GH = TFile(datadir + "EfficienciesAndSF_ISO_GH.root") + file_sf_muon_tracking = TFile(datadir + "Tracking_EfficienciesAndSF_BCDEFGH.root") + + h_sf_muon_id_BCDEF = file_sf_muon_id_BCDEF.Get("MC_NUM_TightID_DEN_genTracks_PAR_pt_eta/abseta_pt_ratio") + h_sf_muon_iso_BCDEF = file_sf_muon_iso_BCDEF.Get("TightISO_TightID_pt_eta/abseta_pt_ratio") + h_sf_muon_id_GH = file_sf_muon_id_GH.Get("MC_NUM_TightID_DEN_genTracks_PAR_pt_eta/abseta_pt_ratio") + h_sf_muon_iso_GH = file_sf_muon_iso_GH.Get("TightISO_TightID_pt_eta/abseta_pt_ratio") + g_sf_muon_tracking = file_sf_muon_tracking.Get("ratio_eff_aeta_dr030e030_corr") + + sf_muon_id_BCDEF = scaleFactor(h_sf_muon_id_BCDEF) + sf_muon_iso_BCDEF = scaleFactor(h_sf_muon_iso_BCDEF) + sf_muon_id_GH = scaleFactor(h_sf_muon_id_GH) + sf_muon_iso_GH = scaleFactor(h_sf_muon_iso_GH) + sf_muon_tracking = scaleFactor_tracking(g_sf_muon_tracking, "sf_muon_tracking") + # electron scalefactor hists from https://twiki.cern.ch/twiki/bin/view/CMS/EgammaIDRecipesRun2 + # -> Electron cut-based 80XID WPs. Scale factors for 80X -> Tight cut-based ID WP scale factor (root file) + # id and iso is combined in one file, there is no tracking for the electron + file_sf_el_id = TFile(datadir + "egammaEffi.txt_EGM2D.root") + h_sf_el_id = file_sf_el_id.Get("EGamma_SF2D") + sf_el_id = scaleFactor(h_sf_el_id) + + muonhandle = Handle("vector") +#muonhandle = Handle('edm::OwnVector >') + muonlabel = ("goodMuons") +#muonlabel = ("GoodMuon") + electronhandle = Handle("vector") +#electronhandle = Handle('edm::OwnVector >') + electronlabel = ("goodElectrons") +#electronlabel = ("GoodElectron") + jethandle = Handle('vector') + #jetlabel = ('slimmedJets') + jetlabel = ("goodJets") + ### LHE weights + LHEweighthandle = Handle('LHEEventProduct') + LHEweightlabel = ("externalLHEProducer") + + ## Pileup reweighting + puphandle = Handle('vector') + puplabel = ("slimmedAddPileupInfo") + + triggerhandle = Handle("edm::TriggerResults") + triggerlabel = ("TriggerResults::HLT") + + hs_ll_pt = THStack("hs_ll_pt", "ll pt") + hs_ll_eta = THStack("hs_ll_eta", "ll eta") + hs_tl_pt = THStack("hs_tl_pt", "tl pt") + hs_tl_eta = THStack("hs_tl_eta", "tl eta") + hs_j_n = THStack("hs_njets", "jet number") + + hs_list = hs_ll_pt, hs_ll_eta, hs_tl_pt, hs_tl_eta, hs_j_n + + collection = [] + hists_data = None + + tsize = 20 + tfont = 43 + + def __init__(self, name, isData=False, color=1): + self.name = name + self.isData = isData + self.h_ll_pt = TH1F(self.name + "_ll_pt", "pt of leading lepton in " + self.name, 10, array('d', [0, 25, 30, 35, 40, 45, 50, 60, 70, 100, 200])) + self.h_ll_eta = TH1F(self.name + "_ll_eta", "eta of the leading leption in " + self.name, 20, -2.4, 2.4) + self.h_tl_pt = TH1F(self.name + "_tl_pt", "pt of leading lepton in " + self.name, 11, array('d', [0, 20, 25, 30, 35, 40, 45, 50, 60, 70, 100, 200])) + self.h_tl_eta = TH1F(self.name + "_tl_eta", "eta of the leading leption in " + self.name, 20, -2.4, 2.4) + self.h_njets = TH1F(self.name + "_njets", "number of jets in " +self.name, 10, -0.5, 9.5) + + #self.h_ll_pt.Sumw2() + #self.h_ll_eta.Sumw2() + #self.h_tl_pt.Sumw2() + #self.h_tl_eta.Sumw2() + #self.h_njets.Sumw2() + + self.hists = self.h_ll_pt, self.h_ll_eta, self.h_tl_pt, self.h_tl_eta, self.h_njets + + self.integral = 0. + + #remove status boxes + for hist in self.hists: + hist.SetStats(False) + + if self.isData==False: + process.collection.append(self) + for hist in self.hists: + hist.SetFillColor(color) + else: + process.hists_data = self + + self.h_ll_pt.SetXTitle("leading lept. p_{T} (GeV)") + self.h_ll_eta.SetXTitle("leading lept. eta") + self.h_tl_pt.SetXTitle("trailing lept. p_{T} (GeV)") + self.h_tl_eta.SetXTitle("trailing lept. eta") + self.h_njets.SetXTitle("N_{Jet}") + + self.h_ll_pt.SetYTitle("Events/GeV") + self.h_ll_pt.GetYaxis().SetTitleOffset(1.2) + self.h_ll_eta.SetYTitle("Events") + self.h_tl_pt.SetYTitle("Events/GeV") + self.h_tl_eta.SetYTitle("Events") + self.h_njets.SetYTitle("Events") + + self.h_njets.GetXaxis().SetNdivisions(10,0,0) + + def fillHists(self, infile, islist=True, useLHEWeights=False, isData=False, eventIDs=[]): + print("start to fill hist: " + str(datetime.datetime.now())) + nevents = 0 + directory = '' + + if islist is True: + print("loop over filelist " + str(infile)) + filelist = open(infile, 'r') + ifile = filelist.readline() + ifile = ifile[:-1] + directory = infile[:-len(infile.split('/')[-1])] + else: + ifile = infile + + while ifile != '': + print("process file ", ifile) + events = Events(directory + ifile) + + ### EVENT LOOP + for event in events: + eventID = event.eventAuxiliary().event() + + if isData: # avoid double counting + if eventID in eventIDs: + # print("this event was already processed") + continue + eventIDs = np.append(eventIDs, eventID) + + nevents += 1 + leadingPt = 10 + leadingEta = 0 + trailingPt = 10 + trailingEta = 0 + leadingMuonPt = 10 + leadingMuonEta = 0 + leadingElPt = 10 + leadingElEta = 0 + leadingElSuEta = 0 + weight = 1. + weight_BCDEF = 1. + weight_GH = 1. + + lumi_BCDEF = 5.404+2.396+4.243+4.054+3.105 + lumi_G = 7.544 + lumi_H = 8.453 + + + + # Fill muon hists + event.getByLabel(process.muonlabel, process.muonhandle) + muons = process.muonhandle.product() # get the product + + for muon in muons: + if muon.pt() > leadingMuonPt: + leadingMuonPt = muon.pt() + leadingMuonEta = muon.eta() + + # Fill electron hists + event.getByLabel(process.electronlabel, process.electronhandle) + electrons = process.electronhandle.product() + for electron in electrons: + + if electron.pt() > leadingElPt: + leadingElPt = electron.pt() + leadingElEta = electron.eta() + leadingElSuEta = electron.superCluster().eta() + + if leadingMuonPt > leadingElPt: + leadingPt = leadingMuonPt + leadingEta = leadingMuonEta + trailingPt = leadingElPt + trailingEta = leadingElEta + else: + leadingPt = leadingElPt + leadingEta = leadingElEta + trailingPt = leadingMuonPt + trailingEta = leadingMuonEta + + + + # Fill jet hists + event.getByLabel(process.jetlabel, process.jethandle) + jets = process.jethandle.product() + + numJets = len(jets) + + # Compute Event weights + if useLHEWeights == True: + event.getByLabel(process.LHEweightlabel, process.LHEweighthandle) + weights = process.LHEweighthandle.product() + weight = weights.weights()[0].wgt / abs(weights.weights()[0].wgt) + + if isData == False: + event.getByLabel(process.puplabel, process.puphandle) + pupInfos = process.puphandle.product() + for pupInfo in pupInfos: + if (pupInfo.getBunchCrossing() == 0): + weight *= process.pupw.getPileupWeight(int(pupInfo.getTrueNumInteractions())) + + weight *= process.sf_emu_trigger.getScalefactor(abs(leadingElEta),abs(leadingMuonEta)) + weight *= process.sf_muon_tracking.getScalefactor(abs(leadingMuonEta)) + weight_BCDEF *= process.sf_muon_id_BCDEF.getScalefactor(abs(leadingMuonEta), leadingMuonPt) + weight_BCDEF *= process.sf_muon_iso_BCDEF.getScalefactor(abs(leadingMuonEta), leadingMuonPt) + weight_GH *= process.sf_muon_id_GH.getScalefactor(abs(leadingMuonEta), leadingMuonPt) + weight_GH *= process.sf_muon_iso_GH.getScalefactor(abs(leadingMuonEta), leadingMuonPt) + weight *= process.sf_el_id.getScalefactor(abs(leadingElSuEta), leadingElPt) + + ### For 2016H there was used another emu trigger (DZ version) then for B to G, + # compute weights for events that pass the "2016H trigger" and for those who pass the "B to G trigger" + trigger_BCDEFG = 0 + trigger_H = 0 + event.getByLabel(process.triggerlabel, process.triggerhandle) + trigger = process.triggerhandle.product() + for i in range(trigger.size()): + #print(i, events.object().triggerNames(trigger).triggerName(i), trigger.accept(i)) + if events.object().triggerNames(trigger).triggerName(i)[0:49] == "HLT_Mu23_TrkIsoVVL_Ele12_CaloIdL_TrackIdL_IsoVL_v" or \ + events.object().triggerNames(trigger).triggerName(i)[0:48] == "HLT_Mu8_TrkIsoVVL_Ele23_CaloIdL_TrackIdL_IsoVL_v": + + if trigger.accept(i): + trigger_BCDEFG = 1 + + if events.object().triggerNames(trigger).triggerName(i)[0:52] == "HLT_Mu23_TrkIsoVVL_Ele12_CaloIdL_TrackIdL_IsoVL_DZ_v" or \ + events.object().triggerNames(trigger).triggerName(i)[0:51] == "HLT_Mu8_TrkIsoVVL_Ele23_CaloIdL_TrackIdL_IsoVL_DZ_v": + if trigger.accept(i): + trigger_H = 1 + + weight = weight * (weight_BCDEF * lumi_BCDEF * trigger_BCDEFG + + weight_GH * lumi_G * trigger_BCDEFG + + weight_GH * lumi_H * trigger_H) + + + + + self.h_ll_pt.Fill(leadingPt, weight) + self.h_ll_eta.Fill(leadingEta, weight) + self.h_tl_pt.Fill(trailingPt, weight) + self.h_tl_eta.Fill(trailingEta, weight) + self.h_njets.Fill(numJets, weight) + + if islist is True: + ifile = filelist.readline() + ifile = ifile[:-1] + else: + break + + if islist is True: + print("close filelist") + filelist.close() + + print("processed events ", nevents) + self.save_hists() + + def scale_hists(self,sigma=1.,lumi=1.,eff=1.): + for hist in self.hists: + hist.Scale(sigma*lumi*eff) + self.integral = self.hists[0].Integral() + + self.divideByBinWidth(self.hists[0]) + self.divideByBinWidth(self.hists[2]) + #self.hists[1].Scale(4) + #self.hists[3].Scale(4) + + def divideByBinWidth(self,hist): + for ibin in range(0, len(hist)): + content = hist.GetBinContent(ibin) + content /= hist.GetBinWidth(ibin) + hist.SetBinContent(ibin, content) + + def save_hists(self): + f1 = TFile(self.name + ".root", "RECREATE") + for hist in self.hists: + hist.Write() + f1.Close() + + def load_hists(self): + print "load " + self.name + + if(os.path.isfile(self.name + ".root")): + f1 = TFile(self.name + ".root", "READ") + self.h_ll_pt.Add(f1.Get(self.name + "_ll_pt")) + self.h_ll_eta.Add(f1.Get(self.name + "_ll_eta")) + self.h_tl_pt.Add(f1.Get(self.name + "_tl_pt")) + self.h_tl_eta.Add(f1.Get(self.name + "_tl_eta")) + self.h_njets.Add(f1.Get(self.name + "_njets")) + + f1.Close() + else: + print("no file for this process found") + + def draw_hists(self): + canvas = TCanvas(self.name, "", 800, 600) + for hist in self.hists: + canvas.Clear() + hist.Draw("HIST") + canvas.Print(hist.GetName()+".png") + + def setMinMax(self, y_min_ll_pt, y_max_ll_pt, y_min_ll_eta, y_max_ll_eta, y_min_tl_pt, y_max_tl_pt, y_min_tl_eta, y_max_tl_eta, y_min_njets, y_max_njets): + self.hists[0].SetMinimum(y_min_ll_pt) + self.hists[0].SetMaximum(y_max_ll_pt) + self.hists[1].SetMinimum(y_min_ll_eta) + self.hists[1].SetMaximum(y_max_ll_eta) + self.hists[2].SetMinimum(y_min_tl_pt) + self.hists[2].SetMaximum(y_max_tl_pt) + self.hists[3].SetMinimum(y_min_tl_eta) + self.hists[3].SetMaximum(y_max_tl_eta) + self.hists[4].SetMinimum(y_min_njets) + self.hists[4].SetMaximum(y_max_njets) + + @classmethod + def make_stacks(cls): + for proc in cls.collection: + for i,hist in enumerate(proc.hists): + cls.hs_list[i].Add(hist) + + + @classmethod + def print_stacks(cls): + cls.make_stacks() + + leg = TLegend(0.59, 0.54, 0.89, 0.84) + leg.SetBorderSize(0) + leg.SetTextFont(cls.tfont) + leg.AddEntry(cls.hists_data.hists[0], "data", "ep") + leg.AddEntry(cls.collection[8].hists[0], "t#bar{t}", "f") + leg.AddEntry(cls.collection[6].hists[0], "Wt/W#bar{t}", "f") + leg.AddEntry(cls.collection[5].hists[0], "W+Jets", "f") + leg.AddEntry(cls.collection[2].hists[0], "VV", "f") + leg.AddEntry(cls.collection[0].hists[0], "DY", "f") + + leg21 = TLegend(0.19, 0.69, 0.39, 0.84) + leg21.SetBorderSize(0) + leg21.SetTextFont(cls.tfont) + leg21.AddEntry(cls.hists_data.hists[0], "data", "ep") + leg21.AddEntry(cls.collection[8].hists[0], "t#bar{t}", "f") + leg21.AddEntry(cls.collection[6].hists[0], "Wt/W#bar{t}", "f") + + leg22 = TLegend(0.69, 0.69, 0.89, 0.84) + leg22.SetBorderSize(0) + leg22.SetTextFont(cls.tfont) + leg22.AddEntry(cls.collection[5].hists[0], "W+Jets", "f") + leg22.AddEntry(cls.collection[2].hists[0], "VV", "f") + leg22.AddEntry(cls.collection[0].hists[0], "DY", "f") + + leg.SetFillStyle(0) # make transparent + leg21.SetFillStyle(0) + leg22.SetFillStyle(0) + + canvas = TCanvas("c1","Example 2 pads (20,80)",200,10,800,600) + #canvas.SetBottomMargin(0.25) + pad1 = TPad("pad1", "The pad 80% of the height", 0, 0.2, 1, 1.0) + pad2 = TPad("pad2", "The pad 20% of the height", 0, 0.05, 1, 0.2) + + pad1.SetLeftMargin(0.15) + pad1.SetRightMargin(0.1) + pad2.SetLeftMargin(0.15) + pad2.SetRightMargin(0.1) + + pad1.SetBottomMargin(0.02) + pad2.SetTopMargin(0.02) + pad2.SetBottomMargin(0.25) + + pad1.Draw() + pad2.Draw() + for i, h_data in enumerate(cls.hists_data.hists): + + h_data.SetTitle("") + h_data.GetXaxis().SetLabelFont(cls.tfont) + h_data.GetXaxis().SetLabelSize(0) + h_data.GetXaxis().SetTitleSize(0) + h_data.GetYaxis().SetLabelFont(cls.tfont) + h_data.GetYaxis().SetLabelSize(cls.tsize) + h_data.GetYaxis().SetTitleFont(cls.tfont) + h_data.GetYaxis().SetTitleSize(cls.tsize) + h_data.GetYaxis().SetTitleOffset(1.5) + + h_data.SetMarkerStyle(20) + #pad1.Clear() + #pad2.Clear() + pad1.cd() + + h_data.Draw("E1") + cls.hs_list[i].Draw("HIST same") + h_data.Draw("E1 same") + gPad.RedrawAxis() # draw axis in foreground + + if i%2 == 0: + leg.Draw("same") + else: + leg21.Draw("same") + leg22.Draw("same") + + pad2.cd() + h_ratio = TH1F(cls.hs_list[i].GetStack().Last()) + + h_ratio.SetTitle("") + h_ratio.GetXaxis().SetLabelFont(cls.tfont) + h_ratio.GetXaxis().SetLabelSize(cls.tsize) + h_ratio.GetYaxis().SetLabelFont(cls.tfont) + h_ratio.GetYaxis().SetLabelSize(cls.tsize) + h_ratio.GetXaxis().SetTitleFont(cls.tfont) + h_ratio.GetXaxis().SetTitleSize(cls.tsize) + h_ratio.GetYaxis().SetTitleFont(cls.tfont) + h_ratio.GetYaxis().SetTitleSize(cls.tsize) + h_ratio.GetXaxis().SetTitleOffset(7.) + h_ratio.GetYaxis().SetTitleOffset(1.5) + h_ratio.GetYaxis().SetNdivisions(403) + h_ratio.GetXaxis().SetTickLength(0.15) + + if i ==4: + h_ratio.GetXaxis().SetNdivisions(10, 0, 0) + + h_ratio.SetMarkerStyle(20) + h_ratio.SetXTitle(h_data.GetXaxis().GetTitle()) + + h_ratio.Divide(h_data) + h_ratio.GetYaxis().SetTitle("mc/data") + + h_ratio.SetMinimum(0.75) + h_ratio.SetMaximum(1.25) + + + middleLine = TLine(h_ratio.GetXaxis().GetXmin(),1. ,h_ratio.GetXaxis().GetXmax(), 1.) + + h_ratio.Draw("e1") + middleLine.Draw("same") + gStyle.SetErrorX(0.) + canvas.Print("stacks_" + str(i) + ".png") + + + +dy10to50 = process("dy10to50", color=593) +dy50 = process("dy50", color=593) +ww = process("ww", color=393) +wz = process("wz", color=393) +zz = process("zz", color=393) +wjets = process("wjets", color=409) +wt = process("wt", color=609) +wantit = process("wantit", color=609) +tt = process("tt", color=625) +data = process("data",isData=True) + +#dy10to50.fillHists("dy10to50.txt",useLHEWeights=True) +#dy50.fillHists("dy50.txt",useLHEWeights=True) +#ww.fillHists("ww.txt") +#wz.fillHists("wz.txt") +#zz.fillHists("zz.txt") +#wjets.fillHists("wjets.txt",useLHEWeights=True) +#wt.fillHists("wt.txt") +#wantit.fillHists("wantit.txt") +#tt.fillHists("tt.txt") +#data.fillHists("data.txt",isData=True) + +for ip in process.collection: + ip.load_hists() +process.hists_data.load_hists() + +#lumi_total = 35.9 +#lumi_data = 7.544 # (7.544 + 8.453) + +sigma_tt = 831760 +sigma_dy50 = 6025200 +sigma_dy10to50 = 22635100 +sigma_ww = 118700 +sigma_wz = 44900 +sigma_zz = 15400 +sigma_wantit = 35600 +sigma_wt = 35600 +sigma_wjets = 61526700 + +# Efficiency = 1/(number of events before cuts), in case of lheweights its the +# effective number of events (in this case: number of events with positive weights - number of events with negative weights) +eff_tt = 1./(77081156 + 77867738) +#eff_tt = 1./(77013872) +eff_dy50 = 1./81781052 # effective number of events +#eff_dy50 = 1./122055388 # true number of events +eff_dy10to50 = 1./47946519 # effective number of events +#eff_dy10to50 = 1./65888233 # true number of events +eff_ww = 1./994012 +eff_wz = 1./1000000 +eff_zz = 1./998034 +eff_wantit = 1./6933094 +eff_wt = 1./6952830 +eff_wjets = 1./16497031 # effective number of events +#eff_wjets = 1./24120319 # true number of events + +data.scale_hists() +tt.scale_hists(sigma = sigma_tt, eff = eff_tt) +dy50.scale_hists(sigma = sigma_dy50, eff = eff_dy50) +dy10to50.scale_hists(sigma = sigma_dy10to50, eff = eff_dy10to50) +ww.scale_hists(sigma = sigma_ww, eff = eff_ww) +wz.scale_hists(sigma = sigma_wz, eff = eff_wz) +zz.scale_hists(sigma = sigma_zz, eff = eff_zz) +wantit.scale_hists(sigma = sigma_wantit, eff = eff_wantit) +wt.scale_hists(sigma = sigma_wt, eff = eff_wt) +wjets.scale_hists(sigma = sigma_wjets, eff = eff_wjets) + +data.setMinMax(0,8000,0,40000,0,17000,0,40000,0,140000) + +title='plots_180514_full2016' +directory = os.path.dirname('./'+title+'/') +# make a canvas, draw, and save it +if os.path.exists(directory): + print("ERROR: output directory exists, please choose another name or delite the old one") + import sys + sys.exit() +else: + os.makedirs(directory) + os.chdir(directory) + +integral_mc = 0. +integral_data = process.hists_data.integral +for ip in process.collection: + integral_mc += ip.integral +print "integral of mc events is " + str(integral_mc) +print "integral of data events is " + str(integral_data) +print "data/mc = " +str(integral_data/integral_mc) + +#for ip in process.collection: +# ip.draw_hists() +#process.hists_data.draw_hists() + +process.print_stacks() \ No newline at end of file diff --git a/DeepNtuplizer/plots/plot_event_ttsemilep.py b/DeepNtuplizer/plots/plot_event_ttsemilep.py new file mode 100644 index 00000000000..db531c729f6 --- /dev/null +++ b/DeepNtuplizer/plots/plot_event_ttsemilep.py @@ -0,0 +1,613 @@ +from __future__ import print_function +from __future__ import division + + +import pdb +from ROOT import TH1F, TFile, TCanvas, gROOT, gStyle, THStack, TLegend, TGaxis, gPad, TPad, TGraph, TLine +import numpy as np +import os +from array import array +import datetime +from DataFormats.FWLite import Events, Handle + +gROOT.SetBatch() # don't pop up canvases +gROOT.SetStyle('Plain') # white background +gStyle.SetFillStyle(0) +TGaxis.SetMaxDigits(4) # Force scientific notation for numbers with more than 4 digits +#gStyle.SetBorderSize(0) +#gStyle.SetTextFont(43) +#gStyle.SetTextSize(14) + + +class pileupweight: + """ DOC """ + + def __init__(self, file_pileup_MC, file_pileup_Data): + hist_mc = file_pileup_MC.Get("pileup") + hist_data = file_pileup_Data.Get("pileup") + self.pupWeights = [0] + hist_mc.Scale(1. / hist_mc.Integral()) + hist_data.Scale(1. / hist_data.Integral()) + for bin in range(1, hist_data.GetNbinsX() + 2): + if hist_mc.GetBinContent(bin) != 0: + w = hist_data.GetBinContent(bin) / hist_mc.GetBinContent(bin) + self.pupWeights.append(w) + else: + self.pupWeights.append(1) + + def getPileupWeight(self,nInteractions): + if nInteractions >= len(self.pupWeights): + return 0 + return self.pupWeights[nInteractions] + +class scaleFactor(): + """ + class for isolation and id scale factors + hist is the 2d histogram with the scale factor with abs(eta) on x and pt on y + """ + def __init__(self,hist): + self.hist = hist + self.xaxis = self.hist.GetXaxis() + self.yaxis = self.hist.GetYaxis() + self.binx = 0 + self.biny = 0 + self.scalefactor = 1. + self.maxBinX = self.xaxis.GetLast() + self.maxBinY = self.yaxis.GetLast() + self.minX = self.xaxis.GetXmin() + self.maxX = self.xaxis.GetXmax() + self.minY = self.yaxis.GetXmin() + self.maxY = self.yaxis.GetXmax() + + def getScalefactor(self,x,y): + + self.binx = self.xaxis.FindBin(x) + self.biny = self.yaxis.FindBin(y) + if x >= self.maxX: self.binx = self.maxBinX + if y >= self.maxY: self.biny = self.maxBinY + self.scalefactor = self.hist.GetBinContent(self.binx, self.biny) + return self.scalefactor + + +class scaleFactor_tracking(): + """ + class for tracking scale factors + the tracking scale factors are saved in a TGraphAsymmErrors object which has to be written in a histogram first + """ + def __init__(self,graph, name): + self.graph = graph + self.npoints = self.graph.GetN() + self.x_centers = self.graph.GetX() + self.y_centers = self.graph.GetY() + self.x_lows = np.zeros(self.npoints) + self.x_highs = np.zeros(self.npoints) + self.y_lows = np.zeros(self.npoints) + self.y_highs = np.zeros(self.npoints) + self.x_edges = np.zeros(self.npoints +1) + + for i in range(0,self.npoints): + self.x_lows[i] = self.graph.GetErrorXlow(i) + self.x_highs[i] = self.graph.GetErrorXhigh(i) + self.y_lows[i] = self.graph.GetErrorYlow(i) + self.y_lows[i] = self.graph.GetErrorYhigh(i) + + self.x_edges[i] = self.x_centers[i] - self.x_lows[i] + self.x_edges[self.npoints] = self.x_centers[self.npoints - 1] + self.x_highs[self.npoints-1] + self.hist = TH1F(name, name,self.npoints,self.x_edges) + + for i in range(0,self.npoints): + self.hist.SetBinContent(i+1, self.y_centers[i]) + self.hist.SetBinError(i+1, max(self.y_lows[i], self.y_highs[i])) + + self.bin = 0 + self.axis = self.hist.GetXaxis() + self.scalefactor = 1. + + def getScalefactor(self,x): + self.bin = self.axis.FindBin(abs(x)) + self.scalefactor = self.hist.GetBinContent(self.bin) + return self.scalefactor + + +class process: + """ DOC """ + datadir = "/afs/desy.de/user/d/dwalter/CMSSW_8_0_29/src/DeepNTuples/DeepNtuplizer/data/" + pupw = pileupweight( + TFile(datadir + "MyMCPileupHist.root"), + TFile(datadir + "MyDataPileupHist.root") + ) + + lumi_BCDEF = 5.404 + 2.396 + 4.256 + 4.054 + 3.105 + lumi_GH = 7.179 + 8.746 + + # muon scalefactor hists from https://twiki.cern.ch/twiki/bin/viewauth/CMS/MuonWorkInProgressAndPagResults + file_sf_muon_trigger_GH = TFile(datadir + "EfficienciesAndSF_Period4.root") + file_sf_muon_id_GH = TFile(datadir + "EfficienciesAndSF_ID_GH.root") + file_sf_muon_iso_GH = TFile(datadir + "EfficienciesAndSF_ISO_GH.root") + file_sf_muon_trigger_BCDEF = TFile(datadir + "EfficienciesAndSF_RunBtoF.root") + file_sf_muon_id_BCDEF = TFile(datadir + "EfficienciesAndSF_ID_BCDEF.root") + file_sf_muon_iso_BCDEF = TFile(datadir + "EfficienciesAndSF_ISO_BCDEF.root") + file_sf_muon_tracking = TFile(datadir + "Tracking_EfficienciesAndSF_BCDEFGH.root") + + h_sf_muon_trigger_BCDEF = file_sf_muon_trigger_BCDEF.Get("IsoMu24_OR_IsoTkMu24_PtEtaBins/abseta_pt_ratio") + h_sf_muon_id_BCDEF = file_sf_muon_id_BCDEF.Get("MC_NUM_TightID_DEN_genTracks_PAR_pt_eta/abseta_pt_ratio") + h_sf_muon_iso_BCDEF = file_sf_muon_iso_BCDEF.Get("TightISO_TightID_pt_eta/abseta_pt_ratio") + h_sf_muon_trigger_GH = file_sf_muon_trigger_GH.Get("IsoMu24_OR_IsoTkMu24_PtEtaBins/abseta_pt_ratio") + h_sf_muon_id_GH = file_sf_muon_id_GH.Get("MC_NUM_TightID_DEN_genTracks_PAR_pt_eta/abseta_pt_ratio") + h_sf_muon_iso_GH = file_sf_muon_iso_GH.Get("TightISO_TightID_pt_eta/abseta_pt_ratio") + g_sf_muon_tracking = file_sf_muon_tracking.Get("ratio_eff_aeta_dr030e030_corr") + + sf_muon_trigger_BCDEF = scaleFactor(h_sf_muon_trigger_BCDEF) + sf_muon_id_BCDEF = scaleFactor(h_sf_muon_id_BCDEF) + sf_muon_iso_BCDEF = scaleFactor(h_sf_muon_iso_BCDEF) + sf_muon_trigger_GH = scaleFactor(h_sf_muon_trigger_GH) + sf_muon_id_GH = scaleFactor(h_sf_muon_id_GH) + sf_muon_iso_GH = scaleFactor(h_sf_muon_iso_GH) + sf_muon_tracking = scaleFactor_tracking(g_sf_muon_tracking,"sf_muon_tracking") + + muonhandle = Handle("vector") + muonlabel = ("goodMuons") + + methandle = Handle("vector") + metlabel = ("slimmedMETs") + jethandle = Handle('vector') + jetlabel = ('slimmedJets') +#jetlabel = ("GoodJets") + ### LHE weights + LHEweighthandle = Handle('LHEEventProduct') + LHEweightlabel = ("externalLHEProducer") + + ## Pileup reweighting + puphandle = Handle('vector') + puplabel = ("slimmedAddPileupInfo") + + hs_muon_pt = THStack("hs_muon_pt", "muon pt") + hs_muon_eta = THStack("hs_muon_eta", "muon eta") + hs_mtW = THStack("hs_mtW", "mtW") + hs_ptmiss = THStack("hs_ptmiss", "missing pt") + hs_j_n = THStack("hs_njets", "jet number") + + hs_list = hs_muon_pt, hs_muon_eta, hs_mtW, hs_ptmiss, hs_j_n + + collection = [] + hists_data = None + + tsize = 20 + tfont = 43 + + def __init__(self, name, isData=False, color=1): + self.name = name + self.isData = isData + self.h_muon_pt = TH1F(self.name + "_muon_pt", "pt of muon in " + self.name, 9, array('d', [0, 30, 35, 40, 45, 50, 60, 70, 100, 200])) + self.h_muon_eta = TH1F(self.name + "_muon_eta", "eta of the muon in " + self.name, 20, -2.4, 2.4) + self.h_mtW = TH1F(self.name + "_mtW", "reconstructed transverse mass of W " + self.name, 10, array('d', [0, 50, 60, 70, 80, 90, 100, 110, 120, 150, 200])) + self.h_ptmiss = TH1F(self.name + "_ptmiss", "missing transverse momentum in " + self.name, 11, array('d', [0, 10, 20, 30, 40, 50, 60, 70, 100, 120, 150, 200])) + self.h_njets = TH1F(self.name + "_njets", "number of jets in " +self.name, 20, -0.5, 19.5) + + self.hists = self.h_muon_pt, self.h_muon_eta, self.h_mtW, self.h_ptmiss, self.h_njets + + self.integral = 0. + + #remove status boxes + for hist in self.hists: + hist.SetStats(False) + + if self.isData==False: + process.collection.append(self) + for hist in self.hists: + hist.SetFillColor(color) + else: + process.hists_data = self + + self.h_muon_pt.SetXTitle("muon p_{T} (GeV)") + self.h_muon_eta.SetXTitle("muon eta") + self.h_mtW.SetXTitle("m_{T}^{W} (GeV)") + self.h_ptmiss.SetXTitle("p_{t}^{miss} (GeV)") + self.h_njets.SetXTitle("N_{Jet}") + + self.h_muon_pt.SetYTitle("Events/GeV") + self.h_muon_pt.GetYaxis().SetTitleOffset(1.2) + self.h_muon_eta.SetYTitle("Events/0.25") + self.h_mtW.SetYTitle("Events/GeV") + self.h_ptmiss.SetYTitle("Events/GeV") + self.h_njets.SetYTitle("Events") + + self.h_njets.GetXaxis().SetNdivisions(10,0,0) + + def fillHists(self, infile, islist=True, useLHEWeights=False, isData=False, eventIDs=[]): + print("start to fill hist: " + str(datetime.datetime.now())) + nevents = 0 + njets = 0 + directory = '' + + if islist is True: + print("loop over filelist " + str(infile)) + filelist = open(infile, 'r') + ifile = filelist.readline() + ifile = ifile[:-1] + directory = infile[:-len(infile.split('/')[-1])] + else: + ifile = infile + + while ifile != '': + print("process file ", ifile) + events = Events(directory + ifile) + + ### EVENT LOOP + for event in events: + eventID = event.eventAuxiliary().event() + + #if isData: # avoid double counting + # if eventID in eventIDs: + # # print("this event was already processed") + # continue + # eventIDs = np.append(eventIDs, eventID) + + nevents += 1 + njets += 4 + + + weight = 1. + weight_BCDEF = 1. + weight_GH = 1. + + # Fill muon hists + event.getByLabel(process.muonlabel, process.muonhandle) + muons = process.muonhandle.product() # get the product + muon = muons[0] + + muonPt = muon.pt() + muonEta = muon.eta() + + + # Fill mtw hists + event.getByLabel(process.metlabel, process.methandle) + mets = process.methandle.product() + met = mets[0] + + mtW = np.sqrt(2*(muon.pt()*met.pt() - muon.px()*met.px() - muon.py()*met.py())) + + # Fill jet hists + event.getByLabel(process.jetlabel, process.jethandle) + jets = process.jethandle.product() + + numJets = len(jets) + + # Compute Event weights + if useLHEWeights == True: + event.getByLabel(process.LHEweightlabel, process.LHEweighthandle) + weights = process.LHEweighthandle.product() + weight = weights.weights()[0].wgt / abs(weights.weights()[0].wgt) + + if isData == False: + event.getByLabel(process.puplabel, process.puphandle) + pupInfos = process.puphandle.product() + for pupInfo in pupInfos: + if (pupInfo.getBunchCrossing() == 0): + weight *= process.pupw.getPileupWeight(int(pupInfo.getTrueNumInteractions())) + + weight_BCDEF *= process.sf_muon_trigger_BCDEF.getScalefactor(abs(muonEta), muonPt) + weight_BCDEF *= process.sf_muon_id_BCDEF.getScalefactor(abs(muonEta), muonPt) + weight_BCDEF *= process.sf_muon_iso_BCDEF.getScalefactor(abs(muonEta), muonPt) + weight_GH *= process.sf_muon_trigger_GH.getScalefactor(abs(muonEta), muonPt) + weight_GH *= process.sf_muon_id_GH.getScalefactor(abs(muonEta), muonPt) + weight_GH *= process.sf_muon_iso_GH.getScalefactor(abs(muonEta), muonPt) + weight *= process.sf_muon_tracking.getScalefactor(abs(muonEta)) + + weight = weight * (weight_BCDEF * process.lumi_BCDEF + weight_GH * process.lumi_GH) + + + self.h_muon_pt.Fill(muonPt, weight) + self.h_muon_eta.Fill(muonEta, weight) + self.h_mtW.Fill(mtW, weight) + self.h_ptmiss.Fill(met.pt(), weight) + self.h_njets.Fill(numJets, weight) + + if islist is True: + ifile = filelist.readline() + ifile = ifile[:-1] + else: + break + + if islist is True: + print("close filelist") + filelist.close() + + print("processed events: ", nevents) + print("number of good jets: ", njets) + self.save_hists() + + def scale_hists(self,sigma=1.,lumi=1.,eff=1.): + for hist in self.hists: + hist.Scale(sigma*lumi*eff) + self.integral = self.hists[0].Integral() + + self.divideByBinWidth(self.hists[0]) + self.divideByBinWidth(self.hists[2]) + self.divideByBinWidth(self.hists[3]) + self.hists[1].Scale(4) + + def divideByBinWidth(self,hist): + for ibin in range(0, len(hist)): + content = hist.GetBinContent(ibin) + content /= hist.GetBinWidth(ibin) + hist.SetBinContent(ibin, content) + + def save_hists(self): + f1 = TFile(self.name + ".root", "RECREATE") + for hist in self.hists: + hist.Write() + f1.Close() + + def load_hists(self): + print("load " + self.name) + + if(os.path.isfile(self.name + ".root")): + f1 = TFile(self.name + ".root", "READ") + self.h_muon_pt.Add(f1.Get(self.name + "_muon_pt")) + self.h_muon_eta.Add(f1.Get(self.name + "_muon_eta")) + self.h_mtW.Add(f1.Get(self.name + "_mtW")) + self.h_ptmiss.Add(f1.Get(self.name + "_ptmiss")) + self.h_njets.Add(f1.Get(self.name + "_njets")) + + f1.Close() + else: + print("no file for this process") + + def draw_hists(self): + canvas = TCanvas(self.name, "", 800, 600) + for hist in self.hists: + canvas.Clear() + hist.Draw("HIST") + canvas.Print(hist.GetName()+".png") + + def setMinMax(self, y_min_muon_pt, y_max_muon_pt, y_min_muon_eta, y_max_muon_eta, y_min_mtW, y_max_mtW, y_min_ptmiss, y_max_ptmiss, y_min_njets, y_max_njets): + self.hists[0].SetMinimum(y_min_muon_pt) + self.hists[0].SetMaximum(y_max_muon_pt) + self.hists[1].SetMinimum(y_min_muon_eta) + self.hists[1].SetMaximum(y_max_muon_eta) + self.hists[2].SetMinimum(y_min_mtW) + self.hists[2].SetMaximum(y_max_mtW) + self.hists[3].SetMinimum(y_min_ptmiss) + self.hists[3].SetMaximum(y_max_ptmiss) + self.hists[4].SetMinimum(y_min_njets) + self.hists[4].SetMaximum(y_max_njets) + + @classmethod + def make_stacks(cls): + for proc in cls.collection: + for i,hist in enumerate(proc.hists): + cls.hs_list[i].Add(hist) + + + @classmethod + def print_stacks(cls): + cls.make_stacks() + + leg = TLegend(0.59, 0.54, 0.89, 0.84) + leg.SetBorderSize(0) + leg.SetTextFont(cls.tfont) + leg.AddEntry(cls.hists_data.hists[0], "data", "ep") + leg.AddEntry(cls.collection[8].hists[0], "t#bar{t}", "f") + leg.AddEntry(cls.collection[4].hists[0], "single t", "f") + leg.AddEntry(cls.collection[2].hists[0], "W+Jets", "f") + #leg.AddEntry(cls.collection[2].hists[0], "VV", "f") + leg.AddEntry(cls.collection[0].hists[0], "DY", "f") + + leg21 = TLegend(0.19, 0.69, 0.39, 0.84) + leg21.SetBorderSize(0) + leg21.SetTextFont(cls.tfont) + leg21.AddEntry(cls.hists_data.hists[0], "data", "ep") + leg21.AddEntry(cls.collection[8].hists[0], "t#bar{t}", "f") + leg21.AddEntry(cls.collection[4].hists[0], "single t", "f") + + leg22 = TLegend(0.69, 0.69, 0.89, 0.84) + leg22.SetBorderSize(0) + leg22.SetTextFont(cls.tfont) + leg22.AddEntry(cls.collection[2].hists[0], "W+Jets", "f") + #leg22.AddEntry(cls.collection[2].hists[0], "VV", "f") + leg22.AddEntry(cls.collection[0].hists[0], "DY", "f") + + leg.SetFillStyle(0) # make transparent + leg21.SetFillStyle(0) + leg22.SetFillStyle(0) + + canvas = TCanvas("c1","Example 2 pads (20,80)",200,10,800,600) + #canvas.SetBottomMargin(0.25) + pad1 = TPad("pad1", "The pad 80% of the height", 0, 0.2, 1, 1.0) + pad2 = TPad("pad2", "The pad 20% of the height", 0, 0.05, 1, 0.2) + + pad1.SetLeftMargin(0.15) + pad1.SetRightMargin(0.1) + pad2.SetLeftMargin(0.15) + pad2.SetRightMargin(0.1) + + pad1.SetBottomMargin(0.02) + pad2.SetTopMargin(0.02) + pad2.SetBottomMargin(0.25) + + pad1.Draw() + pad2.Draw() + for i, h_data in enumerate(cls.hists_data.hists): + + h_data.SetTitle("") + h_data.GetXaxis().SetLabelFont(cls.tfont) + h_data.GetXaxis().SetLabelSize(0) + h_data.GetXaxis().SetTitleSize(0) + h_data.GetYaxis().SetLabelFont(cls.tfont) + h_data.GetYaxis().SetLabelSize(cls.tsize) + h_data.GetYaxis().SetTitleFont(cls.tfont) + h_data.GetYaxis().SetTitleSize(cls.tsize) + h_data.GetYaxis().SetTitleOffset(1.5) + + h_data.SetMarkerStyle(20) + #pad1.Clear() + #pad2.Clear() + pad1.cd() + + h_data.Draw("E1") + cls.hs_list[i].Draw("HIST same") + h_data.Draw("E1 same") + gPad.RedrawAxis() # draw axis in foreground + + if i != 1: + leg.Draw("same") + else: + leg21.Draw("same") + leg22.Draw("same") + + pad2.cd() + h_ratio = TH1F(cls.hs_list[i].GetStack().Last()) + + h_ratio.SetTitle("") + h_ratio.GetXaxis().SetLabelFont(cls.tfont) + h_ratio.GetXaxis().SetLabelSize(cls.tsize) + h_ratio.GetYaxis().SetLabelFont(cls.tfont) + h_ratio.GetYaxis().SetLabelSize(cls.tsize) + h_ratio.GetXaxis().SetTitleFont(cls.tfont) + h_ratio.GetXaxis().SetTitleSize(cls.tsize) + h_ratio.GetYaxis().SetTitleFont(cls.tfont) + h_ratio.GetYaxis().SetTitleSize(cls.tsize) + h_ratio.GetXaxis().SetTitleOffset(7.) + h_ratio.GetYaxis().SetTitleOffset(1.5) + h_ratio.GetYaxis().SetNdivisions(403) + h_ratio.GetXaxis().SetTickLength(0.15) + + if i ==4: + h_ratio.GetXaxis().SetNdivisions(10, 0, 0) + + h_ratio.SetMarkerStyle(20) + h_ratio.SetXTitle(h_data.GetXaxis().GetTitle()) + + h_ratio.Divide(h_data) + h_ratio.GetYaxis().SetTitle("mc/data") + + h_ratio.SetMinimum(0.75) + h_ratio.SetMaximum(1.25) + + + middleLine = TLine(h_ratio.GetXaxis().GetXmin(),1. ,h_ratio.GetXaxis().GetXmax(), 1.) + + h_ratio.Draw("e1") + middleLine.Draw("same") + gStyle.SetErrorX(0.) + canvas.Print("stacks_" + str(i) + ".png") + + + +dy10to50 = process("dy10to50", color=593) +dy50 = process("dy50", color=593) +#ww = process("ww", color=393) +#wz = process("wz", color=393) +#zz = process("zz", color=393) +#wjets = process("wjets", color=409) +#w1jets = process("w1jets",color=409) +#w2jets = process("w2jets",color=409) +w3jets = process("w3jets",color=409) +w4jets = process("w4jets",color=409) +wt = process("wt", color=609) +wantit = process("wantit", color=609) +t_tchannel = process("t_tchannel", color=609) +antit_tchannel = process("antit_tchannel", color=609) +tt = process("tt", color=625) +data = process("data",isData=True) + +#for ip in process.collection: +# ip.load_hists() +#process.hists_data.load_hists() + +dy10to50.fillHists("dy10to50.txt",useLHEWeights=True) +dy50.fillHists("dy50.txt",useLHEWeights=True) +#ww.fillHists("ww.txt") +#wz.fillHists("wz.txt") +#zz.fillHists("zz.txt") +#wjets.fillHists("wjets.txt",useLHEWeights=True) +#w1jets.fillHists("w1jets.txt") +#w2jets.fillHists("w2jets.txt") +w3jets.fillHists("w3jets.txt") +w4jets.fillHists("w4jets.txt") +wt.fillHists("wt.txt") +wantit.fillHists("wantit.txt") +t_tchannel.fillHists("t_tchannel.txt") +antit_tchannel.fillHists("antit_tchannel.txt") +tt.fillHists("tt.txt") +data.fillHists("data.txt",isData=True) + + + +sigma_tt = 831760 +sigma_dy50 = 6025200 +sigma_dy10to50 = 22635100 +sigma_ww = 118700 +sigma_wz = 44900 +sigma_zz = 15400 +sigma_wantit = 35600 +sigma_wt = 35600 +sigma_t_tchannel = 136020 +sigma_antit_tchannel = 80950 +sigma_wjets = 61526700 +sigma_w1jets = 9625000 +sigma_w2jets = 3161000 +sigma_w3jets = 948200 #+-423 computed with https://twiki.cern.ch/twiki/bin/view/CMS/HowToGenXSecAnalyzer#Automated_scripts_to_compute_the +sigma_w4jets = 494600 + +# Efficiency = 1/(number of events before cuts), in case of lheweights its the +# effective number of events (in this case: number of events with positive weights - number of events with negative weights) +eff_tt = 1./(77081156 + 77867738) +eff_dy50 = 1./81781052 # effective number of events +#eff_dy50 = 1./122055388 # true number of events +eff_dy10to50 = 1./47946519 # effective number of events +#eff_dy10to50 = 1./65888233 # true number of events +eff_ww = 1./994012 +eff_wz = 1./1000000 +eff_zz = 1./998034 +eff_wantit = 1./6933094 +eff_wt = 1./6952830 +eff_t_tchannel = 1./67240808 +eff_antit_tchannel = 1./38811017 +eff_wjets = 1./16497031 # effective number of events +eff_w1jets = 1./45367044 +eff_w2jets = 1./29878415 +eff_w3jets = 1./19798117 +eff_w4jets = 1./9170576 +#eff_wjets = 1./24120319 # true number of events + +data.scale_hists() +tt.scale_hists(sigma = sigma_tt, eff = eff_tt) +dy50.scale_hists(sigma = sigma_dy50, eff = eff_dy50) +dy10to50.scale_hists(sigma = sigma_dy10to50, eff = eff_dy10to50) +#ww.scale_hists(sigma = sigma_ww, eff = eff_ww) +#wz.scale_hists(sigma = sigma_wz, eff = eff_wz) +#zz.scale_hists(sigma = sigma_zz, eff = eff_zz) +wantit.scale_hists(sigma = sigma_wantit, eff = eff_wantit) +wt.scale_hists(sigma = sigma_wt, eff = eff_wt) +t_tchannel.scale_hists(sigma = sigma_t_tchannel, eff = eff_t_tchannel) +antit_tchannel.scale_hists(sigma = sigma_antit_tchannel, eff = eff_antit_tchannel) +#wjets.scale_hists(sigma = sigma_wjets, eff = eff_wjets) +#w1jets.scale_hists(sigma = sigma_w1jets, eff = eff_w1jets) +#w2jets.scale_hists(sigma = sigma_w2jets, eff = eff_w2jets) +w3jets.scale_hists(sigma = sigma_w3jets, eff = eff_w3jets) +w4jets.scale_hists(sigma = sigma_w4jets, eff = eff_w4jets) + +data.setMinMax(0,20000,0,300000,0,18500,0,15000,0,140000) +title='plots_180515_mainContributions' + +directory = os.path.dirname('./'+title+'/') +# make a canvas, draw, and save it +if not os.path.exists(directory): + os.makedirs(directory) +os.chdir(directory) + +integral_mc = 0. +integral_data = process.hists_data.integral +for ip in process.collection: + integral_mc += ip.integral +print("integral of mc events is " + str(integral_mc)) +print("integral of data events is " + str(integral_data)) +print("data/mc = " +str(integral_data/integral_mc)) + +#for ip in process.collection: +# ip.draw_hists() +#process.hists_data.draw_hists() + +process.print_stacks() \ No newline at end of file diff --git a/DeepNtuplizer/plots/plot_features.py b/DeepNtuplizer/plots/plot_features.py new file mode 100644 index 00000000000..1a77acb38e4 --- /dev/null +++ b/DeepNtuplizer/plots/plot_features.py @@ -0,0 +1,275 @@ +####### +# plot distributions of different featues of deepntuples root tree +# add root files with ntuples to files_mc or files_data +# itinitialize feature to plot and bin range +# for twodimensional features you have to choose 'all' if you want to plot all flattened or 'max' to plot the maximum value of each tuple +# +####### + +import pdb +import ROOT +import os +from array import array + +import matplotlib.pyplot as plt +import matplotlib.patches as mpatches +import root_numpy as rn +import numpy as np + +ROOT.gROOT.SetBatch() # don't pop up canvases +ROOT.gROOT.SetStyle('Plain') # white background +ROOT.gStyle.SetFillStyle(0) + +def isInt(a): + try: + int(a) + return True + except: + return False + +def hist_bin_uncertainty(data, weights, bin_edges): + """ + The statistical uncertainity per bin of the binned data. + If there are weights then the uncertainity will be the root of the + sum of the weights squared. + If there are no weights (weights = 1) this reduces to the root of + the number of events. + Args: + data: `array`, the data being histogrammed. + weights: `array`, the associated weights of the `data`. + bin_edges: `array`, the edges of the bins of the histogram. + Returns: + bin_uncertainties: `array`, the statistical uncertainity on the bins. + """ + # Bound the data and weights to be within the bin edges + in_range_index = [idx for idx in range(len(data)) + if data[idx] > min(bin_edges) and data[idx] < max(bin_edges)] + in_range_data = np.asarray([data[idx] for idx in in_range_index]) + + if weights is None or np.array_equal(weights, np.ones(len(weights))): + # Default to weights of 1 and thus uncertainty = sqrt(N) + in_range_weights = np.ones(len(in_range_data)) + else: + in_range_weights = np.asarray([weights[idx] for idx in in_range_index]) + + # Bin the weights with the same binning as the data + bin_index = np.digitize(in_range_data, bin_edges) + # N.B.: range(1, bin_edges.size) is used instead of set(bin_index) as if + # there is a gap in the data such that a bin is skipped no index would appear + # for it in the set + binned_weights = np.asarray( + [in_range_weights[np.where(bin_index == idx)[0]] for idx in range(1, len(bin_edges))]) + bin_uncertainties = np.asarray( + [np.sqrt(np.sum(np.square(w))) for w in binned_weights]) + return bin_uncertainties + + +class feature: + """ DOC """ + + + def __init__(self, name, bins,ymax=None,fillopt='', xlabel='',description=''): + self.name = name + self.bins = bins + self.fillopt=fillopt + self.xlabel=xlabel + self.description=description + self.ymax=ymax + + +class datacollection: + """ DOC """ + + + def __init__(self, filenames_mc, filenames_data = None): + + self.features = [] + + self.scalefactor_data = 1. + self.filenames_mc = filenames_mc + self.filenames_data = filenames_data + + + #pdb.set_trace() + + + def fill(self): + + print 'load samples' + branches = [] + for ifeature in self.features: + branches.append(ifeature.name) + self.collection_mc = rn.root2array(self.filenames_mc, "deepntuplizer/tree", list(set(branches))) + self.weight_mc = rn.root2array(self.filenames_mc, "deepntuplizer/tree", "event_weight") + if self.filenames_data != None: + self.collection_data = rn.root2array(self.filenames_data, "deepntuplizer/tree", list(set(branches))) + self.weight_data = rn.root2array(self.filenames_data, "deepntuplizer/tree", "event_weight") + + #pdb.set_trace() + + + def add_feature(self,fname, bins, ymax=None, fillopt='',xlabel='',description=''): + self.features.append(feature(fname,bins,ymax,fillopt,xlabel,description)) + + + def scale_data(self,sf): + self.scalefactor_data *= sf + + + def plot_features(self): + for ifeature in self.features: + print 'make plot for feature '+ifeature.name + + if 'all' in ifeature.fillopt: + flat_col_mc = np.concatenate(self.collection_mc[ifeature.name], axis=-1) + freqlist = [len(j) for j in self.collection_mc[ifeature.name]] + flat_weight_mc = np.repeat(self.weight_mc, freqlist) + + flat_col_data = np.concatenate(self.collection_data[ifeature.name], axis=-1) + freqlist = [len(j) for j in self.collection_data[ifeature.name]] + flat_weight_data = np.repeat(self.weight_data, freqlist) + + entries_mc, binEdges = np.histogram(flat_col_mc, ifeature.bins, weights=flat_weight_mc) + entries_data, _ = np.histogram(flat_col_data, ifeature.bins, weights=flat_weight_data) + err_mc = hist_bin_uncertainty(flat_col_mc, flat_weight_mc, binEdges) + err_data = hist_bin_uncertainty(flat_col_data, flat_weight_data, binEdges) + + if 'max' in ifeature.fillopt: + max_col_mc = np.empty((self.collection_mc[ifeature.name].shape)) + max_col_data = np.empty((self.collection_data[ifeature.name].shape)) + for j,arr in enumerate(self.collection_mc[ifeature.name]): + if len(arr) > 0: + max_col_mc[j] = np.amax(arr) + else: + max_col_mc[j] = 0. + + for j, arr in enumerate(self.collection_data[ifeature.name]): + if len(arr) > 0: + max_col_data[j] = np.amax(arr) + else: + max_col_data[j] = 0. + entries_mc, binEdges = np.histogram(max_col_mc, ifeature.bins, weights=self.weight_mc) + entries_data, _ = np.histogram(max_col_data, ifeature.bins, weights=self.weight_data) + err_mc = hist_bin_uncertainty(max_col_mc, self.weight_mc, binEdges) + err_data = hist_bin_uncertainty(max_col_data, self.weight_data, binEdges) + + if isInt(ifeature.fillopt): + ientry = int(ifeature.fillopt) + col_mc = np.empty((self.collection_mc[ifeature.name].shape)) + col_data = np.empty((self.collection_data[ifeature.name].shape)) + for j, arr in enumerate(self.collection_mc[ifeature.name]): + if len(arr) > ientry: + col_mc[j] = arr[ientry] + else: + col_mc[j] = -100 + + for j, arr in enumerate(self.collection_data[ifeature.name]): + if len(arr) > ientry: + col_data[j] = arr[ientry] + else: + col_data[j] = -100 + entries_mc, binEdges = np.histogram(col_mc, ifeature.bins, weights=self.weight_mc) + entries_data, _ = np.histogram(col_data, ifeature.bins, weights=self.weight_data) + err_mc = hist_bin_uncertainty(col_mc, self.weight_mc, binEdges) + err_data = hist_bin_uncertainty(col_data, self.weight_data, binEdges) + + if ifeature.fillopt == '': + entries_mc, binEdges = np.histogram(self.collection_mc[ifeature.name], ifeature.bins, weights=self.weight_mc) + entries_data, _ = np.histogram(self.collection_data[ifeature.name], ifeature.bins, weights=self.weight_data) + err_mc = hist_bin_uncertainty(self.collection_mc[ifeature.name],self.weight_mc,binEdges) + err_data = hist_bin_uncertainty(self.collection_data[ifeature.name], self.weight_data, binEdges) + + plt.cla() + fig, ax = plt.subplots() + + bincenters = 0.5 * (binEdges[1:] + binEdges[:-1]) + binWidth = binEdges[1:] - binEdges[:-1] + ax.bar(bincenters, entries_mc, yerr = err_mc, width = binWidth, align='center', color='w', edgecolor='k',ecolor='k', label='mc') + ax.errorbar(bincenters,entries_data*self.scalefactor_data, yerr=err_data*self.scalefactor_data, fmt='k.',ecolor='k', label='data') + ax.legend(loc='best',numpoints=1, fontsize=20) + ax.set_xlim(binEdges[0], binEdges[-1]) + ax.set_ylabel('weighted frequency', fontsize=20) + ax.set_xlabel(ifeature.xlabel, fontsize=20) + ax.set_title(ifeature.description) + ax.set_ylim(bottom=0) + if ifeature.ymax != None: + ax.set_ylim(top=ifeature.ymax) + x_min, x_max = ax.get_xlim() + y_min, y_max = ax.get_ylim() + #ax.text(x_min + 0.01*(x_max-x_min),y_max*0.95,ifeature.description) + + ax.ticklabel_format(style='sci', scilimits=(-3, 4), axis='both') + plt.gcf().subplots_adjust(bottom=0.15) + + plt.savefig(ifeature.name+"_"+ifeature.fillopt+".png") + + + + +#files_mc = ["tt_1_0.root","tt_2_0.root","tt_3_0.root","tt_4_0.root","tt_5_0.root", +# "dy50_1_0.root","dy10to50_1_0.root","wantit_1_0.root","wt_1_0.root", +# "ww_1_0.root","wz_1_0.root","zz_1_0.root","wjets_1_0.root"] +files_mc = ["tt.root","tt_backup.root","tt_backup_2.root", + "dy50.root","dy10to50.root","wantit.root","wt.root", + "ww.root","wz.root","zz.root","wjets.root"] +files_data = "muonEG_H_0.root" + +dc = datacollection(files_mc,files_data) + +#dc.add_feature("jet_pt", np.arange(0,390,15), xlabel='$p_\mathrm{T}(j)$ (GeV)', description='') +#dc.add_feature("jet_eta", np.arange(-2.4,2.6,0.2), xlabel='$\eta(j)$', description='') +#dc.add_feature("nCpfcand", np.arange(-0.5,30.5,1), xlabel='$N_\mathrm{cPF}$', description='') +#dc.add_feature("nNpfcand", np.arange(-0.5,30.5,1), xlabel='$N_\mathrm{nPF}$', description='') +#dc.add_feature("npv", np.arange(-0.5,50.5,1), xlabel='$N_\mathrm{PV}$', description='') +#dc.add_feature("nsv", np.arange(-0.5,5.5,1), xlabel='$N_\mathrm{SV}$', description='') +#dc.add_feature("TagVarCSV_trackSumJetEtRatio", np.arange(0.,1.04,0.04), xlabel='trackSumJetEtRatio', description='') +#dc.add_feature("TagVarCSV_trackSumJetDeltaR", np.arange(0.,0.26,0.01), xlabel='trackSumJetDeltaR', description='') +#dc.add_feature("TagVarCSV_vertexCategory", np.arange(-0.5,3.5,1), xlabel='vertexCategory', description='') +#dc.add_feature("TagVarCSV_trackSip2dValAboveCharm", np.array((-1.01,-0.99,-0.1,-0.08,-0.06,-0.04,-0.02,0.,0.02,0.04,0.06,0.08,0.1,0.12,0.14,0.16,0.18,0.2,0.22)), xlabel='trackSip2dValAboveCharm', description='') +#dc.add_feature("TagVarCSV_trackSip2dSigAboveCharm", np.arange(-1,21,1), xlabel='trackSip2dSigAboveCharm', description='') +#dc.add_feature("TagVarCSV_trackSip3dValAboveCharm", np.array((-1.01,-0.99,-0.1,-0.08,-0.06,-0.04,-0.02,0.,0.02,0.04,0.06,0.08,0.1,0.12,0.14,0.16,0.18,0.2,0.22)), xlabel='trackSip3dValAboveCharm', description='') +#dc.add_feature("TagVarCSV_trackSip3dSigAboveCharm", np.arange(-1,21,1), xlabel='trackSip3dSigAboveCharm', description='') +#dc.add_feature("TagVarCSV_jetNSelectedTracks", np.arange(-0.5,15.5,1), xlabel='jetNSelectedTracks', description='') +#dc.add_feature("TagVarCSV_jetNTracksEtaRel", np.arange(-0.5,15.5,1), xlabel='jetNTracksEtaRel', description='') + +dc.add_feature("pfDeepCSVJetTags_probb", np.arange(0, 1.05, 0.05), xlabel='pfDeepCSVJetTags_probb', description='') +dc.add_feature("pfDeepCSVJetTags_probbb", np.arange(0, 1.05, 0.05), xlabel='pfDeepCSVJetTags_probbb', description='') +dc.add_feature("pfDeepCSVJetTags_probc", np.arange(0, 1.05, 0.05), xlabel='pfDeepCSVJetTags_probc', description='') +dc.add_feature("pfDeepCSVJetTags_probcc", np.arange(0, 1.05, 0.05), xlabel='pfDeepCSVJetTags_probcc', description='') +dc.add_feature("pfDeepCSVJetTags_probudsg", np.arange(0, 1.05, 0.05), xlabel='pfDeepCSVJetTags_probudsg', description='') + + +#dc.add_feature("Cpfcan_pt", np.arange(0,105,5),fillopt='max',ymax=225000, xlabel='$p_\mathrm{T}(j)$ (GeV)', description='charged particle flow candidate') +#dc.add_feature("Npfcan_pt", np.arange(0,105,5),fillopt='max',xlabel='$p_\mathrm{T}$ (GeV)', description='neutral particle flow candidate') +#dc.add_feature("sv_pt", np.arange(0,110,10),fillopt='max', xlabel='$p_\mathrm{T}$ (GeV)', description='secondary vertice') + +#dc.add_feature("Cpfcan_ptrel", np.arange(-1,-0.616,0.016),fillopt='0',xlabel=r'$\frac{p_\mathrm{T}(cPF)}{p_\mathrm{T}(j)}$', description='first entry') +#dc.add_feature("Npfcan_deltaR", np.arange(-0.6,0.024,0.024),fillopt='0',xlabel='$\Delta R_m(nPF,SV)$', description='first entry') +#dc.add_feature("Cpfcan_drminsv", np.arange(-0.5,0.02,0.02),fillopt='0',xlabel='$\Delta R_m(cPF,SV)$', description='first entry') +#dc.add_feature("sv_d3dsig", np.arange(0,102,2),fillopt='0',xlabel='$S_{3D}(SV)$', description='first entry') + + +#dc.add_feature("Cpfcan_pt", np.arange(0,55,2.5),fillopt='all',ymax=5500000, xlabel='$p_\mathrm{T}$ (GeV)', description='charged particle flow candidate') +#dc.add_feature("Npfcan_pt", np.arange(0,55,2.5),fillopt='all',ymax=6500000, xlabel='$p_\mathrm{T}$ (GeV)', description='neutral particle flow candidate') +#dc.add_feature("sv_pt", np.arange(0,110,10),fillopt='all',ymax=190000, xlabel='$p_\mathrm{T}$ (GeV)', description='secondary vertice') +#dc.add_feature("sv_deltaR", np.arange(-0.5,-0.116,0.016),fillopt='all', xlabel='$\Delta R$', description='secondary vertice') + + + + + +dc.fill() + + +dc.scale_data(1./0.8701594758966494)#*35.9/8.651) + + +title='features' + +directory = os.path.dirname('./plots_'+title+'/') +# make a canvas, draw, and save it +if not os.path.exists(directory): + os.makedirs(directory) +os.chdir(directory) + +dc.plot_features() \ No newline at end of file diff --git a/DeepNtuplizer/plugins/DeepNtuplizer.cc b/DeepNtuplizer/plugins/DeepNtuplizer.cc index 0b1c4296688..6890dd87eae 100644 --- a/DeepNtuplizer/plugins/DeepNtuplizer.cc +++ b/DeepNtuplizer/plugins/DeepNtuplizer.cc @@ -13,6 +13,8 @@ #include "../interface/ntuple_bTagVars.h" #include "../interface/ntuple_FatJetInfo.h" #include "../interface/ntuple_DeepVertex.h" +#include "../interface/ntuple_eventInfo.h" + //ROOT includes #include "TTree.h" @@ -52,7 +54,7 @@ #include "TrackingTools/IPTools/interface/IPTools.h" #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h" - +#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" #if defined( __GXX_EXPERIMENTAL_CXX0X__) #include "TrackingTools/TransientTrack/interface/TransientTrack.h" @@ -108,10 +110,10 @@ DeepNtuplizer::DeepNtuplizer(const edm::ParameterSet& iConfig): vtxToken_(consumes(iConfig.getParameter("vertices"))), svToken_(consumes( iConfig.getParameter("secVertices"))), - jetToken_(consumes >(iConfig.getParameter("jets"))), - puToken_(consumes>(iConfig.getParameter("pupInfo"))), - rhoToken_(consumes(iConfig.getParameter("rhoInfo"))), - t_qgtagger(iConfig.getParameter("qgtagger")) + jetToken_(consumes >(iConfig.getParameter("jets"))), + puToken_(consumes>(iConfig.getParameter("pupInfo"))), + rhoToken_(consumes(iConfig.getParameter("rhoInfo"))), + t_qgtagger(iConfig.getParameter("qgtagger")) { /* * Initialise the modules here @@ -127,13 +129,17 @@ DeepNtuplizer::DeepNtuplizer(const edm::ParameterSet& iConfig): const bool useHerwigCompatibleMatching=iConfig.getParameter("useHerwigCompatible"); const bool isHerwig=iConfig.getParameter("isHerwig"); + const bool isData=iConfig.getParameter("isData"); + ntuple_content::useoffsets = iConfig.getParameter("useOffsets"); applySelection_=iConfig.getParameter("applySelection"); + ntuple_SV* svmodule=new ntuple_SV("", jetR); addModule(svmodule); + //Loose IVF vertices //ntuple_SV* svmodule_LooseIVF=new ntuple_SV("LooseIVF_", jetR); //svmodule_LooseIVF->setSVToken( @@ -156,12 +162,15 @@ DeepNtuplizer::DeepNtuplizer(const edm::ParameterSet& iConfig): jetinfo->setUseHerwigCompatibleMatching(useHerwigCompatibleMatching); jetinfo->setIsHerwig(isHerwig); - jetinfo->setGenJetMatchReclusterToken( - consumes >( - iConfig.getParameter( "genJetMatchRecluster" ))); - jetinfo->setGenJetMatchWithNuToken( - consumes >( - iConfig.getParameter( "genJetMatchWithNu" ))); + if(!isData){ + jetinfo->setGenJetMatchReclusterToken( + consumes >( + iConfig.getParameter( "genJetMatchRecluster" ))); + jetinfo->setGenJetMatchWithNuToken( + consumes >( + iConfig.getParameter( "genJetMatchWithNu" ))); + } + jetinfo->setGenParticlesToken( consumes( @@ -182,7 +191,18 @@ DeepNtuplizer::DeepNtuplizer(const edm::ParameterSet& iConfig): addModule(pfcands); - addModule(new ntuple_bTagVars()); + ntuple_bTagVars * btagvars = new ntuple_bTagVars(); + addModule(btagvars); + + ntuple_eventInfo *evweight = new ntuple_eventInfo(); + if(!isData){ + evweight->setLHEToken(consumes(iConfig.getParameter("lheInfo"))); + evweight->setMuonsToken(consumes >(iConfig.getParameter("sfMuons"))); + evweight->setElectronsToken(consumes >(iConfig.getParameter("sfElectrons"))); + evweight->setTriggerToken(consumes(iConfig.getParameter("triggerToken"))); + } + addModule(evweight); + if(runFatJets_){ auto *fatjetinfo = new ntuple_FatJetInfo(jetR); @@ -196,9 +216,9 @@ DeepNtuplizer::DeepNtuplizer(const edm::ParameterSet& iConfig): * * parse the input parameters (if any) */ - for(auto& m: modules_) + for(auto& m: modules_){ m->getInput(iConfig); - + } } @@ -214,7 +234,6 @@ DeepNtuplizer::~DeepNtuplizer() void DeepNtuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) { - //global info edm::Handle vertices; @@ -225,7 +244,8 @@ DeepNtuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) iEvent.getByToken(svToken_, secvertices); edm::Handle > pupInfo; - iEvent.getByToken(puToken_, pupInfo); + if(!iEvent.isRealData()) + iEvent.getByToken(puToken_, pupInfo); edm::Handle rhoInfo; iEvent.getByToken(rhoToken_,rhoInfo); @@ -236,7 +256,8 @@ DeepNtuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) for(auto& m:modules_){ m->setPrimaryVertices(vertices.product()); m->setSecVertices(secvertices.product()); - m->setPuInfo(pupInfo.product()); + if(!iEvent.isRealData()) + m->setPuInfo(pupInfo.product()); m->setRhoInfo(rhoInfo.product()); m->readSetup(iSetup); m->readEvent(iEvent); @@ -251,13 +272,16 @@ DeepNtuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) edm::View::const_iterator jetIter; // loop over the jets + //for (edm::View::const_iterator jetIter = jets->begin(); jetIter != jets->end(); ++jetIter) { for(size_t j=0;jbegin()+jetidx; const pat::Jet& jet = *jetIter; + if(jet.genJet()) njetswithgenjet_++; @@ -273,6 +297,7 @@ DeepNtuplizer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) njetsselected_++; } } // end of looping over the jets + } diff --git a/DeepNtuplizer/plugins/ElectronIdAdder.cc b/DeepNtuplizer/plugins/ElectronIdAdder.cc new file mode 100644 index 00000000000..2a701e2b6c5 --- /dev/null +++ b/DeepNtuplizer/plugins/ElectronIdAdder.cc @@ -0,0 +1,95 @@ +#include + +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" + +#include "DataFormats/PatCandidates/interface/Electron.h" +#include "DataFormats/VertexReco/interface/Vertex.h" + +class ElectronIdAdder : public edm::EDProducer { + public: + explicit ElectronIdAdder(const edm::ParameterSet& iConfig); + ~ElectronIdAdder(){} + + void produce(edm::Event & iEvent, const edm::EventSetup& iSetup) override; + + private: + edm::EDGetTokenT > src_; + edm::EDGetTokenT > idMap_; + edm::EDGetTokenT > vSrc_; + + edm::Handle > electrons; + edm::Handle > map; + edm::Handle > vertex; + + double minPt_; + double maxAbsEta_; +}; + +ElectronIdAdder::ElectronIdAdder(const edm::ParameterSet& iConfig): + src_(consumes >(iConfig.getParameter("src"))), + idMap_(consumes >(iConfig.getParameter("idMap"))), + vSrc_(consumes > (iConfig.getParameter("vSrc"))){ + produces >(); + + minPt_=(iConfig.getParameter("minPt")); + maxAbsEta_=(iConfig.getParameter("maxAbsEta")); + + } + +void ElectronIdAdder::produce(edm::Event & iEvent, const edm::EventSetup & es) { + + iEvent.getByToken(src_, electrons); + iEvent.getByToken(idMap_, map); + iEvent.getByToken(vSrc_, vertex); + + std::vector * out = new std::vector; + out->reserve(electrons->size()); + + for (size_t i = 0; i < electrons->size(); ++i){ + const auto iElectron = electrons->ptrAt(i); + bool idresult = (*map)[iElectron]; + bool notInCrack = false; + double absd0 = 1.; + double absdz = 1.; + bool inAbsD0 = false; + bool inAbsDz = false; + double global_default_value = 1.5; + double SCeta = (iElectron->superCluster().isAvailable()) ? iElectron->superCluster()->position().eta() : global_default_value; + double absSCeta = fabs(SCeta); + + + if( iElectron->superCluster().isAvailable() ){ + notInCrack = ( absSCeta<1.4442 || absSCeta>1.5660 ); + } + if( iElectron->gsfTrack().isAvailable() ){ + absd0 = fabs(iElectron->gsfTrack()->dxy(vertex->ptrAt(0)->position())); + absdz = fabs(iElectron->gsfTrack()->dz(vertex->ptrAt(0)->position())); + } + if( iElectron->superCluster().isAvailable() ){ + if( absSCeta < 1.4442){ + inAbsD0 = absd0 < 0.05; + inAbsDz = absdz < 0.1; + } + if( absSCeta > 1.5660){ + inAbsD0 = absd0 < 0.1; + inAbsDz = absdz < 0.2; + } + } + + pat::Electron newele(*iElectron); + + if(idresult && notInCrack && inAbsD0 && inAbsDz && iElectron->pt() > minPt_ && std::abs(iElectron->eta())push_back(newele); + } + } + std::unique_ptr > ptr(out); + iEvent.put(std::move(ptr)); + +} + +//define this as a plug-in +DEFINE_FWK_MODULE(ElectronIdAdder); \ No newline at end of file diff --git a/DeepNtuplizer/plugins/MuonIdAdder.cc b/DeepNtuplizer/plugins/MuonIdAdder.cc new file mode 100644 index 00000000000..1c7cae7af94 --- /dev/null +++ b/DeepNtuplizer/plugins/MuonIdAdder.cc @@ -0,0 +1,66 @@ +#include + +#include "FWCore/Framework/interface/EDProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" + +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/VertexReco/interface/Vertex.h" + +class MuonIdAdder : public edm::EDProducer { + public: + explicit MuonIdAdder(const edm::ParameterSet&); + ~MuonIdAdder(){} + void produce(edm::Event & iEvent, const edm::EventSetup& iSetup) override; + + private: + edm::EDGetTokenT > src_; + edm::EDGetTokenT > vSrc_; + + edm::Handle > muons; + edm::Handle > vertex; + + double minPt_; + double maxAbsEta_; + double maxRMI_; +}; + +MuonIdAdder::MuonIdAdder(const edm::ParameterSet& iConfig): + src_(consumes > (iConfig.getParameter("src"))), + vSrc_(consumes > (iConfig.getParameter("vSrc"))){ + produces >(); + + minPt_=(iConfig.getParameter("minPt")); + maxAbsEta_=(iConfig.getParameter("maxAbsEta")); + maxRMI_ = (iConfig.getParameter("maxRMI")); + } + +void MuonIdAdder::produce(edm::Event & iEvent, const edm::EventSetup & iSetup) { + + + iEvent.getByToken(src_, muons); + iEvent.getByToken(vSrc_, vertex); + + std::vector * out = new std::vector; + out->reserve(muons->size()); + + for (size_t i = 0; i < muons->size(); ++i){ + const auto muon = muons->ptrAt(i); + bool idresult = muon->isTightMuon(*(vertex->ptrAt(0))); + double rmi = (muon->pfIsolationR04().sumChargedHadronPt +std::max(0.0, (muon->pfIsolationR04().sumNeutralHadronEt +muon->pfIsolationR04().sumPhotonEt -0.5*muon->pfIsolationR04().sumPUPt)))/muon->pt(); + + pat::Muon newmuon(*muon); + + if(idresult && rmi < maxRMI_ && muon->pt() > minPt_ && std::abs(muon->eta())push_back(newmuon); + } + + } + std::unique_ptr > ptr(out); + iEvent.put(std::move(ptr)); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(MuonIdAdder); \ No newline at end of file diff --git a/DeepNtuplizer/plugins/globalInfo.cc b/DeepNtuplizer/plugins/globalInfo.cc new file mode 100644 index 00000000000..1572d72fe98 --- /dev/null +++ b/DeepNtuplizer/plugins/globalInfo.cc @@ -0,0 +1,157 @@ +// -*- C++ -*- +// +// Package: DeepNTuples/DeepNtuplizer +// Class: globalInfo +// +/**\class globalInfo globalInfo.cc DeepNTuples/DeepNtuplizer/plugins/globalInfo.cc + Description: This analyzer is used to count the (effective) number of initial events. + In the LHE case, the events with negative lhe weights have to be substracted from these with positive lhe weights +*/ +// +// Original Author: David Walter +// Created: Tue, 20 Feb 2018 17:18:02 GMT +// +// + +// system include files +#include +#include +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" +#include "TTree.h" + + + +class globalInfo : public edm::one::EDAnalyzer { + public: + globalInfo(const edm::ParameterSet&); + + ~globalInfo(); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + + private: + virtual void beginJob() override; + virtual void analyze(const edm::Event&, const edm::EventSetup&) override; + virtual void endJob() override; + + edm::Service fs; + edm::EDGetTokenT t_LHEInfo; + edm::Handle h_LHEWeights; + + int nEvents_ = 0; + int nNegLHEEvents_ = 0; + int nEffEvents_ = 0; + bool useLHEWeights_ = true; + + TTree * tree_ = new TTree(); + + + // ----------member data --------------------------- +}; + +// +// constants, enums and typedefs +// + +// +// static data member definitions +// + +// +// constructors and destructor +// +globalInfo::globalInfo(const edm::ParameterSet& iConfig){ + + useLHEWeights_ = (iConfig.getParameter("useLHEWeights")); + if(useLHEWeights_){ + t_LHEInfo = consumes (iConfig.getParameter("lheInfo")); + } + +} + + +globalInfo::~globalInfo() +{ + + // do anything here that needs to be done at desctruction time + // (e.g. close files, deallocate resources etc.) + +} + + +// +// member functions +// + +// ------------ method called for each event ------------ +void +globalInfo::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + nEvents_++; + + if(useLHEWeights_){ + iEvent.getByToken(t_LHEInfo, h_LHEWeights); + double lheWeight = h_LHEWeights->weights()[0].wgt/std::abs(h_LHEWeights->weights()[0].wgt); + if(lheWeight < 0.){ + nNegLHEEvents_++; + } + } +} + + +// ------------ method called once each job just before starting event loop ------------ +void +globalInfo::beginJob() +{ +} + +// ------------ method called once each job just after ending the event loop ------------ +void +globalInfo::endJob() +{ + nEffEvents_ = nEvents_ - 2 * nNegLHEEvents_; + std::cout<<"total number of initial events is: "<make("tree" ,"tree" )); + + tree_->Branch("nEffEvents",&nEffEvents_,"nEffEvents/I"); + + tree_->Fill(); + + +} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void +globalInfo::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + //The following says we do not know what parameters are allowed so do no validation + // Please change this to state exactly what you do use, even if it is no parameters + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(globalInfo); \ No newline at end of file diff --git a/DeepNtuplizer/plugins/lheAnalyzer.cc b/DeepNtuplizer/plugins/lheAnalyzer.cc new file mode 100644 index 00000000000..83c17970604 --- /dev/null +++ b/DeepNtuplizer/plugins/lheAnalyzer.cc @@ -0,0 +1,141 @@ +// -*- C++ -*- +// +// Package: DeepNTuples/DeepNtuplizer +// Class: puAnalyzer +// +/**\class lheAnalyzer lheAnalyzer.cc DeepNTuples/DeepNtuplizer/plugins/lheAnalyzer.cc + Description: This analyzer is used in LHEcounter.py to produce a histogram with the lhe weights of the input files. + This is needed to calculate the effective event number of mc events which is the total number of events minus the once with negative lhe weights + Implementation: + [Notes on implementation] +*/ +// +// Original Author: David Walter +// Created: Wed, 22 Nov 2017 09:40:19 GMT +// +// + + +// system include files +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" + +#include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" +// +// class declaration +// + +// If the analyzer does not use TFileService, please remove +// the template argument to the base class so the class inherits +// from edm::one::EDAnalyzer<> and also remove the line from +// constructor "usesResource("TFileService");" +// This will improve performance in multithreaded jobs. + +class lheAnalyzer : public edm::one::EDAnalyzer { + public: + lheAnalyzer(const edm::ParameterSet&); + + ~lheAnalyzer(); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + + private: + virtual void beginJob() override; + virtual void analyze(const edm::Event&, const edm::EventSetup&) override; + virtual void endJob() override; + + // ----------member data --------------------------- + edm::EDGetTokenT t_LHEInfo; + TH1F * hist_lheWeight_0; +}; + +// +// constants, enums and typedefs +// + +// +// static data member definitions +// + +// +// constructors and destructor +// +lheAnalyzer::lheAnalyzer(const edm::ParameterSet& iConfig): + t_LHEInfo(consumes (iConfig.getParameter("lheInfo"))){ + + edm::Service fs; + hist_lheWeight_0 = fs->make( "lheweight" , "lhe weight [0]", 10, -2., 2. ); + +} + + +lheAnalyzer::~lheAnalyzer() +{ + + // do anything here that needs to be done at desctruction time + // (e.g. close files, deallocate resources etc.) + +} + + +// +// member functions +// + +// ------------ method called for each event ------------ +void +lheAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + + edm::Handle h_LHEWeights; + + iEvent.getByToken(t_LHEInfo, h_LHEWeights); + + double theWeight; + + theWeight = h_LHEWeights->weights()[0].wgt/std::abs(h_LHEWeights->weights()[0].wgt); + + + + hist_lheWeight_0->Fill(theWeight); + +} + + +// ------------ method called once each job just before starting event loop ------------ +void +lheAnalyzer::beginJob() +{ +} + +// ------------ method called once each job just after ending the event loop ------------ +void +lheAnalyzer::endJob() +{ +} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void +lheAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + //The following says we do not know what parameters are allowed so do no validation + // Please change this to state exactly what you do use, even if it is no parameters + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(lheAnalyzer); \ No newline at end of file diff --git a/DeepNtuplizer/plugins/puAnalyzer.cc b/DeepNtuplizer/plugins/puAnalyzer.cc new file mode 100644 index 00000000000..19969030fe0 --- /dev/null +++ b/DeepNtuplizer/plugins/puAnalyzer.cc @@ -0,0 +1,146 @@ +// -*- C++ -*- +// +// Package: DeepNTuples/DeepNtuplizer +// Class: puAnalyzer +// +/**\class puAnalyzer puAnalyzer.cc DeepNTuples/DeepNtuplizer/plugins/puAnalyzer.cc + Description: This analyzer is used in pupMCHistProd.py to calculate the pileup from the number of interactions per bunch crossing + Implementation: + [Notes on implementation] +*/ +// +// Original Author: David Walter +// Created: Mon, 20 Nov 2017 18:49:19 GMT +// +// + + +// system include files +#include +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h" + +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" +// +// class declaration +// + +// If the analyzer does not use TFileService, please remove +// the template argument to the base class so the class inherits +// from edm::one::EDAnalyzer<> and also remove the line from +// constructor "usesResource("TFileService");" +// This will improve performance in multithreaded jobs. + +class puAnalyzer : public edm::one::EDAnalyzer { + public: + puAnalyzer(const edm::ParameterSet&); + + ~puAnalyzer(); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + + private: + virtual void beginJob() override; + virtual void analyze(const edm::Event&, const edm::EventSetup&) override; + virtual void endJob() override; + + // ----------member data --------------------------- + edm::EDGetTokenT< std::vector > EDMPUInfoToken; + TH1F * h_pileup; +}; + +// +// constants, enums and typedefs +// + +// +// static data member definitions +// + +// +// constructors and destructor +// +puAnalyzer::puAnalyzer(const edm::ParameterSet& iConfig): + EDMPUInfoToken(consumes > (iConfig.getParameter("pileupInfo"))){ + + edm::Service fs; + h_pileup = fs->make( "pileup" , "Number of interactions per bunch crossing", 70, 0., 70 ); + +} + + +puAnalyzer::~puAnalyzer() +{ + + // do anything here that needs to be done at desctruction time + // (e.g. close files, deallocate resources etc.) + +} + + +// +// member functions +// + +// ------------ method called for each event ------------ +void +puAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + //using namespace edm; + + edm::Handle< std::vector > PupInfo; + iEvent.getByToken(EDMPUInfoToken, PupInfo); + + + std::vector::const_iterator PVI; + float Tnpv = -1; + for(PVI = PupInfo->begin(); PVI != PupInfo->end(); ++PVI) { + + int BX = PVI->getBunchCrossing(); + + if(BX == 0) { + Tnpv = PVI->getTrueNumInteractions(); + continue; + } + } + h_pileup->Fill(Tnpv); + +} + + +// ------------ method called once each job just before starting event loop ------------ +void +puAnalyzer::beginJob() +{ +} + +// ------------ method called once each job just after ending the event loop ------------ +void +puAnalyzer::endJob() +{ +} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void +puAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + //The following says we do not know what parameters are allowed so do no validation + // Please change this to state exactly what you do use, even if it is no parameters + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(puAnalyzer); \ No newline at end of file diff --git a/DeepNtuplizer/plugins/ttsemilepFilter.cc b/DeepNtuplizer/plugins/ttsemilepFilter.cc new file mode 100644 index 00000000000..12e3780b8c1 --- /dev/null +++ b/DeepNtuplizer/plugins/ttsemilepFilter.cc @@ -0,0 +1,271 @@ +// -*- C++ -*- +// +// Package: /DeepNTuplizer +// Class: DeepNTuplizer +// +/**\class ttsemilepFilter ttsemilepFilter.cc DeepNTuplizer/plugins/ttsemilepFilter.cc + + Description: This filter only takes the events which fulfill the following criteria: + - exactly one muon in src_muons + - no electron in src_electrons + - exactly four jets in src_jets + - transverse mass of muon and met > minMT_muonMETpair + + Implementation: + [Notes on implementation] +*/ +// +// Original Author: David Walter +// Created: Thu, 12 Apr 2018 17:00:08 GMT +// +// + + +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDFilter.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" + +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/PatCandidates/interface/Electron.h" +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "DataFormats/PatCandidates/interface/MET.h" + + +// +// class declaration +// + +class ttsemilepFilter : public edm::stream::EDFilter<> { + public: + explicit ttsemilepFilter(const edm::ParameterSet&); + ~ttsemilepFilter(); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + private: + virtual void beginStream(edm::StreamID) override; + virtual bool filter(edm::Event&, const edm::EventSetup&) override; + virtual void endStream() override; + + //virtual void beginRun(edm::Run const&, edm::EventSetup const&) override; + //virtual void endRun(edm::Run const&, edm::EventSetup const&) override; + //virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override; + //virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override; + + // ----------member data --------------------------- + + edm::EDGetTokenT > src_muons; + edm::EDGetTokenT > src_electrons; + edm::EDGetTokenT > src_jets; + edm::EDGetTokenT > src_mets; + + + edm::Handle > muons; + edm::Handle > electrons; + edm::Handle > jets; + edm::Handle > mets; + + double minMT_muonMETpair = 0; + unsigned int nElectrons = 0; + unsigned int nMuons = 0; + + /* + double min_chi2 = 0; + double m_w = 80.4; + double m_t = 173.1; + + int first_q = 0; + int second_q = 0; + int first_b = 0; + int second_b = 0; + + const std::string btag = "pfCSV_v2JetTags:probb"; + */ + +}; + +// +// constants, enums and typedefs +// + +// +// static data member definitions +// + +// +// constructors and destructor +// +ttsemilepFilter::ttsemilepFilter(const edm::ParameterSet& iConfig): + src_muons(consumes > (iConfig.getParameter("src_muons"))), + src_electrons(consumes > (iConfig.getParameter("src_electrons"))), + src_jets(consumes > (iConfig.getParameter("src_jets"))), + src_mets(consumes > (iConfig.getParameter("src_mets"))){ + //now do what ever initialization is needed + + //minPt_muon=(iConfig.getParameter("minPt_muon")); + //minPt_jet=(iConfig.getParameter("minPt_jet")); + minMT_muonMETpair=(iConfig.getParameter("cut_minMT_leptonMETpair")); + nElectrons=(iConfig.getParameter("nElectrons")); + nMuons=(iConfig.getParameter("nMuons")); + +} + + +ttsemilepFilter::~ttsemilepFilter() +{ + + // do anything here that needs to be done at destruction time + // (e.g. close files, deallocate resources etc.) + +} + + +// +// member functions +// + +// ------------ method called on each new Event ------------ +bool +ttsemilepFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + + iEvent.getByToken(src_muons, muons); + iEvent.getByToken(src_electrons, electrons); + iEvent.getByToken(src_jets, jets); + iEvent.getByToken(src_mets, mets); + + + if(muons->size()!=nMuons) + return false; + + if(electrons->size()!=nElectrons) + return false; + + if(jets->size()!=4) + return false; + + if(mets->size()==0){ //this should normally not happen, is the file corrupted?? + std::cout<<"there was an event with no met object, strange ..."<ptrAt(0); + double mWt = -1.; + + if(nMuons==1 && nElectrons==0){ + const auto muon = muons->ptrAt(0); + mWt = std::sqrt(2*(muon->pt()*met->pt() - muon->px()*met->px() - muon->py()*met->py())); + } + else if (nElectrons==1 && nMuons==0){ + const auto electron = electrons->ptrAt(0); + mWt = std::sqrt(2*(electron->pt()*met->pt() - electron->px()*met->px() - electron->py()*met->py())); + } + + + + + if(mWt == -1.){ //no W reconstruction + } + else if(mWt < minMT_muonMETpair){ //leptonic decaying W reconstruction + return false; + } + + + + //hadronically decaying t reconstruction + /* + for(size_t i = 0; i < jets->size(); ++i){ + for(size_t j = i+1; j < jets->size(); ++j){ + if(i!=j){ + double m_2jet = (jets->ptrAt(i)->p4() + jets->ptrAt(j)->p4()).M(); + + for(size_t k = 0; k < jets->size(); ++k){ + if(i!= j && i!=k && j!=k){ + double m_3jet = (jets->ptrAt(i)->p4() + jets->ptrAt(j)->p4() + jets->ptrAt(k)->p4()).M(); + + double chi2 = - std::log(std::pow(m_w - m_2jet,2)/m_2jet + std::pow(m_t - m_3jet,2)/m_3jet); + std::cout<<"combination ("< 1: +print ("running over these files:") +print (process.source.fileNames) + +process.globalInfo = cms.EDAnalyzer('globalInfo', + lheInfo = cms.InputTag("externalLHEProducer"), + useLHEWeights=cms.bool(options.lheWeights) + ) + +process.source.skipEvents = cms.untracked.uint32(options.skipEvents) +process.maxEvents = cms.untracked.PSet( + input=cms.untracked.int32(options.maxEvents) +) + +if int(release.replace("_", "")) >= 8400 or int(release.replace("_", "")) == 8029: + bTagInfos = [ + 'pfImpactParameterTagInfos', + 'pfInclusiveSecondaryVertexFinderTagInfos', + 'pfDeepCSVTagInfos'] +else: + bTagInfos = [ + 'pfImpactParameterTagInfos', + 'pfInclusiveSecondaryVertexFinderTagInfos', + 'deepNNTagInfos', + ] + +if int(release.replace("_", "")) >= 8400 or int(release.replace("_", "")) == 8029: + bTagDiscriminators = [ + 'softPFMuonBJetTags', + 'softPFElectronBJetTags', + 'pfJetBProbabilityBJetTags', + 'pfJetProbabilityBJetTags', + 'pfCombinedInclusiveSecondaryVertexV2BJetTags', + 'pfDeepCSVJetTags:probudsg', # to be fixed with new names + 'pfDeepCSVJetTags:probb', + 'pfDeepCSVJetTags:probc', + 'pfDeepCSVJetTags:probbb', + 'pfDeepCSVJetTags:probcc', + ] +else: + bTagDiscriminators = [ + 'softPFMuonBJetTags', + 'softPFElectronBJetTags', + 'pfJetBProbabilityBJetTags', + 'pfJetProbabilityBJetTags', + 'pfCombinedInclusiveSecondaryVertexV2BJetTags', + 'deepFlavourJetTags:probudsg', # to be fixed with new names + 'deepFlavourJetTags:probb', + 'deepFlavourJetTags:probc', + 'deepFlavourJetTags:probbb', + 'deepFlavourJetTags:probcc', + ] + +###### ttbar selection +outFileName = options.outputFile + '_' + str(options.job) + '.root' +print ('Using output file ' + outFileName) + +if options.deepNtuplizer: + process.TFileService = cms.Service("TFileService", + fileName=cms.string(outFileName)) +else: + process.MINIAODSIMEventContent.outputCommands.extend([ + 'keep *_goodElectrons_*_*', + 'keep *_goodMuons_*_*', + 'keep *_goodJets_*_*', + 'keep *_GoodOFLeptonPair_*_*' + ]) + + process.outmod = cms.OutputModule("PoolOutputModule", + process.MINIAODSIMEventContent, + SelectEvents = cms.untracked.PSet( + SelectEvents = cms.vstring('p') + ), + fileName = cms.untracked.string(outFileName) + ) + +from PhysicsTools.SelectorUtils.tools.vid_id_tools import * +dataFormat = DataFormat.MiniAOD +switchOnVIDElectronIdProducer(process, dataFormat) +my_id_modules = [ + 'RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_Summer16_80X_V1_cff' +] +for idmod in my_id_modules: + setupAllVIDIdsInModule(process,idmod,setupVIDElectronSelection) + + +#HighLevelTrigger +if options.dataRun in ['B','C','D','E','F','G']: + HLTlistSM = cms.vstring("HLT_Mu23_TrkIsoVVL_Ele12_CaloIdL_TrackIdL_IsoVL_v*", #Run B-G + "HLT_Mu8_TrkIsoVVL_Ele23_CaloIdL_TrackIdL_IsoVL_v*", #Run B-G + ) + +if options.dataRun in ['H']: + HLTlistSM = cms.vstring("HLT_Mu23_TrkIsoVVL_Ele12_CaloIdL_TrackIdL_IsoVL_DZ_v*", #Run H + "HLT_Mu8_TrkIsoVVL_Ele23_CaloIdL_TrackIdL_IsoVL_DZ_v*", #Run H + ) +if not options.isData: + HLTlistSM = cms.vstring("HLT_Mu23_TrkIsoVVL_Ele12_CaloIdL_TrackIdL_IsoVL_v*", #Run B-G + "HLT_Mu8_TrkIsoVVL_Ele23_CaloIdL_TrackIdL_IsoVL_v*", #Run B-G + "HLT_Mu23_TrkIsoVVL_Ele12_CaloIdL_TrackIdL_IsoVL_DZ_v*", #Run H + "HLT_Mu8_TrkIsoVVL_Ele23_CaloIdL_TrackIdL_IsoVL_DZ_v*", #Run H + ) + +process.TriggerSel = cms.EDFilter("HLTHighLevel", + TriggerResultsTag = cms.InputTag("TriggerResults","","HLT"), + HLTPaths = cms.vstring(HLTlistSM), + eventSetupPathsKey = cms.string(''), # not empty => use read paths from AlCaRecoTriggerBitsRcd via this key + andOr = cms.bool(True), # how to deal with multiple triggers: True (OR) accept if ANY is true, False (AND) accept if ALL are true + throw = cms.bool(True) # throw exception on unknown path names +) +# Electron Selection +process.goodElectrons = cms.EDProducer("ElectronIdAdder", + src=cms.InputTag("slimmedElectrons"), + vSrc=cms.InputTag("offlineSlimmedPrimaryVertices"), + idMap=cms.InputTag("egmGsfElectronIDs:cutBasedElectronID-Summer16-80X-V1-tight"), + minPt=cms.double(20.0), + maxAbsEta=cms.double(2.4), + ) + + +### Muon Selection +process.goodMuons = cms.EDProducer("MuonIdAdder", + src=cms.InputTag("slimmedMuons"), + vSrc=cms.InputTag("offlineSlimmedPrimaryVertices"), + minPt=cms.double(20.0), + maxAbsEta=cms.double(2.4), + maxRMI=cms.double(0.15) #RMI = relative muon isolation + ) + + +### Jets + +process.goodJets = cms.EDProducer("PATJetCleaner", + src = cms.InputTag("slimmedJets"), + preselection = cms.string("pt>30 && abs(eta) < 2.4 && neutralHadronEnergyFraction < 0.99 " + "&& neutralEmEnergyFraction < 0.99 && (chargedMultiplicity+neutralMultiplicity) > 1 " + "&& chargedHadronEnergyFraction > 0.0 " + "&& chargedMultiplicity > 0.0 && chargedEmEnergyFraction < 0.99 "), + checkOverlaps = cms.PSet( + muons = cms.PSet( + src = cms.InputTag("goodMuons"), + algorithm = cms.string("byDeltaR"), + preselection = cms.string(""), + deltaR = cms.double(0.4), + checkRecoComponents = cms.bool(False), # don't check if they share some AOD object ref + pairCut = cms.string(""), + requireNoOverlaps = cms.bool(True), # overlaps cause the jet to be discarded + ), + electrons = cms.PSet( + src = cms.InputTag("goodElectrons"), + algorithm = cms.string("byDeltaR"), + preselection = cms.string(""), + deltaR = cms.double(0.4), + checkRecoComponents = cms.bool(False), # don't check if they share some AOD object ref + pairCut = cms.string(""), + requireNoOverlaps = cms.bool(True), # overlaps cause the jet to be discarded + ) + ), + finalCut = cms.string('') + ) + +process.GoodOFLeptonPair = cms.EDProducer("CandViewShallowCloneCombiner", + decay = cms.string("goodMuons@+ goodElectrons@-"), + cut = cms.string("20.0 < mass" + "&& (daughter(0).pt > 25 || daughter(1).pt > 25)") + ) + + +process.FinalSel = cms.EDFilter("CandViewCountFilter", + src = cms.InputTag("GoodOFLeptonPair"), + minNumber = cms.uint32(1), + ) +### end selection + +# Jet Energy Corrections +if options.isData: + jetCorrections = cms.vstring(['L1FastJet', 'L2Relative', 'L3Absolute', 'L2L3Residual']) +else: + jetCorrections = cms.vstring(['L1FastJet', 'L2Relative', 'L3Absolute']) + +from PhysicsTools.PatAlgos.tools.jetTools import updateJetCollection +updateJetCollection( + process, + labelName="DeepFlavour", + # jetSource=cms.InputTag('slimmedJetsAK8PFPuppiSoftDropPacked', 'SubJets'), # 'subjets from AK8' + jetSource=cms.InputTag('goodJets'), # 'ak4Jets' + jetCorrections=('AK4PFchs', jetCorrections, 'None'), + pfCandidates=cms.InputTag('packedPFCandidates'), + pvSource=cms.InputTag("offlineSlimmedPrimaryVertices"), + svSource=cms.InputTag('slimmedSecondaryVertices'), + muSource=cms.InputTag('slimmedMuons'), + elSource=cms.InputTag('slimmedElectrons'), + btagInfos=bTagInfos, + btagDiscriminators=bTagDiscriminators, + explicitJTA=False + ) + +if hasattr(process, 'updatedPatJetsTransientCorrectedDeepFlavour'): + process.updatedPatJetsTransientCorrectedDeepFlavour.addTagInfos = cms.bool(True) + process.updatedPatJetsTransientCorrectedDeepFlavour.addBTagInfo = cms.bool(True) +else: + raise ValueError( + 'I could not find updatedPatJetsTransientCorrectedDeepFlavour to embed the tagInfos, please check the cfg') + +# QGLikelihood +process.load("DeepNTuples.DeepNtuplizer.QGLikelihood_cfi") +process.es_prefer_jec = cms.ESPrefer("PoolDBESSource", "QGPoolDBESSource") +process.load('RecoJets.JetProducers.QGTagger_cfi') +process.QGTagger.srcJets = cms.InputTag("selectedUpdatedPatJetsDeepFlavour") +process.QGTagger.jetsLabel = cms.string('QGL_AK4PFchs') + +from RecoJets.JetProducers.ak4GenJets_cfi import ak4GenJets + +process.ak4GenJetsWithNu = ak4GenJets.clone(src='packedGenParticles') + +## Filter out neutrinos from packed GenParticles +process.packedGenParticlesForJetsNoNu = cms.EDFilter("CandPtrSelector", src=cms.InputTag("packedGenParticles"), + cut=cms.string( + "abs(pdgId) != 12 && abs(pdgId) != 14 && abs(pdgId) != 16")) +## Define GenJets +process.ak4GenJetsRecluster = ak4GenJets.clone(src='packedGenParticlesForJetsNoNu') + +process.patGenJetMatchWithNu = cms.EDProducer("GenJetMatcher", # cut on deltaR; pick best by deltaR + src=cms.InputTag("selectedUpdatedPatJetsDeepFlavour"), + # RECO jets (any View is ok) + matched=cms.InputTag("ak4GenJetsWithNu"), + # GEN jets (must be GenJetCollection) + mcPdgId=cms.vint32(), # n/a + mcStatus=cms.vint32(), # n/a + checkCharge=cms.bool(False), # n/a + maxDeltaR=cms.double(0.4), # Minimum deltaR for the match + # maxDPtRel = cms.double(3.0), # Minimum deltaPt/Pt for the match (not used in GenJetMatcher) + resolveAmbiguities=cms.bool(True), + # Forbid two RECO objects to match to the same GEN object + resolveByMatchQuality=cms.bool(False), + # False = just match input in order; True = pick lowest deltaR pair first + ) + +process.patGenJetMatchRecluster = cms.EDProducer("GenJetMatcher", # cut on deltaR; pick best by deltaR + src=cms.InputTag("selectedUpdatedPatJetsDeepFlavour"), + # RECO jets (any View is ok) + matched=cms.InputTag("ak4GenJetsRecluster"), + # GEN jets (must be GenJetCollection) + mcPdgId=cms.vint32(), # n/a + mcStatus=cms.vint32(), # n/a + checkCharge=cms.bool(False), # n/a + maxDeltaR=cms.double(0.4), # Minimum deltaR for the match + # maxDPtRel = cms.double(3.0), # Minimum deltaPt/Pt for the match (not used in GenJetMatcher) + resolveAmbiguities=cms.bool(True), + # Forbid two RECO objects to match to the same GEN object + resolveByMatchQuality=cms.bool(False), + # False = just match input in order; True = pick lowest deltaR pair first + ) + +process.genJetSequence = cms.Sequence(process.packedGenParticlesForJetsNoNu + * process.ak4GenJetsWithNu + * process.ak4GenJetsRecluster + * process.patGenJetMatchWithNu + * process.patGenJetMatchRecluster) + +# Very Loose IVF SV collection +from PhysicsTools.PatAlgos.tools.helpers import loadWithPrefix + +loadWithPrefix(process, 'RecoVertex.AdaptiveVertexFinder.inclusiveVertexing_cff', "looseIVF") +process.looseIVFinclusiveCandidateVertexFinder.primaryVertices = cms.InputTag("offlineSlimmedPrimaryVertices") +process.looseIVFinclusiveCandidateVertexFinder.tracks = cms.InputTag("packedPFCandidates") +process.looseIVFinclusiveCandidateVertexFinder.vertexMinDLen2DSig = cms.double(0.) +process.looseIVFinclusiveCandidateVertexFinder.vertexMinDLenSig = cms.double(0.) +process.looseIVFinclusiveCandidateVertexFinder.fitterSigmacut = 20 + +process.looseIVFcandidateVertexArbitrator.primaryVertices = cms.InputTag("offlineSlimmedPrimaryVertices") +process.looseIVFcandidateVertexArbitrator.tracks = cms.InputTag("packedPFCandidates") +process.looseIVFcandidateVertexArbitrator.secondaryVertices = cms.InputTag("looseIVFcandidateVertexMerger") +process.looseIVFcandidateVertexArbitrator.fitterSigmacut = 20 + +#datapath=os.environ['CMSSW_BASE']+'/src/DeepNTuples/DeepNtuplizer/data/' +datapath='' + +# DeepNtuplizer configurations +process.load("DeepNTuples.DeepNtuplizer.DeepNtuplizer_cfi") +process.deepntuplizer.jets = cms.InputTag('selectedUpdatedPatJetsDeepFlavour') +process.deepntuplizer.bDiscriminators = bTagDiscriminators +process.deepntuplizer.bDiscriminators.append('pfCombinedMVAV2BJetTags') +process.deepntuplizer.LooseSVs = cms.InputTag("looseIVFinclusiveCandidateSecondaryVertices") +process.deepntuplizer.applySelection = cms.bool(options.selectJets) + +process.deepntuplizer.removeUndefined = cms.bool(False) +process.deepntuplizer.isData = cms.bool(options.isData) +process.deepntuplizer.useLHEWeights = cms.bool(options.lheWeights) + +process.deepntuplizer.pileupData=cms.string(datapath+"pileup_data_2016.root") +process.deepntuplizer.pileupMC=cms.string(datapath+"pileup_MC_2016.root") + +process.deepntuplizer.sfMuons = cms.InputTag("goodMuons") +process.deepntuplizer.sfElectrons=cms.InputTag("goodElectrons") + +process.deepntuplizer.periods=cms.vstring("2016BtoF","2016G","2016H") +process.deepntuplizer.lumis=cms.vdouble(5.404+2.396+4.243+4.054+3.105,7.544,8.453) + +process.deepntuplizer.triggers=cms.vstring("HLT_Mu23_TrkIsoVVL_Ele12_CaloIdL_TrackIdL_IsoVL_v*||HLT_Mu8_TrkIsoVVL_Ele23_CaloIdL_TrackIdL_IsoVL_v*", #Run B-G + "HLT_Mu23_TrkIsoVVL_Ele12_CaloIdL_TrackIdL_IsoVL_v*||HLT_Mu8_TrkIsoVVL_Ele23_CaloIdL_TrackIdL_IsoVL_v*", #Run B-G + "HLT_Mu23_TrkIsoVVL_Ele12_CaloIdL_TrackIdL_IsoVL_DZ_v*||HLT_Mu8_TrkIsoVVL_Ele23_CaloIdL_TrackIdL_IsoVL_DZ_v*") #Run H + +#emu trigger scalefactors from +# https://gitlab.cern.ch/ttH/reference/blob/955cff0b09f2c95bc480ae0bf3145aab9ce08fcc/definitions/Moriond17.md#64-trigger-scale-factors +# https://gitlab.cern.ch/ttH/reference/blob/955cff0b09f2c95bc480ae0bf3145aab9ce08fcc/definitions/Moriond17.md#41-triggers +process.deepntuplizer.sfTrigger_emu=cms.vstring(datapath + "triggerSummary_emu_ReReco2016_ttH.root") +process.deepntuplizer.sfTrigger_emu_Hist = cms.vstring("scalefactor_eta2d_with_syst") + +process.deepntuplizer.sfMuonId = cms.vstring(datapath+"EfficienciesAndSF_ID_BCDEF.root", + datapath + "EfficienciesAndSF_ID_GH.root", + datapath + "EfficienciesAndSF_ID_GH.root") +process.deepntuplizer.sfMuonId_Hist = cms.vstring("MC_NUM_TightID_DEN_genTracks_PAR_pt_eta/abseta_pt_ratio", + "MC_NUM_TightID_DEN_genTracks_PAR_pt_eta/abseta_pt_ratio", + "MC_NUM_TightID_DEN_genTracks_PAR_pt_eta/abseta_pt_ratio") +process.deepntuplizer.sfMuonIso=cms.vstring(datapath + "EfficienciesAndSF_ISO_BCDEF.root", + datapath + "EfficienciesAndSF_ISO_GH.root", + datapath + "EfficienciesAndSF_ISO_GH.root") +process.deepntuplizer.sfMuonIso_Hist=cms.vstring("TightISO_TightID_pt_eta/abseta_pt_ratio", + "TightISO_TightID_pt_eta/abseta_pt_ratio", + "TightISO_TightID_pt_eta/abseta_pt_ratio") +process.deepntuplizer.sfMuonTracking=cms.vstring(datapath+"Tracking_EfficienciesAndSF_BCDEFGH.root") +process.deepntuplizer.sfMuonTracking_Hist=cms.vstring("ratio_eff_aeta_dr030e030_corr") +process.deepntuplizer.sfElIdAndIso=cms.vstring(datapath+"egammaEffi.txt_EGM2D.root") +process.deepntuplizer.sfElIdAndIso_Hist=cms.vstring("EGamma_SF2D") + + +if int(release.replace("_", "")) >= 840: + process.deepntuplizer.tagInfoName = cms.string('pfDeepCSV') + +process.deepntuplizer.gluonReduction = cms.double(options.gluonReduction) + +# 1631 +process.ProfilerService = cms.Service( + "ProfilerService", + firstEvent=cms.untracked.int32(1631), + lastEvent=cms.untracked.int32(1641), + paths=cms.untracked.vstring('p') +) + + +if options.deepNtuplizer: + if options.isData: + process.p = cms.Path(process.globalInfo + process.TriggerSel + process.FinalSel + process.QGTagger + process.deepntuplizer) + else: + process.p = cms.Path(process.globalInfo + process.TriggerSel + process.FinalSel + process.QGTagger + process.genJetSequence * process.deepntuplizer) +else: + process.p = cms.Path(process.TriggerSel + process.FinalSel) + + process.endp = cms.EndPath(process.outmod) + diff --git a/DeepNtuplizer/production/DeepNtuplizer_tt_semilep_e.py b/DeepNtuplizer/production/DeepNtuplizer_tt_semilep_e.py new file mode 100644 index 00000000000..a19689f3bc5 --- /dev/null +++ b/DeepNtuplizer/production/DeepNtuplizer_tt_semilep_e.py @@ -0,0 +1,418 @@ +# this can be run with CMSSW 8_0_29; in CMSSW 8_0_25 the module 'cutBasedElectronID_Summer16_80X_V1_cff' is missing +#basically deepntuplizer with tt semileptonic single electron selection + +import FWCore.ParameterSet.Config as cms + +import FWCore.ParameterSet.VarParsing as VarParsing +### parsing job options +import sys + +options = VarParsing.VarParsing() + +options.register('inputScript', '', VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.string,"input Script") +options.register('outputFile', 'output', VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.string, "output File (w/o .root)") +options.register('maxEvents', -1, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int,"maximum events") +options.register('skipEvents', 0, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "skip N events") +options.register('job', 0, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int,"job number") +options.register('nJobs', 1, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "total jobs") +options.register('gluonReduction', 0.0, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.float, "gluon reduction") +options.register('selectJets', True, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool,"select jets with good gen match") +options.register('globalTag', '', VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.string,"global tag for jet energy correction") +options.register('isData', False, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "switch off generator jets") +options.register('deepNtuplizer',True, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "run deepNtuplizer or just the ttbar selection") +options.register('lheWeights',False, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "use LHE weights") +options.register('crossSection', 1., VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.float,"cross section") +options.register('nEvents', 1, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int,"number of events to process") + + +import os + +release = os.environ['CMSSW_VERSION'][6:12] +print("Using release " + release) + +options.register( + 'inputFiles', '', + VarParsing.VarParsing.multiplicity.list, + VarParsing.VarParsing.varType.string, + "input files (default is the tt RelVal)" +) + +if hasattr(sys, "argv"): + options.parseArguments() + +process = cms.Process("semilepElectronSelectedDNNFiller") + +process.load("FWCore.MessageService.MessageLogger_cfi") +process.load("Configuration.EventContent.EventContent_cff") +process.load('Configuration.StandardSequences.Services_cff') +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +if options.globalTag == '': + if options.isData == False: + process.GlobalTag.globaltag = '80X_mcRun2_asymptotic_2016_TrancheIV_v8' # For MC Jet Energy correction + if options.isData == True: + process.GlobalTag.globaltag = '80X_dataRun2_2016SeptRepro_v7' # For Data Jet Energy correction +else: + process.GlobalTag.globaltag = options.globalTag + +process.maxEvents = cms.untracked.PSet(input=cms.untracked.int32(-1)) +#process.PoolSource = cms.untracked.PSet( +# firstEvent = cms.untracked.uint32(33285817), +#) + +process.load('FWCore.MessageService.MessageLogger_cfi') +process.MessageLogger.cerr.FwkReport.reportEvery = 10 +if options.inputScript == '': # this is probably for testing + process.MessageLogger.cerr.FwkReport.reportEvery = 100 + +process.options = cms.untracked.PSet( + allowUnscheduled=cms.untracked.bool(True), + wantSummary=cms.untracked.bool(False), + SkipEvent = cms.untracked.vstring('ProductNotFound') +) + + +process.load('DeepNTuples.DeepNtuplizer.samples.TTJetsPhase1_cfg') # default input + +if options.inputFiles: + process.source.fileNames = options.inputFiles + +if options.inputScript != '' and options.inputScript != 'DeepNTuples.DeepNtuplizer.samples.TTJetsPhase1_cfg': + process.load(options.inputScript) + + +#process.source.fileNames=['file:./00CC509E-0C3B-E711-98F2-0242AC130004.root'] #isData=True +#process.source.fileNames=['file:./0693E0E7-97BE-E611-B32F-0CC47A78A3D8.root'] #isData=False +#process.source.fileNames=['file:./EE95DEDC-96BE-E611-B45D-A0000420FE80.root'] #isData=False +#process.source.fileNames=['file:./0693E0E7-97BE-E611-B32F-0CC47A78A3D8.root'] #store/mc/RunIISummer16MiniAODv2/TT_TuneCUETP8M2T4_13TeV-powheg-pythia8/MINIAODSIM/PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v6-v1/50000/ +#process.source.fileNames=['file:./00E02A09-853C-E711-93FF-3417EBE644A7.root'] #store/data/Run2016H/SingleMuon/MINIAOD/18Apr2017-v1/00000/ +#process.source.fileNames=['file:./F21AE451-7EBA-E611-9399-0025905B858E.root'] #store/mc/RunIISummer16MiniAODv2/TT_TuneCUETP8M2T4_13TeV-powheg-pythia8/MINIAODSIM/PUMoriond17_backup_80X_mcRun2_asymptotic_2016_TrancheIV_v6-v1/00000/ + +#process.source.fileNames=['file:./00CC509E-0C3B-E711-98F2-0242AC130004.root'] #Run2016B/SingleMuon/MINIAOD/18Apr2017_ver2-v1/ +numberOfFiles = len(process.source.fileNames) +numberOfJobs = options.nJobs +jobNumber = options.job + +process.source.fileNames = process.source.fileNames[jobNumber:numberOfFiles:numberOfJobs] +#if options.nJobs > 1: +print ("running over these files:") +print (process.source.fileNames) + +process.globalInfo = cms.EDAnalyzer('globalInfo', + lheInfo = cms.InputTag("externalLHEProducer"), + useLHEWeights=cms.bool(options.lheWeights) + ) + +process.source.skipEvents = cms.untracked.uint32(options.skipEvents) +process.maxEvents = cms.untracked.PSet( + input=cms.untracked.int32(options.maxEvents) +) + +if int(release.replace("_", "")) >= 8400 or int(release.replace("_", "")) == 8029: + bTagInfos = [ + 'pfImpactParameterTagInfos', + 'pfInclusiveSecondaryVertexFinderTagInfos', + 'pfDeepCSVTagInfos'] +else: + bTagInfos = [ + 'pfImpactParameterTagInfos', + 'pfInclusiveSecondaryVertexFinderTagInfos', + 'deepNNTagInfos', + ] + +if int(release.replace("_", "")) >= 8400 or int(release.replace("_", "")) == 8029: + bTagDiscriminators = [ + 'softPFMuonBJetTags', + 'softPFElectronBJetTags', + 'pfJetBProbabilityBJetTags', + 'pfJetProbabilityBJetTags', + 'pfCombinedInclusiveSecondaryVertexV2BJetTags', + 'pfDeepCSVJetTags:probudsg', # to be fixed with new names + 'pfDeepCSVJetTags:probb', + 'pfDeepCSVJetTags:probc', + 'pfDeepCSVJetTags:probbb', + 'pfDeepCSVJetTags:probcc', + ] +else: + bTagDiscriminators = [ + 'softPFMuonBJetTags', + 'softPFElectronBJetTags', + 'pfJetBProbabilityBJetTags', + 'pfJetProbabilityBJetTags', + 'pfCombinedInclusiveSecondaryVertexV2BJetTags', + 'deepFlavourJetTags:probudsg', # to be fixed with new names + 'deepFlavourJetTags:probb', + 'deepFlavourJetTags:probc', + 'deepFlavourJetTags:probbb', + 'deepFlavourJetTags:probcc', + ] + +###### semilep selection +outFileName = options.outputFile + '_' + str(options.job) + '.root' +print ('Using output file ' + outFileName) + +if options.deepNtuplizer: + process.TFileService = cms.Service("TFileService", + fileName=cms.string(outFileName)) +else: + process.MINIAODSIMEventContent.outputCommands.extend([ + 'keep *_goodElectrons_*_*', + 'keep *_goodMuons_*_*', + 'keep *_goodJets_*_*', + ]) + + process.outmod = cms.OutputModule("PoolOutputModule", + process.MINIAODSIMEventContent, + SelectEvents = cms.untracked.PSet( + SelectEvents = cms.vstring('p') + ), + fileName = cms.untracked.string(outFileName) + ) + +from PhysicsTools.SelectorUtils.tools.vid_id_tools import * +dataFormat = DataFormat.MiniAOD +switchOnVIDElectronIdProducer(process, dataFormat) +my_id_modules = [ + 'RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_Summer16_80X_V1_cff' +] +for idmod in my_id_modules: + setupAllVIDIdsInModule(process,idmod,setupVIDElectronSelection) + + +#HighLevelTrigger +HLTlistSM = cms.vstring("HLT_Ele27_WPTight_Gsf_v*" + ) +process.TriggerSel = cms.EDFilter("HLTHighLevel", + TriggerResultsTag = cms.InputTag("TriggerResults","","HLT"), + HLTPaths = cms.vstring(HLTlistSM), + eventSetupPathsKey = cms.string(''), # not empty => use read paths from AlCaRecoTriggerBitsRcd via this key + andOr = cms.bool(True), # how to deal with multiple triggers: True (OR) accept if ANY is true, False (AND) accept if ALL are true + throw = cms.bool(True) # throw exception on unknown path names +) +### Electron Selection +process.goodElectrons = cms.EDProducer("ElectronIdAdder", + src=cms.InputTag("slimmedElectrons"), + vSrc=cms.InputTag("offlineSlimmedPrimaryVertices"), + idMap=cms.InputTag("egmGsfElectronIDs:cutBasedElectronID-Summer16-80X-V1-tight"), + minPt=cms.double(20.0), + maxAbsEta=cms.double(2.4), + ) + +### Muon Selection +process.goodMuons = cms.EDProducer("MuonIdAdder", + src=cms.InputTag("slimmedMuons"), + vSrc=cms.InputTag("offlineSlimmedPrimaryVertices"), + minPt=cms.double(30.0), + maxAbsEta=cms.double(2.4), + maxRMI=cms.double(0.15) #RMI = relative muon isolation + ) + +### Jet Selection +process.goodJets = cms.EDProducer("PATJetCleaner", + src = cms.InputTag("slimmedJets"), + preselection = cms.string("pt>30 && abs(eta) < 2.4 && neutralHadronEnergyFraction < 0.99 " + "&& neutralEmEnergyFraction < 0.99 && (chargedMultiplicity+neutralMultiplicity) > 1 " + "&& chargedHadronEnergyFraction > 0.0 " + "&& chargedMultiplicity > 0.0 && chargedEmEnergyFraction < 0.99 "), + checkOverlaps = cms.PSet( + muons = cms.PSet( + src = cms.InputTag("goodMuons"), + algorithm = cms.string("byDeltaR"), + preselection = cms.string(""), + deltaR = cms.double(0.4), + checkRecoComponents = cms.bool(False), # don't check if they share some AOD object ref + pairCut = cms.string(""), + requireNoOverlaps = cms.bool(True), # overlaps cause the jet to be discarded + ), + electrons = cms.PSet( + src = cms.InputTag("goodElectrons"), + algorithm = cms.string("byDeltaR"), + preselection = cms.string(""), + deltaR = cms.double(0.4), + checkRecoComponents = cms.bool(False), # don't check if they share some AOD object ref + pairCut = cms.string(""), + requireNoOverlaps = cms.bool(True), # overlaps cause the jet to be discarded + ) + ), + finalCut = cms.string('') + ) + + +### Event Filter +process.ttsemilepFilter = cms.EDFilter("ttsemilepFilter", + src_muons = cms.InputTag("goodMuons"), + src_electrons = cms.InputTag("goodElectrons"), + src_jets = cms.InputTag("goodJets"), + src_mets = cms.InputTag("slimmedMETs"), + cut_minMT_leptonMETpair = cms.double(50.0), + nElectrons=cms.uint32(1), + nMuons=cms.uint32(0) + ) + + + +### end selection + +# Jet Energy Corrections +if options.isData: + jetCorrections = cms.vstring(['L1FastJet', 'L2Relative', 'L3Absolute', 'L2L3Residual']) +else: + jetCorrections = cms.vstring(['L1FastJet', 'L2Relative', 'L3Absolute']) + +from PhysicsTools.PatAlgos.tools.jetTools import updateJetCollection +updateJetCollection( + process, + labelName="DeepFlavour", + # jetSource=cms.InputTag('slimmedJetsAK8PFPuppiSoftDropPacked', 'SubJets'), # 'subjets from AK8' + jetSource=cms.InputTag('goodJets'), # 'ak4Jets' + jetCorrections=('AK4PFchs', jetCorrections, 'None'), + pfCandidates=cms.InputTag('packedPFCandidates'), + pvSource=cms.InputTag("offlineSlimmedPrimaryVertices"), + svSource=cms.InputTag('slimmedSecondaryVertices'), + muSource=cms.InputTag('slimmedMuons'), + elSource=cms.InputTag('slimmedElectrons'), + btagInfos=bTagInfos, + btagDiscriminators=bTagDiscriminators, + explicitJTA=False + ) + +if hasattr(process, 'updatedPatJetsTransientCorrectedDeepFlavour'): + process.updatedPatJetsTransientCorrectedDeepFlavour.addTagInfos = cms.bool(True) + process.updatedPatJetsTransientCorrectedDeepFlavour.addBTagInfo = cms.bool(True) +else: + raise ValueError( + 'I could not find updatedPatJetsTransientCorrectedDeepFlavour to embed the tagInfos, please check the cfg') + +# QGLikelihood +process.load("DeepNTuples.DeepNtuplizer.QGLikelihood_cfi") +process.es_prefer_jec = cms.ESPrefer("PoolDBESSource", "QGPoolDBESSource") +process.load('RecoJets.JetProducers.QGTagger_cfi') +process.QGTagger.srcJets = cms.InputTag("selectedUpdatedPatJetsDeepFlavour") +process.QGTagger.jetsLabel = cms.string('QGL_AK4PFchs') + +from RecoJets.JetProducers.ak4GenJets_cfi import ak4GenJets + +process.ak4GenJetsWithNu = ak4GenJets.clone(src='packedGenParticles') + +## Filter out neutrinos from packed GenParticles +process.packedGenParticlesForJetsNoNu = cms.EDFilter("CandPtrSelector", src=cms.InputTag("packedGenParticles"), + cut=cms.string( + "abs(pdgId) != 12 && abs(pdgId) != 14 && abs(pdgId) != 16")) +## Define GenJets +process.ak4GenJetsRecluster = ak4GenJets.clone(src='packedGenParticlesForJetsNoNu') + +process.patGenJetMatchWithNu = cms.EDProducer("GenJetMatcher", # cut on deltaR; pick best by deltaR + src=cms.InputTag("selectedUpdatedPatJetsDeepFlavour"), + # RECO jets (any View is ok) + matched=cms.InputTag("ak4GenJetsWithNu"), + # GEN jets (must be GenJetCollection) + mcPdgId=cms.vint32(), # n/a + mcStatus=cms.vint32(), # n/a + checkCharge=cms.bool(False), # n/a + maxDeltaR=cms.double(0.4), # Minimum deltaR for the match + # maxDPtRel = cms.double(3.0), # Minimum deltaPt/Pt for the match (not used in GenJetMatcher) + resolveAmbiguities=cms.bool(True), + # Forbid two RECO objects to match to the same GEN object + resolveByMatchQuality=cms.bool(False), + # False = just match input in order; True = pick lowest deltaR pair first + ) + +process.patGenJetMatchRecluster = cms.EDProducer("GenJetMatcher", # cut on deltaR; pick best by deltaR + src=cms.InputTag("selectedUpdatedPatJetsDeepFlavour"), + # RECO jets (any View is ok) + matched=cms.InputTag("ak4GenJetsRecluster"), + # GEN jets (must be GenJetCollection) + mcPdgId=cms.vint32(), # n/a + mcStatus=cms.vint32(), # n/a + checkCharge=cms.bool(False), # n/a + maxDeltaR=cms.double(0.4), # Minimum deltaR for the match + # maxDPtRel = cms.double(3.0), # Minimum deltaPt/Pt for the match (not used in GenJetMatcher) + resolveAmbiguities=cms.bool(True), + # Forbid two RECO objects to match to the same GEN object + resolveByMatchQuality=cms.bool(False), + # False = just match input in order; True = pick lowest deltaR pair first + ) + +process.genJetSequence = cms.Sequence(process.packedGenParticlesForJetsNoNu + * process.ak4GenJetsWithNu + * process.ak4GenJetsRecluster + * process.patGenJetMatchWithNu + * process.patGenJetMatchRecluster) + +# Very Loose IVF SV collection +from PhysicsTools.PatAlgos.tools.helpers import loadWithPrefix + +loadWithPrefix(process, 'RecoVertex.AdaptiveVertexFinder.inclusiveVertexing_cff', "looseIVF") +process.looseIVFinclusiveCandidateVertexFinder.primaryVertices = cms.InputTag("offlineSlimmedPrimaryVertices") +process.looseIVFinclusiveCandidateVertexFinder.tracks = cms.InputTag("packedPFCandidates") +process.looseIVFinclusiveCandidateVertexFinder.vertexMinDLen2DSig = cms.double(0.) +process.looseIVFinclusiveCandidateVertexFinder.vertexMinDLenSig = cms.double(0.) +process.looseIVFinclusiveCandidateVertexFinder.fitterSigmacut = 20 + +process.looseIVFcandidateVertexArbitrator.primaryVertices = cms.InputTag("offlineSlimmedPrimaryVertices") +process.looseIVFcandidateVertexArbitrator.tracks = cms.InputTag("packedPFCandidates") +process.looseIVFcandidateVertexArbitrator.secondaryVertices = cms.InputTag("looseIVFcandidateVertexMerger") +process.looseIVFcandidateVertexArbitrator.fitterSigmacut = 20 + +#datapath=os.environ['CMSSW_BASE']+'/src/DeepNTuples/DeepNtuplizer/data/' +datapath='' + +# DeepNtuplizer configurations +process.load("DeepNTuples.DeepNtuplizer.DeepNtuplizer_cfi") +process.deepntuplizer.jets = cms.InputTag('selectedUpdatedPatJetsDeepFlavour') +process.deepntuplizer.bDiscriminators = bTagDiscriminators +process.deepntuplizer.bDiscriminators.append('pfCombinedMVAV2BJetTags') +process.deepntuplizer.LooseSVs = cms.InputTag("looseIVFinclusiveCandidateSecondaryVertices") +process.deepntuplizer.applySelection = cms.bool(options.selectJets) + +process.deepntuplizer.isData = cms.bool(options.isData) + +process.deepntuplizer.removeUndefined = cms.bool(False) + +if not options.isData: + + process.deepntuplizer.useLHEWeights = cms.bool(options.lheWeights) + + process.deepntuplizer.pileupData=cms.string(datapath+"pileup_data_2016.root") + process.deepntuplizer.pileupMC=cms.string(datapath+"pileup_MC_2016.root") + + process.deepntuplizer.sfElectrons = cms.InputTag("goodElectrons") + + process.deepntuplizer.periods=cms.vstring("2016All",) + process.deepntuplizer.lumis=cms.vdouble(1.,) + process.deepntuplizer.crossSection=cms.double(options.crossSection) + process.deepntuplizer.nEvents=cms.uint32(options.nEvents) + + process.deepntuplizer.sfTrigger_e=cms.vstring(datapath + "TriggerSF_Run2016All_v1.root",) + process.deepntuplizer.sfTrigger_e_Hist=cms.vstring("Ele27_WPTight_Gsf",) + + process.deepntuplizer.sfElIdAndIso=cms.vstring(datapath+"egammaEffi.txt_EGM2D.root",) + process.deepntuplizer.sfElIdAndIso_Hist=cms.vstring("EGamma_SF2D",) + + process.deepntuplizer.gluonReduction = cms.double(options.gluonReduction) + +if int(release.replace("_", "")) >= 840: + process.deepntuplizer.tagInfoName = cms.string('pfDeepCSV') + + +# 1631 +#process.ProfilerService = cms.Service( +# "ProfilerService", +# firstEvent=cms.untracked.int32(1631), +# lastEvent=cms.untracked.int32(1641), +# paths=cms.untracked.vstring('p') +#) + + +if options.deepNtuplizer: + if options.isData: + process.p = cms.Path(process.globalInfo + process.TriggerSel + process.ttsemilepFilter + process.QGTagger + process.deepntuplizer) + else: + process.p = cms.Path(process.globalInfo + process.TriggerSel + process.ttsemilepFilter + process.QGTagger + process.genJetSequence * process.deepntuplizer) +else: + process.p = cms.Path(process.TriggerSel + process.ttsemilepFilter) + + process.endp = cms.EndPath(process.outmod) + diff --git a/DeepNtuplizer/production/DeepNtuplizer_tt_semilep_mu.py b/DeepNtuplizer/production/DeepNtuplizer_tt_semilep_mu.py new file mode 100644 index 00000000000..df624a611fb --- /dev/null +++ b/DeepNtuplizer/production/DeepNtuplizer_tt_semilep_mu.py @@ -0,0 +1,428 @@ +# this can be run with CMSSW 8_0_29; in CMSSW 8_0_25 the module 'cutBasedElectronID_Summer16_80X_V1_cff' is missing +#basically deepntuplizer with tt semileptonic selection + +import FWCore.ParameterSet.Config as cms + +import FWCore.ParameterSet.VarParsing as VarParsing +### parsing job options +import sys + +options = VarParsing.VarParsing() + +options.register('inputScript', '', VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.string,"input Script") +options.register('outputFile', 'output', VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.string, "output File (w/o .root)") +options.register('maxEvents', -1, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int,"maximum events") +options.register('skipEvents', 0, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "skip N events") +options.register('job', 0, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int,"job number") +options.register('nJobs', 1, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "total jobs") +options.register('gluonReduction', 0.0, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.float, "gluon reduction") +options.register('selectJets', True, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool,"select jets with good gen match") +options.register('globalTag', '', VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.string,"global tag for jet energy correction") +options.register('isData', False, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "switch off generator jets") +options.register('deepNtuplizer',True, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "run deepNtuplizer or just the ttbar selection") +options.register('lheWeights',False, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "use LHE weights") +options.register('crossSection', 1., VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.float,"cross section") +options.register('nEvents', 1, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int,"number of events to process") + +import os + +release = os.environ['CMSSW_VERSION'][6:12] +print("Using release " + release) + +options.register( + 'inputFiles', '', + VarParsing.VarParsing.multiplicity.list, + VarParsing.VarParsing.varType.string, + "input files (default is the tt RelVal)" +) + +if hasattr(sys, "argv"): + options.parseArguments() + +process = cms.Process("semilepMuonSelectedDNNFiller") + +process.load("FWCore.MessageService.MessageLogger_cfi") +process.load("Configuration.EventContent.EventContent_cff") +process.load('Configuration.StandardSequences.Services_cff') +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') + +if options.globalTag == '': + if options.isData == False: + process.GlobalTag.globaltag = '80X_mcRun2_asymptotic_2016_TrancheIV_v8' # For MC Jet Energy correction + if options.isData == True: + process.GlobalTag.globaltag = '80X_dataRun2_2016SeptRepro_v7' # For Data Jet Energy correction +else: + process.GlobalTag.globaltag = options.globalTag + +process.maxEvents = cms.untracked.PSet(input=cms.untracked.int32(-1)) +#process.PoolSource = cms.untracked.PSet( +# firstEvent = cms.untracked.uint32(33285817), +#) + +process.load('FWCore.MessageService.MessageLogger_cfi') +process.MessageLogger.cerr.FwkReport.reportEvery = 10 +if options.inputScript == '': # this is probably for testing + process.MessageLogger.cerr.FwkReport.reportEvery = 100 + +process.options = cms.untracked.PSet( + allowUnscheduled=cms.untracked.bool(True), + wantSummary=cms.untracked.bool(False), + SkipEvent = cms.untracked.vstring('ProductNotFound') +) + + +process.load('DeepNTuples.DeepNtuplizer.samples.TTJetsPhase1_cfg') # default input + +if options.inputFiles: + process.source.fileNames = options.inputFiles + +if options.inputScript != '' and options.inputScript != 'DeepNTuples.DeepNtuplizer.samples.TTJetsPhase1_cfg': + process.load(options.inputScript) + + +#process.source.fileNames=['file:./00E02A09-853C-E711-93FF-3417EBE644A7.root'] #store/data/Run2016H/SingleMuon/MINIAOD/18Apr2017-v1/00000/00E02A09-853C-E711-93FF-3417EBE644A7.root +#process.source.fileNames=['file:./02AAA49F-C83A-E711-B51F-0CC47AA53D82.root'] #isData=True +#process.source.fileNames=['file:./0693E0E7-97BE-E611-B32F-0CC47A78A3D8.root'] #isData=False +#process.source.fileNames=['file:./EE95DEDC-96BE-E611-B45D-A0000420FE80.root'] #isData=False +#process.source.fileNames=['file:./0693E0E7-97BE-E611-B32F-0CC47A78A3D8.root'] #store/mc/RunIISummer16MiniAODv2/TT_TuneCUETP8M2T4_13TeV-powheg-pythia8/MINIAODSIM/PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v6-v1/50000/ +#process.source.fileNames=['file:./00E02A09-853C-E711-93FF-3417EBE644A7.root'] #store/data/Run2016H/SingleMuon/MINIAOD/18Apr2017-v1/00000/ +#process.source.fileNames=['file:./F21AE451-7EBA-E611-9399-0025905B858E.root'] #store/mc/RunIISummer16MiniAODv2/TT_TuneCUETP8M2T4_13TeV-powheg-pythia8/MINIAODSIM/PUMoriond17_backup_80X_mcRun2_asymptotic_2016_TrancheIV_v6-v1/00000/ + +#process.source.fileNames=['file:./00CC509E-0C3B-E711-98F2-0242AC130004.root'] #Run2016B/SingleMuon/MINIAOD/18Apr2017_ver2-v1/ +numberOfFiles = len(process.source.fileNames) +numberOfJobs = options.nJobs +jobNumber = options.job + +process.source.fileNames = process.source.fileNames[jobNumber:numberOfFiles:numberOfJobs] +#if options.nJobs > 1: +print ("running over these files:") +print (process.source.fileNames) + +process.globalInfo = cms.EDAnalyzer('globalInfo', + lheInfo = cms.InputTag("externalLHEProducer"), + useLHEWeights=cms.bool(options.lheWeights) + ) + +process.source.skipEvents = cms.untracked.uint32(options.skipEvents) +process.maxEvents = cms.untracked.PSet( + input=cms.untracked.int32(options.maxEvents) +) + +if int(release.replace("_", "")) >= 8400 or int(release.replace("_", "")) == 8029: + bTagInfos = [ + 'pfImpactParameterTagInfos', + 'pfInclusiveSecondaryVertexFinderTagInfos', + 'pfDeepCSVTagInfos'] +else: + bTagInfos = [ + 'pfImpactParameterTagInfos', + 'pfInclusiveSecondaryVertexFinderTagInfos', + 'deepNNTagInfos', + ] + +if int(release.replace("_", "")) >= 8400 or int(release.replace("_", "")) == 8029: + bTagDiscriminators = [ + 'softPFMuonBJetTags', + 'softPFElectronBJetTags', + 'pfJetBProbabilityBJetTags', + 'pfJetProbabilityBJetTags', + 'pfCombinedInclusiveSecondaryVertexV2BJetTags', + 'pfDeepCSVJetTags:probudsg', # to be fixed with new names + 'pfDeepCSVJetTags:probb', + 'pfDeepCSVJetTags:probc', + 'pfDeepCSVJetTags:probbb', + 'pfDeepCSVJetTags:probcc', + ] +else: + bTagDiscriminators = [ + 'softPFMuonBJetTags', + 'softPFElectronBJetTags', + 'pfJetBProbabilityBJetTags', + 'pfJetProbabilityBJetTags', + 'pfCombinedInclusiveSecondaryVertexV2BJetTags', + 'deepFlavourJetTags:probudsg', # to be fixed with new names + 'deepFlavourJetTags:probb', + 'deepFlavourJetTags:probc', + 'deepFlavourJetTags:probbb', + 'deepFlavourJetTags:probcc', + ] + +###### semilep selection +outFileName = options.outputFile + '_' + str(options.job) + '.root' +print ('Using output file ' + outFileName) + +if options.deepNtuplizer: + process.TFileService = cms.Service("TFileService", + fileName=cms.string(outFileName)) +else: + process.MINIAODSIMEventContent.outputCommands.extend([ + 'keep *_goodElectrons_*_*', + 'keep *_goodMuons_*_*', + 'keep *_goodJets_*_*', + ]) + + process.outmod = cms.OutputModule("PoolOutputModule", + process.MINIAODSIMEventContent, + SelectEvents = cms.untracked.PSet( + SelectEvents = cms.vstring('p') + ), + fileName = cms.untracked.string(outFileName) + ) + +from PhysicsTools.SelectorUtils.tools.vid_id_tools import * +dataFormat = DataFormat.MiniAOD +switchOnVIDElectronIdProducer(process, dataFormat) +my_id_modules = [ + 'RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_Summer16_80X_V1_cff' +] +for idmod in my_id_modules: + setupAllVIDIdsInModule(process,idmod,setupVIDElectronSelection) + + +#HighLevelTrigger +HLTlistSM = cms.vstring("HLT_IsoTkMu24_v*", #Run B-H + "HLT_IsoMu24_v*" #Run B-H + ) +process.TriggerSel = cms.EDFilter("HLTHighLevel", + TriggerResultsTag = cms.InputTag("TriggerResults","","HLT"), + HLTPaths = cms.vstring(HLTlistSM), + eventSetupPathsKey = cms.string(''), # not empty => use read paths from AlCaRecoTriggerBitsRcd via this key + andOr = cms.bool(True), # how to deal with multiple triggers: True (OR) accept if ANY is true, False (AND) accept if ALL are true + throw = cms.bool(True) # throw exception on unknown path names +) +### Electron Selection +process.goodElectrons = cms.EDProducer("ElectronIdAdder", + src=cms.InputTag("slimmedElectrons"), + vSrc=cms.InputTag("offlineSlimmedPrimaryVertices"), + idMap=cms.InputTag("egmGsfElectronIDs:cutBasedElectronID-Summer16-80X-V1-tight"), + minPt=cms.double(20.0), + maxAbsEta=cms.double(2.4), + ) + +### Muon Selection +process.goodMuons = cms.EDProducer("MuonIdAdder", + src=cms.InputTag("slimmedMuons"), + vSrc=cms.InputTag("offlineSlimmedPrimaryVertices"), + minPt=cms.double(30.0), + maxAbsEta=cms.double(2.4), + maxRMI=cms.double(0.15) #RMI = relative muon isolation + ) + +### Jet Selection +process.goodJets = cms.EDProducer("PATJetCleaner", + src = cms.InputTag("slimmedJets"), + preselection = cms.string("pt>30 && abs(eta) < 2.4 && neutralHadronEnergyFraction < 0.99 " + "&& neutralEmEnergyFraction < 0.99 && (chargedMultiplicity+neutralMultiplicity) > 1 " + "&& chargedHadronEnergyFraction > 0.0 " + "&& chargedMultiplicity > 0.0 && chargedEmEnergyFraction < 0.99 "), + checkOverlaps = cms.PSet( + muons = cms.PSet( + src = cms.InputTag("goodMuons"), + algorithm = cms.string("byDeltaR"), + preselection = cms.string(""), + deltaR = cms.double(0.4), + checkRecoComponents = cms.bool(False), # don't check if they share some AOD object ref + pairCut = cms.string(""), + requireNoOverlaps = cms.bool(True), # overlaps cause the jet to be discarded + ), + electrons = cms.PSet( + src = cms.InputTag("goodElectrons"), + algorithm = cms.string("byDeltaR"), + preselection = cms.string(""), + deltaR = cms.double(0.4), + checkRecoComponents = cms.bool(False), # don't check if they share some AOD object ref + pairCut = cms.string(""), + requireNoOverlaps = cms.bool(True), # overlaps cause the jet to be discarded + ) + ), + finalCut = cms.string('') + ) + + +### Event Filter +process.ttsemilepFilter = cms.EDFilter("ttsemilepFilter", + src_muons = cms.InputTag("goodMuons"), + src_electrons = cms.InputTag("goodElectrons"), + src_jets = cms.InputTag("goodJets"), + src_mets = cms.InputTag("slimmedMETs"), + cut_minMT_leptonMETpair = cms.double(50.0), + nElectrons=cms.uint32(0), + nMuons=cms.uint32(1) + ) + + +### end selection + +# Jet Energy Corrections +if options.isData: + jetCorrections = cms.vstring(['L1FastJet', 'L2Relative', 'L3Absolute', 'L2L3Residual']) +else: + jetCorrections = cms.vstring(['L1FastJet', 'L2Relative', 'L3Absolute']) + +from PhysicsTools.PatAlgos.tools.jetTools import updateJetCollection +updateJetCollection( + process, + labelName="DeepFlavour", + # jetSource=cms.InputTag('slimmedJetsAK8PFPuppiSoftDropPacked', 'SubJets'), # 'subjets from AK8' + jetSource=cms.InputTag('goodJets'), # 'ak4Jets' + jetCorrections=('AK4PFchs', jetCorrections, 'None'), + pfCandidates=cms.InputTag('packedPFCandidates'), + pvSource=cms.InputTag("offlineSlimmedPrimaryVertices"), + svSource=cms.InputTag('slimmedSecondaryVertices'), + muSource=cms.InputTag('slimmedMuons'), + elSource=cms.InputTag('slimmedElectrons'), + btagInfos=bTagInfos, + btagDiscriminators=bTagDiscriminators, + explicitJTA=False + ) + +if hasattr(process, 'updatedPatJetsTransientCorrectedDeepFlavour'): + process.updatedPatJetsTransientCorrectedDeepFlavour.addTagInfos = cms.bool(True) + process.updatedPatJetsTransientCorrectedDeepFlavour.addBTagInfo = cms.bool(True) +else: + raise ValueError( + 'I could not find updatedPatJetsTransientCorrectedDeepFlavour to embed the tagInfos, please check the cfg') + +# QGLikelihood +process.load("DeepNTuples.DeepNtuplizer.QGLikelihood_cfi") +process.es_prefer_jec = cms.ESPrefer("PoolDBESSource", "QGPoolDBESSource") +process.load('RecoJets.JetProducers.QGTagger_cfi') +process.QGTagger.srcJets = cms.InputTag("selectedUpdatedPatJetsDeepFlavour") +process.QGTagger.jetsLabel = cms.string('QGL_AK4PFchs') + +from RecoJets.JetProducers.ak4GenJets_cfi import ak4GenJets + +process.ak4GenJetsWithNu = ak4GenJets.clone(src='packedGenParticles') + +## Filter out neutrinos from packed GenParticles +process.packedGenParticlesForJetsNoNu = cms.EDFilter("CandPtrSelector", src=cms.InputTag("packedGenParticles"), + cut=cms.string( + "abs(pdgId) != 12 && abs(pdgId) != 14 && abs(pdgId) != 16")) +## Define GenJets +process.ak4GenJetsRecluster = ak4GenJets.clone(src='packedGenParticlesForJetsNoNu') + +process.patGenJetMatchWithNu = cms.EDProducer("GenJetMatcher", # cut on deltaR; pick best by deltaR + src=cms.InputTag("selectedUpdatedPatJetsDeepFlavour"), + # RECO jets (any View is ok) + matched=cms.InputTag("ak4GenJetsWithNu"), + # GEN jets (must be GenJetCollection) + mcPdgId=cms.vint32(), # n/a + mcStatus=cms.vint32(), # n/a + checkCharge=cms.bool(False), # n/a + maxDeltaR=cms.double(0.4), # Minimum deltaR for the match + # maxDPtRel = cms.double(3.0), # Minimum deltaPt/Pt for the match (not used in GenJetMatcher) + resolveAmbiguities=cms.bool(True), + # Forbid two RECO objects to match to the same GEN object + resolveByMatchQuality=cms.bool(False), + # False = just match input in order; True = pick lowest deltaR pair first + ) + +process.patGenJetMatchRecluster = cms.EDProducer("GenJetMatcher", # cut on deltaR; pick best by deltaR + src=cms.InputTag("selectedUpdatedPatJetsDeepFlavour"), + # RECO jets (any View is ok) + matched=cms.InputTag("ak4GenJetsRecluster"), + # GEN jets (must be GenJetCollection) + mcPdgId=cms.vint32(), # n/a + mcStatus=cms.vint32(), # n/a + checkCharge=cms.bool(False), # n/a + maxDeltaR=cms.double(0.4), # Minimum deltaR for the match + # maxDPtRel = cms.double(3.0), # Minimum deltaPt/Pt for the match (not used in GenJetMatcher) + resolveAmbiguities=cms.bool(True), + # Forbid two RECO objects to match to the same GEN object + resolveByMatchQuality=cms.bool(False), + # False = just match input in order; True = pick lowest deltaR pair first + ) + +process.genJetSequence = cms.Sequence(process.packedGenParticlesForJetsNoNu + * process.ak4GenJetsWithNu + * process.ak4GenJetsRecluster + * process.patGenJetMatchWithNu + * process.patGenJetMatchRecluster) + +# Very Loose IVF SV collection +from PhysicsTools.PatAlgos.tools.helpers import loadWithPrefix + +loadWithPrefix(process, 'RecoVertex.AdaptiveVertexFinder.inclusiveVertexing_cff', "looseIVF") +process.looseIVFinclusiveCandidateVertexFinder.primaryVertices = cms.InputTag("offlineSlimmedPrimaryVertices") +process.looseIVFinclusiveCandidateVertexFinder.tracks = cms.InputTag("packedPFCandidates") +process.looseIVFinclusiveCandidateVertexFinder.vertexMinDLen2DSig = cms.double(0.) +process.looseIVFinclusiveCandidateVertexFinder.vertexMinDLenSig = cms.double(0.) +process.looseIVFinclusiveCandidateVertexFinder.fitterSigmacut = 20 + +process.looseIVFcandidateVertexArbitrator.primaryVertices = cms.InputTag("offlineSlimmedPrimaryVertices") +process.looseIVFcandidateVertexArbitrator.tracks = cms.InputTag("packedPFCandidates") +process.looseIVFcandidateVertexArbitrator.secondaryVertices = cms.InputTag("looseIVFcandidateVertexMerger") +process.looseIVFcandidateVertexArbitrator.fitterSigmacut = 20 + +#datapath=os.environ['CMSSW_BASE']+'/src/DeepNTuples/DeepNtuplizer/data/' +datapath='' + +# DeepNtuplizer configurations +process.load("DeepNTuples.DeepNtuplizer.DeepNtuplizer_cfi") +process.deepntuplizer.jets = cms.InputTag('selectedUpdatedPatJetsDeepFlavour') +process.deepntuplizer.bDiscriminators = bTagDiscriminators +process.deepntuplizer.bDiscriminators.append('pfCombinedMVAV2BJetTags') +process.deepntuplizer.LooseSVs = cms.InputTag("looseIVFinclusiveCandidateSecondaryVertices") +process.deepntuplizer.applySelection = cms.bool(options.selectJets) + +process.deepntuplizer.removeUndefined = cms.bool(False) +process.deepntuplizer.isData = cms.bool(options.isData) + + +process.deepntuplizer.sfMuons = cms.InputTag("goodMuons") + +if not options.isData: + process.deepntuplizer.useLHEWeights = cms.bool(options.lheWeights) + process.deepntuplizer.pileupData = cms.string(datapath + "pileup_data_2016.root") + process.deepntuplizer.pileupMC = cms.string(datapath + "pileup_MC_2016.root") + + process.deepntuplizer.periods=cms.vstring("2016BtoF","2016GH") + process.deepntuplizer.lumis=cms.vdouble(5.404+2.396+4.256+4.054+3.105, 7.179+8.746) + process.deepntuplizer.crossSection = cms.double(options.crossSection) + process.deepntuplizer.nEvents = cms.uint32(options.nEvents) + + process.deepntuplizer.sfTrigger_mu=cms.vstring(datapath + "EfficienciesAndSF_RunBtoF.root", + datapath + "EfficienciesAndSF_Period4.root") + process.deepntuplizer.sfTrigger_mu_Hist=cms.vstring("IsoMu24_OR_IsoTkMu24_PtEtaBins/abseta_pt_ratio", + "IsoMu24_OR_IsoTkMu24_PtEtaBins/abseta_pt_ratio") + + process.deepntuplizer.sfMuonId = cms.vstring(datapath+"EfficienciesAndSF_ID_BCDEF.root", + datapath+"EfficienciesAndSF_ID_GH.root") + process.deepntuplizer.sfMuonId_Hist = cms.vstring("MC_NUM_TightID_DEN_genTracks_PAR_pt_eta/abseta_pt_ratio", + "MC_NUM_TightID_DEN_genTracks_PAR_pt_eta/abseta_pt_ratio") + + process.deepntuplizer.sfMuonIso=cms.vstring(datapath+"EfficienciesAndSF_ISO_BCDEF.root", + datapath+"EfficienciesAndSF_ISO_GH.root") + process.deepntuplizer.sfMuonIso_Hist=cms.vstring("TightISO_TightID_pt_eta/abseta_pt_ratio", + "TightISO_TightID_pt_eta/abseta_pt_ratio") + + process.deepntuplizer.sfMuonTracking=cms.vstring(datapath+"Tracking_EfficienciesAndSF_BCDEFGH.root") + process.deepntuplizer.sfMuonTracking_Hist=cms.vstring("ratio_eff_aeta_dr030e030_corr") + + process.deepntuplizer.gluonReduction = cms.double(options.gluonReduction) + +if int(release.replace("_", "")) >= 840: + process.deepntuplizer.tagInfoName = cms.string('pfDeepCSV') + + +# 1631 +#process.ProfilerService = cms.Service( +# "ProfilerService", +# firstEvent=cms.untracked.int32(1631), +# lastEvent=cms.untracked.int32(1641), +# paths=cms.untracked.vstring('p') +#) + + +if options.deepNtuplizer: + if options.isData: + process.p = cms.Path(process.globalInfo + process.TriggerSel + process.ttsemilepFilter + process.QGTagger + process.deepntuplizer) + else: + process.p = cms.Path(process.globalInfo + process.TriggerSel + process.ttsemilepFilter + process.QGTagger + process.genJetSequence * process.deepntuplizer) +else: + process.p = cms.Path(process.TriggerSel + process.ttsemilepFilter) + + process.endp = cms.EndPath(process.outmod) + diff --git a/DeepNtuplizer/python/DeepNtuplizer_cfi.py b/DeepNtuplizer/python/DeepNtuplizer_cfi.py index 7994705291b..8a88a7048e1 100644 --- a/DeepNtuplizer/python/DeepNtuplizer_cfi.py +++ b/DeepNtuplizer/python/DeepNtuplizer_cfi.py @@ -1,4 +1,6 @@ import FWCore.ParameterSet.Config as cms +import os + deepntuplizer = cms.EDAnalyzer('DeepNtuplizer', vertices = cms.InputTag("offlineSlimmedPrimaryVertices"), @@ -7,6 +9,7 @@ jetR = cms.double(0.4), runFatJet = cms.bool(False), pupInfo = cms.InputTag("slimmedAddPileupInfo"), + lheInfo = cms.InputTag("externalLHEProducer"), rhoInfo = cms.InputTag("fixedGridRhoFastjetAll"), SVs = cms.InputTag("slimmedSecondaryVertices"), LooseSVs = cms.InputTag("inclusiveCandidateSecondaryVertices"), @@ -31,5 +34,47 @@ useHerwigCompatible=cms.bool(False), isHerwig=cms.bool(False), useOffsets=cms.bool(True), - applySelection=cms.bool(True) + applySelection=cms.bool(True), + + removeUndefined=cms.bool(True), + isData=cms.bool(False), + + #All following is for event weight computation + useLHEWeights=cms.bool(True), + + #leave empty string if you don't want to use pileup weights + pileupData=cms.string(""), + pileupMC=cms.string(""), + + periods=cms.vstring(), + lumis=cms.vdouble(), #weights are weighted with the luminosity of the period + crossSection=cms.double(1.0), + nEvents=cms.uint32(1), + + #to use different triggers for the periods + triggerToken=cms.InputTag("TriggerResults::HLT"), + triggers=cms.vstring(), + + #scalefactor information + sfMuons = cms.InputTag("slimmedMuons"), + sfElectrons=cms.InputTag("slimmedElectrons"), + + # The scalefactors can be computed on different ways: + # - empty vector to not use a specific scalefactor, + # - one string to use the same hist for all periods + # - number of strings equal to number of periods to use a different hist for each period + sfTrigger_mu=cms.vstring(), + sfTrigger_mu_Hist = cms.vstring(), + sfTrigger_e=cms.vstring(), + sfTrigger_e_Hist=cms.vstring(), + sfTrigger_emu=cms.vstring(), + sfTrigger_emu_Hist=cms.vstring(), + sfMuonId = cms.vstring(), + sfMuonId_Hist = cms.vstring(), + sfMuonIso=cms.vstring(), + sfMuonIso_Hist=cms.vstring(), + sfMuonTracking=cms.vstring(), + sfMuonTracking_Hist=cms.vstring(), + sfElIdAndIso=cms.vstring(), + sfElIdAndIso_Hist=cms.vstring(), ) diff --git a/DeepNtuplizer/python/LHEcounter.py b/DeepNtuplizer/python/LHEcounter.py new file mode 100644 index 00000000000..5aabb11fa17 --- /dev/null +++ b/DeepNtuplizer/python/LHEcounter.py @@ -0,0 +1,93 @@ +### +# Program to produce a histogram of the lhe weights to compute the effective event number of a sample with lhe weights +### +import FWCore.ParameterSet.Config as cms + +import FWCore.ParameterSet.VarParsing as VarParsing +### parsing job options +import sys + + +options = VarParsing.VarParsing() + +options.register('inputScript','',VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.string,"input Script") +options.register('outputFile','output',VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.string,"output File (w/o .root)") +options.register('maxEvents',-1,VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.int,"maximum events") +options.register('skipEvents', 0, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "skip N events") +options.register('job', 0, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "job number") +options.register('nJobs', 1, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "total jobs") + + +import os +release=os.environ['CMSSW_VERSION'][6:11] +print("Using release "+release) + +options.register( + 'inputFiles','', + VarParsing.VarParsing.multiplicity.list, + VarParsing.VarParsing.varType.string, + "input files " + ) + +if hasattr(sys, "argv"): + options.parseArguments() + + +process = cms.Process("lheWeightCounter") + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) ) + +process.load('FWCore.MessageService.MessageLogger_cfi') +process.MessageLogger.cerr.FwkReport.reportEvery = 10 +if options.inputScript == '': #this is probably for testing + process.MessageLogger.cerr.FwkReport.reportEvery = 100 + +process.options = cms.untracked.PSet( + allowUnscheduled = cms.untracked.bool(True), + wantSummary=cms.untracked.bool(True) +) + + +sampleListFile = 'DeepNTuples.DeepNtuplizer.samples.singleMuon_2016_cfg' +process.load(sampleListFile) #default input + +if options.inputFiles: + process.source.fileNames = options.inputFiles + +if options.inputScript != '' and options.inputScript != sampleListFile: + process.load(options.inputScript) + +#process.source.fileNames=['file:/afs/cern.ch/work/d/dwalter/data/ttbar/data.root'] #store/data/Run2016H/SingleMuon/MINIAOD/18Apr2017-v1/00000/00E02A09-853C-E711-93FF-3417EBE644A7.root + + + +numberOfFiles = len(process.source.fileNames) +numberOfJobs = options.nJobs +jobNumber = options.job + +process.source.fileNames = process.source.fileNames[jobNumber:numberOfFiles:numberOfJobs] +if options.nJobs > 1: + print ("running over these files:") + print (process.source.fileNames) + +process.source.skipEvents = cms.untracked.uint32(options.skipEvents) +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32 (options.maxEvents) +) + +process.options = cms.untracked.PSet( + SkipEvent = cms.untracked.vstring('ProductNotFound') +) + + +process.TFileService = cms.Service("TFileService", + fileName=cms.string("lheWeights_"+str(options.job) +".root"), + closeFileFast=cms.untracked.bool(True) + ) + +process.lheAnalyzer = cms.EDAnalyzer('lheAnalyzer', + lheInfo = cms.InputTag("externalLHEProducer") + ) + + +process.endp = cms.EndPath(process.lheAnalyzer) \ No newline at end of file diff --git a/DeepNtuplizer/python/PupMCHistProd.py b/DeepNtuplizer/python/PupMCHistProd.py new file mode 100644 index 00000000000..5b86114116b --- /dev/null +++ b/DeepNtuplizer/python/PupMCHistProd.py @@ -0,0 +1,88 @@ +### +# Program to produce a histogram of the number of interactions per bunch crossing of a monte carlo sample, uses puAnalyzer.cc +# This histogram is needed to calculate the pileup reweighting +### +import FWCore.ParameterSet.Config as cms + +import FWCore.ParameterSet.VarParsing as VarParsing +### parsing job options +import sys + + +options = VarParsing.VarParsing() + +options.register('inputScript','',VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.string,"input Script") +options.register('outputFile','output',VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.string,"output File (w/o .root)") +options.register('maxEvents',-1,VarParsing.VarParsing.multiplicity.singleton,VarParsing.VarParsing.varType.int,"maximum events") +options.register('skipEvents', 0, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "skip N events") +options.register('job', 0, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "job number") +options.register('nJobs', 1, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.int, "total jobs") +options.register('isData', True, VarParsing.VarParsing.multiplicity.singleton, VarParsing.VarParsing.varType.bool, "switch off generator jets") + +import os +release=os.environ['CMSSW_VERSION'][6:11] +print("Using release "+release) + +options.register( + 'inputFiles','', + VarParsing.VarParsing.multiplicity.list, + VarParsing.VarParsing.varType.string, + "input files (default is the tt RelVal)" + ) + +if hasattr(sys, "argv"): + options.parseArguments() + + +process = cms.Process("pupMCHistProd") + + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) ) + +process.load('FWCore.MessageService.MessageLogger_cfi') +process.MessageLogger.cerr.FwkReport.reportEvery = 10 +if options.inputScript == '': #this is probably for testing + process.MessageLogger.cerr.FwkReport.reportEvery = 100 + +process.options = cms.untracked.PSet( + allowUnscheduled = cms.untracked.bool(True), + wantSummary=cms.untracked.bool(True) +) + + +sampleListFile = 'DeepNTuples.DeepNtuplizer.samples.singleMuon_2016_cfg' +process.load(sampleListFile) #default input + +if options.inputFiles: + process.source.fileNames = options.inputFiles + +if options.inputScript != '' and options.inputScript != sampleListFile: + process.load(options.inputScript) + +process.source.fileNames=['file:/afs/cern.ch/work/d/dwalter/data/ttbar/TT/output_0_1.root'] #store/data/Run2016H/SingleMuon/MINIAOD/18Apr2017-v1/00000/00E02A09-853C-E711-93FF-3417EBE644A7.root + +numberOfFiles = len(process.source.fileNames) +numberOfJobs = options.nJobs +jobNumber = options.job + +process.source.fileNames = process.source.fileNames[jobNumber:numberOfFiles:numberOfJobs] +if options.nJobs > 1: + print ("running over these files:") + print (process.source.fileNames) + + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32 (options.maxEvents) +) + + +process.TFileService = cms.Service("TFileService", + fileName=cms.string("MyMCPileupHistogram_"+str(options.job) +".root"), + closeFileFast=cms.untracked.bool(True) + ) +process.puAnalyzer = cms.EDAnalyzer("puAnalyzer", + pileupInfo=cms.InputTag("slimmedAddPileupInfo") + ) + + +process.endp = cms.EndPath( process.puAnalyzer) \ No newline at end of file diff --git a/DeepNtuplizer/python/PupWeightMaker.py b/DeepNtuplizer/python/PupWeightMaker.py new file mode 100644 index 00000000000..b25be612e1b --- /dev/null +++ b/DeepNtuplizer/python/PupWeightMaker.py @@ -0,0 +1,28 @@ +###### +# takes data and MC histograms of pileup distribution and computes the pileup weights +###### + +import ROOT +import pdb + +file_pileup_MC = ROOT.TFile( + "/afs/cern.ch/work/d/dwalter/DeepNTuples/CMSSW_8_0_29/src/DeepNTuples/DeepNtuplizer/data/MyMCPileupHist.root") +file_pileup_Data = ROOT.TFile( + "/afs/cern.ch/work/d/dwalter/DeepNTuples/CMSSW_8_0_29/src/DeepNTuples//DeepNtuplizer/data/MyDataPileupHist.root") + +hist_pileup_MC = file_pileup_MC.Get("pileup") +hist_pileup_Data = file_pileup_Data.Get("pileup") + +def makePileupWeights(hist_data, hist_mc): + pupWeights = [1.0] + hist_mc.Scale(1. / hist_mc.Integral()) + hist_data.Scale(1. / hist_data.Integral()) + for bin in range(1,hist_data.GetNbinsX()+2): + if hist_mc.GetBinContent(bin) != 0: + w = hist_data.GetBinContent(bin)/hist_mc.GetBinContent(bin) + pupWeights.append(w) + else: + pupWeights.append(1) + return pupWeights + +print makePileupWeights(hist_pileup_Data,hist_pileup_MC) \ No newline at end of file diff --git a/DeepNtuplizer/python/QGLikelihood_cfi.py b/DeepNtuplizer/python/QGLikelihood_cfi.py index 65caa455bbe..ffd7f4f022b 100644 --- a/DeepNtuplizer/python/QGLikelihood_cfi.py +++ b/DeepNtuplizer/python/QGLikelihood_cfi.py @@ -3,7 +3,8 @@ qgDatabaseVersion = 'cmssw8020_v2' -databasepath=os.environ['CMSSW_BASE']+'/src/DeepNTuples/DeepNtuplizer/data/QGL_cmssw8020_v2.db' +#databasepath=os.environ['CMSSW_BASE']+'/src/DeepNTuples/DeepNtuplizer/data/QGL_cmssw8020_v2.db' +databasepath='QGL_cmssw8020_v2.db' from CondCore.CondDB.CondDB_cfi import * QGPoolDBESSource = cms.ESSource("PoolDBESSource", diff --git a/DeepNtuplizer/src/mergeDescriptor.cc b/DeepNtuplizer/src/mergeDescriptor.cc index d7a2058e639..c5f0401bafd 100644 --- a/DeepNtuplizer/src/mergeDescriptor.cc +++ b/DeepNtuplizer/src/mergeDescriptor.cc @@ -25,6 +25,7 @@ #include "DeepNTuples/DeepNtuplizer/interface/ntuple_pfCands.h" #include "DeepNTuples/DeepNtuplizer/interface/ntuple_SV.h" #include "DeepNTuples/DeepNtuplizer/interface/ntuple_DeepVertex.h" +#include "DeepNTuples/DeepNtuplizer/interface/ntuple_eventInfo.h" static bool debug=true; @@ -133,6 +134,7 @@ std::vector mergeDescriptor::createChains( branchinfos.push_back(new ntuple_bTagVars()); branchinfos.push_back(new ntuple_pfCands()); // branchinfos.push_back(new ntuple_DeepVertex()); + branchinfos.push_back(new ntuple_eventInfo()); std::vector chains; for(size_t i=0;i("jetAbsEtaMin")); jetAbsEtaMax_=(iConfig.getParameter("jetAbsEtaMax")); + removeUndefined_=(iConfig.getParameter("removeUndefined")); + vector disc_names = iConfig.getParameter >("bDiscriminators"); for(auto& name : disc_names) { discriminators_[name] = 0.; @@ -57,6 +59,7 @@ void ntuple_JetInfo::initBranches(TTree* tree){ addBranch(tree,"isG",&isG_, "isG_/i"); addBranch(tree,"isUndefined",&isUndefined_, "isUndefined_/i"); addBranch(tree,"genDecay",&genDecay_, "genDecay_/f"); + addBranch(tree,"isRealData",&isRealData_, "isRealData_/i"); //truth labeling with fallback to physics definition for light/gluon/undefined of standard flavor definition addBranch(tree,"isPhysB",&isPhysB_, "isPhysB_/i"); @@ -137,15 +140,20 @@ void ntuple_JetInfo::initBranches(TTree* tree){ void ntuple_JetInfo::readEvent(const edm::Event& iEvent){ + isData_ = iEvent.isRealData(); + iEvent.getByToken(qglToken_, qglHandle); iEvent.getByToken(ptDToken_, ptDHandle); iEvent.getByToken(axis2Token_, axis2Handle); iEvent.getByToken(multToken_, multHandle); - iEvent.getByToken(genJetMatchReclusterToken_, genJetMatchRecluster); - iEvent.getByToken(genJetMatchWithNuToken_, genJetMatchWithNu); + if(!isData_){ + iEvent.getByToken(genJetMatchReclusterToken_, genJetMatchRecluster); + iEvent.getByToken(genJetMatchWithNuToken_, genJetMatchWithNu); - iEvent.getByToken(genParticlesToken_, genParticlesHandle); + iEvent.getByToken(genParticlesToken_, genParticlesHandle); + + } iEvent.getByToken(muonsToken_, muonsHandle); @@ -165,88 +173,92 @@ void ntuple_JetInfo::readEvent(const edm::Event& iEvent){ Bhadron_daughter_.clear(); //std::cout << " start search for a b in this event "<(genC); - - if((abs(gen.pdgId())>500&&abs(gen.pdgId())<600)||(abs(gen.pdgId())>5000&&abs(gen.pdgId())<6000)) { - - //std::cout<0){ - - if( (abs(gen.daughter(0)->pdgId())>500&&abs(gen.daughter(0)->pdgId())<600)||(abs(gen.daughter(0)->pdgId())>5000&&abs(gen.daughter(0)->pdgId())<6000)) - { - if(gen.daughter(0)->numberOfDaughters()>0) - { - - const reco::GenParticle &daughter_ = static_cast< const reco::GenParticle &>(*(gen.daughter(0)->daughter(0))); - - if(daughter_.vx()!=gen.vx()) - { - Bhadron_daughter_.push_back(daughter_); - } - // else { - // std::cout << "only b daughters " << endl; - // } - } - else Bhadron_daughter_.push_back(gen); - - } - else{ - // std::cout<vx()<< " oh " <pdgId() <(*gen.daughter(0)); - Bhadron_daughter_.push_back(daughter_); - } - - }// if daughter is there - else { - - //std::cout << " lonly B hadron, has NO daughter??? "<(genC); - if(abs(gen.pdgId())==12||abs(gen.pdgId())==14||abs(gen.pdgId())==16) { - const reco::GenParticle* mother = static_cast< const reco::GenParticle*> (gen.mother()); - if(mother!=NULL) { - if((abs(mother->pdgId())>500&&abs(mother->pdgId())<600)||(abs(mother->pdgId())>5000&&abs(mother->pdgId())<6000)) { - neutrinosLepB.emplace_back(gen); + if(!isData_){ + for (const reco::Candidate &genC : *genParticlesHandle){ + + const reco::GenParticle &gen = static_cast< const reco::GenParticle &>(genC); + + if((abs(gen.pdgId())>500&&abs(gen.pdgId())<600)||(abs(gen.pdgId())>5000&&abs(gen.pdgId())<6000)) { + + //std::cout<0){ + + if( (abs(gen.daughter(0)->pdgId())>500&&abs(gen.daughter(0)->pdgId())<600)||(abs(gen.daughter(0)->pdgId())>5000&&abs(gen.daughter(0)->pdgId())<6000)) + { + if(gen.daughter(0)->numberOfDaughters()>0) + { + + const reco::GenParticle &daughter_ = static_cast< const reco::GenParticle &>(*(gen.daughter(0)->daughter(0))); + + if(daughter_.vx()!=gen.vx()) + { + Bhadron_daughter_.push_back(daughter_); + } + else { + // std::cout << "only b daughters " << endl; + Bhadron_daughter_.push_back(gen); + } + } + else Bhadron_daughter_.push_back(gen); + + } + else{ + // std::cout<vx()<< " oh " <pdgId() <(*gen.daughter(0)); + Bhadron_daughter_.push_back(daughter_); + } + + }// if daughter is there + else { + + //std::cout << " lonly B hadron, has NO daughter??? "<pdgId())>400&&abs(mother->pdgId())<500)||(abs(mother->pdgId())>4000&&abs(mother->pdgId())<5000)) { - neutrinosLepB_C.emplace_back(gen); + } + } + + + for (const reco::Candidate &genC : *genParticlesHandle) { + const reco::GenParticle &gen = static_cast< const reco::GenParticle &>(genC); + if(abs(gen.pdgId())==12||abs(gen.pdgId())==14||abs(gen.pdgId())==16) { + const reco::GenParticle* mother = static_cast< const reco::GenParticle*> (gen.mother()); + if(mother!=NULL) { + if((abs(mother->pdgId())>500&&abs(mother->pdgId())<600)||(abs(mother->pdgId())>5000&&abs(mother->pdgId())<6000)) { + neutrinosLepB.emplace_back(gen); + } + if((abs(mother->pdgId())>400&&abs(mother->pdgId())<500)||(abs(mother->pdgId())>4000&&abs(mother->pdgId())<5000)) { + neutrinosLepB_C.emplace_back(gen); + } + } + else { + std::cout << "No mother" << std::endl; } } - else { - std::cout << "No mother" << std::endl; + + int id(std::abs(gen.pdgId())); + int status(gen.status()); + + if (id == 21 && status >= 21 && status <= 59) { //// Pythia8 hard scatter, ISR, or FSR + if ( gen.numberOfDaughters() == 2 ) { + const reco::Candidate* d0 = gen.daughter(0); + const reco::Candidate* d1 = gen.daughter(1); + if ( std::abs(d0->pdgId()) == 5 && std::abs(d1->pdgId()) == 5 + && d0->pdgId()*d1->pdgId() < 0 && reco::deltaR(*d0, *d1) < 0.4) gToBB.push_back(gen) ; + if ( std::abs(d0->pdgId()) == 4 && std::abs(d1->pdgId()) == 4 + && d0->pdgId()*d1->pdgId() < 0 && reco::deltaR(*d0, *d1) < 0.4) gToCC.push_back(gen) ; + } } - } - int id(std::abs(gen.pdgId())); - int status(gen.status()); - - if (id == 21 && status >= 21 && status <= 59) { //// Pythia8 hard scatter, ISR, or FSR - if ( gen.numberOfDaughters() == 2 ) { - const reco::Candidate* d0 = gen.daughter(0); - const reco::Candidate* d1 = gen.daughter(1); - if ( std::abs(d0->pdgId()) == 5 && std::abs(d1->pdgId()) == 5 - && d0->pdgId()*d1->pdgId() < 0 && reco::deltaR(*d0, *d1) < 0.4) gToBB.push_back(gen) ; - if ( std::abs(d0->pdgId()) == 4 && std::abs(d1->pdgId()) == 4 - && d0->pdgId()*d1->pdgId() < 0 && reco::deltaR(*d0, *d1) < 0.4) gToCC.push_back(gen) ; + if(id == 15 && false){ + alltaus_.push_back(gen); } - } - if(id == 15 && false){ - alltaus_.push_back(gen); } - + //technically a branch fill but per event, therefore here } - //technically a branch fill but per event, therefore here } //use either of these functions @@ -263,11 +275,12 @@ bool ntuple_JetInfo::fillBranches(const pat::Jet & jet, const size_t& jetidx, co if ( fabs(jet.eta()) < jetAbsEtaMin_ || fabs(jet.eta()) > jetAbsEtaMax_ ) returnval=false; // apply jet eta cut - // often we have way to many gluons that we do not need. This randomply reduces the gluons - if (gluonReduction_>0 && jet.partonFlavour()==21) - if(TRandom_.Uniform()>gluonReduction_) returnval=false; + if(!isData_){ + // often we have way to many gluons that we do not need. This randomply reduces the gluons + if (gluonReduction_>0 && jet.partonFlavour()==21) + if(TRandom_.Uniform()>gluonReduction_) returnval=false; + } - if(jet.genJet()==NULL)returnval=false; //branch fills @@ -277,12 +290,15 @@ bool ntuple_JetInfo::fillBranches(const pat::Jet & jet, const size_t& jetidx, co npv_ = vertices()->size(); - for (auto const& v : *pupInfo()) { - int bx = v.getBunchCrossing(); - if (bx == 0) { - ntrueInt_ = v.getTrueNumInteractions(); + if(!isData_){ + for (auto const& v : *pupInfo()) { + int bx = v.getBunchCrossing(); + if (bx == 0) { + ntrueInt_ = v.getTrueNumInteractions(); + } } } + rho_ = rhoInfo()[0]; @@ -299,7 +315,7 @@ bool ntuple_JetInfo::fillBranches(const pat::Jet & jet, const size_t& jetidx, co //std::vector > p= coll->ptrs(); isB_=0; isGBB_=0; isBB_=0; isC_=0; isGCC_=0; isCC_=0; isUD_=0;isTau_=0; - isS_=0; isG_=0, isLeptonicB_=0, isLeptonicB_C_=0, isUndefined_=0; + isS_=0; isG_=0, isLeptonicB_=0, isLeptonicB_C_=0, isUndefined_=0, isRealData_=0; auto muIds = deep_ntuples::jet_muonsIds(jet,*muonsHandle); auto elecIds = deep_ntuples::jet_electronsIds(jet,*electronsHandle); @@ -372,11 +388,15 @@ bool ntuple_JetInfo::fillBranches(const pat::Jet & jet, const size_t& jetidx, co } } - if(!jet.genJet()){//for data - isUndefined_=1;isPhysUndefined_=1; + if(!jet.genJet() || jet.genJet()==NULL){ + isUndefined_=1; + isPhysUndefined_=1; + if(removeUndefined_) returnval=false; } - if(isUndefined_ && isPhysUndefined_) returnval=false; //skip event, if neither standard flavor definition nor physics definition fallback define a "proper flavor" + if(isData_){ + isRealData_=1; + } pat::JetCollection h; @@ -470,9 +490,6 @@ bool ntuple_JetInfo::fillBranches(const pat::Jet & jet, const size_t& jetidx, co y_axis2_ = std::get<5>(qgtuple); y_pt_dr_log_=std::get<6>(qgtuple); - - - return returnval; } diff --git a/DeepNtuplizer/src/ntuple_eventInfo.cc b/DeepNtuplizer/src/ntuple_eventInfo.cc new file mode 100644 index 00000000000..df93764e3dc --- /dev/null +++ b/DeepNtuplizer/src/ntuple_eventInfo.cc @@ -0,0 +1,310 @@ +/* + * ntuple_eventInfo.cc + * + * Created on: 20 Feb 2018 + * Author: dwalter + */ + +#include "../interface/ntuple_eventInfo.h" + + + +void ntuple_eventInfo::getInput(const edm::ParameterSet& iConfig){ + + useLHEWeights_ = (iConfig.getParameter("useLHEWeights")); + + periods = (iConfig.getParameter>("periods")); + lumis = (iConfig.getParameter>("lumis")); + sigma = (iConfig.getParameter("crossSection")); + nEvents = (iConfig.getParameter("nEvents")); + + pupDataDir_ = (iConfig.getParameter("pileupData")); + pupMCDir_ = (iConfig.getParameter("pileupMC")); + + sfTrigger_mu_Dir_ = (iConfig.getParameter>("sfTrigger_mu")); + sfTrigger_mu_Name_ = (iConfig.getParameter>("sfTrigger_mu_Hist")); + sfTrigger_e_Dir_ = (iConfig.getParameter>("sfTrigger_e")); + sfTrigger_e_Name_ = (iConfig.getParameter>("sfTrigger_e_Hist")); + sfTrigger_emu_Dir_ = (iConfig.getParameter>("sfTrigger_emu")); + sfTrigger_emu_Name_ = (iConfig.getParameter>("sfTrigger_emu_Hist")); + sfMuonId_Dir_ = (iConfig.getParameter>("sfMuonId")); + sfMuonId_Name_ = (iConfig.getParameter>("sfMuonId_Hist")); + sfMuonIso_Dir_ = (iConfig.getParameter>("sfMuonIso")); + sfMuonIso_Name_ = (iConfig.getParameter>("sfMuonIso_Hist")); + sfElIdAndIso_Dir_ = (iConfig.getParameter>("sfElIdAndIso")); + sfElIdAndIso_Name_ = (iConfig.getParameter>("sfElIdAndIso_Hist")); + + sfMuonTracking_Dir_ = (iConfig.getParameter>("sfMuonTracking")); + sfMuonTracking_Name_ = (iConfig.getParameter>("sfMuonTracking_Hist")); + + triggers = (iConfig.getParameter>("triggers")); + + + if(pupDataDir_=="" || pupMCDir_=="") + std::cout<<"no pileup histograms, proceed without pileup reweighting. \n"; + else{ + TFile *pupMCFile = new TFile(pupMCDir_.c_str()); + TFile *pupDataFile = new TFile(pupDataDir_.c_str()); + + TH1F *pupMCHist = (TH1F*)pupMCFile->Get("pileup"); + TH1F *pupDataHist = (TH1F*)pupDataFile->Get("pileup"); + + pupMCHist->Scale(1./pupMCHist->Integral()); + pupDataHist->Scale(1./pupDataHist->Integral()); + + pupWeights.push_back(1.0); + for(int bin = 1; bin < pupMCHist->GetNbinsX() + 1; bin++){ + pupWeights.push_back(pupDataHist->GetBinContent(bin)/pupMCHist->GetBinContent(bin)); + } + } + + initializeScalefactor(sfTrigger_mu_Dir_, sfTrigger_mu_Name_, &sfTrigger_mu_Hist, periods); + initializeScalefactor(sfTrigger_e_Dir_, sfTrigger_e_Name_, &sfTrigger_e_Hist, periods); + initializeScalefactor(sfTrigger_emu_Dir_, sfTrigger_emu_Name_, &sfTrigger_emu_Hist, periods); + initializeScalefactor(sfMuonId_Dir_, sfMuonId_Name_, &sfMuonId_Hist, periods); + initializeScalefactor(sfMuonIso_Dir_, sfMuonIso_Name_, &sfMuonIso_Hist, periods); + initializeScalefactor(sfElIdAndIso_Dir_, sfElIdAndIso_Name_, &sfElIdAndIso_Hist, periods); + initializeScalefactor(sfMuonTracking_Dir_, sfMuonTracking_Name_, &sfMuonTracking_Hist, periods); + +} + + +void ntuple_eventInfo::initBranches(TTree* tree){ + + addBranch(tree,"event_weight", &event_weight_,"event_weight_/f" ); +} + + +void ntuple_eventInfo::readEvent(const edm::Event& iEvent){ + + if(!iEvent.isRealData()){ + TriggerInfo triggerInfo(iEvent,triggerToken_); + + + iEvent.getByToken(lheToken_, lheInfo); + iEvent.getByToken(muonToken_, muons); + iEvent.getByToken(electronToken_, electrons); + + ///// pileup weights + for (auto const& v : *pupInfo()) { + int bx = v.getBunchCrossing(); + if (bx == 0) { + ntrueInt_ = v.getTrueNumInteractions(); + } + } + double pupWeight = 1.; + if(ntrueInt_ < pupWeights.size()){ + pupWeight = pupWeights.at(ntrueInt_); + } + + ///// lhe weights + double lheWeight = 1.; + if(useLHEWeights_){ + lheWeight = lheInfo->weights()[0].wgt/std::abs(lheInfo->weights()[0].wgt); + } + + ///// scalefactors + + double leadingMuon_pt = 0.; + double leadingMuon_eta = 0.; + for (size_t i = 0; i < muons->size(); ++i){ + const auto & muon = (*muons).at(i); + if(muon.pt() > leadingMuon_pt){ + leadingMuon_pt = muon.pt(); + leadingMuon_eta = muon.eta(); + } + } + double leadingElectron_pt = 0.; + double leadingElectron_sueta = 0.; + double leadingElectron_eta = 0.; + for (size_t i = 0; i < electrons->size(); ++i){ + const auto & electron = (*electrons).at(i); + if(electron.pt() > leadingElectron_pt){ + leadingElectron_pt = electron.pt(); + leadingElectron_sueta = electron.superCluster()->eta(); + leadingElectron_eta = electron.eta(); + } + } + + + event_weight_ = 0.; + for(unsigned int i=0; i< periods.size(); i++){ + + double isf = 1.; + bool itriggered = 1; + + isf *= getScalefactor(std::abs(leadingMuon_eta), leadingMuon_pt, sfTrigger_mu_Hist, i); + isf *= getScalefactor(leadingElectron_pt, std::abs(leadingElectron_sueta),sfTrigger_e_Hist, i, true); + isf *= getScalefactor(std::abs(leadingElectron_eta), std::abs(leadingMuon_eta), sfTrigger_emu_Hist, i); + isf *= getScalefactor(std::abs(leadingMuon_eta), leadingMuon_pt, sfMuonId_Hist, i); + isf *= getScalefactor(std::abs(leadingMuon_eta), leadingMuon_pt, sfMuonIso_Hist, i); + isf *= getScalefactor(std::abs(leadingElectron_sueta), leadingElectron_pt, sfElIdAndIso_Hist, i); + isf *= getScalefactor(std::abs(leadingMuon_eta), sfMuonTracking_Hist, i); + + if(triggers.size() != 0){ + itriggered = triggerInfo.IsAnyTriggered(triggers[i]); + //std::cout<<"istriggered "< hists, unsigned int period){ + + if(hists.size()==0) return 1.; + else if(hists.size()==1){ + int bin = hists[0]->GetXaxis()->FindBin(x); + if(x > hists[0]->GetXaxis()->GetXmax()) //dont take overflow bins, but the last ones + bin -= 1; + return(hists[0]->GetBinContent(bin)); + } + else if(hists.size() > period){ + int bin = hists[period]->GetXaxis()->FindBin(x); + if(x > hists[period]->GetXaxis()->GetXmax()) //dont take overflow bins, but the last ones + bin -= 1; + return(hists[period]->GetBinContent(bin)); + } + else std::cout<<"ERROR: something went wrong in getScalefactor function"< dirs, std::vector names, std::vector *sfHists, std::vector periods){ + if(dirs.size()==0) + std::cout<<"no scalefactor "<push_back((TH2F*)file->Get(names[i].c_str())); + } + } + else if(dirs.size()==1){ + std::cout<<"initialize one scalefactor histogram for all periods " <push_back((TH2F*)file->Get(names[0].c_str())); + } + else{ + std::cout<<"ERROR: something went wrong by initialization of a scalefactor " < dirs, std::vector names, std::vector *sfHists, std::vector periods){ + if(dirs.size()==0) + std::cout<<"no scalefactor "<Get(names[i].c_str()); + readHistoFromGraph(sfTGraph, &hist, "sfMuonTrackingHist"); + sfHists->push_back(hist); + } + } + else if(dirs.size()==1){ + std::cout<<"initialize one scalefactor histogram for all periods " <Get(names[0].c_str()); + readHistoFromGraph(sfTGraph, &hist, "sfMuonTrackingHist"); + sfHists->push_back(hist); + } + else{ + std::cout<<"ERROR: something went wrong by initialization of a scalefactor " <GetN(); + const double* x_centers = graph->GetX(); + const double* y_centers = graph->GetY(); + double x_lows[npoints]; + double x_highs[npoints]; + double y_lows[npoints]; + double y_highs[npoints]; + for(int i=0; iGetErrorXlow(i); + x_highs[i] = graph->GetErrorXhigh(i); + y_lows[i] = graph->GetErrorYlow(i); + y_highs[i] = graph->GetErrorYhigh(i); + } + + double x_edges[npoints+1]; + for(int i=0; iSetDirectory(0); // without this the histo will get deleted when a currently open TFile is closed + + for(int i=0; iSetBinContent(i+1, y_centers[i]); + (*h)->SetBinError(i+1, std::max(y_lows[i], y_highs[i])); +// cout << (*h)->GetBinError(i) << endl; + } +} diff --git a/DeepNtuplizer/src/ntuple_pfCands.cc b/DeepNtuplizer/src/ntuple_pfCands.cc index 88333895b1f..a391e0dbde1 100644 --- a/DeepNtuplizer/src/ntuple_pfCands.cc +++ b/DeepNtuplizer/src/ntuple_pfCands.cc @@ -273,25 +273,25 @@ bool ntuple_pfCands::fillBranches(const pat::Jet & jet, const size_t& jetidx, co -mindrsvpfcand(PackedCandidate), PackedCandidate->pt()/jet_uncorr_pt)); } else{ - sortedneutrals.push_back(sorting::sortingClass - (i, -1, -mindrsvpfcand(PackedCandidate), PackedCandidate->pt()/jet_uncorr_pt)); + sorting::sortingClass ineutral = sorting::sortingClass(i, -1, -mindrsvpfcand(PackedCandidate), PackedCandidate->pt()/jet_uncorr_pt); + sortedneutrals.push_back(ineutral); } } } + std::sort(sortedcharged.begin(),sortedcharged.end(),sorting::sortingClass::compareByABCInv); std::sort(sortedneutrals.begin(),sortedneutrals.end(),sorting::sortingClass::compareByABCInv); + // counts neutral and charged candicates n_Cpfcand_ = std::min(sortedcharged.size(),max_pfcand_); n_Npfcand_ = std::min(sortedneutrals.size(),max_pfcand_); std::vector sortedchargedindices,sortedneutralsindices; - sortedchargedindices=sorting::invertSortingVector(sortedcharged); - sortedneutralsindices=sorting::invertSortingVector(sortedneutrals); - - + sortedchargedindices=sorting::invertSortingVector(sortedcharged); + sortedneutralsindices=sorting::invertSortingVector(sortedneutrals); @@ -459,6 +459,11 @@ bool ntuple_pfCands::fillBranches(const pat::Jet & jet, const size_t& jetidx, co nCpfcand_=n_Cpfcand_; nNpfcand_=n_Npfcand_; + if(jet.eta() == -1.039){ + std::cout<<" end pfCands ..."<& triggerBitsToken + ){ + + edm::Handle h_triggerBits; + iEvent.getByToken(triggerBitsToken, h_triggerBits); + // todo : use trigger objects + // edm::Handle h_triggerObjects; + // iEvent.getByToken(triggerObjectsToken, h_triggerObjects); + + const edm::TriggerNames &names = iEvent.triggerNames(*h_triggerBits); + + for (unsigned int i = 0; i < h_triggerBits->size(); ++i) { + string name=names.triggerName(i); +// cout << "Trigger #" << i << " " << name << endl; + triggernames.push_back(name); + triggers[name]=h_triggerBits->accept(i); + } +} + +bool TriggerInfo::Exists(std::string triggername) const { + if(triggers.count(triggername)==0){ + return false; + } + return true; +} + + +bool TriggerInfo::IsTriggered(std::string triggername) const { + if(triggername=="any"||triggername=="Any"||triggername=="none"||triggername=="None") return true; + int asterix = triggername.find("*"); + if(asterix==int(std::string::npos)){ + return Exists(triggername)&&triggers.at(triggername); + } + else{ + if(asterix==int(triggername.length()-1)){ + triggername.pop_back(); + auto it=triggers.lower_bound(triggername); + while(it!=triggers.end()&&it->first.find(triggername) == 0){ + if(it->second) return true; + it++; + } + return false; + } + else{ + cerr << "* wildcard only allowed at the end" << endl; + return false; + } + } + return false; + +} + +bool TriggerInfo::IsAnyTriggered(std::vector< std::string > triggerNames) const { + for(auto name=triggerNames.begin(); name!=triggerNames.end();name++){ + if(IsTriggered(*name)) return true; + } + return false; +} + +bool TriggerInfo::IsAnyTriggered(std::string triggerstring) const{ + std::vector triggernames; + + int posOr; + while(true){ + posOr = triggerstring.find("||"); + if(posOr==int(std::string::npos)){ + triggernames.push_back(triggerstring); + return IsAnyTriggered(triggernames); + } + triggernames.push_back(triggerstring.substr(0,posOr)); + triggerstring=triggerstring.substr(posOr+2,triggerstring.length()); + //std::cout<<"triggerstring = "< TriggerInfo::GetTriggers() const{ + return triggers; +} + + +void TriggerInfo::ListTriggers() { + int counter=0; + for(unsigned int i=0; i