From c23a321638b99e6fb0371378dfdc0aa43c8f7cc5 Mon Sep 17 00:00:00 2001 From: Dylan Teague Date: Mon, 11 Nov 2019 05:15:43 -0500 Subject: [PATCH] Big commit of 3top changes (other changes removed for PR) commit before testing new branch stuff off shoring selection to analysis, fixing up fake rate code (need to commit more often...) rm unnecessary TTTSelector and fixed addFile funcitons changed logging to Fireworks version and made changes for PR fixing fake rate things and changing logging a bit (less messages) Add fillHist (run analysis as script) so can rebase --- Templates/config.dteague | 1 + Utilities/python/CombineCardTools.py | 7 +- Utilities/python/ConfigureJobs.py | 58 +- Utilities/python/SelectorTools.py | 14 +- Utilities/python/UserInput.py | 2 +- Utilities/scripts/addUnrolledHistsToFile.py | 4 +- Utilities/scripts/makeFakeRates.py | 137 +++-- Utilities/scripts/prepareCombine.py | 4 +- Utilities/scripts/printEventDetails.py | 2 +- fillHists.sh | 20 +- interface/BranchManager.h | 23 +- interface/FakeRateSelector.h | 62 ++- interface/HelperFunctions.h | 42 +- interface/SelectorBase.h | 122 ++--- interface/TTTSelector.h | 149 ----- interface/ThreeLepSelector.h | 174 +++--- interface/WZSelector.h | 4 +- interface/WZSelectorBase.h | 6 + interface/dynEnum.h | 20 + src/BuildFile.xml | 1 + src/FakeRateSelector.cc | 158 ++++-- src/SelectorBase.cc | 120 ++-- src/TTTSelector.cc | 484 ----------------- src/ThreeLepSelector.cc | 572 +++++++++++++------- src/WZSelector.cc | 12 +- src/WZSelectorBase.cc | 17 +- src/classes.h | 2 - src/classes_def.xml | 1 - 28 files changed, 986 insertions(+), 1232 deletions(-) delete mode 100644 interface/TTTSelector.h create mode 100644 interface/dynEnum.h delete mode 100644 src/TTTSelector.cc diff --git a/Templates/config.dteague b/Templates/config.dteague index ed12f3d8..7a616df6 100644 --- a/Templates/config.dteague +++ b/Templates/config.dteague @@ -2,6 +2,7 @@ user = dteague analysis = 3top dataset_manager_path = /afs/cern.ch/work/d/%(user)s/ +dataset_manager_name = AnalysisDatasetManager combine_path = %(dataset_manager_path)s/HiggsCombine/CMSSW_8_1_0/src/HiggsAnalysis/CombinedLimit output_path = /afs/cern.ch/work/d/%(user)s/testHold fakeRate_output = %(output_path)s/FakeRates diff --git a/Utilities/python/CombineCardTools.py b/Utilities/python/CombineCardTools.py index c5bb1029..2f909fa7 100644 --- a/Utilities/python/CombineCardTools.py +++ b/Utilities/python/CombineCardTools.py @@ -80,7 +80,8 @@ def getVariations(self): def getVariationsForProcess(self, process): if process not in self.variations.keys(): - raise ValueError("Variations not defined for process %s" % process) + + raise ValueError("Variations not defined for process %s" % process) return self.variations[process] def setVariationsByProcess(self, process, variations): @@ -204,7 +205,7 @@ def setCombineChannels(self, groups): self.channelsToCombine = groups def combineChannels(self, group, processName, central=True): - variations = self.variations[group.GetName()][:] + variations = self.variations[group.GetName()][:] if group.GetName() in self.variations else [] fitVariable = self.getFitVariable(group.GetName()) if central: variations.append("") @@ -263,7 +264,7 @@ def loadHistsForProcess(self, processName, scaleNorm=1, expandedTheory=True): #TODO:Make optional processedHists = [] for chan in self.channels: - histName = "_".join([fitVariable, chan]) if chan != "all" else fitVariable + histName = "_".join([fitVariable, chan]) hist = group.FindObject(histName) if not hist: logging.warning("Failed to produce hist %s for process %s" % (histName, processName)) diff --git a/Utilities/python/ConfigureJobs.py b/Utilities/python/ConfigureJobs.py index 3c0c4a34..88e6eadf 100644 --- a/Utilities/python/ConfigureJobs.py +++ b/Utilities/python/ConfigureJobs.py @@ -181,17 +181,15 @@ def getListOfHDFSFiles(file_path): # TODO: Would be good to switch the order of the last two arguments # completely deprecate manager_path without breaking things -def getListOfFiles(filelist, selection, manager_path="", analysis=""): - if manager_path is "": - manager_path = getManagerPath() +def getListOfFiles(filelist, analysis="", input_tier=""): + manager_path = getManagerPath() data_path = "%s/%s/FileInfo" % (manager_path, getManagerName()) - group_path = "%s/AnalysisDatasetManager/PlotGroups" % manager_path + group_path = "%s/%s/PlotGroups" % (manager_path, getManagerName()) data_info = UserInput.readAllInfo("/".join([data_path, "data/*"])) mc_info = UserInput.readAllInfo("/".join([data_path, "montecarlo/*"])) - analysis_info = UserInput.readInfo("/".join([data_path, analysis, selection])) \ - if analysis != "" else [] - valid_names = (data_info.keys() + mc_info.keys()) if not analysis_info else analysis_info.keys() - group_names = UserInput.readAllInfo("%s/%s.py" %(group_path, analysis)) if analysis else dict() + + valid_names = UserInput.readInfo("/".join([data_path, analysis, input_tier])) if input_tier and analysis else (data_info.keys() + mc_info.keys()) + group_names = UserInput.readAllInfo("%s/%s.py" %(group_path, analysis)) names = [] for name in filelist: if ".root" in name: @@ -199,17 +197,8 @@ def getListOfFiles(filelist, selection, manager_path="", analysis=""): # Allow negative contributions elif name[0] == "-" : names.append(name) - elif "WZxsec2016" in name: - dataset_file = manager_path + \ - "%s/FileInfo/WZxsec2016/%s.json" % (getManagerPath(), selection) - allnames = json.load(open(dataset_file)).keys() - if "nodata" in name: - nodata = [x for x in allnames if "data" not in x] - names += nodata - elif "data" in name: - names += [x for x in allnames if "data" in x] - else: - names += allnames + elif "*" in name: + names += fnmatch.filter(valid_names, name) elif name in group_names: names += group_names[name]['Members'] elif "*" in name: @@ -247,11 +236,10 @@ def fillTemplatedFile(template_file_name, out_file_name, template_dict): with open(out_file_name, "w") as outFile: outFile.write(result) -def getListOfFilesWithXSec(filelist, manager_path="", selection="ntuples"): - if manager_path is "": - manager_path = getManagerPath() +def getListOfFilesWithXSec(filelist, analysis="", input_tier=""): + manager_path = getManagerPath() data_path = "%s/%s/FileInfo" % (manager_path, getManagerName()) - files = getListOfFiles(filelist, selection, manager_path) + files = getListOfFiles(filelist, analysis, input_tier) mc_info = UserInput.readAllInfo("/".join([data_path, "montecarlo/*"])) info = {} for file_name in files: @@ -269,21 +257,20 @@ def getListOfFilesWithXSec(filelist, manager_path="", selection="ntuples"): info.update({file_name : file_info["cross_section"]*kfac}) return info -def getListOfFilesWithPath(filelist, analysis, selection, das=True, manager_path=""): - if manager_path is "": - manager_path = getManagerPath() +def getListOfFilesWithPath(filelist, analysis, input_tier, das=True): + manager_path = getManagerPath() data_path = "%s/%s/FileInfo" % (manager_path, getManagerName()) - files = getListOfFiles(filelist, selection, manager_path, analysis) - selection_info = UserInput.readInfo("/".join([data_path, analysis, selection])) + files = getListOfFiles(filelist, analysis, input_tier) + input_tier_info = UserInput.readInfo("/".join([data_path, analysis, input_tier])) info = {} for file_name in files: - if das and "DAS" not in selection_info[file_name].keys(): - logging.error("DAS path not defined for file %s in analysis %s/%s" % (file_name, analysis, selection)) + if das and "DAS" not in input_tier_info[file_name].keys(): + logging.error("DAS path not defined for file %s in analysis %s/%s" % (file_name, analysis, input_tier)) continue - elif not das and "file_path" not in selection_info[file_name].keys(): - logging.error("File_path not defined for file %s in analysis %s/%s" % (file_name, analysis, selection)) + elif not das and "file_path" not in input_tier_info[file_name].keys(): + logging.error("File_path not defined for file %s in analysis %s/%s" % (file_name, analysis, input_tier)) continue - info.update({file_name : selection_info[file_name]["DAS" if das else "file_path"]}) + info.update({file_name : input_tier_info[file_name]["DAS" if das else "file_path"]}) return info def getPreviousStep(selection, analysis): @@ -331,9 +318,8 @@ def getConfigFileName(config_file_name): raise ValueError("Invalid configuration file. Tried to read %s which does not exist" % \ config_file_name) -def getInputFilesPath(sample_name, selection, analysis, manager_path=""): - if manager_path is "": - manager_path = getManagerPath() +def getInputFilesPath(sample_name, selection, analysis): + manager_path = getManagerPath() if ".root" in sample_name: logging.info("Using simple file %s" % sample_name) return sample_name diff --git a/Utilities/python/SelectorTools.py b/Utilities/python/SelectorTools.py index 778b3ce6..130405ab 100755 --- a/Utilities/python/SelectorTools.py +++ b/Utilities/python/SelectorTools.py @@ -21,10 +21,10 @@ def __init__(self, analysis, selection, input_tier, year): "ZZGen" : "ZZGenSelector", "WGen" : "WGenSelector", "ZGen" : "ZGenSelector", - "TTT" : "TTTSelector", "ThreeLep" : "ThreeLepSelector", "Eff" : "Efficiency", "Efficiency" : "Efficiency", + "FR" : "FakeRateSelector", } self.subanalysis = None @@ -95,7 +95,9 @@ def setInputs(self, inputs): self.addTNamed("ntupleType", self.ntupleType) self.addTNamed("selection", self.selection) self.addTNamed("year", self.year) - + self.addTNamed("logLevel", logging.getLevelName(logging.root.level)) + print(logging.getLevelName(logging.root.level)) + def setSelection(self, selection): self.selection = selection @@ -158,7 +160,7 @@ def clearDatasets(self): def setDatasets(self, datalist): analysis = self.subanalysis if self.subanalysis else self.analysis - datasets = ConfigureJobs.getListOfFiles(datalist, self.input_tier, analysis=analysis) + datasets = ConfigureJobs.getListOfFiles(datalist, analysis, self.input_tier) for dataset in datasets: if "@" in dataset: @@ -189,8 +191,10 @@ def isBackground(self): self.selector_name = self.selector_name.replace("Selector", "BackgroundSelector") def processDataset(self, dataset, file_path, chan): + logging.info("Processing dataset %s" % dataset) select = getattr(ROOT, self.selector_name)() + select.SetInputList(self.inputs) self.addTNamed("name", dataset) if dataset in self.regions: @@ -215,8 +219,8 @@ def processDataset(self, dataset, file_path, chan): sumweights_hist = ROOT.TH1D("sumweights", "sumweights", 1000, 0, 1000) sumweights_hist.SetDirectory(ROOT.gROOT) self.current_file = ROOT.TFile.Open(currfile_name, "update") + self.processLocalFiles(select, file_path, addSumweights, chan) - output_list = select.GetOutputList() processes = [dataset] + (self.regions[dataset] if dataset in self.regions else []) self.writeOutput(output_list, chan, processes, dataset, addSumweights) @@ -296,7 +300,7 @@ def processParallelByDataset(self, datasets, chan): # Store arrays in temp files, since it can get way too big to keep around in memory tempfiles = [self.tempfileName(d) for d in datasets] self.combineParallelFiles(tempfiles, chan) - + # Pool.map can only take in one argument, so expand the array def processDatasetHelper(self, args): self.processDataset(*args) diff --git a/Utilities/python/UserInput.py b/Utilities/python/UserInput.py index bf74fff2..98ad52a8 100644 --- a/Utilities/python/UserInput.py +++ b/Utilities/python/UserInput.py @@ -16,7 +16,7 @@ def getDefaultParser(allow_from_file=True): required=False, help="Name of selection to make, " " as defined in Cuts//.json") parser.add_argument("-a", "--analysis", type=str, - required=False, default="WZxsec2016", + required=False, help="Analysis name, used in selecting the cut json") parser.add_argument("--selectorArgs", nargs='+', type=str, help="List of additional configurations to send to selector") diff --git a/Utilities/scripts/addUnrolledHistsToFile.py b/Utilities/scripts/addUnrolledHistsToFile.py index a192de84..67f09df7 100644 --- a/Utilities/scripts/addUnrolledHistsToFile.py +++ b/Utilities/scripts/addUnrolledHistsToFile.py @@ -23,7 +23,7 @@ ['wzjj-vbfnlo-sf', 'wzjj-vbfnlo-of', ] + \ ['wz3lnu-mg5amcnlo','wz3lnu-powheg', 'zz4l-mg5amcnlo'] + \ ['AllData', 'WZxsec2016data', 'DataEWKCorrected'], - 'Wselection'), + input_tier='Wselection'), ["mjj_etajj_2D_%s" % c for c in ConfigureJobs.getChannels()] + \ ["mjj_etajj_2D_Fakes_%s" % c for c in ConfigureJobs.getChannels()] + \ ["mjj_etajj_2D_%s_%s" % (var, c) for c in ConfigureJobs.getChannels() @@ -39,7 +39,7 @@ ['wzjj-vbfnlo-sf', 'wzjj-vbfnlo-of', ] + \ ['wz3lnu-mg5amcnlo','wz3lnu-powheg', 'zz4l-mg5amcnlo'] + \ ['AllData', 'WZxsec2016data', 'DataEWKCorrected'], - 'Wselection'), + input_tier='Wselection'), ["mjj_dRjj_2D_%s" % c for c in ConfigureJobs.getChannels()] + \ ["mjj_dRjj_2D_Fakes_%s" % c for c in ConfigureJobs.getChannels()] + \ ["mjj_dRjj_2D_%s_%s" % (var, c) for c in ConfigureJobs.getChannels() diff --git a/Utilities/scripts/makeFakeRates.py b/Utilities/scripts/makeFakeRates.py index 54d2b6ee..c152f911 100755 --- a/Utilities/scripts/makeFakeRates.py +++ b/Utilities/scripts/makeFakeRates.py @@ -5,10 +5,12 @@ from python import UserInput,OutputTools from python import ConfigureJobs from python import SelectorTools,HistTools +import sys +import logging ROOT.gROOT.SetBatch(True) -channels = ["eee", "eem", "emm", "mmm"] + def getComLineArgs(): parser = UserInput.getDefaultParser() parser.add_argument("--proof", "-p", @@ -17,6 +19,16 @@ def getComLineArgs(): default=35.87, help="luminosity value (in fb-1)") parser.add_argument("--output_file", "-o", type=str, default="", help="Output file name") + parser.add_argument("-j", "--numCores", type=int, default=1, + help="Number of cores to use (parallelize by dataset)") + parser.add_argument("--input_tier", type=str, + default="", help="Selection stage of input files") + parser.add_argument("--year", type=str, + default="default", help="Year of Analysis") + parser.add_argument("--maxEntries", "-m", type=int, + default=-1, help="Max entries to process") + parser.add_argument("--debug", action='store_true', + help="Print verbose info") return vars(parser.parse_args()) def getHistNames(channels): @@ -26,36 +38,35 @@ def getHistNames(channels): return base_hists return [x+"_"+y for x in base_hists for y in channels] +def addOverflow(hist): + nxbins = hist.GetNbinsX() + nybins = hist.GetNbinsY() + for i in range(1, nybins+1): + setbin = hist.GetBin(nxbins, i) + obin = hist.GetBin(nxbins+1, i) + hist.SetBinContent(setbin, hist.GetBinContent(obin)+hist.GetBinContent(setbin)) + for i in range(1,nxbins): + setbin = hist.GetBin(i, nybins) + obin = hist.GetBin(i, nybins+1) + hist.SetBinContent(setbin, hist.GetBinContent(obin)+hist.GetBinContent(setbin)) + + return hist + # Turn off overflow for FR hists (> 50 is pretty much all EWK anyway) -def makeCompositeHists(name, members, addRatios=True, overflow=False): +def makeCompositeHists(infile, name, members, channels, addRatios=True, overflow=False): composite = ROOT.TList() composite.SetName(name) for directory in [str(i) for i in members.keys()]: - for histname in getHistNames(["eee", "eem", "emm", "mmm"]): - hist = fOut.Get("/".join([directory, histname])) + for histname in getHistNames(channels): + hist = infile.Get("/".join([directory, histname])) if hist: sumhist = composite.FindObject(hist.GetName()) if "data" not in directory and hist.GetEntries() > 0: - sumweights_hist = fOut.Get("/".join([directory, "sumweights"])) + sumweights_hist = infile.Get("/".join([directory, "sumweights"])) sumweights = sumweights_hist.Integral() hist.Scale(members[directory]*1000*args['lumi']/sumweights) if overflow and isinstance(hist, ROOT.TH1): - xbins = hist.GetNbinsX() - ybins = hist.GetNbinsY() - for i in range(1,xbins): - setbin = hist.GetBin(i, ybins) - obin = hist.GetBin(i, ybins+1) - hist.SetBinContent(setbin, - hist.GetBinContent(obin)+hist.GetBinContent(setbin)) - for i in range(1, ybins): - setbin = hist.GetBin(xbins, i) - obin = hist.GetBin(xbins+1, i) - hist.SetBinContent(setbin, - hist.GetBinContent(obin)+hist.GetBinContent(setbin)) - setbin = hist.GetBin(xbins, ybins) - obin = hist.GetBin(xbins+1, ybins+1) - hist.SetBinContent(setbin, - hist.GetBinContent(obin)+hist.GetBinContent(setbin)) + addOverflow(hist) else: raise RuntimeError("hist %s was not produced for " "dataset %s!" % (histname, directory)) @@ -64,15 +75,7 @@ def makeCompositeHists(name, members, addRatios=True, overflow=False): composite.Add(sumhist) else: sumhist.Add(hist) - for hist_name in getHistNames([]): - etot = composite.FindObject(hist_name+"_eee").Clone() - etot.SetName(hist_name+"_allE") - etot.Add(composite.FindObject(hist_name+"_emm")) - composite.Add(etot) - mtot = composite.FindObject(hist_name+"_mmm").Clone() - mtot.SetName(hist_name+"_allMu") - mtot.Add(composite.FindObject(hist_name+"_eem")) - composite.Add(mtot) + if addRatios: ratios = getRatios(composite) for ratio in ratios: @@ -92,16 +95,29 @@ def getRatios(hists): ratios.append(ratio) return ratios + + +runAnalysis = True +channels = ["e", "m"] + + +manager_path = ConfigureJobs.getManagerPath() +if manager_path not in sys.path: + sys.path.insert(0, "/".join([manager_path, + "AnalysisDatasetManager", "Utilities/python"])) + args = getComLineArgs() +logging.basicConfig(level=(logging.DEBUG if args['debug'] else logging.INFO)) proof = 0 if args['proof']: ROOT.TProof.Open("workers=12") proof = ROOT.gProof today = datetime.date.today().strftime("%d%b%Y") -fileName = "data/fakeRate%s-%s.root" % (today, args['selection']) if args['output_file'] == "" \ - else args['output_file'] -fOut = ROOT.TFile(fileName, "recreate") +fileName = args['output_file'] +if not fileName: + fileName = "data/fakeRate%s-%s.root" % (today, args['selection']) +# scale factors fScales = ROOT.TFile('data/scaleFactors.root') muonIsoSF = fScales.Get('muonIsoSF') muonIdSF = fScales.Get('muonTightIdSF') @@ -110,18 +126,43 @@ def getRatios(hists): pileupSF = fScales.Get('pileupSF') sf_inputs = [electronTightIdSF, electronGsfSF, muonIsoSF, muonIdSF, pileupSF] -SelectorTools.applySelector(args["filenames"], - "FakeRateSelector", args['selection'], fOut, - extra_inputs=sf_inputs, proof=args['proof'], - addSumweights=True) - -alldata = makeCompositeHists("AllData", ConfigureJobs.getListOfFilesWithXSec(["WZxsec2016data"])) -OutputTools.writeOutputListItem(alldata, fOut) -allewk = makeCompositeHists("AllEWK", ConfigureJobs.getListOfFilesWithXSec( - ConfigureJobs.getListOfEWKFilenames()), False) -OutputTools.writeOutputListItem(allewk, fOut) -allnonprompt = makeCompositeHists("NonpromptMC", ConfigureJobs.getListOfFilesWithXSec( - ConfigureJobs.getListOfNonpromptFilenames())) -OutputTools.writeOutputListItem(allnonprompt, fOut) -final = HistTools.getDifference(fOut, "DataEWKCorrected", "AllData", "AllEWK", getRatios) -OutputTools.writeOutputListItem(final, fOut) +if runAnalysis: + # Setup and run actual analysis + selector = SelectorTools.SelectorDriver(args['analysis'], args['selection'], args['input_tier'], args['year']) + selector.setNumCores(args['numCores']) + selector.setOutputfile(fileName) + selector.setInputs(sf_inputs) + selector.setMaxEntries(args['maxEntries']) + selector.setNtupeType("NanoAOD") + if args['filenames']: + selector.setDatasets(args['filenames']) + else: + selector.setFileList(*args['inputs_from_file']) + + + mc = selector.applySelector() + + selector.outputFile().Close() + #OutputTools.addMetaInfo(ROOT.TFile(fileName, "UPDATE")) + + +# Create end files +fOut = ROOT.TFile(fileName, "UPDATE") + +combinedHist = makeCompositeHists(fOut, "AllData", ConfigureJobs.getListOfFilesWithXSec(args["filenames"], args['analysis']), channels) +for i in combinedHist: + i.Print() +OutputTools.writeOutputListItem(combinedHist, fOut) +combinedHist.Clear() +# alldata = makeCompositeHists("AllData", ConfigureJobs.getListOfFilesWithXSec(["WZxsec2016data"]), channels) +# OutputTools.writeOutputListItem(alldata, fOut) +# allewk = makeCompositeHists("AllEWK", ConfigureJobs.getListOfFilesWithXSec( +# ConfigureJobs.getListOfEWKFilenames()), channels, False) +# OutputTools.writeOutputListItem(allewk, fOut) +# allnonprompt = makeCompositeHists("NonpromptMC", ConfigureJobs.getListOfFilesWithXSec( +# ConfigureJobs.getListOfNonpromptFilenames()), channels) +# OutputTools.writeOutputListItem(allnonprompt, fOut) +# final = HistTools.getDifference(fOut, "DataEWKCorrected", "AllData", "AllEWK", getRatios) +# OutputTools.writeOutputListItem(final, fOut) + + diff --git a/Utilities/scripts/prepareCombine.py b/Utilities/scripts/prepareCombine.py index e24a215b..5fbaf07d 100755 --- a/Utilities/scripts/prepareCombine.py +++ b/Utilities/scripts/prepareCombine.py @@ -212,7 +212,7 @@ def addInterference(fOut, variable, addControlRegion): if variable == "MTWZ": rebin = array.array('d', ConfigureJobs.getBinning(isVBS=isVBS, isHiggs=args['higgs'])) alldata = HistTools.makeCompositeHists(fIn, "AllData", - ConfigureJobs.getListOfFilesWithXSec(["WZxsec2016data"], manager_path), args['lumi'], + ConfigureJobs.getListOfFilesWithXSec(["WZxsec2016data"]), args['lumi'], [variable +"_"+ c for c in chans], rebin=rebin) for chan in chans: hist = alldata.FindObject(variable+"_"+chan) @@ -274,7 +274,7 @@ def addInterference(fOut, variable, addControlRegion): plots += ["backgroundControlYield_lheWeights_" + c for c in chans] group = HistTools.makeCompositeHists(fIn, plot_group, ConfigureJobs.getListOfFilesWithXSec( - config_factory.getPlotGroupMembers(plot_group), manager_path), args['lumi'], plots, rebin=rebin) + config_factory.getPlotGroupMembers(plot_group)), args['lumi'], plots, rebin=rebin) if scaleWZ and plot_group in wz_scalefacs.keys(): for h in group: if h.InheritsFrom("TH1"): diff --git a/Utilities/scripts/printEventDetails.py b/Utilities/scripts/printEventDetails.py index 711c626e..73afe754 100644 --- a/Utilities/scripts/printEventDetails.py +++ b/Utilities/scripts/printEventDetails.py @@ -41,7 +41,7 @@ def getEventSelectionExpr(path, comparison, channel): isfile = any(os.path.isfile(name) or os.path.exists(name.rstrip("/*")) for name in args.filelist) -filelist = ConfigureJobs.getListOfFiles(args.filelist, args.selection) if \ +filelist = ConfigureJobs.getListOfFiles(args.filelist, input_tier=args.selection) if \ not isfile else args.filelist print filelist states = [x.strip() for x in args.channels.split(",")] diff --git a/fillHists.sh b/fillHists.sh index 6f6922ae..154dff91 100755 --- a/fillHists.sh +++ b/fillHists.sh @@ -17,12 +17,21 @@ if [ "$modTime" -gt "$runTime" ]; then fi EXE=./Utilities/scripts/makeHistFile.py +# EXE=./Utilities/scripts/makeFakeRates.py ANA=ThreeLep YEAR=2016 -INP=TwoLep_Met40 -# SEL=FourTopCutBasedEl -SEL=FourTopMVAEl +# INP=TwoLep_Met40 +INP=NoSkim_MultYear +# INP=TwoLep_Met25 + +# # SEL=FourTopCutBasedEl +# SEL=FourTopMVAEl +SEL=MVAStudy +# ANA=FR +# INP=FakeRateEOS +# YEAR=2016 +# SEL=FakeRate # EXE=./Utilities/scripts/makeHistFile.py # ANA=Efficiency:ThreeLep @@ -33,9 +42,6 @@ SEL=FourTopMVAEl # SEL=test - - - $EXE -a $ANA --input_tier $INP -s $SEL \ - --year $YEAR --test ${@} + --year $YEAR ${@} diff --git a/interface/BranchManager.h b/interface/BranchManager.h index 0cc22825..b99e4c80 100644 --- a/interface/BranchManager.h +++ b/interface/BranchManager.h @@ -35,11 +35,26 @@ struct BranchManager { for(auto& it: branchHolder) { it->GetEntry(entry); } - } - void SetSpecificEntry(int entry, std::string name) { - specificBranch[name]->GetEntry(entry); - } + template + void SetSpecificBranch(std::string name, T& holder) { + specificBranch[name] = {}; + fChain->SetBranchAddress(name.c_str(), &holder, &specificBranch[name]); + } + + void SetEntry(int entry) { + for(auto& it: branchHolder) { + it->GetEntry(entry); + } + } + + bool branchExists(std::string name) { + return bool(fChain->FindBranch(name.c_str())); + } + + void SetSpecificEntry(int entry, std::string name) { + specificBranch[name]->GetEntry(entry); + } void CleanUp() { branchHolder.clear(); diff --git a/interface/FakeRateSelector.h b/interface/FakeRateSelector.h index d25cadf9..28bf90d4 100644 --- a/interface/FakeRateSelector.h +++ b/interface/FakeRateSelector.h @@ -6,38 +6,62 @@ #include #include #include +#include "Math/GenVector/VectorUtil.h" #include #include +#include +#include // Headers needed by this particular selector -#include "Analysis/VVAnalysis/interface/WZSelectorBase.h" -#include +#include "Analysis/VVAnalysis/interface/SelectorBase.h" +#include "Analysis/VVAnalysis/interface/ThreeLepSelector.h" +#include "Analysis/VVAnalysis/interface/ScaleFactor.h" +#include "Analysis/VVAnalysis/interface/BranchManager.h" +#include "Analysis/VVAnalysis/interface/GoodParticle.h" -class FakeRateSelector : public WZSelectorBase { +class FakeRateSelector : public SelectorBase { public : - TH2D* passingTight2D_; - TH1D* passingTight1DPt_; - TH1D* passingTight1DEta_; - TH2D* passingLoose2D_; - TH1D* passingLoose1DPt_; - TH1D* passingLoose1DEta_; - TH2D* ratio2D_; - TH1D* ratio1DPt_; - TH1D* ratio1DEta_; - - Float_t type1_pfMETEt; - UInt_t nCBVIDVetoElec; - UInt_t nWZLooseMuon; + int sel_MVAStudy, sel_FakeRate; + + + - TBranch* b_type1_pfMETEt; - TBranch* b_nCBVIDVetoElec; - TBranch* b_nWZLooseMuon; + std::unordered_map passingTight2D_map; + std::unordered_map passingTight1DPt_map; + std::unordered_map passingTight1DEta_map; + std::unordered_map passingLoose2D_map; + std::unordered_map passingLoose1DPt_map; + std::unordered_map passingLoose1DEta_map; + + TTreeReader fReader; + TTreeReaderValue HLT_Mu8 = {fReader, "HLT_Mu8"}; + TTreeReaderValue HLT_Ele8_CaloIdM_TrackIdM_PFJet30 = {fReader, "HLT_Ele8_CaloIdM_TrackIdM_PFJet30"}; + + double weight = 1.0; + + const float FR_MAX_PT_ = 50; + const float FR_MAX_ETA_ = 2.5; + // Readers to access the data (delete the ones you do not need). virtual void SetupNewDirectory() override; virtual void FillHistograms(Long64_t entry, std::pair variation) override; + void clearValues(); + void fillReco(std::vector& genList, const std::vector& recoList); + + // overloaded or necessary functions + virtual void SetBranchesNanoAOD() override; + void LoadBranchesNanoAOD(Long64_t entry, std::pair variation) override; + // Readers to access the data (delete the ones you do not need). + virtual void SlaveBegin(TTree *tree) override {return;} + virtual void Init(TTree *tree) override; //{return;} + ///ignore + void LoadBranchesUWVV(Long64_t entry, std::pair variation) override {return;} + virtual void SetBranchesUWVV() override {return;} + + ThreeLepSelector Analyzer; ClassDefOverride(FakeRateSelector,0); }; diff --git a/interface/HelperFunctions.h b/interface/HelperFunctions.h index e676c090..58a74bae 100644 --- a/interface/HelperFunctions.h +++ b/interface/HelperFunctions.h @@ -4,13 +4,15 @@ #include #include #include - +#include #include "Analysis/VVAnalysis/interface/GoodParticle.h" +/* #include */ double JetSphericity(std::vector& jets) { TMatrixDSym sphere(3); double pSum = 0; for(auto jet: jets) { + if(!jet.passedJetSel()) continue; for(int i = 1; i <= 3; i++) { sphere(i-1,i-1) += pow(jet[i], 2); pSum += pow(jet[i], 2); @@ -29,10 +31,48 @@ double JetSphericity(std::vector& jets) { double JetCentrality(std::vector& jets, double HT) { double eTot = 0; for(auto jet: jets) { + if(!jet.passedJetSel()) continue; eTot += jet.E(); } return HT/eTot; } +std::pair EventShape(std::vector& jets, std::vector& leps, double Met2, double MetPhi) { + TMatrixDSym shape(3); + double pSum = 0; + for(auto jet: jets) { + for(int i = 1; i <= 3; i++) { + shape(i-1,i-1) += pow(jet[i], 2); + pSum += pow(jet[i], 2); + for(int j = i+1; j <= 3; j++) { + shape(i-1,j-1) += jet[i]*jet[j]; + shape(j-1,i-1) += jet[i]*jet[j]; + } + } + } + for(auto lep: leps) { + for(int i = 1; i <= 3; i++) { + shape(i-1,i-1) += pow(lep[i], 2); + pSum += pow(lep[i], 2); + for(int j = i+1; j <= 3; j++) { + shape(i-1,j-1) += lep[i]*lep[j]; + shape(j-1,i-1) += lep[i]*lep[j]; + } + } + } + pSum += Met2; + shape(1, 1 ) += Met2*pow(cos(MetPhi), 2); + shape(2, 2 ) += Met2*pow(sin(MetPhi), 2); + shape(1, 2 ) += Met2*sin(MetPhi/2)/2; + shape(2, 1 ) += Met2*sin(MetPhi/2)/2; + shape *= 1/pSum; + TMatrixDSymEigen eigenGetter(shape); + TVectorD eigen = eigenGetter.GetEigenValues(); + return std::make_pair(3.*(eigen(1)*eigen(2) + eigen(1)*eigen(0) + eigen(0)*eigen(2)), 27*eigen(0)*eigen(1)*eigen(2)); + +} + +#include "Math/GenVector/VectorUtil.h" + #endif /* HELPERFUNCTIONS_H */ diff --git a/interface/SelectorBase.h b/interface/SelectorBase.h index b14c2942..32ed22e6 100644 --- a/interface/SelectorBase.h +++ b/interface/SelectorBase.h @@ -17,16 +17,18 @@ #include #include #include "Analysis/VVAnalysis/interface/ScaleFactor.h" +#include "Analysis/VVAnalysis/interface/dynEnum.h" +#include "Fireworks/Core/interface/fwLog.h" -#define PAIR(NAME_) {#NAME_, NAME_} enum Channel { - e, m, - ep, en, mp, mn, - ee, em, mm, - eee, eem, emm, mmm, - eeee, eemm, mmee, mmmm, - Inclusive, Unknown, lll, all, + e, m, + ep, en, mp, mn, + ee, em, mm, + eee, eem, emm, mmm, + eeee, eemm, mmee, mmmm, + Inclusive, Unknown, lll, all, + SS, OS, mult, one, }; enum Systematic { @@ -108,96 +110,30 @@ class SelectorBase : public TSelector { /* | |___| |\ | |_| | | | \__ \ */ /* |_____|_| \_|\___/|_| |_|___/ */ /*********************************/ - - enum NtupleType { - NanoAOD, - UWVV, - Bacon, - }; - - enum Selection { - Default, None, - tightleptons, ZZGenFiducial, - Wselection, Zselection, - WselectionAntiIso, - Wselection_Full, FakeRateSelectionLoose, - FakeRateSelectionTight, VBSselection_Loose, - VBSselection_NoZeppenfeld, VBSselection_Tight, - VBSselection_Loose_Full, VBSselection_NoZeppenfeld_Full, - VBSselection_Tight_Full, VBSBackgroundControl, - VBSBackgroundControlATLAS, VBSBackgroundControl_Full, - VBSBackgroundControlLoose, VBSBackgroundControlLoose_Full, - Inclusive2Jet, Inclusive2Jet_Full, - TightWithLooseVeto, FourTopPlots, - FourTopCutBasedEl, FourTopMVAEl, - BEfficiency, test, - }; - enum Year { - yrdefault, yr2016, yr2017, yr2018 + enum NtupleType { + NanoAOD, + UWVV, + Bacon, }; typedef std::pair ChannelPair; - + typedef std::unordered_map HistMap1D; typedef std::unordered_map HistMap2D; typedef std::unordered_map HistMap3D; typedef std::pair SystPair; typedef std::map SystMap; - /****************************/ - /* __ __ */ - /* | \/ | __ _ _ __ ___ */ - /* | |\/| |/ _` | '_ \/ __| */ - /* | | | | (_| | |_) \__ \ */ - /* |_| |_|\__,_| .__/|___/ */ - /* |_| */ - /****************************/ - - std::map selectionMap_ = { - {"Default", Default}, - {"None", None}, - {"tightleptons", tightleptons}, - {"ZZGenFiducial", ZZGenFiducial}, - {"Wselection", Wselection}, - {"WselectionAntiIso", WselectionAntiIso}, - {"Zselection", Zselection}, - {"Wselection_Full", Wselection_Full}, - {"FakeRateSelectionLoose", FakeRateSelectionLoose}, - {"FakeRateSelectionTight", FakeRateSelectionTight}, - {"VBSselection_Loose", VBSselection_Loose}, - {"VBSselection_NoZeppenfeld", VBSselection_NoZeppenfeld}, - {"VBSselection_Tight", VBSselection_Tight}, - {"VBSselection_Loose_Full", VBSselection_Loose_Full}, - {"VBSselection_NoZeppenfeld_Full", VBSselection_NoZeppenfeld_Full}, - {"VBSselection_Tight_Full", VBSselection_Tight_Full}, - {"VBSBackgroundControl", VBSBackgroundControl}, - {"VBSBackgroundControlATLAS", VBSBackgroundControlATLAS}, - {"VBSBackgroundControl_Full", VBSBackgroundControl_Full}, - {"VBSBackgroundControlLoose", VBSBackgroundControlLoose}, - {"VBSBackgroundControlLoose_Full", VBSBackgroundControlLoose_Full}, - {"Inclusive2Jet", Inclusive2Jet}, - {"Inclusive2Jet_Full", Inclusive2Jet_Full}, - {"TightWithLooseVeto", TightWithLooseVeto}, - {"FourTopPlots", FourTopPlots}, - {"FourTopCutBasedEl", FourTopCutBasedEl}, - {"FourTopMVAEl", FourTopMVAEl}, - }; - std::map yearMap_ = { - {"default", yrdefault}, - {"2016", yr2016}, - {"2017", yr2017}, - {"2018", yr2018}, - }; - std::map channelMap_ = { - {"e", e}, {"m", m}, - {"ep", ep}, {"mp", mp}, {"ep", ep}, {"mn", mn}, - {"ee", ee}, {"em", em}, {"mm", mm}, - {"eee", eee}, {"eem", eem}, {"emm", emm}, {"mmm", mmm}, - {"eeee", eeee}, {"eemm", eemm}, {"mmee", mmee}, {"mmmm", mmmm}, - {"Inclusive", Inclusive}, {"lll", lll}, {"all", all}, + {"e", e}, {"m", m}, + {"ep", ep}, {"mp", mp}, {"ep", ep}, {"mn", mn}, + {"ee", ee}, {"em", em}, {"mm", mm}, + {"eee", eee}, {"eem", eem}, {"emm", emm}, {"mmm", mmm}, + {"eeee", eeee}, {"eemm", eemm}, {"mmee", mmee}, {"mmmm", mmmm}, + {"Inclusive", Inclusive}, {"lll", lll}, {"all", all}, + {"SS", SS}, {"OS", OS}, {"mult", mult}, {"one", one}, }; @@ -301,12 +237,20 @@ class SelectorBase : public TSelector { std::string channelName_ = "Unnamed"; Channel channel_ = Unknown; NtupleType ntupleType_ = NanoAOD; - std::string selectionName_ = "tightleptons"; - Selection selection_ = tightleptons; - Year year_ = yrdefault; - bool isMC_; + // Enums + DynEnum enumFactory; + std::unordered_map selectionMap_; + std::unordered_map addSelection_; + std::map yearMap_; + int selection_; + int year_; + + // default enum values + int yrdefault; + bool isMC_; + float GetPrefiringEfficiencyWeight(std::vector* jetPt, std::vector* jetEta); virtual std::string GetNameFromFile() { return ""; } void InitializeHistogramFromConfig(std::string name, ChannelPair channel, std::vector& histData); diff --git a/interface/TTTSelector.h b/interface/TTTSelector.h deleted file mode 100644 index 7dc41340..00000000 --- a/interface/TTTSelector.h +++ /dev/null @@ -1,149 +0,0 @@ -#ifndef TTTSelector_h -#define TTTSelector_h - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -// Headers needed by this particular selector -#include -#include "Analysis/VVAnalysis/interface/ScaleFactor.h" -#include "Analysis/VVAnalysis/interface/SelectorBase.h" -#include "Analysis/VVAnalysis/interface/BranchManager.h" -#include "Analysis/VVAnalysis/interface/GoodParticle.h" - -typedef ROOT::Math::LorentzVector> LorentzVector; - -class TTTSelector : public SelectorBase { - public : - /*****************************************/ - /* ____ ____ ___ __ __ ___ __ __ */ - /* || )) || \\ // \\ ||\ || // || || */ - /* ||=) ||_// ||=|| ||\\|| (( ||==|| */ - /* ||_)) || \\ || || || \|| \\__ || || */ - /*****************************************/ - - ScaleFactor* pileupSF_; - ScaleFactor* muonSF_; - ScaleFactor* eIdSF_ ; - ScaleFactor* eGsfSF_; - ScaleFactor* mIdSF_; - ScaleFactor* mIsoSF_; - - // Common variables - Float_t genWeight; - Float_t MET; - Float_t type1_pfMETPhi; - - //NanoAOD variables - static const unsigned int N_KEEP_MU_E_ = 15; - static const unsigned int N_KEEP_JET_ = 35; - - UInt_t nElectron; - Float_t Electron_pt[N_KEEP_MU_E_]; - Float_t Electron_eta[N_KEEP_MU_E_]; - Float_t Electron_phi[N_KEEP_MU_E_]; - Float_t Electron_mass[N_KEEP_MU_E_]; - Int_t Electron_cutBased[N_KEEP_MU_E_]; - Int_t Electron_charge[N_KEEP_MU_E_]; - Float_t Electron_MVA[N_KEEP_MU_E_]; - Float_t Electron_miniPFRelIso_all[N_KEEP_MU_E_]; - Float_t Electron_dxy[N_KEEP_MU_E_]; - Float_t Electron_dz[N_KEEP_MU_E_]; - Float_t Electron_sip3d[N_KEEP_MU_E_]; - - - UInt_t nMuon; - Float_t Muon_pt[N_KEEP_MU_E_]; - Float_t Muon_eta[N_KEEP_MU_E_]; - Float_t Muon_phi[N_KEEP_MU_E_]; - Float_t Muon_mass[N_KEEP_MU_E_]; - Int_t Muon_charge[N_KEEP_MU_E_]; - Bool_t Muon_tightId[N_KEEP_MU_E_]; - Bool_t Muon_mediumId[N_KEEP_MU_E_]; - UChar_t Muon_tkIsoId[N_KEEP_MU_E_]; - Float_t Muon_pfRelIso04_all[N_KEEP_MU_E_]; - Float_t Muon_miniPFRelIso_all[N_KEEP_MU_E_]; - Float_t Muon_dxy[N_KEEP_MU_E_]; - Float_t Muon_dz[N_KEEP_MU_E_]; - Float_t Muon_sip3d[N_KEEP_MU_E_]; - - Int_t numPU; - - UInt_t nJet; - Float_t Jet_btagCSVV2[N_KEEP_JET_]; - Float_t Jet_btagDeepB[N_KEEP_JET_]; - Float_t Jet_eta[N_KEEP_JET_]; - Float_t Jet_phi[N_KEEP_JET_]; - Float_t Jet_pt[N_KEEP_JET_]; - Float_t Jet_mass[N_KEEP_JET_]; - Float_t Jet_neHEF[N_KEEP_JET_]; - Float_t Jet_neEmEF[N_KEEP_JET_]; - Int_t Jet_nConstituents[N_KEEP_JET_]; - Float_t Jet_chHEF[N_KEEP_JET_]; - Float_t Jet_chEmEF[N_KEEP_JET_]; - - ClassDefOverride(TTTSelector,0); - - /*******************************************************/ - /* __ __ ___ ____ __ ___ ____ __ ____ __ */ - /* || || // \\ || \\ || // \\ || )) || || (( \ */ - /* \\ // ||=|| ||_// || ||=|| ||=) || ||== \\ */ - /* \V/ || || || \\ || || || ||_)) ||__| ||___ \_)) */ - /*******************************************************/ - - Float_t weight; - BranchManager b; - std::vector goodParts; - std::vector charge; - double HT; - int nTightJet, nBJets; - - /************************************************************/ - /* _____ __ __ __ __ ___ ______ __ ___ __ __ __ */ - /* || || || ||\ || // | || | || // \\ ||\ || (( \ */ - /* ||== || || ||\\|| (( || || (( )) ||\\|| \\ */ - /* || \\_// || \|| \\__ || || \\_// || \|| \_)) */ - /************************************************************/ - - /// Captial I for particle definition - bool IsGoodMuon(size_t); - bool IsGoodCBElectron(size_t); - bool IsGoodJet(size_t); - bool IsGoodBJet(size_t); - bool IsGoodMVAElectron2016(size_t); - bool IsGoodMVAElectron2017(size_t); - - // lowercase I for helper function (kinda shitty convention, may change) - bool isLooseJetId(size_t); - bool isTightJetId(size_t); - bool isOverlap(size_t); - bool passFullIso(LorentzVector&, int, int); - - //// General Functions - int getSRBin(); - void clearValues(); - void ApplyScaleFactors(); - - // Overloaded or necesary functions - virtual void SetBranchesNanoAOD() override; - void LoadBranchesNanoAOD(Long64_t entry, std::pair variation) override; - void FillHistograms(Long64_t entry, std::pair variation) override; - virtual void SetupNewDirectory() override; - // Readers to access the data (delete the ones you do not need). - virtual void SlaveBegin(TTree *tree) override; - virtual void Init(TTree *tree) override; - - ///ignore - void LoadBranchesUWVV(Long64_t entry, std::pair variation) override; - virtual void SetBranchesUWVV() override; -}; - -#endif - diff --git a/interface/ThreeLepSelector.h b/interface/ThreeLepSelector.h index 3d8850dc..c6a07ab8 100644 --- a/interface/ThreeLepSelector.h +++ b/interface/ThreeLepSelector.h @@ -7,12 +7,16 @@ #include #include #include +#include #include #include #include #include #include #include +#include +#include +#include // Headers needed by this particular selector #include "Analysis/VVAnalysis/interface/ScaleFactor.h" @@ -28,7 +32,7 @@ enum PID {PID_MUON = 13, PID_ELECTRON = 11, PID_BJET = 5, PID_CJET = 4, PID_JET} class ThreeLepSelector : public SelectorBase { public : - #include "Analysis/VVAnalysis/interface/FourTopScales.h" +#include "Analysis/VVAnalysis/interface/FourTopScales.h" /*****************************************/ /* ____ ____ ___ __ __ ___ __ __ */ @@ -49,86 +53,73 @@ public : Float_t MET; Float_t type1_pfMETPhi; - //NanoAOD variables - static const unsigned int N_KEEP_MU_E_ = 15; - static const unsigned int N_KEEP_JET_ = 35; - static const unsigned int N_KEEP_GEN_ = 300; + int MVAStudy, FourTopMVAEl, FourTopCutBasedEl, FakeRate; + int yr2016, yr2017, yr2018; + TTreeReader fReader, tmpReader; - UInt_t nElectron; - Float_t Electron_pt[N_KEEP_MU_E_]; - Float_t Electron_eta[N_KEEP_MU_E_]; - Float_t Electron_phi[N_KEEP_MU_E_]; - Float_t Electron_mass[N_KEEP_MU_E_]; - Int_t Electron_cutBased[N_KEEP_MU_E_]; - Int_t Electron_charge[N_KEEP_MU_E_]; - Float_t Electron_MVA[N_KEEP_MU_E_]; - Float_t Electron_miniPFRelIso_all[N_KEEP_MU_E_]; - Float_t Electron_dxy[N_KEEP_MU_E_]; - Float_t Electron_dz[N_KEEP_MU_E_]; - Float_t Electron_sip3d[N_KEEP_MU_E_]; - Bool_t Electron_convVeto[N_KEEP_MU_E_]; - UChar_t Electron_lostHits[N_KEEP_MU_E_]; - Int_t Electron_tightCharge[N_KEEP_MU_E_]; - Float_t Electron_sieie[N_KEEP_MU_E_]; - Float_t Electron_hoe[N_KEEP_MU_E_]; - Float_t Electron_deltaEtaSC[N_KEEP_MU_E_]; - Float_t Electron_eInvMinusPInv[N_KEEP_MU_E_]; - Float_t Electron_dr03EcalRecHitSumEt[N_KEEP_MU_E_]; - Float_t Electron_dr03HcalDepth1TowerSumEt[N_KEEP_MU_E_]; - Float_t Electron_dr03TkSumPt[N_KEEP_MU_E_]; - Int_t Electron_vidBitmap[N_KEEP_MU_E_]; - Float_t Electron_jetRelIso[N_KEEP_MU_E_]; - Float_t Electron_eCorr[N_KEEP_MU_E_]; - - UInt_t nMuon; - Float_t Muon_pt[N_KEEP_MU_E_]; - Float_t Muon_eta[N_KEEP_MU_E_]; - Float_t Muon_phi[N_KEEP_MU_E_]; - Float_t Muon_mass[N_KEEP_MU_E_]; - Int_t Muon_charge[N_KEEP_MU_E_]; - Bool_t Muon_tightId[N_KEEP_MU_E_]; - Bool_t Muon_mediumId[N_KEEP_MU_E_]; - UChar_t Muon_tkIsoId[N_KEEP_MU_E_]; - Float_t Muon_pfRelIso04_all[N_KEEP_MU_E_]; - Float_t Muon_miniPFRelIso_all[N_KEEP_MU_E_]; - Float_t Muon_dxy[N_KEEP_MU_E_]; - Float_t Muon_dz[N_KEEP_MU_E_]; - Float_t Muon_sip3d[N_KEEP_MU_E_]; - Bool_t Muon_isGlobal[N_KEEP_MU_E_]; - Bool_t Muon_isTracker[N_KEEP_MU_E_]; - Bool_t Muon_isPFcand[N_KEEP_MU_E_]; - Int_t Muon_tightCharge[N_KEEP_MU_E_]; - Float_t Muon_jetPtRelv2[N_KEEP_MU_E_]; - Float_t Muon_jetRelIso[N_KEEP_MU_E_]; + + TTreeReaderArray Electron_pt = {fReader, "Electron_pt"}; + TTreeReaderArray Electron_eta = {fReader, "Electron_eta"}; + TTreeReaderArray Electron_phi = {fReader, "Electron_phi"}; + TTreeReaderArray Electron_mass = {fReader, "Electron_mass"}; + TTreeReaderArray Electron_cutBased = {fReader, "Electron_cutBased"}; + TTreeReaderArray Electron_charge = {fReader, "Electron_charge"}; + TTreeReaderArray Electron_MVA = {fReader, "Electron_MVA"}; + TTreeReaderArray Electron_miniPFRelIso_all = {fReader, "Electron_miniPFRelIso_all"}; + TTreeReaderArray Electron_dxy = {fReader, "Electron_dxy"}; + TTreeReaderArray Electron_dz = {fReader, "Electron_dz"}; + TTreeReaderArray Electron_sip3d = {fReader, "Electron_sip3d"}; + TTreeReaderArray Electron_convVeto = {fReader, "Electron_convVeto"}; + TTreeReaderArray Electron_lostHits = {fReader, "Electron_lostHits"}; + TTreeReaderArray Electron_tightCharge = {fReader, "Electron_tightCharge"}; + TTreeReaderArray Electron_sieie = {fReader, "Electron_sieie"}; + TTreeReaderArray Electron_hoe = {fReader, "Electron_hoe"}; + TTreeReaderArray Electron_deltaEtaSC = {fReader, "Electron_deltaEtaSC"}; + TTreeReaderArray Electron_eInvMinusPInv = {fReader, "Electron_eInvMinusPInv"}; + TTreeReaderArray Electron_dr03EcalRecHitSumEt = {fReader, "Electron_dr03EcalRecHitSumEt"}; + TTreeReaderArray Electron_dr03HcalDepth1TowerSumEt = {fReader, "Electron_dr03HcalDepth1TowerSumEt"}; + TTreeReaderArray Electron_dr03TkSumPt = {fReader, "Electron_dr03TkSumPt"}; + //TTreeReaderArray Electron_vidBitmap = {fReader, "Electron_vidBitmap"}; + TTreeReaderArray Electron_jetRelIso = {fReader, "Electron_jetRelIso"}; + TTreeReaderArray Electron_eCorr = {fReader, "Electron_eCorr"}; + + TTreeReaderArray Muon_pt = {fReader, "Muon_pt"}; + TTreeReaderArray Muon_eta = {fReader, "Muon_eta"}; + TTreeReaderArray Muon_phi = {fReader, "Muon_phi"}; + TTreeReaderArray Muon_mass = {fReader, "Muon_mass"}; + TTreeReaderArray Muon_charge = {fReader, "Muon_charge"}; + TTreeReaderArray Muon_tightId = {fReader, "Muon_tightId"}; + TTreeReaderArray Muon_mediumId = {fReader, "Muon_mediumId"}; + TTreeReaderArray Muon_tkIsoId = {fReader, "Muon_tkIsoId"}; + TTreeReaderArray Muon_pfRelIso04_all = {fReader, "Muon_pfRelIso04_all"}; + TTreeReaderArray Muon_miniPFRelIso_all = {fReader, "Muon_miniPFRelIso_all"}; + TTreeReaderArray Muon_dxy = {fReader, "Muon_dxy"}; + TTreeReaderArray Muon_dz = {fReader, "Muon_dz"}; + TTreeReaderArray Muon_sip3d = {fReader, "Muon_sip3d"}; + TTreeReaderArray Muon_isGlobal = {fReader, "Muon_isGlobal"}; + TTreeReaderArray Muon_isTracker = {fReader, "Muon_isTracker"}; + TTreeReaderArray Muon_isPFcand = {fReader, "Muon_isPFcand"}; + TTreeReaderArray Muon_tightCharge = {fReader, "Muon_tightCharge"}; + + TTreeReaderArray Jet_btagCSVV2 = {fReader, "Jet_btagCSVV2"}; + TTreeReaderArray Jet_btagDeepB = {fReader, "Jet_btagDeepB"}; + TTreeReaderArray Jet_eta = {fReader, "Jet_eta"}; + TTreeReaderArray Jet_phi = {fReader, "Jet_phi"}; + TTreeReaderArray Jet_pt = {fReader, "Jet_pt"}; + TTreeReaderArray Jet_mass = {fReader, "Jet_mass"}; + TTreeReaderArray Jet_nConstituents = {fReader, "Jet_nConstituents"}; + TTreeReaderArray Jet_jetId = {fReader, "Jet_jetId"}; + TTreeReaderArray Jet_hadronFlavour = {fReader, "Jet_hadronFlavour"}; + TTreeReaderArray Jet_rawFactor = {fReader, "Jet_rawFactor"}; - Int_t numPU; Float_t Pileup_nTrueInt; - + UInt_t nElectron; + UInt_t nMuon; UInt_t nJet; - Float_t Jet_btagCSVV2[N_KEEP_JET_]; - Float_t Jet_btagDeepB[N_KEEP_JET_]; - Float_t Jet_eta[N_KEEP_JET_]; - Float_t Jet_phi[N_KEEP_JET_]; - Float_t Jet_pt[N_KEEP_JET_]; - Float_t Jet_mass[N_KEEP_JET_]; - Float_t Jet_neHEF[N_KEEP_JET_]; - Float_t Jet_neEmEF[N_KEEP_JET_]; - Int_t Jet_nConstituents[N_KEEP_JET_]; - Float_t Jet_chHEF[N_KEEP_JET_]; - Float_t Jet_chEmEF[N_KEEP_JET_]; - Int_t Jet_jetId[N_KEEP_JET_]; - Int_t Jet_hadronFlavour[N_KEEP_JET_]; - Float_t Jet_rawFactor[N_KEEP_JET_]; - Float_t Jet_L1[N_KEEP_JET_]; - Float_t Jet_L2L3[N_KEEP_JET_]; - - Int_t GenPart_pdgId[N_KEEP_GEN_]; - Int_t GenPart_genPartIdxMother[N_KEEP_GEN_]; - Int_t GenPart_status[N_KEEP_GEN_]; UInt_t nGenPart; - + Bool_t Flag_goodVertices; Bool_t Flag_globalSuperTightHalo2016Filter; Bool_t Flag_HBHENoiseFilter; @@ -137,6 +128,19 @@ public : Bool_t Flag_BadPFMuonFilter; Bool_t Flag_ecalBadCalibFilter; + Bool_t HLT_DoubleMu8_Mass8_PFHT300; + Bool_t HLT_Mu8_Ele8_CaloIdM_TrackIdM_Mass8_PFHT300; + Bool_t HLT_DoubleEle8_CaloIdM_TrackIdM_Mass8_PFHT300; + Bool_t HLT_AK8PFJet450; + Bool_t HLT_PFJet450; + + /// Can be excluded + TTreeReaderArray Jet_L1 = {tmpReader, "Jet_L1"}; + TTreeReaderArray Jet_L2L3 = {tmpReader, "Jet_L2L3"}; + TTreeReaderArray GenPart_pdgId = {tmpReader, "GenPart_pdgId"}; + TTreeReaderArray GenPart_genPartIdxMother = {tmpReader, "GenPart_genPartIdxMother"}; + TTreeReaderArray GenPart_status = {tmpReader, "GenPart_status"}; + ClassDefOverride(ThreeLepSelector,0); @@ -152,6 +156,9 @@ public : std::vector goodLeptons; std::vector looseLeptons; std::vector goodJets; + std::vector jetList; + std::vector bjetList; + double HT; int nJets, nBJets; bool passZVeto; @@ -166,20 +173,6 @@ public : TH2D* h_btag_eff_c; TH2D* h_btag_eff_udsg; - - std::map eventVec; - std::map info; - - Bool_t HLT_DoubleMu8_Mass8_PFHT300; - Bool_t HLT_Mu8_Ele8_CaloIdM_TrackIdM_Mass8_PFHT300; - Bool_t HLT_DoubleEle8_CaloIdM_TrackIdM_Mass8_PFHT300; - Bool_t HLT_AK8PFJet450; - Bool_t HLT_PFJet450; - - Float_t PV_x; - Float_t PV_y; - Float_t PV_z; - /************************************************************/ /* _____ __ __ __ __ ___ ______ __ ___ __ __ __ */ /* || || || ||\ || // | || | || // \\ ||\ || (( \ */ @@ -218,6 +211,10 @@ public : bool MetFilter(); float getBtagEffFromFile(double, double, int); double getWDecayScaleFactor(); + std::vector::iterator findJet(std::vector::iterator&, int); + std::map treeMap = {{"tree", nullptr}}; + float bNJets, bnBJets, bHT, bMET, bl1Pt, bl2Pt, blMass, bsphere, bCentral, bShape1, bShape2, bnLeps, bDilepCharge, bnlBJets, bntBJets, bnlLeps; + float bjMass, bjdr, bj1Pt, bj2Pt, bj3Pt, bj4Pt, bj5Pt, bj6Pt, bj7Pt, bj8Pt, bb1Pt, bb2Pt, bb3Pt, bb4Pt; //// General Functions int getSRBin() const; @@ -239,4 +236,3 @@ public : }; #endif - diff --git a/interface/WZSelector.h b/interface/WZSelector.h index be83480f..ab65db92 100644 --- a/interface/WZSelector.h +++ b/interface/WZSelector.h @@ -12,6 +12,8 @@ public : bool isaQGC_ = false; bool doaQGC_ = false; + int Inclusive2Jet_Full, Wselection_Full, Inclusive2Jet, FakeRateSelectionTight, FakeRateSelectionLoose; + std::vector* scaleWeights = NULL; std::vector* pdfWeights = NULL; std::vector lheWeights; @@ -129,7 +131,7 @@ public : void FillHistograms(Long64_t entry, std::pair variation) override; void FillVBSHistograms(float weight, bool noBlind, std::pair variation); - bool PassesBaseSelection(Long64_t entry, bool tightLeps, Selection selection); + bool PassesBaseSelection(Long64_t entry, bool tightLeps, int selection); bool PassesVBSSelection(bool noBlind); bool PassesVBSBackgroundControlSelection(); bool PassesFullWZSelection(Long64_t entry); diff --git a/interface/WZSelectorBase.h b/interface/WZSelectorBase.h index 9f3f4c4c..b9bdfe9d 100644 --- a/interface/WZSelectorBase.h +++ b/interface/WZSelectorBase.h @@ -27,6 +27,12 @@ public : ScaleFactor* mIdSF_; ScaleFactor* mIsoSF_; + int VBSBackgroundControlLoose_Full, VBSBackgroundControl_Full, VBSBackgroundControl, VBSBackgroundControlLoose, VBSBackgroundControlATLAS, VBSselection_Loose_Full, VBSselection_Tight_Full, VBSselection_Tight, VBSselection_Loose, VBSselection_NoZeppenfeld, VBSselection_NoZeppenfeld_Full; + + + + + // Derived variables bool isVBS_; bool passesLeptonVeto; diff --git a/interface/dynEnum.h b/interface/dynEnum.h new file mode 100644 index 00000000..f97df394 --- /dev/null +++ b/interface/dynEnum.h @@ -0,0 +1,20 @@ +#include +#include + +class DynEnum +{ + private: + int count = 1; + std::unordered_map map; + public: + int addEnum(std::string name) { + map[name] = count; + return count++; + } + int getEnum(std::string name) { + return map.at(name); + } + bool foundEnum(std::string name) { + return map.find(name) != map.end(); + } +}; diff --git a/src/BuildFile.xml b/src/BuildFile.xml index 94ff9a4a..a9032981 100644 --- a/src/BuildFile.xml +++ b/src/BuildFile.xml @@ -2,6 +2,7 @@ + diff --git a/src/FakeRateSelector.cc b/src/FakeRateSelector.cc index 5be36887..65e8b506 100644 --- a/src/FakeRateSelector.cc +++ b/src/FakeRateSelector.cc @@ -1,53 +1,129 @@ #include "Analysis/VVAnalysis/interface/FakeRateSelector.h" #include +#include "TLorentzVector.h" -void FakeRateSelector::FillHistograms(Long64_t entry, std::pair variation) { - if (!passesLeptonVeto) - return; - if (l1Pt < 25 || l2Pt < 15) - return; - if (ZMass > 101.1876 || ZMass < 81.1876) - return; - if (type1_pfMETEt > 25) - return; - if (l3MtToMET > 30) - return; - if (!tightZLeptons()) - return; - - if (!IsGenMatched3l()) - return; - - float pt_fillval = l3Pt; - float eta_fillval = std::abs(l3Eta); - - float loose_weight = weight; - if (channel_ == eee || channel_ == emm) { - loose_weight /= eIdSF_->Evaluate2D(std::abs(l3Eta), l3Pt); + +#define Fill1D(NAME, VALUE_) HistFullFill(histMap1D_, NAME, variation.first, VALUE_, weight); +#define Fill2D(NAME, VALUE1_, VALUE2_) HistFullFill(histMap2D_, NAME, variation.first, VALUE1_, VALUE2_, weight); + +typedef ROOT::Math::LorentzVector>LorentzVector; +namespace VectorUtil=ROOT::Math::VectorUtil; + +void FakeRateSelector::Init(TTree *tree) { + allChannels_ = {{e, "e"}, {m, "m"}}; + + SelectorBase::Init(tree); + + TList* slaveClassList = (TList*)GetInputList()->Clone(); + slaveClassList->Add(new TNamed("isSlaveClass", "isSlaveClass")); + Analyzer.SetInputList(slaveClassList); + Analyzer.Init(tree); + Analyzer.SetBranches(); + fReader.SetTree(tree); + +} + +void FakeRateSelector::LoadBranchesNanoAOD(Long64_t entry, std::pair variation) { + Analyzer.LoadBranchesNanoAOD(entry,variation); + fReader.SetLocalEntry(entry); + + channelName_ = "all"; + if(Analyzer.looseLeptons.size() == 1) { + if(Analyzer.looseLeptons[0].Id() == PID_MUON) channelName_ = "m"; + else channelName_ = "e"; + } + + channel_ = channelMap_[channelName_]; +} + +void FakeRateSelector::SetBranchesNanoAOD() {} + +void FakeRateSelector::FillHistograms(Long64_t entry, std::pair variation) { + // need 3 loose leptons + size_t N_LEPS_FAKE = 1; + + //if (Analyzer.looseLeptons.size() != N_LEPS_FAKE) return; + //std::cout << Analyzer.looseLeptons.size() << "\n"; + int tagLepIdx = -1; + for(size_t i=0; i < Analyzer.looseLeptons.size(); ++i) { + auto lep = Analyzer.looseLeptons[i]; + if(Analyzer.passFakeableCuts(lep)) { + if(tagLepIdx >= 0) return; + if((*HLT_Mu8 && channel_ == m) || (*HLT_Ele8_CaloIdM_TrackIdM_PFJet30 && channel_ == e)) { + tagLepIdx = i; + } else { + return; + } + } } - else if (channel_ == eem || channel_ == mmm) { - loose_weight /= mIsoSF_->Evaluate2D(std::abs(l3Eta), l3Pt); + if(tagLepIdx < 0) return; + LorentzVector met(Analyzer.MET, 0, Analyzer.type1_pfMETPhi, Analyzer.MET); + LorentzVector tagLep = Analyzer.looseLeptons[tagLepIdx].v; + double mt = std::sqrt(2.*met.Et()*tagLep.pt()*(1.-cos(VectorUtil::DeltaPhi(tagLep, met)))); + + if (mt > 20) return; + if (met.Pt() > 20) return; // need to fix with skim + bool hasJet = false; + for(auto jet: Analyzer.goodJets) { + if(reco::deltaR(jet.v, tagLep) > 1.) { + hasJet = true; + break; + } } - passingLoose2D_->Fill(pt_fillval, eta_fillval, loose_weight); - passingLoose1DPt_->Fill(pt_fillval, loose_weight); - passingLoose1DEta_->Fill(eta_fillval, loose_weight); - if (lepton3IsTight()) { - passingTight2D_->Fill(pt_fillval, eta_fillval, weight); - passingTight1DPt_->Fill(pt_fillval, weight); - passingTight1DEta_->Fill(eta_fillval, weight); + if(!hasJet) return; + + float tagPt = tagLep.Pt(); + float tagEta = std::abs(tagLep.Eta()); + + // int closeIdx = Analyzer.getCloseJetIndex(tagLep); + // LorentzVector closeJet = Analyzer.get4Vector(PID_JET, closeIdx); + // if(Analyzer.LepRelPt(tagLep, closeJet) > 7.2) + // tagPt *= (1 + std::max(0, Im-I1)); + // else + // tagPt = std::max(tagPt, closeJet.Pt()*I2); + + float loose_weight = 1.; + // float loose_weight = weight; + // if (channel_ == eee || channel_ == emm) { + // loose_weight /= eIdSF_->Evaluate2D(std::abs(l3Eta), l3Pt); + // }p + // else if (channel_ == eem || channel_ == mmm) { + // loose_weight /= mIsoSF_->Evaluate2D(std::abs(l3Eta), l3Pt); + // } + + passingLoose2D_map[channel_]->Fill(tagPt, tagEta, loose_weight); + passingLoose1DPt_map[channel_]->Fill(tagPt, loose_weight); + passingLoose1DEta_map[channel_]->Fill(tagEta, loose_weight); + if (Analyzer.goodLeptons.size() == N_LEPS_FAKE) { + passingTight2D_map[channel_]->Fill(tagPt, tagEta, weight); + passingTight1DPt_map[channel_]->Fill(tagPt, weight); + passingTight1DEta_map[channel_]->Fill(tagEta, weight); } } void FakeRateSelector::SetupNewDirectory() { - WZSelectorBase::SetupNewDirectory(); + SelectorBase::SetupNewDirectory(); - const int nvarbins = 3; - double variable_pt_bins[nvarbins+1] = {10, 20, 30, FR_MAX_PT_}; - AddObject(passingTight2D_, ("passingTight2D_"+channelName_).c_str(), "#eta; p_{T} [GeV]", nvarbins, variable_pt_bins, 3, 0, 2.5); - AddObject(passingTight1DPt_, ("passingTight1DPt_"+channelName_).c_str(), "Tight leptons; p_{T} [GeV]", nvarbins, variable_pt_bins); - AddObject(passingTight1DEta_, ("passingTight1DEta_"+channelName_).c_str(), "Tight leptons; #eta", 3, 0, 2.5); + const int nPtBins = 5; + const int nEtaBins = 3; + double pt_bins[nPtBins+1] = {10, 15, 25, 35, 50, 70}; + double eta_El_bins[nEtaBins+1] = {0, 0.8, 1.479, 2.5}; + double eta_Mu_bins[nEtaBins+1] = {0, 1.2, 2.1, 2.4}; + std::unordered_map eta_bins= {{e, eta_El_bins}, {m, eta_Mu_bins}}; + for(auto pair : allChannels_) { + passingTight2D_map[pair.first] = nullptr; + passingTight1DPt_map[pair.first] = nullptr; + passingTight1DEta_map[pair.first] = nullptr; + passingLoose2D_map[pair.first] = nullptr; + passingLoose1DPt_map[pair.first] = nullptr; + passingLoose1DEta_map[pair.first] = nullptr; + + AddObject(passingTight2D_map[pair.first], ("passingTight2D_"+pair.second).c_str(), "#eta; p_{T} [GeV]", nPtBins, pt_bins, nEtaBins, eta_bins[pair.first]); + AddObject(passingTight1DPt_map[pair.first], ("passingTight1DPt_"+pair.second).c_str(), "Tight leptons; p_{T} [GeV]", nPtBins, pt_bins); + AddObject(passingTight1DEta_map[pair.first], ("passingTight1DEta_"+pair.second).c_str(), "Tight leptons; #eta", nEtaBins, eta_bins[pair.first]); - AddObject(passingLoose2D_, ("passingLoose2D_"+channelName_).c_str(), "#eta; p_{T} [GeV]", nvarbins, variable_pt_bins, 3, 0, 2.5); - AddObject(passingLoose1DPt_, ("passingLoose1DPt_"+channelName_).c_str(), "Loose leptons; p_{T} [GeV]", nvarbins, variable_pt_bins); - AddObject(passingLoose1DEta_, ("passingLoose1DEta_"+channelName_).c_str(), "Loose leptons; #eta", 3, 0, 2.5); + AddObject(passingLoose2D_map[pair.first], ("passingLoose2D_"+pair.second).c_str(), "#eta; p_{T} [GeV]", nPtBins, pt_bins, nEtaBins, eta_bins[pair.first]); + AddObject(passingLoose1DPt_map[pair.first], ("passingLoose1DPt_"+pair.second).c_str(), "Loose leptons; p_{T} [GeV]", nPtBins, pt_bins); + AddObject(passingLoose1DEta_map[pair.first], ("passingLoose1DEta_"+pair.second).c_str(), "Loose leptons; #eta", nEtaBins, eta_bins[pair.first]); + } } diff --git a/src/SelectorBase.cc b/src/SelectorBase.cc index ef578b8b..adaec8ef 100644 --- a/src/SelectorBase.cc +++ b/src/SelectorBase.cc @@ -25,78 +25,110 @@ void SelectorBase::Init(TTree *tree) fChain = tree; TString option = GetOption(); - + std::string selectionName_ = "Unknown"; + std::string yearName_ = "yrdefault"; + fwlog::LogLevel fwlogLevel = fwlog::kError; + + yrdefault = enumFactory.addEnum("yrdefault"); + for(auto pair: yearMap_) + pair.second = enumFactory.addEnum(pair.first); + selectionMap_.insert(addSelection_.begin(), addSelection_.end()); + for(auto pair: selectionMap_) + pair.second = enumFactory.addEnum(pair.first); + if (GetInputList() != nullptr) { TNamed* ntupleType = (TNamed *) GetInputList()->FindObject("ntupleType"); - TNamed* name = (TNamed *) GetInputList()->FindObject("name"); - TNamed* chan = (TNamed *) GetInputList()->FindObject("channel"); - TNamed* selection = (TNamed *) GetInputList()->FindObject("selection"); + TNamed* name = (TNamed *) GetInputList()->FindObject("name"); + TNamed* chan = (TNamed *) GetInputList()->FindObject("channel"); + TNamed* selection = (TNamed *) GetInputList()->FindObject("selection"); TNamed* year = (TNamed *) GetInputList()->FindObject("year"); + TNamed* logLevel = (TNamed *) GetInputList()->FindObject("logLevel"); if (ntupleType != nullptr) { std::string ntupleName = ntupleType->GetTitle(); - if (ntupleName == "NanoAOD") - ntupleType_ = NanoAOD; - else if (ntupleName == "UWVV") - ntupleType_ = UWVV; - else if (ntupleName == "Bacon") - ntupleType_ = Bacon; - else + if (ntupleName == "NanoAOD") ntupleType_ = NanoAOD; + else if (ntupleName == "UWVV") ntupleType_ = UWVV; + else if (ntupleName == "Bacon") ntupleType_ = Bacon; + else throw std::invalid_argument("Unsupported ntuple type!"); + } else { + std::cerr << "INFO: Assuming NanoAOD ntuples" << std::endl; + ntupleType_ = NanoAOD; + } + + if (name != nullptr) { + name_ = name->GetTitle(); + } else { + name_ = GetNameFromFile(); } - else { - std::cerr << "INFO: Assuming NanoAOD ntuples" << std::endl; - ntupleType_ = NanoAOD; - } - - if (name != nullptr) { - name_ = name->GetTitle(); - } - else { - name_ = GetNameFromFile(); - } - if (name_.empty()){ - std::cerr << "INFO: Using default name \"Unknown\" for file" << std::endl; - name_ = "Unknown"; - } - if(year != nullptr) { - year_ = yearMap_[year->GetTitle()]; + + if(year != nullptr ) { + yearName_ = year->GetTitle(); } if (chan != nullptr) { channelName_ = chan->GetTitle(); + } else if (ntupleType_ == UWVV) { + channelName_ = fChain->GetTree()->GetDirectory()->GetName(); } - else if (ntupleType_ == UWVV) - channelName_ = fChain->GetTree()->GetDirectory()->GetName(); - if (selection != nullptr) { + + if (selection != nullptr) { selectionName_ = selection->GetTitle(); } + + if(logLevel != nullptr) { + std::string logName = logLevel->GetTitle(); + if( logName == "DEBUG") fwlogLevel = fwlog::kDebug; + else if( logName == "INFO") fwlogLevel = fwlog::kInfo; + else if( logName == "WARNING") fwlogLevel = fwlog::kWarning; + else fwlogLevel = fwlog::kError; + } + } + fwlog::setPresentLogLevel(fwlogLevel); + + // Setup Name + if (name_.empty()){ + std::cerr << "INFO: Using default name \"Unknown\" for file" << std::endl; + name_ = "Unknown"; } - if (selectionMap_.find(selectionName_) != selectionMap_.end()) { - selection_ = selectionMap_[selectionName_]; + // Setup Selection + if (enumFactory.foundEnum(selectionName_)) { + selection_ = enumFactory.getEnum(selectionName_); + } else { + selection_ = enumFactory.addEnum(selectionName_); + fwLog(fwlog::kDebug) << "Selection ("<< selectionName_ << ") not found.\n"; } - else - throw std::invalid_argument(selectionName_ + " is not a valid selection!"); + + // Setup Year + if (enumFactory.foundEnum(yearName_)) { + year_ = enumFactory.getEnum(yearName_); + } else { + year_ = enumFactory.addEnum(yearName_); + fwLog(fwlog::kDebug)<< "Year ("<< yearName_ << ") not found.\n"; + } + + fwlog::setPresentLogLevel(fwlogLevel); + + isMC_ = false; if (name_.find("data") == std::string::npos){ - isMC_ = true; + isMC_ = true; } if (doSystematics_ && isMC_ && !isNonprompt_) - variations_.insert(systematics_.begin(), systematics_.end()); - + variations_.insert(systematics_.begin(), systematics_.end()); + if (channelMap_.find(channelName_) != channelMap_.end()) - channel_ = channelMap_[channelName_]; + channel_ = channelMap_[channelName_]; else { - std::string message = "Invalid channel choice! "; - message += "Choice was " + channelName_ + "\n"; - message += "Valid choices are: "; - for (const auto& chan : channelMap_) + std::string message = "Invalid channel choice! "; + message += "Choice was " + channelName_ + "\n"; + message += "Valid choices are: "; + for (const auto& chan : channelMap_) message += chan.first + ", " ; throw std::invalid_argument(message); } - // only make the directory iff class isn't being run as a slave class ///// TNamed* isSlaveClass = (TNamed *) GetInputList()->FindObject("isSlaveClass"); if(isSlaveClass != nullptr) return; diff --git a/src/TTTSelector.cc b/src/TTTSelector.cc deleted file mode 100644 index 0c9462e0..00000000 --- a/src/TTTSelector.cc +++ /dev/null @@ -1,484 +0,0 @@ -#include "Analysis/VVAnalysis/interface/TTTSelector.h" - -#include -#include -#include "TParameter.h" - -#define Fill1D(NAME, VALUE_) HistFullFill(histMap1D_, NAME, variation.first, VALUE_, weight); -#define Fill2D(NAME, VALUE1_, VALUE2_) HistFullFill(histMap2D_, NAME, variation.first, VALUE1_, VALUE2_, weight); - - - -enum Lepton {Muon=13, Electron=11}; - -typedef ROOT::Math::LorentzVector> LorentzVector; - -void TTTSelector::SlaveBegin(TTree *) { - return; - pileupSF_ = (ScaleFactor *) GetInputList()->FindObject("pileupSF"); - if (pileupSF_ == nullptr ) - Abort("Must pass pileup weights SF"); - eIdSF_ = (ScaleFactor *) GetInputList()->FindObject("electronTightIdSF"); - if (eIdSF_ == nullptr ) - Abort("Must pass electron ID SF"); - eGsfSF_ = (ScaleFactor *) GetInputList()->FindObject("electronGsfSF"); - if (eGsfSF_ == nullptr ) - Abort("Must pass electron GSF SF"); - mIdSF_ = (ScaleFactor *) GetInputList()->FindObject("muonTightIdSF"); - if (mIdSF_ == nullptr ) - Abort("Must pass muon ID SF"); - mIsoSF_ = (ScaleFactor *) GetInputList()->FindObject("muonIsoSF"); - if (mIsoSF_ == nullptr ) - Abort("Must pass muon Iso SF"); - - prefireEff_ = (TEfficiency*) GetInputList()->FindObject("prefireEfficiencyMap"); - if (prefireEff_ == nullptr ) - Abort("Must pass prefiring efficiency map"); - -} - -void TTTSelector::Init(TTree *tree){ - b.SetTree(tree); - - allChannels_ = {{ee, "ee"}, {mm, "mm"}, {em, "em"}, {all, "all"}}; - - hists1D_ = {"CutFlow", "ZMass", "ptl1", "etal1", "ptl2", "etal2", "SR", "bjetpt", "jetpt", "nbjet", "njet", "jetphi","bjetphi", "phil1", "phil2","HT","MET", "nelec", "nmuon", "lept_charge", "before", "after", "nGoodParts"}; - hists2D_ = {"bJetvsJets"}; - - SelectorBase::Init(tree); - -} - -void TTTSelector::SetBranchesNanoAOD() { - // NECESSARY!!!! - b.CleanUp(); - - b.SetBranch("nElectron", nElectron); - b.SetBranch("Electron_pt", Electron_pt); - b.SetBranch("Electron_eta", Electron_eta); - b.SetBranch("Electron_phi", Electron_phi); - b.SetBranch("Electron_charge", Electron_charge); - b.SetBranch("Electron_mass", Electron_mass); - b.SetBranch("Electron_mass", Electron_miniPFRelIso_all); - b.SetBranch("Electron_miniPFRelIso_all", Electron_miniPFRelIso_all); - b.SetBranch("Electron_dxy", Electron_dxy); - b.SetBranch("Electron_dz", Electron_dz); - b.SetBranch("Electron_sip3d", Electron_sip3d); - if(year_ == yr2018) { - b.SetBranch("Electron_mvaFall17V2noIso", Electron_MVA); - b.SetBranch("Electron_cutBased", Electron_cutBased); - } else if(year_ == yr2017) { - b.SetBranch("Electron_mvaFall17V1noIso", Electron_MVA); - b.SetBranch("Electron_cutBased_Fall17_V1", Electron_cutBased); - } else if(year_ == yr2016 || year_ == yrdefault) { - b.SetBranch("Electron_mvaSpring16GP", Electron_MVA); - b.SetBranch("Electron_cutBased_Sum16", Electron_cutBased); - } - - b.SetBranch("nMuon", nMuon); - b.SetBranch("Muon_pt", Muon_pt); - b.SetBranch("Muon_eta", Muon_eta); - b.SetBranch("Muon_phi", Muon_phi); - b.SetBranch("Muon_tightId", Muon_tightId); - b.SetBranch("Muon_mediumId", Muon_mediumId); - b.SetBranch("Muon_tkIsoId", Muon_tkIsoId); - b.SetBranch("Muon_pfRelIso04_all", Muon_pfRelIso04_all); - b.SetBranch("Muon_miniPFRelIso_all", Muon_miniPFRelIso_all); - b.SetBranch("Muon_charge", Muon_charge); - b.SetBranch("Muon_mass", Muon_mass); - b.SetBranch("Muon_dxy", Electron_dxy); - b.SetBranch("Muon_dz", Electron_dz); - b.SetBranch("Muon_sip3d", Electron_sip3d); - - b.SetBranch("nJet", nJet); - b.SetBranch("Jet_btagCSVV2", Jet_btagCSVV2); - b.SetBranch("Jet_btagDeepB", Jet_btagDeepB); - b.SetBranch("Jet_eta", Jet_eta); - b.SetBranch("Jet_phi", Jet_phi); - b.SetBranch("Jet_pt", Jet_pt); - b.SetBranch("Jet_mass", Jet_mass); - - b.SetBranch("Jet_neHEF", Jet_neHEF); - b.SetBranch("Jet_neEmEF", Jet_neEmEF); - b.SetBranch("Jet_nConstituents", Jet_nConstituents); - b.SetBranch("Jet_chHEF", Jet_chHEF); - b.SetBranch("Jet_chEmEF", Jet_chEmEF); - - b.SetBranch("MET_pt", MET); - b.SetBranch("MET_phi", type1_pfMETPhi); - - if (isMC_) { - b.SetBranch("genWeight", genWeight); - b.SetBranch("Pileup_nPU", numPU); - } - -} - -void TTTSelector::SetBranchesUWVV() { - return; -} - -void TTTSelector::LoadBranchesUWVV(Long64_t entry, std::pair variation) { - return; -} - -/// Make to seperate fuctionality -void TTTSelector::clearValues() { - weight = 1; - HT = 0; - nTightJet = 0; - nBJets = 0; - goodParts.clear(); - charge.clear(); -} - -void TTTSelector::LoadBranchesNanoAOD(Long64_t entry, std::pair variation) { - clearValues(); - b.SetEntry(entry); - - if (nElectron > N_KEEP_MU_E_ || nMuon > N_KEEP_MU_E_) { - std::string message = "Found more electrons or muons than max read number.\n Found "; - message += std::to_string(nElectron); - message += " electrons.\n Found "; - message += std::to_string(nMuon); - message += " Muons\n --> Max read number was "; - message += std::to_string(N_KEEP_MU_E_); - message += "\nExiting because this can cause problems. Increase N_KEEP_MU_E_ to avoid this error.\n"; - throw std::domain_error(message); - } - - ///////////////////// - // Setup Electrons // - ///////////////////// - - for (size_t i = 0; i < nElectron; i++) { - // if(IsGoodMVAElectron2016(i)) { - if(IsGoodCBElectron(i)) { - // // Extra Iso requirement - // LorentzVector lep(Electron_pt[i], Electron_eta[i], Electron_phi[i], Electron_mass[i]); - // if(!passFullIso(lep, 0.8, 7.2)) continue; - - // Setup goodPart - goodParts.push_back(GoodPart(Electron_pt[i], Electron_eta[i], Electron_phi[i], Electron_mass[i])); - goodParts.back().SetPdgId(Electron*Electron_charge[i]); - charge.push_back(Electron_charge[i]); - } - } - - ///////////////// - // Setup Muons // - ///////////////// - - for (size_t i = 0; i < nMuon; i++) { - if(IsGoodMuon(i)) { - // // Extra Iso requirement - // LorentzVector lep(Muon_pt[i], Muon_eta[i], Muon_phi[i], Muon_mass[i]); - // if(!passFullIso(lep, 0.76, 7.2)) continue; - - // Setup goodPart - goodParts.push_back(GoodPart(Muon_pt[i], Muon_eta[i], Muon_phi[i], Muon_mass[i])); - goodParts.back().SetPdgId(Muon*Muon_charge[i]); - charge.push_back(Muon_charge[i]); - - } - } - - //////////////// - // Setup Jets // - //////////////// - - for(size_t i = 0; i < nJet; i++) { - if(goodParts.size() != 2) break; - - //common to both jets, may put more here - if(!isLooseJetId(i)) continue; - - // regular jets - if(IsGoodJet(i)) { - nTightJet++; - HT += Jet_pt[i]; - } - - // bjet - if(IsGoodBJet(i)) nBJets++; - } - - channel_ = channelMap_[channelName_]; - if(goodParts.size() != 2) { - channel_ = Unknown; - channelName_ = "Unknown"; - } else if(goodParts[0].Id() == Muon && goodParts[1].Id() == Muon) { - channel_ = mm; - channelName_ = "mm"; - } else if(goodParts[0].Id() == Electron && goodParts[1].Id() == Electron) { - channel_ = ee; - channelName_ = "ee"; - } else { - channel_ = em; - channelName_ ="em"; - /// fix order of leptons by pt - if(goodParts[0].Pt() < goodParts[1].Pt()) { - std::swap(goodParts[0], goodParts[1]); - } - } - - if (isMC_) { - ApplyScaleFactors(); - } -} - - -void TTTSelector::ApplyScaleFactors() { - weight *= (genWeight > 0) ? 1.0 : -1.0; - // weight = 1; - // skipping for now - return; - //// to add!! -} - -//////////////////////////////////////////////////// -// Functions for defining particles: used in main // -// loop of leptons/jets for getting multiplicity // -//////////////////////////////////////////////////// - - - -bool TTTSelector::IsGoodMuon(size_t index) { - return ( (Muon_pt[index] > 20) && - (abs(Muon_eta[index]) < 2.4) && - Muon_mediumId[index] && - (Muon_miniPFRelIso_all[index] < 0.16) && - (Muon_dz[index] < 0.1) && - (Muon_dxy[index] < 0.05) && - (Muon_sip3d[index] < 4)); -} - -bool TTTSelector::IsGoodCBElectron(size_t index) { - return ((Electron_pt[index] > 20) && - (abs(Electron_eta[index]) < 2.5) && - (Electron_miniPFRelIso_all[index] < 0.12) && - (Electron_cutBased[index] == 4) && - (Electron_dz[index] < 0.1) && - (Electron_dxy[index] < 0.05) && - (Electron_sip3d[index] < 4)); -} - -bool TTTSelector::IsGoodMVAElectron2016(size_t index) { - bool mvaRec = false; - if(abs(Electron_eta[index]) < 0.8) - mvaRec = Electron_MVA[index] > std::max(0.52, 0.77-0.025*(Electron_pt[index]-15)); - else if(abs(Electron_eta[index]) < 1.479) - mvaRec = Electron_MVA[index] > std::max(0.11, 0.56-0.045*(Electron_pt[index]-15)); - else if(abs(Electron_eta[index]) < 2.5) - mvaRec = Electron_MVA[index] > std::max(-0.01, 0.48-0.049*(Electron_pt[index]-15)); - - return ((Electron_pt[index] > 20) && - (Electron_miniPFRelIso_all[index] < 0.12) && - (mvaRec) && - (Electron_dz[index] < 0.1) && - (Electron_dxy[index] < 0.05) && - (Electron_sip3d[index] < 4) ); -} - -bool TTTSelector::IsGoodMVAElectron2017(size_t index) { - bool mvaRec = false; - if(abs(Electron_eta[index]) < 0.8) - mvaRec = std::max(0.52, 0.77-0.025*(Electron_pt[index]-15)); - else if(abs(Electron_eta[index]) < 1.479) - mvaRec = std::max(0.11, 0.56-0.045*(Electron_pt[index]-15)); - else if(abs(Electron_eta[index]) < 2.5) - mvaRec = std::max(-0.01, 0.48-0.049*(Electron_pt[index]-15)); - - return ((Electron_pt[index] > 20) && - mvaRec); -} - -//40 -bool TTTSelector::IsGoodJet(size_t index) { - return ((Jet_pt[index] > 40) && - (abs(Jet_eta[index]) < 2.4) && - isOverlap(index) - ); -} - -/// TODO: add toggle for different btag stuff - -//25 -bool TTTSelector::IsGoodBJet(size_t index) { - return ((Jet_pt[index] > 25) && - (abs(Jet_eta[index]) < 2.4) && - (Jet_btagCSVV2[index] > 0.8484) && - //(Jet_btagDeepB[index] > 0.6324) && - isOverlap(index) - ); -} - -bool TTTSelector::isLooseJetId(size_t index) { - return (Jet_neHEF[index] < 0.99 && - Jet_neEmEF[index] < 0.99 && - Jet_nConstituents[index] > 1 && - Jet_chHEF[index] > 0 && - Jet_chEmEF[index] < 0.99 - ); -} - -bool TTTSelector::isTightJetId(size_t index) { - return ((Jet_neHEF[index] < 0.9) && - (Jet_neEmEF[index] < 0.9) && - (Jet_nConstituents[index] > 1) && - (Jet_chHEF[index] > 0) - ); -} - -bool TTTSelector::passFullIso(LorentzVector& lep, int I2, int I3) { - LorentzVector closeJet; - double minDR = 10; - for(size_t index = 0; index < nJet; index++) { - LorentzVector jet(Jet_pt[index], Jet_eta[index], Jet_phi[index], Jet_mass[index]); - double dr = reco::deltaR(jet, lep); - if(minDR > dr) { - closeJet = jet; - minDR = dr; - } - } - - if(lep.Pt()/closeJet.Pt() > I2 ) return true; - - auto diff = closeJet.Vect() - lep.Vect(); - auto cross = diff.Cross(lep.Vect()); - return (cross.Mag2()/diff.Mag2() > I3*I3); -} - -bool TTTSelector::isOverlap(size_t index) { - LorentzVector tmp(Jet_pt[index], Jet_eta[index], Jet_phi[index], Jet_mass[index]); - double dR = 0.4; - bool passOverlap = true; - for(auto lep: goodParts) { - passOverlap = (reco::deltaR(tmp, lep.v) > dR); - if(!passOverlap) return false; - } - return true; -} - -void TTTSelector::FillHistograms(Long64_t entry, std::pair variation) { - - int step = 0; - Fill1D("CutFlow", 0); - - Fill1D("nGoodParts", goodParts.size()); - - /// 2 good leptons - if(goodParts.size() != 2) return; - Fill1D("CutFlow", ++step); - - - // std::cout<<"charge array : " <= 0) - Fill1D("before", 0); - - if(goodParts[0].Charge() * goodParts[1].Charge() <= 0) - Fill1D("before", 1); - - // first lep requirement - if(goodParts[0].Pt() < 25) return; - Fill1D("CutFlow", ++step); - - // same sign requirement - // if(goodParts[0].Charge() * goodParts[1].Charge() <= 0) return; - Fill1D("CutFlow", ++step); - - // met cut 50 - if (MET < 50) return; - Fill1D("CutFlow", ++step); - - // ht cut - if(HT <300 ) return; - Fill1D("CutFlow", ++step); - - // jet cut 4 - if(nTightJet < 4) return; - Fill1D("CutFlow", ++step); - - // bjet cut 2 - if(nBJets < 2) return; - Fill1D("CutFlow", ++step); - - if(goodParts[0].Charge() * goodParts[1].Charge() >= 0) - Fill1D("after", 0); - if(goodParts[0].Charge() * goodParts[1].Charge() <=0) - Fill1D("after", 1); - - // // veto cut - // if(!passesLeptonVeto) - // Fill1D("CutFlow", ++step); - // in SR stuff - - Fill1D("ptl1", goodParts[0].Pt()); - Fill1D("ptl2", goodParts[1].Pt()); - Fill1D("SR", getSRBin()); - Fill1D("njet", nTightJet); - Fill1D("nbjet", nBJets); - Fill1D("phil1", goodParts[0].Phi()); - Fill1D("phil2", goodParts[1].Phi()); - Fill1D("HT", HT ); - Fill1D("MET", MET); - Fill1D("nelec", nElectron); - Fill1D("nmuon", nMuon); - Fill1D("lept_charge", charge[0]); - Fill1D("lept_charge", charge[1]); - - Fill2D("bJetvsJets", nTightJet, nBJets); - - - for(size_t i = 0; i < nJet; i++) { - if(IsGoodJet(i)) { - Fill1D("jetpt", Jet_pt[i]); - Fill1D("jetphi",Jet_phi[i]); - } - if(IsGoodBJet(i)) { - Fill1D("bjetpt", Jet_pt[i]); - Fill1D("bjetphi",Jet_phi[i]); - } - } - } - -void TTTSelector::SetupNewDirectory() { - SelectorBase::SetupNewDirectory(); - - InitializeHistogramsFromConfig(); -} - -int TTTSelector::getSRBin() { - if(nBJets == 2) { - if(nTightJet <= 5) return 0; - if(nTightJet == 6) return 1; - else if(nTightJet == 7) return 2; - else if(nTightJet >= 8) return 3; - } else if(nBJets == 3) { - if(nTightJet == 5) return 4; - else if(nTightJet == 6) return 4; - else if(nTightJet == 7) return 5; - else if(nTightJet >= 8) return 5; - } else if(nBJets >= 4) { - if(nTightJet >= 5) return 6; - } - return -1; - - // if(nBJets == 2) { - // if(nTightJet <= 5) return 0; - // if(nTightJet == 6) return 1; - // else if(nTightJet == 7) return 2; - // else if(nTightJet >= 8) return 3; - // } else if(nBJets == 3) { - // if(nTightJet == 5) return 4; - // else if(nTightJet == 6) return 5; - // else if(nTightJet == 7) return 6; - // else if(nTightJet >= 8) return 7; - // } else if(nBJets >= 4) { - // if(nTightJet >= 5) return 8; - // } - // return -1; - - - } - diff --git a/src/ThreeLepSelector.cc b/src/ThreeLepSelector.cc index df297789..f3660034 100644 --- a/src/ThreeLepSelector.cc +++ b/src/ThreeLepSelector.cc @@ -12,7 +12,11 @@ #define READFROM(data, index, size) (((data) & GETMASK((index), (size))) >> (index)) #define CHGPT(index) (Electron_eCorr[index]) -#define CLOSEJET_REWEIGHT +// #define CLOSEJET_REWEIGHT +// #define SELECTION 2 +// #define TRIGGER +//#define OS +// #define USETREE typedef std::bitset IntBits; @@ -53,91 +57,108 @@ void ThreeLepSelector::SetScaleFactors() { } void ThreeLepSelector::Init(TTree *tree) { + // Setup maps + selectionMap_ = {{"MVAStudy", MVAStudy}, + {"FourTopMVAEl", FourTopMVAEl}, + {"FourTopCutBasedEl", FourTopCutBasedEl}, + {"FakeRate", FakeRate},}; + yearMap_ = {{"2016", yr2016}, {"2017", yr2017}, {"2018", yr2018},}; + + // Continue b.SetTree(tree); - - allChannels_ = {{mm, "mm"}, {ee, "ee"}, {em, "em"}, {all, "all"}, {lll, "lll"}}; - + + + //allChannels_ = {{mm, "mm"}, {ee, "ee"}, {em, "em"}, {all, "all"}, {lll, "lll"}}; + allChannels_ = {{SS, "SS"}, {OS, "OS"}, {mult, "mult"}, {all, "all"}, {one, "one"}}; + hists1D_ = { - "CutFlow", "ZMass", "ptl1", "etal1", "ptl2", "etal2", "SR", - "bjetpt", "jetpt", "nbjet", "njet", "nleps", "CRW_nbjet", "CRW_njet", - "CRZ_nbjet", "CRZ_njet", "CRZ_HT", "CRZ_Met", "Met", "HT", "weight", - "CRW_HT", "CRW_Met", "CRZ_ptl3", "sphericity", "centrality", + "CutFlow", "ptl1", "etal1", "ptl2", "etal2", "SR", + "bjetpt", "jetpt", "nbjet", "njet", "nleps", + //"CRZ_nbjet", "CRZ_njet", "CRZ_HT", "CRZ_Met", + "Met", "HT", "weight","sphericity", "centrality", + //"CRW_HT", "CRW_Met", "CRZ_ptl3", "CRW_nbjet", "CRW_njet", + "ptj1", "ptj2", "ptj3", "ptj1OverHT", + "etaj1", "etaj2","etaj3", "etab1", "etab2","etab3", "dphi_l1j1","dphi_l1j2","dphi_l1j3", + "ptb1", "ptb2", "ptb3", "ptb1OverHT", + "dilepMass", "dilepCharge", "DRLep", "DRjet", "dijetMass", + "Shape1", "Shape2", "LepCos", "JLep1Cos", "JLep2Cos", "JBCos", "DRjb", "etaj", "etab", + "ntightbjet", "nloosebjet","nlooseleps", }; hists2D_ = {"bJetvsJets", "Beff_b_btag", "Beff_j_btag", "Beff_b", "Beff_j"}; SelectorBase::Init(tree); + TNamed* name = (TNamed *) GetInputList()->FindObject("name"); + std::string name_tmp = name->GetTitle(); + if(name_tmp.find("2016") != std::string::npos) year_ = yr2016; + else if(name_tmp.find("2017") != std::string::npos) year_ = yr2017; + else if(name_tmp.find("2018") != std::string::npos) year_ = yr2018; + else year_ = yr2016; + + fReader.SetTree(tree); + +#ifdef USETREE + AddObject(treeMap["tree"], "testTree", "testTree"); + treeMap["tree"]->Branch("NJets", &bNJets); + treeMap["tree"]->Branch("NBJets", &bnBJets); + treeMap["tree"]->Branch("NlooseBJets", &bnlBJets); + treeMap["tree"]->Branch("NtightBJets", &bntBJets); + treeMap["tree"]->Branch("NLeps", &bnLeps); + treeMap["tree"]->Branch("NlooseLeps", &bnlLeps); + treeMap["tree"]->Branch("DilepCharge", &bDilepCharge); + treeMap["tree"]->Branch("HT", &bHT); + treeMap["tree"]->Branch("MET", &bMET); + treeMap["tree"]->Branch("l1Pt", &bl1Pt); + treeMap["tree"]->Branch("l2Pt", &bl2Pt); + treeMap["tree"]->Branch("lepMass", &blMass); + treeMap["tree"]->Branch("jetMass", &bjMass); + treeMap["tree"]->Branch("jetDR", &bjdr); + treeMap["tree"]->Branch("sphericity",&bsphere); + treeMap["tree"]->Branch("centrality", &bCentral); + treeMap["tree"]->Branch("j1Pt", &bj1Pt); + treeMap["tree"]->Branch("j2Pt", &bj2Pt); + treeMap["tree"]->Branch("j3Pt", &bj3Pt); + treeMap["tree"]->Branch("j4Pt", &bj4Pt); + treeMap["tree"]->Branch("j5Pt", &bj5Pt); + treeMap["tree"]->Branch("j6Pt", &bj6Pt); + treeMap["tree"]->Branch("j7Pt", &bj7Pt); + treeMap["tree"]->Branch("j8Pt", &bj8Pt); + treeMap["tree"]->Branch("b1Pt",&bb1Pt); + treeMap["tree"]->Branch("b2Pt",&bb2Pt); + treeMap["tree"]->Branch("b3Pt",&bb3Pt); + treeMap["tree"]->Branch("b4Pt",&bb4Pt); + treeMap["tree"]->Branch("weight",&weight); + treeMap["tree"]->Branch("Shape1",&bShape1); + treeMap["tree"]->Branch("Shape2",&bShape2); + std::cout << "here" << "\n"; +#endif + } void ThreeLepSelector::SetBranchesNanoAOD() { // NECESSARY!!!! b.CleanUp(); + fReader.Restart(); + if(year_ == yr2018 || year_ == yrdefault) { + Electron_MVA = {fReader, "Electron_mvaFall17V2noIso"}; + } else if(year_ == yr2017) { + Electron_MVA = {fReader, "Electron_mvaFall17V1noIso"}; + } else if(year_ == yr2016 || year_ == yrdefault) { + Electron_MVA = {fReader, "Electron_mvaSpring16GP"}; + Electron_cutBased = {fReader, "Electron_cutBased_Sum16"}; + } +#ifdef TRIGGER b.SetBranch("HLT_DoubleMu8_Mass8_PFHT300", HLT_DoubleMu8_Mass8_PFHT300); b.SetBranch("HLT_Mu8_Ele8_CaloIdM_TrackIdM_Mass8_PFHT300", HLT_Mu8_Ele8_CaloIdM_TrackIdM_Mass8_PFHT300); b.SetBranch("HLT_DoubleEle8_CaloIdM_TrackIdM_Mass8_PFHT300", HLT_DoubleEle8_CaloIdM_TrackIdM_Mass8_PFHT300); b.SetBranch("HLT_AK8PFJet450", HLT_AK8PFJet450); b.SetBranch("HLT_PFJet450", HLT_PFJet450); - - b.SetBranch("nElectron", nElectron); - b.SetBranch("Electron_pt", Electron_pt); - b.SetBranch("Electron_eta", Electron_eta); - b.SetBranch("Electron_phi", Electron_phi); - b.SetBranch("Electron_charge", Electron_charge); - b.SetBranch("Electron_mass", Electron_mass); - b.SetBranch("Electron_miniPFRelIso_all", Electron_miniPFRelIso_all); - b.SetBranch("Electron_dxy", Electron_dxy); - b.SetBranch("Electron_dz", Electron_dz); - b.SetBranch("Electron_sip3d", Electron_sip3d); - b.SetBranch("Electron_lostHits", Electron_lostHits); - b.SetBranch("Electron_convVeto", Electron_convVeto); - if(year_ == yr2018) { - b.SetBranch("Electron_mvaFall17V2noIso", Electron_MVA); - b.SetBranch("Electron_cutBased", Electron_cutBased); - } else if(year_ == yr2017) { - b.SetBranch("Electron_mvaFall17noIso", Electron_MVA); - b.SetBranch("Electron_cutBased", Electron_cutBased); - } else if(year_ == yr2016 || year_ == yrdefault) { - b.SetBranch("Electron_mvaSpring16GP", Electron_MVA); - b.SetBranch("Electron_cutBased_Sum16", Electron_cutBased); - } - b.SetBranch("Electron_tightCharge", Electron_tightCharge); - b.SetBranch("Electron_sieie", Electron_sieie); - b.SetBranch("Electron_hoe", Electron_hoe); - b.SetBranch("Electron_deltaEtaSC", Electron_deltaEtaSC); - b.SetBranch("Electron_eInvMinusPInv", Electron_eInvMinusPInv); - b.SetBranch("Electron_dr03EcalRecHitSumEt", Electron_dr03EcalRecHitSumEt); - b.SetBranch("Electron_dr03HcalDepth1TowerSumEt", Electron_dr03HcalDepth1TowerSumEt); - b.SetBranch("Electron_dr03TkSumPt", Electron_dr03TkSumPt); - b.SetBranch("Electron_vidNestedWPBitmapSum16", Electron_vidBitmap); - b.SetBranch("Electron_eCorr", Electron_eCorr); - - b.SetBranch("nMuon", nMuon); - b.SetBranch("Muon_pt", Muon_pt); - b.SetBranch("Muon_eta", Muon_eta); - b.SetBranch("Muon_phi", Muon_phi); - b.SetBranch("Muon_mass", Muon_mass); - b.SetBranch("Muon_charge", Muon_charge); - b.SetBranch("Muon_mediumId", Muon_mediumId); - b.SetBranch("Muon_miniPFRelIso_all", Muon_miniPFRelIso_all); - b.SetBranch("Muon_dxy", Muon_dxy); - b.SetBranch("Muon_dz", Muon_dz); - b.SetBranch("Muon_sip3d", Muon_sip3d); - b.SetBranch("Muon_isGlobal", Muon_isGlobal); - b.SetBranch("Muon_isTracker", Muon_isTracker); - b.SetBranch("Muon_isPFcand", Muon_isPFcand); - b.SetBranch("Muon_tightCharge", Muon_tightCharge); +#endif // TRIGGER + - b.SetBranch("nJet", nJet); - b.SetBranch("Jet_btagCSVV2", Jet_btagCSVV2); - b.SetBranch("Jet_btagDeepB", Jet_btagDeepB); - b.SetBranch("Jet_eta", Jet_eta); - b.SetBranch("Jet_phi", Jet_phi); - b.SetBranch("Jet_pt", Jet_pt); - b.SetBranch("Jet_mass", Jet_mass); - b.SetBranch("Jet_jetId", Jet_jetId); - b.SetBranch("Jet_hadronFlavour", Jet_hadronFlavour); - b.SetBranch("Jet_rawFactor", Jet_rawFactor); - + b.SetBranch("nElectron", nElectron); + b.SetBranch("nMuon", nMuon); + b.SetBranch("nJet", nJet); b.SetBranch("MET_pt", MET); b.SetBranch("MET_phi", type1_pfMETPhi); @@ -145,10 +166,10 @@ void ThreeLepSelector::SetBranchesNanoAOD() { b.SetBranch("luminosityBlock", lumi); if(applyScaleFactors_) { - b.SetBranch("GenPart_pdgId",GenPart_pdgId); - b.SetBranch("GenPart_genPartIdxMother", GenPart_genPartIdxMother); - b.SetBranch("GenPart_status", GenPart_status); b.SetBranch("nGenPart", nGenPart); + GenPart_pdgId = {tmpReader, "GenPart_pdgId"}; + GenPart_genPartIdxMother = {tmpReader, "GenPart_genPartIdxMother"}; + GenPart_status = {tmpReader, "GenPart_status"}; } b.SetBranch("Flag_goodVertices", Flag_goodVertices); @@ -158,10 +179,11 @@ void ThreeLepSelector::SetBranchesNanoAOD() { b.SetBranch("Flag_EcalDeadCellTriggerPrimitiveFilter", Flag_EcalDeadCellTriggerPrimitiveFilter); b.SetBranch("Flag_BadPFMuonFilter", Flag_BadPFMuonFilter); b.SetBranch("Flag_ecalBadCalibFilter", Flag_ecalBadCalibFilter); - + #ifdef CLOSEJET_REWEIGHT - b.SetBranch("Jet_L1", Jet_L1); - b.SetBranch("Jet_L2L3", Jet_L2L3); + Jet_L1 = {fReader, "Jet_L1"}; + Jet_L2L3 = {fReader, "Jet_L2L3"}; + #endif // CLOSEJET_REWEIGHT if (isMC_) { @@ -169,18 +191,34 @@ void ThreeLepSelector::SetBranchesNanoAOD() { b.SetBranch("Pileup_nPU", numPU); b.SetBranch("Pileup_nTrueInt", Pileup_nTrueInt); } - - mvaValues[CBID_LOOSE] = {{-0.46, -0.48, 0, -0.85}, - {-0.03, -0.67, 0, -0.91}, - {0.06, -0.49, 0, -0.83}}; - mvaValues[CBID_TIGHT] = {{10, 0.77, 0, 0.52}, - {10, 0.56, 0, 0.11}, - {10, 0.48, 0, -0.01}}; - /// Setup interpolation used for 25>pt>15 - for(auto& pair: mvaValues) { - for(auto& vec: pair.second) { - vec[2] = (vec[3]-vec[1])/10; + + if(year_ == yr2016) { + mvaValues[CBID_LOOSE] = {{-0.46, -0.48, 0, -0.85}, + {-0.03, -0.67, 0, -0.91}, + {0.06, -0.49, 0, -0.83}}; + mvaValues[CBID_TIGHT] = {{10, 0.77, 0, 0.52}, + {10, 0.56, 0, 0.11}, + {10, 0.48, 0, -0.01}}; + /// Setup interpolation used for 25>pt>15 + for(auto& pair: mvaValues) { + for(auto& vec: pair.second) { + vec[2] = (vec[3]-vec[1])/10; + } } + } else if(year_ == yr2017) { + mvaValues[CBID_LOOSE] = {{0.488, -0.738667, 0.00986667, -0.64}, + {-0.045, -0.825, 0.005, -0.775}, + {0.176, -0.784333, 0.00513333, -0.733}}; + mvaValues[CBID_TIGHT] = {{10, 0.36, 0.032, 0.68}, + {10, 0.225, 0.025, 0.475}, + {10, 0.04, 0.028, 0.32}}; + } else if(year_ == yr2018) { + mvaValues[CBID_LOOSE] = {{1.32, 0.544, 0.066, 1.204}, + {0.192, -0.246, 0.033, 0.084}, + {0.362, -0.653, 0.053, -0.123}}; + mvaValues[CBID_TIGHT] = {{100, 3.157, 0.112, 4.277}, + {100, 2.552, 0.06, 3.152}, + {100, 1.489, 0.087, 2.359}}; } } @@ -194,30 +232,16 @@ void ThreeLepSelector::clearValues() { goodLeptons.clear(); looseLeptons.clear(); goodJets.clear(); + jetList.clear(); + bjetList.clear(); } void ThreeLepSelector::LoadBranchesNanoAOD(Long64_t entry, std::pair variation) { clearValues(); - + fReader.SetLocalEntry(entry); + b.SetEntry(entry); - if (nElectron > N_KEEP_MU_E_ || nMuon > N_KEEP_MU_E_ || nGenPart > N_KEEP_GEN_) { - std::string message = "Found more electrons or muons than max read number.\n Found "; - message += std::to_string(nElectron); - message += " electrons.\n Found "; - message += std::to_string(nMuon); - message += " Muons\n --> Max read number was "; - message += std::to_string(N_KEEP_MU_E_); - message += "\nExiting because this can cause problems. Increase N_KEEP_MU_E_ to avoid this error.\n"; - message += std::to_string(nGenPart); - message += " Gens\n --> Max read number was "; - message += std::to_string(N_KEEP_GEN_); - message += "\nExiting because this can cause problems. Increase N_KEEP_GEN_ to avoid this error.\n"; - throw std::domain_error(message); - } - - if(eventVec[event]) return; - /// basic setups setupElectrons(); setupMuons(); @@ -228,11 +252,9 @@ void ThreeLepSelector::LoadBranchesNanoAOD(Long64_t entry, std::pair= 2 && looseLeptons.size() >= 3) { - for(auto lep : goodLeptons) { - passZVeto = doesPassZVeto(lep, looseLeptons); - if(!passZVeto) break; - } + for(auto lep : goodLeptons) { + passZVeto = doesPassZVeto(lep, looseLeptons); + if(!passZVeto) break; } } @@ -243,7 +265,8 @@ void ThreeLepSelector::setupMuons() { goodLeptons.back().index = i; looseLeptons.push_back(goodLeptons.back()); - if(!passFullIso(goodLeptons.back().v, 0.76, 7.2)) { // Extra Iso requirement + if((year_ == yr2016 && !passFullIso(goodLeptons.back().v, 0.76, 7.2)) || // Extra Iso requirement + ((year_ == yr2017 || year_ == yr2018) && !passFullIso(goodLeptons.back().v, 0.74, 6.8))) { goodLeptons.pop_back(); } } @@ -261,9 +284,10 @@ void ThreeLepSelector::setupElectrons() { goodLeptons.push_back(GoodPart(get4Vector(PID_ELECTRON, i), PID_ELECTRON * Electron_charge[i])); goodLeptons.back().index = i; looseLeptons.push_back(goodLeptons.back()); - - if(!passFullIso(goodLeptons.back().v, 0.8, 7.2)) { // Extra Iso requirement - goodLeptons.pop_back(); + + if((year_ == yr2016 && !passFullIso(goodLeptons.back().v, 0.8, 7.2)) || // Extra Iso requirement + ((year_ == yr2017 || year_ == yr2018) && !passFullIso(goodLeptons.back().v, 0.78, 8.0))) { + goodLeptons.pop_back(); } } else if(isLooseElectron(i)) { @@ -280,7 +304,6 @@ void ThreeLepSelector::setupJets() { } for(size_t i = 0; i < nJet; ++i) { - // if(goodLeptons.size() < 2) break; // only try to find jets if have leptons if(std::find(closeJet.begin(), closeJet.end(), i) != closeJet.end()) continue; /// jet bool passedGoodJet = false; @@ -291,6 +314,7 @@ void ThreeLepSelector::setupJets() { goodJets.back().index = i; nJets++; HT += Jet_pt[i]; + jetList.push_back(goodJets.size() - 1); } // bjet if(isGoodBJet(i)) { @@ -300,41 +324,39 @@ void ThreeLepSelector::setupJets() { } goodJets.back().setPassBJetSel(); nBJets++; + bjetList.push_back(goodJets.size() - 1); } } } void ThreeLepSelector::setupChannel() { - if(goodLeptons.size() > 2) { - channelName_ = "lll"; - if(goodLeptons[1].Charge() * goodLeptons[2].Charge() > 0) { - std::swap(goodLeptons[0], goodLeptons[2]); - } - else if(goodLeptons[0].Charge() * goodLeptons[2].Charge() > 0) { - std::swap(goodLeptons[1], goodLeptons[2]); - } - /// PT swap + if(goodLeptons.size() == 0) { + channelName_ = "Unknown"; + } else if(goodLeptons.size() == 1) { + channelName_ = "one"; + } else if(goodLeptons.size() == 2) { if(goodLeptons[0].Pt() < goodLeptons[1].Pt()) { std::swap(goodLeptons[0], goodLeptons[1]); } - } - else if(goodLeptons.size() < 2) { - channelName_ = "Unknown"; - } - else if(goodLeptons[0].Id() == PID_MUON && goodLeptons[1].Id() == PID_MUON) { - channelName_ = "mm"; - } - else if(goodLeptons[0].Id() == PID_ELECTRON && goodLeptons[1].Id() == PID_ELECTRON) { - channelName_ = "ee"; - } - else { - channelName_ = "em"; + if(goodLeptons[0].Charge()*goodLeptons[1].Charge() > 0) { + channelName_ = "SS"; + } else { + channelName_ = "OS"; + } + } else { + channelName_ = "mult"; if(goodLeptons[0].Pt() < goodLeptons[1].Pt()) { std::swap(goodLeptons[0], goodLeptons[1]); } + if(goodLeptons[0].Pt() < goodLeptons[2].Pt()) { + std::swap(goodLeptons[0], goodLeptons[2]); + } + if(goodLeptons[1].Pt() < goodLeptons[2].Pt()) { + std::swap(goodLeptons[1], goodLeptons[2]); + } } - + channel_ = channelMap_[channelName_]; } @@ -360,7 +382,6 @@ bool ThreeLepSelector::doesPassZVeto(GoodPart& lep, std::vector& loose } void ThreeLepSelector::ApplyScaleFactors() { - if(selection_ == BEfficiency) return; // weight *= (genWeight > 0) ? 1 : -1; weight *= genWeight; @@ -399,7 +420,7 @@ bool ThreeLepSelector::passFakeableCuts(GoodPart& lep) { return (Muon_mediumId[index] && Muon_sip3d[index] < 4 && Muon_tightCharge[index] == 2 - && passFullIso(lep.v, 0.72, 7.2) + //&& passFullIso(lep.v, 0.72, 7.2) ); } @@ -407,7 +428,7 @@ bool ThreeLepSelector::passFakeableCuts(GoodPart& lep) { return (Electron_sip3d[index] < 4 && Electron_tightCharge[index] == 2 && Electron_lostHits[index] == 0 - && passFullIso(lep.v, 0.8, 7.2) + //&& passFullIso(lep.v, 0.8, 7.2) ); } } @@ -417,8 +438,12 @@ bool ThreeLepSelector::isGoodMuon(size_t index) { bool yearCuts = true; if(year_ == yr2016) yearCuts = (Muon_miniPFRelIso_all[index] < 0.16); else yearCuts = (Muon_miniPFRelIso_all[index] < 0.11); + + double ptCut = 20; + if(selection_ == MVAStudy) ptCut = 15; + if(selection_ == FakeRate) ptCut = 10; - return ( (Muon_pt[index] > 20) + return ( (Muon_pt[index] > ptCut) && (Muon_tightCharge[index] == 2) && (abs(Muon_eta[index]) < 2.4) && (Muon_mediumId[index]) @@ -434,7 +459,7 @@ bool ThreeLepSelector::passMVACut(std::vector > mvaCuts, int //// PT Splitting if(Electron_pt[index]/CHGPT(index) < 5) return false; else if(Electron_pt[index]/CHGPT(index) < 10) caseIndex += 0; - else if(Electron_pt[index]/CHGPT(index) < 15) caseIndex += 1; + else if(Electron_pt[index]/CHGPT(index) < 15 && year_ == yr2016) caseIndex += 1; else if(Electron_pt[index]/CHGPT(index) < 25) caseIndex += 2; else caseIndex += 3; //// ETA Splitting @@ -442,8 +467,11 @@ bool ThreeLepSelector::passMVACut(std::vector > mvaCuts, int else if(abs(Electron_eta[index]) < 1.479) caseIndex += 4; else if(abs(Electron_eta[index]) < 2.5) caseIndex += 8; - if(caseIndex % 4 != 2) return Electron_MVA[index] > mvaCuts[caseIndex/4][caseIndex%4]; - else return Electron_MVA[index] > mvaInterpolate(Electron_pt[index]/CHGPT(index), mvaCuts[caseIndex/4]); + double mvaValue = Electron_MVA[index]; + if(year_ == yr2018) mvaValue = atanh(Electron_MVA[index]); + + if(caseIndex % 4 != 2) return mvaValue > mvaCuts[caseIndex/4][caseIndex%4]; + else return mvaValue > mvaInterpolate(Electron_pt[index]/CHGPT(index), mvaCuts[caseIndex/4]); } double ThreeLepSelector::mvaInterpolate(double pt, std::vector cuts) { @@ -454,24 +482,25 @@ double ThreeLepSelector::mvaInterpolate(double pt, std::vector cuts) { bool ThreeLepSelector::isGoodElectron(size_t index) { if(abs(Electron_eta[index]) > 2.5) return false; bool passId = false; + double ptCut = 20; + if(selection_ == MVAStudy) ptCut = 15; + if(selection_ == FakeRate) ptCut = 10; if(selection_ == FourTopMVAEl || selection_ != FourTopCutBasedEl) { if(year_ == yr2016 || year_ == yrdefault) { passId = passMVACut(mvaValues[CBID_TIGHT], index); passId = passId && (Electron_miniPFRelIso_all[index] < 0.12); } - else if(year_ == yr2017) { + else if(year_ == yr2017 || year_ == yr2018) { ///// NEED to fix mva values for 2017 passId = passMVACut(mvaValues[CBID_TIGHT], index); passId = passId && (Electron_miniPFRelIso_all[index] < 0.07); } } else { - passId = (Electron_cutBased[index] >= CBID_MEDIUM); - // if(year_ == yr2016) passId = passId && (Electron_miniPFRelIso_all[index] < 0.12); - // else if(year_ == yr2017) passId = passId && (Electron_miniPFRelIso_all[index] < 0.07); + passId = (Electron_cutBased[index] >= CBID_LOOSE); } - return ((Electron_pt[index]/CHGPT(index) > 20) + return ((Electron_pt[index]/CHGPT(index) > ptCut) && (passId) && (Electron_convVeto[index]) && (Electron_lostHits[index] == 0) @@ -504,7 +533,7 @@ bool ThreeLepSelector::isLooseElectron(size_t index) { passId = passMVACut(mvaValues[CBID_LOOSE], index); } else { - passId = (Electron_cutBased[index] >= CBID_LOOSE); + passId = (Electron_cutBased[index] >= CBID_VETO); } return ((passId)// && (Electron_pt[index]/CHGPT(index) > 7) @@ -520,11 +549,17 @@ bool ThreeLepSelector::isLooseElectron(size_t index) { bool ThreeLepSelector::isGoodJet(size_t index) { bool yearCut = true; double ptCut = 40; + double etaCut = 2.4; + + if(selection_ == MVAStudy) { + ptCut = 25; + etaCut = 4.0; + } + + if(year_ == yr2016) yearCut = IntBits(Jet_jetId[index]).test(0) || IntBits(Jet_jetId[index]).test(1); - if(year_ == yr2016) yearCut = IntBits(Jet_jetId[index]).test(0); - else yearCut = IntBits(Jet_jetId[index]).test(1); return ((Jet_pt[index] > ptCut) && - (abs(Jet_eta[index]) < 2.4) && + (abs(Jet_eta[index]) < etaCut) && (yearCut) ); } @@ -534,12 +569,20 @@ bool ThreeLepSelector::isGoodJet(size_t index) { bool ThreeLepSelector::isGoodBJet(size_t index) { bool yearCut = true; double ptCut = 25; + double etaCut = 2.4; + + if(selection_ == MVAStudy) { + ptCut = 25; + etaCut = 4.0; + } - if(year_ == yr2016) yearCut = (IntBits(Jet_jetId[index]).test(0)) && (Jet_btagDeepB[index] > 0.6324); - else yearCut = (IntBits(Jet_jetId[index]).test(1)) && (Jet_btagDeepB[index] > 0.4941); + yearCut = IntBits(Jet_jetId[index]).test(0) || IntBits(Jet_jetId[index]).test(1); + if(year_ == yr2016) yearCut = yearCut && (Jet_btagDeepB[index] > 0.6324); + else if (year_ == yr2017) yearCut = yearCut && (Jet_btagDeepB[index] > 0.4941); + else if (year_ == yr2018) yearCut = yearCut && (Jet_btagDeepB[index] > 0.4184); return ((Jet_pt[index] > ptCut) - && (abs(Jet_eta[index]) < 2.4) + && (abs(Jet_eta[index]) < etaCut) && (yearCut) ); } @@ -578,14 +621,14 @@ double ThreeLepSelector::LepRelPt(LorentzVector& lep, LorentzVector& closeJet) { // Need to include DeltaPhi Ieta bool ThreeLepSelector::passTriggerEmu(size_t index) { bool etaDepend = false; - int DEtaInCut = READFROM(Electron_vidBitmap[index], 6, 3); - int DPhiInCut = READFROM(Electron_vidBitmap[index], 9, 3); + // int DEtaInCut = READFROM(Electron_vidBitmap[index], 6, 3); + // int DPhiInCut = READFROM(Electron_vidBitmap[index], 9, 3); if(abs(Electron_eta[index]) < 1.479) { - etaDepend = Electron_sieie[index] < 0.011 && (DEtaInCut >= 1) && (DPhiInCut >= 4); + etaDepend = Electron_sieie[index] < 0.011;// && (DEtaInCut >= 1) && (DPhiInCut >= 4); } else { - etaDepend = Electron_sieie[index] < 0.031 && (DEtaInCut >= 1) && (DPhiInCut >= 3); + etaDepend = Electron_sieie[index] < 0.031;// && (DEtaInCut >= 1) && (DPhiInCut >= 3); } return (abs(Electron_eInvMinusPInv[index]) < 0.01 && etaDepend && @@ -617,6 +660,7 @@ void ThreeLepSelector::FillHistograms(Long64_t entry, std::pair 0) - ) - return; - Fill1D("CutFlow", ++step); + if(selection_ != MVAStudy) { + /// 2 good leptons + if(goodLeptons.size() != 2) return; + Fill1D("CutFlow", ++step); + + // first lep requirement + if(goodLeptons[0].Pt() < 25) return; + Fill1D("CutFlow", ++step); + + // same sign requirement + if((goodLeptons.size() == 2 && goodLeptons[0].Charge() * goodLeptons[1].Charge() < 0) || + (goodLeptons.size() == 3 && goodLeptons[0].Charge() * goodLeptons[2].Charge() > 0)) + return; + Fill1D("CutFlow", ++step); + } + + if(channel_ == mult) { + if((goodLeptons[0].Charge() * goodLeptons[1].Charge() > 0) && + (goodLeptons[0].Charge() * goodLeptons[2].Charge() > 0)) return; + } + // jet cut if(nJets < 2) return; Fill1D("CutFlow", ++step); // bjet cut - if(nBJets < 2) return; + if(selection_ == MVAStudy && nBJets < 1) return; + else if(selection_ != MVAStudy && nBJets < 2) return; Fill1D("CutFlow", ++step); // ht cut - if(HT < 300) return; + if(selection_ == MVAStudy && HT < 100) return; + else if(selection_ != MVAStudy && HT < 300) return; Fill1D("CutFlow", ++step); // met cut - if (MET < 50) return; + if(selection_ == MVAStudy && MET < 25) return; + else if(selection_ != MVAStudy && MET < 50) return; Fill1D("CutFlow", ++step); Fill1D("SR", getSRBin()); + if(!passZVeto) return; + Fill1D("CutFlow", ++step); + // if(getSRBin() == -1) { + // return; + // } + // else if(getSRBin() == 0) { + // Fill1D("CRW_Met", MET); + // Fill1D("CRW_HT", HT); + // Fill1D("CRW_njet", nJets); + // Fill1D("CRW_nbjet", nBJets); + // return; + // } + // else if(getSRBin() == 9) { + // Fill1D("CRZ_Met", MET); + // Fill1D("CRZ_HT", HT); + // Fill1D("CRZ_njet", nJets); + // Fill1D("CRZ_nbjet", nBJets); + // Fill1D("CRZ_ptl3", goodLeptons[3].Pt()); + // return; + // } - if(getSRBin() == -1) { - return; - } - else if(getSRBin() == 0) { - Fill1D("CRW_Met", MET); - Fill1D("CRW_HT", HT); - Fill1D("CRW_njet", nJets); - Fill1D("CRW_nbjet", nBJets); - return; - } - else if(getSRBin() == 9) { - Fill1D("CRZ_Met", MET); - Fill1D("CRZ_HT", HT); - Fill1D("CRZ_njet", nJets); - Fill1D("CRZ_nbjet", nBJets); - Fill1D("CRZ_ptl3", goodLeptons[3].Pt()); - return; + HistFullFill(histMap1D_, "weight", variation.first, abs(weight), 1); + + int NlooseBs = 0; + int NtightBs = 0; + double lbjetCut = 0; + double tbjetCut = 0; + if(year_ == yr2016) { + lbjetCut = 0.2219; + tbjetCut = 0.8958; + } else if (year_ == yr2017) { + lbjetCut = 0.1522; + tbjetCut = 0.8001; + } else if (year_ == yr2018) { + lbjetCut = 0.1241; + tbjetCut = 0.7527; + } + for(size_t i=0; i 20) + && (IntBits(Jet_jetId[i]).test(0) || IntBits(Jet_jetId[i]).test(1)) + && Jet_btagDeepB[i] > lbjetCut ) + NlooseBs++; + if((Jet_pt[i] > 20) + && (IntBits(Jet_jetId[i]).test(0) || IntBits(Jet_jetId[i]).test(1)) + && Jet_btagDeepB[i] > tbjetCut ) + NtightBs++; } - Fill1D("CutFlow", ++step); - - HistFullFill(histMap1D_, "weight", variation.first, abs(weight), 1); - Fill1D("Met", MET); - Fill1D("HT", HT); - Fill1D("ptl1", goodLeptons[0].Pt()); - Fill1D("ptl2", goodLeptons[1].Pt()); Fill1D("njet", nJets); Fill1D("nbjet", nBJets); + Fill1D("nloosebjet", NlooseBs); + Fill1D("ntightbjet", NtightBs); Fill2D("bJetvsJets", nJets, nBJets); Fill1D("nleps", goodLeptons.size()); - + Fill1D("nlooseleps", looseLeptons.size()); + Fill1D("nloosevstightleps", looseLeptons.size() - goodLeptons.size()); + Fill1D("dilepCharge", goodLeptons[0].Charge() * goodLeptons[1].Charge() > 0 ? 1 : -1); + Fill1D("HT", HT); + Fill1D("Met", MET); + Fill1D("ptl1", goodLeptons[0].Pt()); + Fill1D("ptl2", goodLeptons[1].Pt()); + Fill1D("dilepMass", (goodLeptons[0].v+goodLeptons[1].v).M()); Fill1D("sphericity", JetSphericity(goodJets)); Fill1D("centrality", JetCentrality(goodJets,HT)); + Fill1D("DRLep", reco::deltaR(goodLeptons[0].v, goodLeptons[1].v)); + Fill1D("DRjet", reco::deltaR(goodJets[0].v, goodJets[1].v)); + Fill1D("dijetMass", (goodJets[0].v, goodJets[1].v).M()); + Fill1D("DRLep", reco::deltaR(goodLeptons[0].v, goodLeptons[1].v)); + + auto event_pair = EventShape(goodJets, goodLeptons, pow(MET, 2), type1_pfMETPhi); + Fill1D("Shape1", event_pair.first); + Fill1D("Shape2", event_pair.second); + Fill1D("LepCos", ROOT::Math::VectorUtil::CosTheta(goodLeptons[0].v, goodLeptons[1].v)); + Fill1D("JLep1Cos", ROOT::Math::VectorUtil::CosTheta(goodLeptons[0].v, goodJets[0].v)); + Fill1D("JLep2Cos", ROOT::Math::VectorUtil::CosTheta(goodLeptons[1].v, goodJets[0].v)); + int goodjet1 = 0; + if(bjetList.at(0) == 0) goodjet1 = 1; + Fill1D("JBCos", ROOT::Math::VectorUtil::CosTheta(goodJets[bjetList.at(0)].v, goodJets[goodjet1].v)); + Fill1D("DRjb", reco::deltaR(goodJets[bjetList.at(0)].v, goodJets[goodjet1].v)); + + for(auto i : jetList) { + Fill1D("jetpt", goodJets[i].Pt()); + Fill1D("etaj", goodJets[i].Eta()); + } + for(auto i : bjetList) { + Fill1D("bjetpt", goodJets[i].Pt()); + Fill1D("etab", goodJets[i].Eta()); + } + int k= 1; + for(auto it: jetList) { + if(k > 3) break; + std::string intStr = std::to_string(k); + LorentzVector jit = goodJets.at(it).v; + Fill1D(("ptj" + intStr).c_str(), jit.Pt()); + Fill1D(("etaj" + intStr).c_str(), jit.Eta()); + Fill1D(("dphi_l1j" + intStr).c_str(), abs(ROOT::Math::VectorUtil::DeltaPhi(goodLeptons[0].v, jit))); + k++; + } + k=1; + for(auto it: bjetList) { + if(k > 3) break; + std::string intStr = std::to_string(k); + LorentzVector jit = goodJets.at(it).v; + Fill1D(("ptb" + intStr).c_str(), jit.Pt()); + Fill1D(("etab" + intStr).c_str(), jit.Eta()); + k++; + } +#ifdef USETREE + bNJets = nJets; + bnBJets = nBJets; + bnlBJets =NlooseBs; + bntBJets =NtightBs; + bnLeps = goodLeptons.size(); + bnlLeps = looseLeptons.size(); + bDilepCharge = goodLeptons[0].Charge() * goodLeptons[1].Charge() > 0 ? 1 : -1; + bHT = HT; + bMET = MET; + bl1Pt = goodLeptons[0].Pt(); + bl2Pt = goodLeptons[1].Pt(); + blMass = (goodLeptons[0].v+goodLeptons[1].v).M(); + bjMass = (goodJets[0].v+goodJets[1].v).M(); + bjdr = reco::deltaR(goodJets[0].v,goodJets[1].v); + bsphere = JetSphericity(goodJets); + bCentral = JetCentrality(goodJets,HT); + bj1Pt = goodJets.at(jetList[0]).Pt(); + bj2Pt = (jetList.size() > 1) ? goodJets.at(jetList[1]).Pt() : 0; + bj3Pt = (jetList.size() > 2) ? goodJets.at(jetList[2]).Pt() : 0; + bj4Pt = (jetList.size() > 3) ? goodJets.at(jetList[3]).Pt() : 0; + bj5Pt = (jetList.size() > 4) ? goodJets.at(jetList[4]).Pt() : 0; + bj6Pt = (jetList.size() > 5) ? goodJets.at(jetList[5]).Pt() : 0; + bj7Pt = (jetList.size() > 6) ? goodJets.at(jetList[6]).Pt() : 0; + bj8Pt = (jetList.size() > 7) ? goodJets.at(jetList[7]).Pt() : 0; + bb1Pt = goodJets.at(bjetList[0]).Pt(); + bb2Pt = (bjetList.size() > 1) ? goodJets.at(bjetList[1]).Pt() : 0; + bb3Pt = (bjetList.size() > 2) ? goodJets.at(bjetList[2]).Pt() : 0; + bb4Pt = (bjetList.size() > 3) ? goodJets.at(bjetList[3]).Pt() : 0; + bShape1 = event_pair.first; + bShape2 = event_pair.second; + treeMap["tree"]->Fill(); +#endif + +} + +std::vector::iterator ThreeLepSelector::findJet(std::vector::iterator& start, int pid) { + while(start != goodJets.end()) { + if(pid == PID_BJET && start->passedBJetSel()) return start; + else if(pid != PID_BJET && start->passedJetSel()) return start; + ++start; + } + return goodJets.end(); } void ThreeLepSelector::SetupNewDirectory() { diff --git a/src/WZSelector.cc b/src/WZSelector.cc index a5e99ec5..3f212107 100644 --- a/src/WZSelector.cc +++ b/src/WZSelector.cc @@ -8,9 +8,15 @@ void WZSelector::Init(TTree *tree) weight_info_ = GetLheWeightInfo(); } + addSelection_ = {{"Inclusive2Jet_Full", Inclusive2Jet_Full}, + {"Wselection_Full", Wselection_Full}, + {"Inclusive2Jet", Inclusive2Jet}, + {"FakeRateSelectionTight", FakeRateSelectionTight}, + {"FakeRateSelectionLoose", FakeRateSelectionLoose}, }; + systematics_ = { - {jetEnergyScaleUp, "CMS_scale_jUp"}, - {jetEnergyScaleDown, "CMS_scale_jDown"}, + {jetEnergyScaleUp, "CMS_scale_jUp"}, + {jetEnergyScaleDown, "CMS_scale_jDown"}, {jetEnergyResolutionUp, "CMS_res_jUp"}, {jetEnergyResolutionDown, "CMS_res_jDown"}, {metUnclusteredEnergyUp, "CMS_scale_unclEnergyUp"}, @@ -505,7 +511,7 @@ bool WZSelector::PassesFullWZSelection(Long64_t entry) { return true; } -bool WZSelector::PassesBaseSelection(Long64_t entry, bool tightLeps, Selection selection) { +bool WZSelector::PassesBaseSelection(Long64_t entry, bool tightLeps, int selection) { //if (!(Flag_BadChargedCandidateFilterPass // && Flag_HBHENoiseFilterPass // && Flag_HBHENoiseIsoFilterPass diff --git a/src/WZSelectorBase.cc b/src/WZSelectorBase.cc index ecc94f85..b1fa11ee 100644 --- a/src/WZSelectorBase.cc +++ b/src/WZSelectorBase.cc @@ -53,9 +53,22 @@ void WZSelectorBase::Init(TTree *tree) allChannels_ = {{eee, "eee"}, {eem, "eem"}, {emm, "emm"}, {mmm, "mmm"}}; + selectionMap_ = {{"VBSBackgroundControlLoose_Full", VBSBackgroundControlLoose_Full}, + {"VBSBackgroundControl_Full", VBSBackgroundControl_Full}, + {"VBSBackgroundControl", VBSBackgroundControl}, + {"VBSBackgroundControlLoose", VBSBackgroundControlLoose}, + {"VBSBackgroundControlATLAS", VBSBackgroundControlATLAS}, + {"VBSselection_Loose_Full", VBSselection_Loose_Full}, + {"VBSselection_Tight_Full", VBSselection_Tight_Full}, + {"VBSselection_Tight", VBSselection_Tight}, + {"VBSselection_Loose", VBSselection_Loose}, + {"VBSselection_NoZeppenfeld", VBSselection_NoZeppenfeld}, + {"VBSselection_NoZeppenfeld_Full", VBSselection_NoZeppenfeld_Full}, + }; + if (isMC_){ - isNonpromptMC_ = false; - isZgamma_ = false; + isNonpromptMC_ = false; + isZgamma_ = false; if (std::find(nonprompt3l_.begin(), nonprompt3l_.end(), name_) != nonprompt3l_.end()) { isNonpromptMC_ = true; } diff --git a/src/classes.h b/src/classes.h index a36d1fb9..5dc4c5a6 100644 --- a/src/classes.h +++ b/src/classes.h @@ -10,7 +10,6 @@ #include "Analysis/VVAnalysis/interface/ZGenSelector.h" #include "Analysis/VVAnalysis/interface/NanoGenSelectorBase.h" #include "Analysis/VVAnalysis/interface/WZSelector.h" -#include "Analysis/VVAnalysis/interface/TTTSelector.h" #include "Analysis/VVAnalysis/interface/ThreeLepSelector.h" #include "Analysis/VVAnalysis/interface/WZSelectorBase.h" #include "Analysis/VVAnalysis/interface/WZBackgroundSelector.h" @@ -33,7 +32,6 @@ namespace{ WZSelector pWZSelector; NanoGenSelectorBase pNanoGenSelectorBase; ZZGenSelector pZZGenSelector; - TTTSelector pTTTSelector; ThreeLepSelector pThreeLepSelector; WGenSelector pWGenSelector; ZGenSelector pZGenSelector; diff --git a/src/classes_def.xml b/src/classes_def.xml index 4bfa0012..c7b463ad 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -1,6 +1,5 @@ -