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 @@ -