From 1126ffbdad75715490352dff98a9b64f4d4ab9cc Mon Sep 17 00:00:00 2001 From: Alexandre Mertens Date: Mon, 29 Feb 2016 18:03:26 +0100 Subject: [PATCH 1/3] Adding MC stat uncertainty and better brazilian band plots --- python/classes.py | 7 +-- python/createCards.py | 113 +++++++++++++++++++++++++++++++++--------- python/makeVBBP.py | 33 +++++++----- 3 files changed, 113 insertions(+), 40 deletions(-) diff --git a/python/classes.py b/python/classes.py index 53f091e..0855653 100644 --- a/python/classes.py +++ b/python/classes.py @@ -38,7 +38,7 @@ def get_sample(self, name, tag): return resultset.one() def prepare_process(self, path, shortname, name, tag): - #sample = self.get_sample(name, tag) + sample = self.get_sample(name, tag) self.name = shortname if 'ZZ' in name: self.type = -1 @@ -54,8 +54,9 @@ def prepare_process(self, path, shortname, name, tag): print self.file os.path.isfile(self.file) print self.file - self.xsection = 1 #sample.source_dataset.xsection - self.sumW = 1 #sample.event_weight_sum + if sample : + self.xsection = sample.source_dataset.xsection + self.sumW = sample.event_weight_sum self.channels = {} return self diff --git a/python/createCards.py b/python/createCards.py index 8b9e729..ae0f0dc 100644 --- a/python/createCards.py +++ b/python/createCards.py @@ -36,7 +36,7 @@ def main(): intL = options.lumi # in pb-1 #tag = 'v1.2.0+7415-19-g7bbca78_ZAAnalysis_1a69757' #path = '/nfs/scratch/fynu/amertens/cmssw/CMSSW_7_4_15/src/cp3_llbb/CommonTools/histFactory/16_01_28_syst/build' - tag = 'v1.1.0+7415-83-g2a9f912_ZAAnalysis_2ff9261' + tag = 'v1.2.0+7415-73-g91bfcf0_ZAAnalysis_4f7ac83' #tag = 'v1.1.0+7415-57-g4bff5ea_ZAAnalysis_b1377a8' path = options.path CHANNEL = options.CHANNEL @@ -47,6 +47,7 @@ def main(): c = ch.CombineHarvester() cats = [(0, "mmbbSR"+cutkey), + #(1, "eebbSR"+cutkey) (1, "mll_mmbbBR"+cutkey), (2, "eebbSR"+cutkey), (3, "mll_eebbBR"+cutkey) @@ -61,32 +62,38 @@ def main(): processes = {} p = Process('data_obs') #DoubleMuon_Run2015D_v1.1.0+7415-57-g4bff5ea_ZAAnalysis_b1377a8_histos.root - p.prepare_process(path, 'data_obs', 'DoubleMuon_DoubleEG_Run2015D', tag) + p.prepare_process(path, 'data_obs', 'DoubleMuon_DoubleEG_Run2015D_v2_2015-12-18', tag) processes['data_obs'] = p if DEBUG: print p # define signal # define backgrounds # zz + + p = Process('zh') + p.prepare_process(path, 'zh', 'ZH_HToBB_ZToLL_M125_13TeV_powheg_pythia8', tag) + processes['zh'] = p + if DEBUG: print p + p = Process('zz') - p.prepare_process(path, 'zz', 'ZZTo2L2Q_13TeV_amcatnloFXFX_madspin_pythia8_MiniAODv2', tag) + p.prepare_process(path, 'zz', 'ZZTo2L2Q_13TeV_amcatnloFXFX_madspin_pythia8_Fall15MiniAODv2', tag) processes['zz'] = p if DEBUG: print p # ttbar p = Process('ttbar') - p.prepare_process(path, 'ttbar', 'TTTo2L2Nu_13TeV-powheg_MiniAODv2', tag) + p.prepare_process(path, 'ttbar', 'TTTo2L2Nu_13TeV-powheg_Fall15MiniAODv2', tag) processes['ttbar'] = p p = Process('ttbar') if DEBUG: print p ''' # drell-yan p = Process('dy1') - p.prepare_process(path, 'dy1', 'DYJetsToLL_M-10to50_TuneCUETP8M1_13TeV-amcatnloFXFX_MiniAODv2', tag) + p.prepare_process(path, 'dy1', 'DYJetsToLL_M-10to50_TuneCUETP8M1_13TeV-amcatnloFXFX_Fall15MiniAODv2', tag) processes['dy1'] = p if DEBUG: print p ''' p = Process('dy2') - p.prepare_process(path, 'dy2', 'DYJetsToLL_M-50_TuneCUETP8M1_13TeV-amcatnloFXFX_MiniAODv2', tag) + p.prepare_process(path, 'dy2', 'DYJetsToLL_M-50_TuneCUETP8M1_13TeV-amcatnloFXFX_Fall15MiniAODv2', tag) processes['dy2'] = p if DEBUG: print p @@ -94,25 +101,25 @@ def main(): c.AddObservations([MASS], [ANALYSIS], [ERA], [CHANNEL], cats) c.AddProcesses([MASS], [ANALYSIS], [ERA], [CHANNEL], ['ZA'], cats, True) - c.AddProcesses([MASS], [ANALYSIS], [ERA], [CHANNEL], ['ttbar', 'dy2','zz'], cats, False) - c.cp().process(['ttbar', 'dy2','ZA']).AddSyst( - c, "lumi", "lnN", ch.SystMap('channel', 'era', 'bin_id') - ([CHANNEL], [ERA], [0,1,2,3], 1.046)) + c.AddProcesses([MASS], [ANALYSIS], [ERA], [CHANNEL], ['ttbar', 'dy2','zz','zh'], cats, False) + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "lumi_13TeV2015", "lnN", ch.SystMap('channel', 'era', 'bin_id') + ([CHANNEL], [ERA], [0,1,2,3], 1.027)) - c.cp().process(['ttbar', 'dy2','ZA']).AddSyst( + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( c, "trig", "lnN", ch.SystMap('channel', 'era', 'bin_id') ([CHANNEL], [ERA], [0,1,2,3], 1.04)) - c.cp().process(['ttbar', 'dy2']).AddSyst( + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( c, "btag", "shape", ch.SystMap()(1.0)) - c.cp().process(['ttbar', 'dy2']).AddSyst( + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( c, "jec", "shape", ch.SystMap()(1.0)) - c.cp().process(['ttbar', 'dy2']).AddSyst( + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( c, "jer", "shape", ch.SystMap()(1.0)) - c.cp().process(['ttbar', 'dy2']).AddSyst( + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( c, "pu", "shape", ch.SystMap()(1.0)) c.cp().process(['ttbar']).AddSyst( @@ -123,11 +130,11 @@ def main(): c.cp().process(['dy2']).AddSyst( c, "DYnorm", "lnN", ch.SystMap('channel', 'era', 'bin_id') - ([CHANNEL], [ERA], [0,1], 1.1)) + ([CHANNEL], [ERA], [0,1], 1.04)) c.cp().process(['ttbar']).AddSyst( c, "TTnorm", "lnN", ch.SystMap('channel', 'era', 'bin_id') - ([CHANNEL], [ERA], [0], 1.1)) + ([CHANNEL], [ERA], [0], 1.08)) nChannels = len(bins) @@ -143,11 +150,11 @@ def main(): f = TFile(outputRoot, "recreate") f.Close() for b in bins: - print b , bins[b] + #print b , bins[b] for p in processes: if p == 'data_obs' : file_in = TFile(processes[p].file,"READ") - print " Getting ", bins[b], " in file ", processes[p].file + #print " Getting ", bins[b], " in file ", processes[p].file h = file_in.Get(bins[b]) h.SetDirectory(0) file_in.Close() @@ -161,23 +168,75 @@ def main(): else : for s1,s2 in systematics.iteritems() : file_in = TFile(processes[p].file,"READ") - print " Getting ", bins[b]+s2, " in file ", processes[p].file + #if s1 == '' : + # print " Getting ", bins[b]+s2, " in file ", processes[p].file h = file_in.Get(bins[b]+s2) + #if s1 == '' : + # print "Val: ", h.GetBinContent(1), "ERROR : ", h.GetBinError(1), 'Entries : ', h.GetEntries() + ''' + if s1 == '' and 'mmbbSR' in bins[b] : + print "Adding MC stat uncertainty in ", p, bins[b] , " of ", + if h.GetBinContent(1) > 0.0 : + MCStatUnc = 1+h.GetBinError(1)/h.GetBinContent(1) + else : + MCStatUnc = 2.0 + print h.GetBinError(1), h.GetBinContent(1), h.GetEntries(), MCStatUnc + + c.cp().process([p]).AddSyst( + c, p+bins[b]+"MCStat", "lnN", ch.SystMap('bin_id') + ([0],MCStatUnc)) + if s1 == '' and 'eebbSR' in bins[b] : + print "Adding MC stat uncertainty in ", p, bins[b] , " of ", + MCStatUnc = 1+h.GetBinError(1)/h.GetBinContent(1) + print h.GetBinError(1), h.GetBinContent(1), h.GetEntries(), MCStatUnc + + c.cp().process([p]).AddSyst( + c, p+bins[b]+"MCStat", "lnN", ch.SystMap('bin_id') + ([2],MCStatUnc)) + ''' + h.SetDirectory(0) file_in.Close() f = TFile(outputRoot, "update") h.SetName("hist_"+bins[b]+"_"+p+s1) - h.Sumw2() + Nbins = h.GetNbinsX() + for it_b in range(0,h.GetNbinsX()+1) : + if h.GetBinContent(it_b) < 0.0 : + print "WARNING : negative bin content in bin ", it_b, " of histo ", "hist_"+bins[b]+"_"+p+s1 + h.SetBinContent(it_b,0.0) + #h.Sumw2() #h.Scale(processes[p].xsection * intL / processes[p].sumW) h.Scale(intL) h.Write() f.Write() f.Close() - + for s1,s2 in systematics.iteritems() : + file_sig_path = "/home/fynu/amertens/scratch/cmssw/CMSSW_7_6_3/src/cp3_llbb/CommonTools/histFactory/Plots_v3/condor/output/HToZATo2L2B_MH-"+str(int(mH))+"_MA-"+str(int(mA))+"_13TeV-madgraph_Fall15MiniAODv1_v1.2.0+7415-73-g91bfcf0_ZAAnalysis_4f7ac83_histos.root" + file_sig = TFile(file_sig_path, "READ") + h = file_sig.Get(bins[b]+s2) + h.SetDirectory(0) + file_sig.Close() + f = TFile(outputRoot, "update") + h.SetName("hist_"+bins[b]+"_"+"ZA"+s1) + + if DEBUG and s1 == '' : + #if s1 == '' : + print file_sig_path + print "saving ", "hist_"+bins[b]+"_"+"ZA"+s1 + print h.GetBinContent(1), "after rescaling", + h.Scale(options.lumifb) + #if DEBUG and s1 == '' : + if s1 == '' : + print h.GetBinContent(1) + h.Write() + f.Write() + f.Close() # Fill signal histograms FIXME: read efficiencies from eff.root - + + ''' + eff_file = TFile("eff.root", "READ") effee_hist = eff_file.Get("effee") eff_ee = effee_hist.Interpolate(mA,mH) @@ -203,7 +262,7 @@ def main(): h4 = TH1F("hist_"+bins['mll_bkgregion_ee']+"_ZA","hist_"+bins['mll_bkgregion_ee']+"_ZA",60,60,120) h4.Write() - + ''' f.Write() f.Close() @@ -212,6 +271,12 @@ def main(): outputRoot, "hist_$BIN_$PROCESS", "hist_$BIN_$PROCESS_$SYSTEMATIC") c.cp().signals().ExtractShapes( outputRoot, "hist_$BIN_$PROCESS", "hist_$BIN_$PROCESS_$SYSTEMATIC") + + bbb = ch.BinByBinFactory().SetAddThreshold(0.05).SetFixNorm(False) + + bbb.AddBinByBin(c.cp().backgrounds(),c) + ch.SetStandardBinNames(c) + writer = ch.CardWriter('$TAG/$MASS/$ANALYSIS_$CHANNEL_$ERA.dat', '$TAG/common/$ANALYSIS_$CHANNEL_$MASS.input_$ERA.root') writer.WriteCards('CARDS/', c) diff --git a/python/makeVBBP.py b/python/makeVBBP.py index d8cffca..b4a0c1c 100644 --- a/python/makeVBBP.py +++ b/python/makeVBBP.py @@ -18,21 +18,21 @@ ### Definitions ### ################### -mH=800 +mH=500 -run_combine = 0 +run_combine = 1 TGAS_1 = TGraphAsymmErrors(0) TGAS_2 = TGraphAsymmErrors(0) -myTG_0 = TGraph(5) -myTG_1 = TGraph(5) -myTG_2 = TGraph(5) -myTG_3 = TGraph(5) -myTG_4 = TGraph(5) -myTG_5 = TGraph(5) +myTG_0 = TGraph(2) +myTG_1 = TGraph(2) +myTG_2 = TGraph(2) +myTG_3 = TGraph(2) +myTG_4 = TGraph(2) +myTG_5 = TGraph(2) -myTG_8TeV = TGraph(5) +myTG_8TeV = TGraph(2) myTGraph = TGraph2D(9) #myTGraph.SetNpx(100) @@ -134,7 +134,7 @@ gStyle.SetOptStat(0) Cname = "C_"+str(mH) -C = TCanvas(Cname,Cname,1200,500) +C = TCanvas(Cname,Cname,1200,1200) C.SetLeftMargin(0.2) C.SetBottomMargin(0.2) C.SetTitle(" ") @@ -176,6 +176,8 @@ TGAS_2.GetXaxis().SetTitle("m_{A} (GeV)") TGAS_2.GetYaxis().SetTitle("#sigma (pp #rightarrow H) #times BR(H #rightarrow Z(ll)A(bb)) (fb)") +TGAS_2.SetMinimum(5) +TGAS_2.SetMaximum(1000) #gPad().SetLogy(1) @@ -189,9 +191,14 @@ leg.AddEntry(myTG_8TeV,"8 TeV","L") leg.Draw() -C.Print("BBP"+str(mH)+".pdf") -C.Print("BBP"+str(mH)+".png") -f = TFile("narrow_BBP"+str(mH)+".root","recreate") +gPad.SetLogy(1) + + + + +C.Print("BBP"+str(mH)+"_v1.pdf") +C.Print("BBP"+str(mH)+"_v1.png") +f = TFile("narrow_BBP"+str(mH)+"_v1.root","recreate") C.Write() f.Close() From 6cb5f159784668d210679238a1526b1721d27655 Mon Sep 17 00:00:00 2001 From: Alexandre Mertens Date: Wed, 9 Mar 2016 17:51:22 +0100 Subject: [PATCH 2/3] adding new BB plots --- python/ZACnC.py | 26 +- python/classes.py | 3 - python/createCards.py | 483 ++++++++++++++++----------------- python/createCards_combined.py | 287 ++++++++++++++++++++ python/createXsecLines.py | 53 ++++ python/drawSyst.py | 2 +- python/getPostfitValues.py | 128 +++++++++ python/makeEfficiencymap.py | 4 +- python/makeVBBP.py | 137 ++++++++-- 9 files changed, 843 insertions(+), 280 deletions(-) create mode 100644 python/createCards_combined.py create mode 100644 python/createXsecLines.py create mode 100644 python/getPostfitValues.py diff --git a/python/ZACnC.py b/python/ZACnC.py index 4c9c883..6bf8b46 100644 --- a/python/ZACnC.py +++ b/python/ZACnC.py @@ -4,14 +4,15 @@ class options_(): # parameters - lumi = 2245.792 #in pb + lumi = 2300.0 #in pb lumifb = lumi/1000.0 #in fb - tag = 'v1.1.0+7415-83-g2a9f912_ZAAnalysis_2ff9261' - path = '/home/fynu/amertens/scratch/cmssw/CMSSW_7_4_15/src/cp3_llbb/CommonTools/histFactory/test_narrow/condor/output/' + tag = 'v1.2.0+7415-73-g91bfcf0_ZAAnalysis_4f7ac83' + path = '/home/fynu/amertens/scratch/cmssw/CMSSW_7_6_3/src/cp3_llbb/CommonTools/histFactory/Plots_v5/condor/output/' CHANNEL = 'mumu' ERA = '13TeV' ANALYSIS = 'HtoZAtoLLBB' + ''' ### Define 2D mapping for the search in the M(bb) - M(llbb) plane ### ''' rangeMassA = [] @@ -32,22 +33,27 @@ class options_(): rangeMassH.append([int(mllbb-dmllbb),int(mllbb+dmllbb),int(mllbb)]) mllbb+=step_mllbb ''' - mbb_list = {50,75,100,125,150,200,225,250,300,350,400,500,600,700} + + #mbb_list = {50,75,100,125,150,200,225,250,300,350,400,500,600,700} #mllbb_list = {150,200,250,300,350,400,450,500,550,650,800,1000} - #mbb_list = [250] + #mbb_list = [50,75,100,125,150,200,225,250,300,350,400] #mbb_list = [50, 75, 100] - mllbb_list = {300, 500, 800} + #mllbb_list = {300, 500, 800} + + mbb_list = [50, 100, 200, 300, 400, 700] + #mbb_list = [100] + mllbb_list = [300, 500, 800] - #rangeMassA = [] + rangeMassA = [] rangeMassA = [] rangeMassH = [] - #for mbb in mbb_list : - # dmbb=0.15*mbb*1.5 - # rangeMassA.append([int(mbb-dmbb),int(mbb+dmbb),int(mbb)]) + for mbb in mbb_list : + dmbb=0.15*mbb*1.5 + rangeMassA.append([int(mbb-dmbb),int(mbb+dmbb),int(mbb)]) for mllbb in mllbb_list : dmllbb=0.15*mllbb*1.5 diff --git a/python/classes.py b/python/classes.py index 0855653..885ac15 100644 --- a/python/classes.py +++ b/python/classes.py @@ -48,12 +48,9 @@ def prepare_process(self, path, shortname, name, tag): self.type = 2 if 'data_obs' in shortname: self.type = 0 - print name + "_" + tag + "_histos.root" self.file = os.path.join(path, name + "_" + tag + "_histos.root") - print self.file os.path.isfile(self.file) - print self.file if sample : self.xsection = sample.source_dataset.xsection self.sumW = sample.event_weight_sum diff --git a/python/createCards.py b/python/createCards.py index ae0f0dc..9baee48 100644 --- a/python/createCards.py +++ b/python/createCards.py @@ -6,6 +6,7 @@ # User imports from classes import * import CombineHarvester.CombineTools.ch as ch +import os # ROOT imports import ROOT from ROOT import gROOT @@ -26,260 +27,248 @@ def main(): mA = float(options.mA_list[cutkey]) print mH, mA - """Main function""" - # start the timer - tstart = datetime.now() - print 'starting...' - # get the options - #options = get_options() - - intL = options.lumi # in pb-1 - #tag = 'v1.2.0+7415-19-g7bbca78_ZAAnalysis_1a69757' - #path = '/nfs/scratch/fynu/amertens/cmssw/CMSSW_7_4_15/src/cp3_llbb/CommonTools/histFactory/16_01_28_syst/build' - tag = 'v1.2.0+7415-73-g91bfcf0_ZAAnalysis_4f7ac83' - #tag = 'v1.1.0+7415-57-g4bff5ea_ZAAnalysis_b1377a8' - path = options.path - CHANNEL = options.CHANNEL - ERA = options.ERA - MASS = str(mH)+"_"+str(mA) - ANALYSIS = options.ANALYSIS - DEBUG = 0 - - c = ch.CombineHarvester() - cats = [(0, "mmbbSR"+cutkey), - #(1, "eebbSR"+cutkey) - (1, "mll_mmbbBR"+cutkey), - (2, "eebbSR"+cutkey), - (3, "mll_eebbBR"+cutkey) - ] - - bins = {} - bins['signalregion_mm'] = "mmbbSR"+cutkey - bins['mll_bkgregion_mm'] = "mll_mmbbBR"+cutkey - bins['signalregion_ee'] = "eebbSR"+cutkey - bins['mll_bkgregion_ee'] = "mll_eebbBR"+cutkey - - processes = {} - p = Process('data_obs') - #DoubleMuon_Run2015D_v1.1.0+7415-57-g4bff5ea_ZAAnalysis_b1377a8_histos.root - p.prepare_process(path, 'data_obs', 'DoubleMuon_DoubleEG_Run2015D_v2_2015-12-18', tag) - processes['data_obs'] = p - if DEBUG: print p - # define signal - # define backgrounds - # zz - - p = Process('zh') - p.prepare_process(path, 'zh', 'ZH_HToBB_ZToLL_M125_13TeV_powheg_pythia8', tag) - processes['zh'] = p - if DEBUG: print p - - p = Process('zz') - p.prepare_process(path, 'zz', 'ZZTo2L2Q_13TeV_amcatnloFXFX_madspin_pythia8_Fall15MiniAODv2', tag) - processes['zz'] = p - if DEBUG: print p - - # ttbar - p = Process('ttbar') - p.prepare_process(path, 'ttbar', 'TTTo2L2Nu_13TeV-powheg_Fall15MiniAODv2', tag) - processes['ttbar'] = p - p = Process('ttbar') - if DEBUG: print p - ''' - # drell-yan - p = Process('dy1') - p.prepare_process(path, 'dy1', 'DYJetsToLL_M-10to50_TuneCUETP8M1_13TeV-amcatnloFXFX_Fall15MiniAODv2', tag) - processes['dy1'] = p - if DEBUG: print p - ''' - p = Process('dy2') - p.prepare_process(path, 'dy2', 'DYJetsToLL_M-50_TuneCUETP8M1_13TeV-amcatnloFXFX_Fall15MiniAODv2', tag) - processes['dy2'] = p - if DEBUG: print p - - - - c.AddObservations([MASS], [ANALYSIS], [ERA], [CHANNEL], cats) - c.AddProcesses([MASS], [ANALYSIS], [ERA], [CHANNEL], ['ZA'], cats, True) - c.AddProcesses([MASS], [ANALYSIS], [ERA], [CHANNEL], ['ttbar', 'dy2','zz','zh'], cats, False) - c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( - c, "lumi_13TeV2015", "lnN", ch.SystMap('channel', 'era', 'bin_id') - ([CHANNEL], [ERA], [0,1,2,3], 1.027)) - - c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( - c, "trig", "lnN", ch.SystMap('channel', 'era', 'bin_id') - ([CHANNEL], [ERA], [0,1,2,3], 1.04)) - - c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( - c, "btag", "shape", ch.SystMap()(1.0)) - - c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( - c, "jec", "shape", ch.SystMap()(1.0)) - - c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( - c, "jer", "shape", ch.SystMap()(1.0)) - - c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( - c, "pu", "shape", ch.SystMap()(1.0)) - - c.cp().process(['ttbar']).AddSyst( - c, "TTpdf", "shape", ch.SystMap()(1.0)) - - c.cp().process(['dy2']).AddSyst( - c, "DYpdf", "shape", ch.SystMap()(1.0)) - - c.cp().process(['dy2']).AddSyst( - c, "DYnorm", "lnN", ch.SystMap('channel', 'era', 'bin_id') - ([CHANNEL], [ERA], [0,1], 1.04)) - - c.cp().process(['ttbar']).AddSyst( - c, "TTnorm", "lnN", ch.SystMap('channel', 'era', 'bin_id') - ([CHANNEL], [ERA], [0], 1.08)) - - - nChannels = len(bins) - nBackgrounds = len([processes[x] for x in processes if processes[x].type > 0]) - nNuisances = 1 - - systematics = {'':'','_btagUp':'__btagup', '_btagDown':'__btagdown', - '_jecUp':'__jecup','_jecDown':'__jecdown', - '_jerUp':'__jerup','_jerDown':'__jerdown', - '_puUp':'__puup', '_puDown':'__pudown', - '_TTpdfUp':'__pdfup','_TTpdfDown':'__pdfdown', '_DYpdfUp':'__pdfup','_DYpdfDown':'__pdfdown'} - outputRoot = "shapes.root" - f = TFile(outputRoot, "recreate") - f.Close() - for b in bins: - #print b , bins[b] - for p in processes: - if p == 'data_obs' : - file_in = TFile(processes[p].file,"READ") - #print " Getting ", bins[b], " in file ", processes[p].file - h = file_in.Get(bins[b]) - h.SetDirectory(0) - file_in.Close() - f = TFile(outputRoot, "update") - h.SetName("hist_"+bins[b]+"_"+p) - h.Write() - f.Write() - f.Close() - - - else : - for s1,s2 in systematics.iteritems() : + file_sig_path = "/home/fynu/amertens/scratch/cmssw/CMSSW_7_6_3/src/cp3_llbb/CommonTools/histFactory/Plots_v3/condor/output/HToZATo2L2B_MH-"+str(int(mH))+"_MA-"+str(int(mA))+"_13TeV-madgraph_Fall15MiniAODv1_v1.2.0+7415-73-g91bfcf0_ZAAnalysis_4f7ac83_histos.root" + if os.path.isfile(file_sig_path) : + + """Main function""" + # start the timer + tstart = datetime.now() + print 'starting...' + # get the options + #options = get_options() + + intL = options.lumi # in pb-1 + #tag = 'v1.2.0+7415-19-g7bbca78_ZAAnalysis_1a69757' + #path = '/nfs/scratch/fynu/amertens/cmssw/CMSSW_7_4_15/src/cp3_llbb/CommonTools/histFactory/16_01_28_syst/build' + tag = 'v1.2.0+7415-73-g91bfcf0_ZAAnalysis_4f7ac83' + #tag = 'v1.1.0+7415-57-g4bff5ea_ZAAnalysis_b1377a8' + path = options.path + CHANNEL = options.CHANNEL + ERA = options.ERA + MASS = str(mH)+"_"+str(mA) + ANALYSIS = options.ANALYSIS + DEBUG = 0 + + c = ch.CombineHarvester() + cats = [(0, "mmbbSR"+cutkey), + #(1, "eebbSR"+cutkey) + (1, "mll_mmbbBR"+cutkey), + (2, "eebbSR"+cutkey), + (3, "mll_eebbBR"+cutkey) + ] + + bins = {} + bins['signalregion_mm'] = "mmbbSR"+cutkey + bins['mll_bkgregion_mm'] = "mll_mmbbBR"+cutkey + bins['signalregion_ee'] = "eebbSR"+cutkey + bins['mll_bkgregion_ee'] = "mll_eebbBR"+cutkey + + processes = {} + p = Process('data_obs') + #DoubleMuon_Run2015D_v1.1.0+7415-57-g4bff5ea_ZAAnalysis_b1377a8_histos.root + p.prepare_process(path, 'data_obs', 'DoubleMuon_DoubleEG_Run2015D_v2_2015-12-18', tag) + processes['data_obs'] = p + if DEBUG: print p + # define signal + # define backgrounds + # zz + + p = Process('zh') + p.prepare_process(path, 'zh', 'ZH_HToBB_ZToLL_M125_13TeV_powheg_pythia8', tag) + processes['zh'] = p + if DEBUG: print p + + p = Process('zz') + p.prepare_process(path, 'zz', 'ZZTo2L2Q_13TeV_amcatnloFXFX_madspin_pythia8_Fall15MiniAODv2', tag) + processes['zz'] = p + if DEBUG: print p + + # ttbar + p = Process('ttbar') + p.prepare_process(path, 'ttbar', 'TTTo2L2Nu_13TeV-powheg_Fall15MiniAODv2', tag) + processes['ttbar'] = p + p = Process('ttbar') + if DEBUG: print p + ''' + # drell-yan + p = Process('dy1') + p.prepare_process(path, 'dy1', 'DYJetsToLL_M-10to50_TuneCUETP8M1_13TeV-amcatnloFXFX_Fall15MiniAODv2', tag) + processes['dy1'] = p + if DEBUG: print p + ''' + p = Process('dy2') + p.prepare_process(path, 'dy2', 'DYJetsToLL_M-50_TuneCUETP8M1_13TeV-amcatnloFXFX_Fall15MiniAODv2', tag) + processes['dy2'] = p + if DEBUG: print p + + + + c.AddObservations([MASS], [ANALYSIS], [ERA], [CHANNEL], cats) + c.AddProcesses([MASS], [ANALYSIS], [ERA], [CHANNEL], ['ZA'], cats, True) + c.AddProcesses([MASS], [ANALYSIS], [ERA], [CHANNEL], ['ttbar', 'dy2','zz','zh'], cats, False) + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "lumi_13TeV2015", "lnN", ch.SystMap('channel', 'era', 'bin_id') + ([CHANNEL], [ERA], [0,1,2,3], 1.027)) + + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "trig", "lnN", ch.SystMap('channel', 'era', 'bin_id') + ([CHANNEL], [ERA], [0,1,2,3], 1.04)) + + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "btag", "shape", ch.SystMap()(1.0)) + + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "jec", "shape", ch.SystMap()(1.0)) + + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "jer", "shape", ch.SystMap()(1.0)) + + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "pu", "shape", ch.SystMap()(1.0)) + + c.cp().process(['ttbar']).AddSyst( + c, "TTpdf", "shape", ch.SystMap()(1.0)) + + c.cp().process(['dy2']).AddSyst( + c, "DYpdf", "shape", ch.SystMap()(1.0)) + + c.cp().process(['dy2']).AddSyst( + c, "DYnorm", "lnN", ch.SystMap('channel', 'era', 'bin_id') + ([CHANNEL], [ERA], [0,1], 1.04)) + + c.cp().process(['ttbar']).AddSyst( + c, "TTnorm", "lnN", ch.SystMap('channel', 'era', 'bin_id') + ([CHANNEL], [ERA], [0], 1.08)) + + + nChannels = len(bins) + nBackgrounds = len([processes[x] for x in processes if processes[x].type > 0]) + + systematics = {'':'','_btagUp':'__btagup', '_btagDown':'__btagdown', + '_jecUp':'__jecup','_jecDown':'__jecdown', + '_jerUp':'__jerup','_jerDown':'__jerdown', + '_puUp':'__puup', '_puDown':'__pudown', + '_TTpdfUp':'__pdfup','_TTpdfDown':'__pdfdown', '_DYpdfUp':'__pdfup','_DYpdfDown':'__pdfdown'} + outputRoot = "shapes.root" + f = TFile(outputRoot, "recreate") + f.Close() + for b in bins: + #print b , bins[b] + for p in processes: + if p == 'data_obs' : file_in = TFile(processes[p].file,"READ") - #if s1 == '' : - # print " Getting ", bins[b]+s2, " in file ", processes[p].file - h = file_in.Get(bins[b]+s2) - #if s1 == '' : - # print "Val: ", h.GetBinContent(1), "ERROR : ", h.GetBinError(1), 'Entries : ', h.GetEntries() - ''' - if s1 == '' and 'mmbbSR' in bins[b] : - print "Adding MC stat uncertainty in ", p, bins[b] , " of ", - if h.GetBinContent(1) > 0.0 : - MCStatUnc = 1+h.GetBinError(1)/h.GetBinContent(1) - else : - MCStatUnc = 2.0 - print h.GetBinError(1), h.GetBinContent(1), h.GetEntries(), MCStatUnc - - c.cp().process([p]).AddSyst( - c, p+bins[b]+"MCStat", "lnN", ch.SystMap('bin_id') - ([0],MCStatUnc)) - if s1 == '' and 'eebbSR' in bins[b] : - print "Adding MC stat uncertainty in ", p, bins[b] , " of ", - MCStatUnc = 1+h.GetBinError(1)/h.GetBinContent(1) - print h.GetBinError(1), h.GetBinContent(1), h.GetEntries(), MCStatUnc - - c.cp().process([p]).AddSyst( - c, p+bins[b]+"MCStat", "lnN", ch.SystMap('bin_id') - ([2],MCStatUnc)) - ''' - + #print " Getting ", bins[b], " in file ", processes[p].file + h = file_in.Get(bins[b]) h.SetDirectory(0) file_in.Close() f = TFile(outputRoot, "update") - h.SetName("hist_"+bins[b]+"_"+p+s1) - Nbins = h.GetNbinsX() - for it_b in range(0,h.GetNbinsX()+1) : - if h.GetBinContent(it_b) < 0.0 : - print "WARNING : negative bin content in bin ", it_b, " of histo ", "hist_"+bins[b]+"_"+p+s1 - h.SetBinContent(it_b,0.0) - #h.Sumw2() - #h.Scale(processes[p].xsection * intL / processes[p].sumW) - h.Scale(intL) - h.Write() - f.Write() - f.Close() - - for s1,s2 in systematics.iteritems() : - file_sig_path = "/home/fynu/amertens/scratch/cmssw/CMSSW_7_6_3/src/cp3_llbb/CommonTools/histFactory/Plots_v3/condor/output/HToZATo2L2B_MH-"+str(int(mH))+"_MA-"+str(int(mA))+"_13TeV-madgraph_Fall15MiniAODv1_v1.2.0+7415-73-g91bfcf0_ZAAnalysis_4f7ac83_histos.root" - file_sig = TFile(file_sig_path, "READ") - h = file_sig.Get(bins[b]+s2) - h.SetDirectory(0) - file_sig.Close() - f = TFile(outputRoot, "update") - h.SetName("hist_"+bins[b]+"_"+"ZA"+s1) - - if DEBUG and s1 == '' : - #if s1 == '' : - print file_sig_path - print "saving ", "hist_"+bins[b]+"_"+"ZA"+s1 - print h.GetBinContent(1), "after rescaling", - h.Scale(options.lumifb) - #if DEBUG and s1 == '' : - if s1 == '' : - print h.GetBinContent(1) + h.SetName("hist_"+bins[b]+"_"+p) h.Write() f.Write() f.Close() - - # Fill signal histograms FIXME: read efficiencies from eff.root - - ''' - - eff_file = TFile("eff.root", "READ") - effee_hist = eff_file.Get("effee") - eff_ee = effee_hist.Interpolate(mA,mH) - effmm_hist = eff_file.Get("effmm") - eff_mm = effmm_hist.Interpolate(mA,mH) - - print "lumi : ", options.lumifb - print "eff at ", mA, mH, ":", eff_ee, eff_mm - print "ZA yields: ", options.lumifb*eff_mm, options.lumifb*eff_ee - - - f = TFile(outputRoot, "update") - h1 = TH1F("hist_"+bins['signalregion_mm']+"_ZA","hist_"+bins['signalregion_mm']+"_ZA",1,0,1) - h1.Fill(0.5,options.lumifb*eff_mm) - h1.Write() - - h2 = TH1F("hist_"+bins['mll_bkgregion_mm']+"_ZA","hist_"+bins['mll_bkgregion_mm']+"_ZA",60,60,120) - h2.Write() - - h3 = TH1F("hist_"+bins['signalregion_ee']+"_ZA","hist_"+bins['signalregion_ee']+"_ZA",1,0,1) - h3.Fill(0.5,options.lumifb*eff_ee) - h3.Write() - - h4 = TH1F("hist_"+bins['mll_bkgregion_ee']+"_ZA","hist_"+bins['mll_bkgregion_ee']+"_ZA",60,60,120) - h4.Write() - ''' - - f.Write() - f.Close() - - c.cp().backgrounds().ExtractShapes( - outputRoot, "hist_$BIN_$PROCESS", "hist_$BIN_$PROCESS_$SYSTEMATIC") - c.cp().signals().ExtractShapes( - outputRoot, "hist_$BIN_$PROCESS", "hist_$BIN_$PROCESS_$SYSTEMATIC") - - bbb = ch.BinByBinFactory().SetAddThreshold(0.05).SetFixNorm(False) - - bbb.AddBinByBin(c.cp().backgrounds(),c) - ch.SetStandardBinNames(c) - - writer = ch.CardWriter('$TAG/$MASS/$ANALYSIS_$CHANNEL_$ERA.dat', - '$TAG/common/$ANALYSIS_$CHANNEL_$MASS.input_$ERA.root') - writer.WriteCards('CARDS/', c) + + + else : + for s1,s2 in systematics.iteritems() : + file_in = TFile(processes[p].file,"READ") + #if s1 == '' : + # print " Getting ", bins[b]+s2, " in file ", processes[p].file + h = file_in.Get(bins[b]+s2) + #if s1 == '' : + # print "Val: ", h.GetBinContent(1), "ERROR : ", h.GetBinError(1), 'Entries : ', h.GetEntries() + h.SetDirectory(0) + file_in.Close() + f = TFile(outputRoot, "update") + h.SetName("hist_"+bins[b]+"_"+p+s1) + Nbins = h.GetNbinsX() + h.Scale(intL) + + for it_b in range(1,h.GetNbinsX()+1) : + """ + if h.GetBinContent(it_b) == 0 : + print "WARNING : no entries in bin ", it_b, " of histo ", "hist_"+bins[b]+"_"+p+s1, "setting high error to 1.8 x 0.73", + h.SetBinContent(it_b,1e-5) + err = 1.8*0.73 + print err + h.SetBinError(it_b,err) + """ + if h.GetBinContent(it_b) <= 0: + print "WARNING : negative bin content in bin ", it_b, " of histo ", "hist_"+bins[b]+"_"+p+s1 + h.SetBinContent(it_b,1e-6) + #h.Sumw2() + #h.Scale(processes[p].xsection * intL / processes[p].sumW) + #h.Scale(intL) + + h.Write() + f.Write() + f.Close() + + for s1,s2 in systematics.iteritems() : + file_sig = TFile(file_sig_path, "READ") + h = file_sig.Get(bins[b]+s2) + h.SetDirectory(0) + file_sig.Close() + f = TFile(outputRoot, "update") + h.SetName("hist_"+bins[b]+"_"+"ZA"+s1) + + if DEBUG and s1 == '' : + print file_sig_path + print "saving ", "hist_"+bins[b]+"_"+"ZA"+s1 + print h.GetBinContent(1), "after rescaling", + h.Scale(options.lumifb) + if DEBUG and s1 == '' : + print h.GetBinContent(1) + h.Write() + f.Write() + f.Close() + + # Fill signal histograms FIXME: read efficiencies from eff.root + + ''' + + eff_file = TFile("eff.root", "READ") + effee_hist = eff_file.Get("effee") + eff_ee = effee_hist.Interpolate(mA,mH) + effmm_hist = eff_file.Get("effmm") + eff_mm = effmm_hist.Interpolate(mA,mH) + + print "lumi : ", options.lumifb + print "eff at ", mA, mH, ":", eff_ee, eff_mm + print "ZA yields: ", options.lumifb*eff_mm, options.lumifb*eff_ee + + + f = TFile(outputRoot, "update") + h1 = TH1F("hist_"+bins['signalregion_mm']+"_ZA","hist_"+bins['signalregion_mm']+"_ZA",1,0,1) + h1.Fill(0.5,options.lumifb*eff_mm) + h1.Write() + + h2 = TH1F("hist_"+bins['mll_bkgregion_mm']+"_ZA","hist_"+bins['mll_bkgregion_mm']+"_ZA",60,60,120) + h2.Write() + + h3 = TH1F("hist_"+bins['signalregion_ee']+"_ZA","hist_"+bins['signalregion_ee']+"_ZA",1,0,1) + h3.Fill(0.5,options.lumifb*eff_ee) + h3.Write() + + h4 = TH1F("hist_"+bins['mll_bkgregion_ee']+"_ZA","hist_"+bins['mll_bkgregion_ee']+"_ZA",60,60,120) + h4.Write() + ''' + + f.Write() + f.Close() + + c.cp().backgrounds().ExtractShapes( + outputRoot, "hist_$BIN_$PROCESS", "hist_$BIN_$PROCESS_$SYSTEMATIC") + c.cp().signals().ExtractShapes( + outputRoot, "hist_$BIN_$PROCESS", "hist_$BIN_$PROCESS_$SYSTEMATIC") + + bbb = ch.BinByBinFactory().SetAddThreshold(0.05).SetFixNorm(False) + + bbb.AddBinByBin(c.cp().backgrounds(),c) + ch.SetStandardBinNames(c) + + writer = ch.CardWriter('$TAG/$MASS/$ANALYSIS_$CHANNEL_$ERA.dat', + '$TAG/common/$ANALYSIS_$CHANNEL_$MASS.input_$ERA.root') + writer.WriteCards('CARDS/', c) # # main diff --git a/python/createCards_combined.py b/python/createCards_combined.py new file mode 100644 index 0000000..4452f8d --- /dev/null +++ b/python/createCards_combined.py @@ -0,0 +1,287 @@ +#! /usr/bin/env python + +# Python imports +#import os, sys, argparse, getpass +from datetime import datetime +# User imports +from classes import * +import CombineHarvester.CombineTools.ch as ch +import os +# ROOT imports +import ROOT +from ROOT import gROOT +from ROOT import TChain, TFile, TCanvas +from ROOT import TH1F +from cp3_llbb.ZATools.ZACnC import * + +#gROOT.Reset() +#gROOT.SetBatch() +#ROOT.PyConfig.IgnoreCommandLineOptions = True + +def main(): + options = options_() + for cutkey in options.cut : + print 'cutkey : ', cutkey + ### get M_A and M_H ### + mH = float(options.mH_list[cutkey]) + mA = float(options.mA_list[cutkey]) + print mH, mA + + file_sig_path = "/home/fynu/amertens/scratch/cmssw/CMSSW_7_6_3/src/cp3_llbb/CommonTools/histFactory/Plots_v5/condor/output/HToZATo2L2B_MH-"+str(int(mH))+"_MA-"+str(int(mA))+"_13TeV-madgraph_Fall15MiniAODv1_v1.2.0+7415-73-g91bfcf0_ZAAnalysis_4f7ac83_histos.root" + if os.path.isfile(file_sig_path) : + + """Main function""" + # start the timer + tstart = datetime.now() + print 'starting...' + # get the options + #options = get_options() + + intL = options.lumi # in pb-1 + #tag = 'v1.2.0+7415-19-g7bbca78_ZAAnalysis_1a69757' + #path = '/nfs/scratch/fynu/amertens/cmssw/CMSSW_7_4_15/src/cp3_llbb/CommonTools/histFactory/16_01_28_syst/build' + tag = 'v1.2.0+7415-73-g91bfcf0_ZAAnalysis_4f7ac83' + #tag = 'v1.1.0+7415-57-g4bff5ea_ZAAnalysis_b1377a8' + path = options.path + CHANNEL = options.CHANNEL + ERA = options.ERA + MASS = str(mH)+"_"+str(mA) + ANALYSIS = options.ANALYSIS + DEBUG = 0 + + c = ch.CombineHarvester() + cats = [(0, "llbbSR"+cutkey), + (1, "mll_llbbBR"+cutkey) + ] + + bins = {} + bins['signalregion_ll'] = "llbbSR"+cutkey + bins['mll_bkgregion_ll'] = "mll_llbbBR"+cutkey + + processes = {} + p = Process('data_obs') + #DoubleMuon_Run2015D_v1.1.0+7415-57-g4bff5ea_ZAAnalysis_b1377a8_histos.root + p.prepare_process(path, 'data_obs', 'DoubleMuon_DoubleEG_Run2015D_v2_2015-12-18', tag) + processes['data_obs'] = p + if DEBUG: print p + # define signal + # define backgrounds + # zz + + p = Process('zh') + p.prepare_process(path, 'zh', 'ZH_HToBB_ZToLL_M125_13TeV_powheg_pythia8', tag) + processes['zh'] = p + if DEBUG: print p + + p = Process('zz') + p.prepare_process(path, 'zz', 'ZZTo2L2Q_13TeV_amcatnloFXFX_madspin_pythia8_Fall15MiniAODv2', tag) + processes['zz'] = p + if DEBUG: print p + + # ttbar + p = Process('ttbar') + p.prepare_process(path, 'ttbar', 'TTTo2L2Nu_13TeV-powheg_Fall15MiniAODv2', tag) + processes['ttbar'] = p + p = Process('ttbar') + if DEBUG: print p + ''' + # drell-yan + p = Process('dy1') + p.prepare_process(path, 'dy1', 'DYJetsToLL_M-10to50_TuneCUETP8M1_13TeV-amcatnloFXFX_Fall15MiniAODv2', tag) + processes['dy1'] = p + if DEBUG: print p + ''' + p = Process('dy2') + p.prepare_process(path, 'dy2', 'DYJetsToLL_M-50_TuneCUETP8M1_13TeV-amcatnloFXFX_Fall15MiniAODv2', tag) + processes['dy2'] = p + if DEBUG: print p + + + + c.AddObservations([MASS], [ANALYSIS], [ERA], [CHANNEL], cats) + c.AddProcesses([MASS], [ANALYSIS], [ERA], [CHANNEL], ['ZA'], cats, True) + c.AddProcesses([MASS], [ANALYSIS], [ERA], [CHANNEL], ['ttbar', 'dy2','zz','zh'], cats, False) + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "lumi_13TeV2015", "lnN", ch.SystMap('channel', 'era', 'bin_id') + ([CHANNEL], [ERA], [0,1,2,3], 1.027)) + + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "trig", "lnN", ch.SystMap('channel', 'era', 'bin_id') + ([CHANNEL], [ERA], [0,1,2,3], 1.04)) + + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "btag", "shape", ch.SystMap()(1.0)) + + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "jec", "shape", ch.SystMap()(1.0)) + + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "jer", "shape", ch.SystMap()(1.0)) + + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "elID", "shape", ch.SystMap()(1.0)) + + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "muID", "shape", ch.SystMap()(1.0)) + + c.cp().process(['ttbar', 'dy2','zz','zh','ZA']).AddSyst( + c, "pu", "shape", ch.SystMap()(1.0)) + + c.cp().process(['ttbar']).AddSyst( + c, "TTpdf", "shape", ch.SystMap()(1.0)) + + c.cp().process(['dy2']).AddSyst( + c, "DYpdf", "shape", ch.SystMap()(1.0)) + + c.cp().process(['ttbar']).AddSyst( + c, "TTscale", "shape", ch.SystMap()(1.0)) + + c.cp().process(['dy2']).AddSyst( + c, "DYscale", "shape", ch.SystMap()(1.0)) + + c.cp().process(['dy2']).AddSyst( + c, "DYnorm", "lnN", ch.SystMap('channel', 'era', 'bin_id') + ([CHANNEL], [ERA], [0,1], 1.04)) + + c.cp().process(['ttbar']).AddSyst( + c, "TTnorm", "lnN", ch.SystMap('channel', 'era', 'bin_id') + ([CHANNEL], [ERA], [0], 1.08)) + + + nChannels = len(bins) + nBackgrounds = len([processes[x] for x in processes if processes[x].type > 0]) + + systematics = {'':'','_btagUp':'__btagup', '_btagDown':'__btagdown', + '_jecUp':'__jecup','_jecDown':'__jecdown', + '_jerUp':'__jerup','_jerDown':'__jerdown', + '_elIDUp':'__elIDup', '_elIDDown':'__elIDdown', + '_muIDUp':'__muIDup', '_muIDDown':'__muIDdown', + '_puUp':'__puup', '_puDown':'__pudown', + '_TTpdfUp':'__pdfup','_TTpdfDown':'__pdfdown', '_DYpdfUp':'__pdfup','_DYpdfDown':'__pdfdown', + '_TTscaleUp':'__scaleup','_TTscaleDown':'__scaledown', '_DYscaleUp':'__scaleup','_DYscaleDown':'__scaledown'} + outputRoot = "shapes.root" + f = TFile(outputRoot, "recreate") + f.Close() + for b in bins: + #print b , bins[b] + for p in processes: + if p == 'data_obs' : + file_in = TFile(processes[p].file,"READ") + #print " Getting ", bins[b], " in file ", processes[p].file + h = file_in.Get(bins[b]) + h.SetDirectory(0) + file_in.Close() + f = TFile(outputRoot, "update") + h.SetName("hist_"+bins[b]+"_"+p) + h.Write() + f.Write() + f.Close() + + + else : + for s1,s2 in systematics.iteritems() : + file_in = TFile(processes[p].file,"READ") + #if s1 == '' : + # print " Getting ", bins[b]+s2, " in file ", processes[p].file + h = file_in.Get(bins[b]+s2) + #if s1 == '' : + # print "Val: ", h.GetBinContent(1), "ERROR : ", h.GetBinError(1), 'Entries : ', h.GetEntries() + h.SetDirectory(0) + file_in.Close() + f = TFile(outputRoot, "update") + h.SetName("hist_"+bins[b]+"_"+p+s1) + Nbins = h.GetNbinsX() + h.Scale(intL) + + for it_b in range(1,h.GetNbinsX()+1) : + """ + if h.GetBinContent(it_b) == 0 : + print "WARNING : no entries in bin ", it_b, " of histo ", "hist_"+bins[b]+"_"+p+s1, "setting high error to 1.8 x 0.73", + h.SetBinContent(it_b,1e-5) + err = 1.8*0.73 + print err + h.SetBinError(it_b,err) + """ + if h.GetBinContent(it_b) <= 0: + print "WARNING : negative bin content in bin ", it_b, " of histo ", "hist_"+bins[b]+"_"+p+s1 + h.SetBinContent(it_b,1e-6) + #h.Sumw2() + #h.Scale(processes[p].xsection * intL / processes[p].sumW) + #h.Scale(intL) + + h.Write() + f.Write() + f.Close() + + for s1,s2 in systematics.iteritems() : + file_sig = TFile(file_sig_path, "READ") + h = file_sig.Get(bins[b]+s2) + h.SetDirectory(0) + file_sig.Close() + f = TFile(outputRoot, "update") + h.SetName("hist_"+bins[b]+"_"+"ZA"+s1) + + if DEBUG and s1 == '' : + print file_sig_path + print "saving ", "hist_"+bins[b]+"_"+"ZA"+s1 + print h.GetBinContent(1), "after rescaling", + h.Scale(options.lumifb) + if DEBUG and s1 == '' : + print h.GetBinContent(1) + h.Write() + f.Write() + f.Close() + + # Fill signal histograms FIXME: read efficiencies from eff.root + + ''' + + eff_file = TFile("eff.root", "READ") + effee_hist = eff_file.Get("effee") + eff_ee = effee_hist.Interpolate(mA,mH) + effmm_hist = eff_file.Get("effmm") + eff_mm = effmm_hist.Interpolate(mA,mH) + + print "lumi : ", options.lumifb + print "eff at ", mA, mH, ":", eff_ee, eff_mm + print "ZA yields: ", options.lumifb*eff_mm, options.lumifb*eff_ee + + + f = TFile(outputRoot, "update") + h1 = TH1F("hist_"+bins['signalregion_mm']+"_ZA","hist_"+bins['signalregion_mm']+"_ZA",1,0,1) + h1.Fill(0.5,options.lumifb*eff_mm) + h1.Write() + + h2 = TH1F("hist_"+bins['mll_bkgregion_mm']+"_ZA","hist_"+bins['mll_bkgregion_mm']+"_ZA",60,60,120) + h2.Write() + + h3 = TH1F("hist_"+bins['signalregion_ee']+"_ZA","hist_"+bins['signalregion_ee']+"_ZA",1,0,1) + h3.Fill(0.5,options.lumifb*eff_ee) + h3.Write() + + h4 = TH1F("hist_"+bins['mll_bkgregion_ee']+"_ZA","hist_"+bins['mll_bkgregion_ee']+"_ZA",60,60,120) + h4.Write() + ''' + + f.Write() + f.Close() + + c.cp().backgrounds().ExtractShapes( + outputRoot, "hist_$BIN_$PROCESS", "hist_$BIN_$PROCESS_$SYSTEMATIC") + c.cp().signals().ExtractShapes( + outputRoot, "hist_$BIN_$PROCESS", "hist_$BIN_$PROCESS_$SYSTEMATIC") + + bbb = ch.BinByBinFactory().SetAddThreshold(0.05).SetFixNorm(False) + + bbb.AddBinByBin(c.cp().backgrounds(),c) + ch.SetStandardBinNames(c) + + writer = ch.CardWriter('$TAG/$MASS/$ANALYSIS_$CHANNEL_$ERA.dat', + '$TAG/common/$ANALYSIS_$CHANNEL_$MASS.input_$ERA.root') + writer.WriteCards('CARDS_combined/', c) + +# +# main +# +if __name__ == '__main__': + main() diff --git a/python/createXsecLines.py b/python/createXsecLines.py new file mode 100644 index 0000000..c8b7f1d --- /dev/null +++ b/python/createXsecLines.py @@ -0,0 +1,53 @@ +#!/usr/bin/python + + +import math +from cp3_llbb.Calculators42HDM.Calc2HDM import * +from ROOT import * + + +mode = 'H' +sqrts = 13000 +type = 2 +tb = 1 +m12 = 0 +mh = 125 +mH = 300 +mA = 50 +mhc = mH + +beta=math.atan(tb) +cba = 0.01 +alpha=math.atan(tb)-math.acos(cba) +sba = math.sin(math.atan(tb)-alpha) + +print 'sba : ', sba + +outputFile = "out.dat" + + +test = Calc2HDM(mode = 'H', sqrts = sqrts, type = type, tb = tb, m12 = m12, mh = mh, mH = mH, mA = mA, mhc = mhc, sba = sba, outputFile = outputFile) +test.computeBR() + +xsec = test.getXsecFromSusHi() + +xsectot_list = [] + + +Xsec_mH500_tb1 = TGraph(50) + +n=0 + +while mA < mH-90: + test.setmA(mA) + test.computeBR() + print "ZA BR", test.HtoZABR + Xsec_mH500_tb1.SetPoint(n,mA,xsec*test.HtoZABR*test.AtobbBR*0.067) + mA+=2 + n+=1 + + + +file_out = TFile("xsec.root","RECREATE") +Xsec_mH500_tb1.Write() + diff --git a/python/drawSyst.py b/python/drawSyst.py index 11e1979..832f1fa 100644 --- a/python/drawSyst.py +++ b/python/drawSyst.py @@ -16,7 +16,7 @@ def main(): gr_statMC = TGraph(4) gr_stat = TGraph(4) - mH_tested = 500.0 + mH_tested = 300.0 syst_list = {"btag":gr_btag, "jec":gr_jec, diff --git a/python/getPostfitValues.py b/python/getPostfitValues.py new file mode 100644 index 0000000..21aea6c --- /dev/null +++ b/python/getPostfitValues.py @@ -0,0 +1,128 @@ +#! /usr/bin/env python + +import os, sys, argparse, math + +# to prevent pyroot to hijack argparse we need to go around +tmpargv = sys.argv[:] +sys.argv = [] +# ROOT imports +from ROOT import gROOT, gSystem, PyConfig, TFile +gROOT.Reset() +gROOT.SetBatch() +PyConfig.IgnoreCommandLineOptions = True +sys.argv = tmpargv + +def mapping(bkg): + if bkg == "ttbar": + return "TTTo2L2Nu_13TeV" + elif bkg == "dy2": + return "DYJetsToLL" + elif bkg == "singletop": + return "ST_" + elif bkg == "ttV": + return "TT(WJets|Z)To" + elif bkg == "VV": + return "(VV|ZZ)To" + +def prettyName(bkg): + if bkg == "ttbar": + return "TT" + elif bkg == "dy2": + return "DY" + elif bkg == "singletop": + return "ST" + elif bkg == "ttV": + return "TT+V" + elif bkg == "VV": + return "VV" + + +parser = argparse.ArgumentParser(description='Compute data/MC scale factors from a MaxLikelihoodFit') + +parser.add_argument('-i', '--input', action='store', type=str, dest='input', help='Path to the mlfit ROOT file created by combine', required=True) + +options = parser.parse_args() + +# Compute scale factors + +mlfit = TFile.Open(options.input) + +prefit_combine_norm = mlfit.Get('norm_prefit') +postfit_combine_norm = mlfit.Get('norm_fit_b') + +channels = [] +prefit_shapes = mlfit.Get('shapes_prefit') +for k in prefit_shapes.GetListOfKeys(): + channels.append(k.GetName()) + +print "Detected channels: ", channels + + +# Construct the list of backgrounds +# Naming is 'category/bkg_name' + +backgrounds = [] + +all_backgrounds = prefit_combine_norm.contentsString().split(',') +for bkg in all_backgrounds: + if not bkg.startswith(channels[0]): + continue + + backgrounds.append(bkg.split('/')[-1]) + +print 'Detected backgrounds: ', backgrounds + +prefit_norm = {} +postfit_norm = {} + +for channel in channels: + for bkg in backgrounds: + name = channel + '/' + bkg + + prefit = prefit_combine_norm[name] + postfit = postfit_combine_norm[name] + + if postfit.getVal() == 0 or prefit.getVal() == 0: + continue + + if not bkg in prefit_norm: + prefit_norm[bkg] = [prefit.getVal(), prefit.getError()**2] + postfit_norm[bkg] = [postfit.getVal(), postfit.getError()**2] + else: + prefit_norm[bkg][0] += prefit.getVal() + postfit_norm[bkg][0] += postfit.getVal() + + prefit_norm[bkg][1] += prefit.getError()**2 + postfit_norm[bkg][1] += postfit.getError()**2 + +scale_factors = {} +print '' +print 'Scale factors: ' +for bkg in prefit_norm.keys(): + scale_factors[bkg] = postfit_norm[bkg][0] / prefit_norm[bkg][0] + print('%s: %f' % (bkg, scale_factors[bkg])) + +print("") + +print 'Scale factors (for plotIt): ' +for bkg, sf in scale_factors.items(): + if sf == 0: + continue + print("%s\n\tscale: %f" % (bkg, sf)) + + +print("") +print("Relative post-fit uncertainties") +for bkg, values in postfit_norm.items(): + print "%s: %.4f" % (bkg, math.sqrt(values[1]) / values[0]) + +print("") +print("Post-fit uncertainties (for plotIt -- YAML)") + +print("systematics:") + +for bkg, values in postfit_norm.items(): + m = mapping(bkg) + if not m: + continue + print(" - %s: {type: const, value: %.4f, on: '%s', pretty-name: '%s postfit'}" % (bkg + "_postfit", 1 + math.sqrt(values[1]) / values[0], m, prettyName(bkg))) diff --git a/python/makeEfficiencymap.py b/python/makeEfficiencymap.py index 355c13b..62717e5 100644 --- a/python/makeEfficiencymap.py +++ b/python/makeEfficiencymap.py @@ -70,12 +70,12 @@ sigfile = TFile(Signal_path+file_name,"READ") tree = sigfile.Get("t") #print options.cut[cutkey] - weight_cut = options.cut[cutkey]+" * mumu_LooseZCandidate_cut" + weight_cut = options.cut[cutkey]+" * mumu_TightZCandidate_cut" tree.Draw("1>>tempHistmm",weight_cut,"") tempHistmm=gDirectory.Get("tempHistmm") eff_mm = tempHistmm.GetEntries()/tree.GetEntriesFast() - weight_cut = options.cut[cutkey]+" * elel_LooseZCandidate_cut" + weight_cut = options.cut[cutkey]+" * elel_TightZCandidate_cut" tree.Draw("1>>tempHistee",weight_cut,"") tempHistee=gDirectory.Get("tempHistee") eff_ee = tempHistee.GetEntries()/tree.GetEntriesFast() diff --git a/python/makeVBBP.py b/python/makeVBBP.py index b4a0c1c..05bc6e3 100644 --- a/python/makeVBBP.py +++ b/python/makeVBBP.py @@ -5,22 +5,62 @@ import math import os import os.path +import numpy as numpy from ROOT import * from ROOT import TMath #as tmath from ZACnC import * + #from ROOT.TMath import LorentzVector #ROOT.gSystem.Load('CMS_label_C') #ROOT.gSystem.Load("DrawCanvas_C") +def drawTheoryGraph(gr, tb): + + x1 = 290 + lenght = 40 + x2 = x1+20 + y1 = gr.Eval(x1) + y2 = gr.Eval(x2) + + pente = math.atan( 420/3.3 *(log10(y2)-log10(y1))/(x2-x1)) + angle = pente*180/math.pi + + x3 = x1 + lenght*math.cos(pente) + y3 = gr.Eval(x3) + + print angle + print x1, y1, x3, y3 + + size = 5 + x_pl = numpy.array([x1-size, x3+size, x3+size, x1-size, x1-size], dtype=float) #[x1+10, x2+10, x2-10, x1-10, x1+10] + y_pl = numpy.array([y1+size, y3+size, y3-size, y1-size, y1+size], dtype=float) + pline = TPolyLine(5,x_pl,y_pl) + pline.SetFillColor(kWhite) + pline.Draw("f") + + l = TLatex(x1,y1,"tan #beta = "+str(tb)) + l.SetTextSize(20) + l.SetTextFont(43) + l.SetTextAlign(12) + l.SetTextAngle(angle) + l.SetTextColor(gr.GetLineColor()) + l.Draw() + SetOwnership( l, False ) + SetOwnership( pline, False ) + + + + ################### ### Definitions ### ################### mH=500 -run_combine = 1 +run_combine = 0 +run_blind = 0 TGAS_1 = TGraphAsymmErrors(0) TGAS_2 = TGraphAsymmErrors(0) @@ -40,8 +80,8 @@ myTGraph.SetName("p-value") #DataCards_path = "../cards/" -DataCards_path = "CARDS/" -RootFiles_path = "../rootfiles/" +DataCards_path = "CARDS_combined/" +RootFiles_path = "../rootfiles_combined/" Signal_path = "/home/fynu/amertens/scratch/cmssw/CMSSW_7_6_3/src/cp3_llbb/ZAAnalysis/" @@ -84,7 +124,10 @@ #combine_cmd = "combine -M ProfileLikelihood --signif -m "+str(int(mbb))+" "+DataCard+" --toysFreq" #combine_cmd = "combine -M ProfileLikelihood --significance --pvalue -m "+str(int(mbb))+" "+DataCard - combine_cmd = "combine -M Asymptotic -m "+str(int(mbb))+" --run=blind "+DataCard + if run_blind == 1 : + combine_cmd = "combine -M Asymptotic -m "+str(int(mbb))+" --run=blind "+DataCard + else : + combine_cmd = "combine -M Asymptotic -m "+str(int(mbb))+" "+DataCard os.system(str(combine_cmd)) mv_cmd = "mv higgsCombineTest.Asymptotic.mH"+str(int(mbb))+".root "+RootFile os.system(str(mv_cmd)) @@ -171,14 +214,72 @@ #TGAS_2.SetMinimum(0) TGAS_2.Draw('A E3') TGAS_1.Draw('E3') -myTG_2.Draw("L") -myTG_8TeV.Draw("L") - -TGAS_2.GetXaxis().SetTitle("m_{A} (GeV)") -TGAS_2.GetYaxis().SetTitle("#sigma (pp #rightarrow H) #times BR(H #rightarrow Z(ll)A(bb)) (fb)") +myTG_2.Draw("L*") +myTG_2.SetMarkerStyle(20) +if run_blind == 0 : + myTG_5.Draw("L*") + myTG_5.SetMarkerStyle(20) +if run_blind == 1 : + myTG_8TeV.Draw("L") + +TGAS_2.GetXaxis().SetTitle("m_{A} [GeV]") +TGAS_2.GetYaxis().SetTitle("#sigma (pp #rightarrow H) #times BR(H #rightarrow Z(ll)A(bb)) [fb]") +TGAS_2.GetYaxis().SetTitleOffset(1.30) TGAS_2.SetMinimum(5) -TGAS_2.SetMaximum(1000) +TGAS_2.SetMaximum(10000) + + +## Adding theory curves +file_theory = TFile("xsec.root") +gr_tb1 = file_theory.Get("Xsec_mH"+str(mH)+"_tb1") +gr_tb1_5 = file_theory.Get("Xsec_mH"+str(mH)+"_tb1_5") + +gr_tb1.SetLineColor(kGray+1) +gr_tb1.SetLineStyle(2) +gr_tb1.SetLineWidth(3) + +gr_tb1_5.SetLineColor(kGray+3) +gr_tb1_5.SetLineStyle(2) +gr_tb1_5.SetLineWidth(3) + +gr_tb1.Draw("C") +gr_tb1_5.Draw("C") + +drawTheoryGraph(gr_tb1, 1) +drawTheoryGraph(gr_tb1_5, 1.5) + +C.Update() + +''' +x1 = 280 +x2 = x1+20 +y1 = gr_tb1.Eval(x1) +y2 = gr_tb1.Eval(x2) + +x3 = (x1+x2)/2 +y3 = gr_tb1.Eval(x3) + +pente = math.atan( 420/3.3 *(log10(y2)-log10(y1))/(x2-x1)) +#pente = 420/3.3 * math.atan((log10(y2)-log10(y1))/(x2-x1)) + +print y1, y2, x1, x2, (log10(y2)-log10(y1))/(x2-x1), pente + +size = 5 +x_pl = numpy.array([x1+size, x2+size, x2-size, x1-size, x1+size], dtype=float) #[x1+10, x2+10, x2-10, x1-10, x1+10] +y_pl = numpy.array([y1+size, y2+size, y2-size, y1-size, y1+size], dtype=float) +pline = TPolyLine(5,x_pl,y_pl) +pline.SetFillColor(kWhite) +pline.Draw("f") + +l = TLatex(x1,y1,"tan #beta = 1") +l.SetTextSize(20) +l.SetTextFont(43) +l.SetTextAlign(12) +l.SetTextAngle(pente*180/math.pi) +l.SetTextColor(kGray+1) +l.Draw() +''' #gPad().SetLogy(1) leg = TLegend(0.59,0.68,0.89,0.88) @@ -188,22 +289,24 @@ leg.AddEntry(TGAS_2,"CL_{s} Expected #pm 2 #sigma","F") leg.AddEntry(myTG_2,"CL_{s} Expected","L") leg.AddEntry(myTG_5,"CL_{s} Observed","L") -leg.AddEntry(myTG_8TeV,"8 TeV","L") +#leg.AddEntry(myTG_8TeV,"8 TeV","L") leg.Draw() gPad.SetLogy(1) +if run_blind : + C.Print("BBP"+str(mH)+"_blind_v2.pdf") + C.Print("BBP"+str(mH)+"_blind_v2.png") + f = TFile("narrow_BBP"+str(mH)+"_blind_v2.root","recreate") +else : + C.Print("BBP"+str(mH)+"_v2.pdf") + C.Print("BBP"+str(mH)+"_v2.png") + f = TFile("narrow_BBP"+str(mH)+"_v2.root","recreate") - -C.Print("BBP"+str(mH)+"_v1.pdf") -C.Print("BBP"+str(mH)+"_v1.png") -f = TFile("narrow_BBP"+str(mH)+"_v1.root","recreate") C.Write() f.Close() - - ########### ### FIN ### ########### From 026b1bfac67d6f91b86e590d343866594b6c1697 Mon Sep 17 00:00:00 2001 From: Alexandre Mertens Date: Thu, 10 Mar 2016 00:13:14 +0100 Subject: [PATCH 3/3] CMS style BBP --- python/makeVBBP.py | 125 +++++++++------ python/makeVBBP_300.py | 345 +++++++++++++++++++++++++++++++++++++++++ python/makeVBBP_800.py | 345 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 769 insertions(+), 46 deletions(-) create mode 100644 python/makeVBBP_300.py create mode 100644 python/makeVBBP_800.py diff --git a/python/makeVBBP.py b/python/makeVBBP.py index 05bc6e3..9ab8a66 100644 --- a/python/makeVBBP.py +++ b/python/makeVBBP.py @@ -10,6 +10,8 @@ from ROOT import TMath #as tmath from ZACnC import * +import CMS_lumi, tdrstyle + #from ROOT.TMath import LorentzVector #ROOT.gSystem.Load('CMS_label_C') @@ -18,13 +20,13 @@ def drawTheoryGraph(gr, tb): - x1 = 290 - lenght = 40 + x1 = 160 + lenght = 20 x2 = x1+20 y1 = gr.Eval(x1) y2 = gr.Eval(x2) - pente = math.atan( 420/3.3 *(log10(y2)-log10(y1))/(x2-x1)) + pente = math.atan( 300/4 *(log10(y2)-log10(y1))/(x2-x1)) angle = pente*180/math.pi x3 = x1 + lenght*math.cos(pente) @@ -33,7 +35,7 @@ def drawTheoryGraph(gr, tb): print angle print x1, y1, x3, y3 - size = 5 + size = 10 x_pl = numpy.array([x1-size, x3+size, x3+size, x1-size, x1-size], dtype=float) #[x1+10, x2+10, x2-10, x1-10, x1+10] y_pl = numpy.array([y1+size, y3+size, y3-size, y1-size, y1+size], dtype=float) pline = TPolyLine(5,x_pl,y_pl) @@ -57,7 +59,7 @@ def drawTheoryGraph(gr, tb): ### Definitions ### ################### -mH=500 +mH=300 run_combine = 0 run_blind = 0 @@ -176,11 +178,53 @@ def drawTheoryGraph(gr, tb): gStyle.SetOptStat(0) +''' Cname = "C_"+str(mH) C = TCanvas(Cname,Cname,1200,1200) C.SetLeftMargin(0.2) C.SetBottomMargin(0.2) C.SetTitle(" ") +''' + +#set the tdr style +tdrstyle.setTDRStyle() + +#change the CMS_lumi variables (see CMS_lumi.py) +CMS_lumi.lumi_13TeV = "2.3 fb^{-1}" +CMS_lumi.writeExtraText = 1 +CMS_lumi.extraText = "Preliminary" +CMS_lumi.lumi_sqrtS = "13 TeV" # used with iPeriod = 0, e.g. for simulation-only plots (default is an empty string) + +iPos = 11 +if( iPos==0 ): CMS_lumi.relPosX = 0.12 + +H_ref = 800; +W_ref = 800; +W = W_ref +H = H_ref + +iPeriod = 4 + +# references for T, B, L, R +T = 0.08*H_ref +B = 0.12*H_ref +L = 0.12*W_ref +R = 0.04*W_ref + +canvas = TCanvas("C_"+str(mH),"C_"+str(mH),50,50,W,H) +canvas.SetFillColor(0) +canvas.SetBorderMode(0) +canvas.SetFrameFillStyle(0) +canvas.SetFrameBorderMode(0) +canvas.SetLeftMargin( L/W ) +canvas.SetRightMargin( R/W ) +canvas.SetTopMargin( T/H ) +canvas.SetBottomMargin( B/H ) +canvas.SetTickx(0) +canvas.SetTicky(0) + +#draw the lumi text on the canvas +CMS_lumi.CMS_lumi(canvas, iPeriod, iPos) TGAS_2.SetFillColor(kYellow) TGAS_2.SetFillStyle(1001) @@ -209,8 +253,8 @@ def drawTheoryGraph(gr, tb): title = "M_{H} = "+str(int(mH)) TGAS_2.SetTitle(title) TGAS_2.GetXaxis().SetTitle("M_{A}") -TGAS_2.GetXaxis().SetTitleSize(0.045) -TGAS_2.GetXaxis().SetTitleOffset(0.8) +#TGAS_2.GetXaxis().SetTitleSize(0.045) +#TGAS_2.GetXaxis().SetTitleOffset(0.8) #TGAS_2.SetMinimum(0) TGAS_2.Draw('A E3') TGAS_1.Draw('E3') @@ -223,17 +267,25 @@ def drawTheoryGraph(gr, tb): myTG_8TeV.Draw("L") TGAS_2.GetXaxis().SetTitle("m_{A} [GeV]") -TGAS_2.GetYaxis().SetTitle("#sigma (pp #rightarrow H) #times BR(H #rightarrow Z(ll)A(bb)) [fb]") +TGAS_2.GetYaxis().SetTitle("#sigma #times BR [fb]") TGAS_2.GetYaxis().SetTitleOffset(1.30) -TGAS_2.SetMinimum(5) -TGAS_2.SetMaximum(10000) +TGAS_2.SetMinimum(10) +TGAS_2.SetMaximum(100000) + + +xAxis = TGAS_2.GetXaxis() +#xAxis.SetNdivisions(6,5,0) + +yAxis = TGAS_2.GetYaxis() +#yAxis.SetNdivisions(6,5,0) +yAxis.SetTitleOffset(1) ## Adding theory curves file_theory = TFile("xsec.root") gr_tb1 = file_theory.Get("Xsec_mH"+str(mH)+"_tb1") -gr_tb1_5 = file_theory.Get("Xsec_mH"+str(mH)+"_tb1_5") +gr_tb1_5 = file_theory.Get("Xsec_mH"+str(mH)+"_tb1.5") gr_tb1.SetLineColor(kGray+1) gr_tb1.SetLineStyle(2) @@ -249,42 +301,14 @@ def drawTheoryGraph(gr, tb): drawTheoryGraph(gr_tb1, 1) drawTheoryGraph(gr_tb1_5, 1.5) -C.Update() +canvas.Update() -''' -x1 = 280 -x2 = x1+20 -y1 = gr_tb1.Eval(x1) -y2 = gr_tb1.Eval(x2) - -x3 = (x1+x2)/2 -y3 = gr_tb1.Eval(x3) - -pente = math.atan( 420/3.3 *(log10(y2)-log10(y1))/(x2-x1)) -#pente = 420/3.3 * math.atan((log10(y2)-log10(y1))/(x2-x1)) - -print y1, y2, x1, x2, (log10(y2)-log10(y1))/(x2-x1), pente - -size = 5 -x_pl = numpy.array([x1+size, x2+size, x2-size, x1-size, x1+size], dtype=float) #[x1+10, x2+10, x2-10, x1-10, x1+10] -y_pl = numpy.array([y1+size, y2+size, y2-size, y1-size, y1+size], dtype=float) -pline = TPolyLine(5,x_pl,y_pl) -pline.SetFillColor(kWhite) -pline.Draw("f") - -l = TLatex(x1,y1,"tan #beta = 1") -l.SetTextSize(20) -l.SetTextFont(43) -l.SetTextAlign(12) -l.SetTextAngle(pente*180/math.pi) -l.SetTextColor(kGray+1) -l.Draw() -''' #gPad().SetLogy(1) leg = TLegend(0.59,0.68,0.89,0.88) leg.SetLineColor(0) leg.SetFillColor(0) +leg.SetShadowColor(0) leg.AddEntry(TGAS_1,"CL_{s} Expected #pm 1 #sigma","F") leg.AddEntry(TGAS_2,"CL_{s} Expected #pm 2 #sigma","F") leg.AddEntry(myTG_2,"CL_{s} Expected","L") @@ -294,17 +318,26 @@ def drawTheoryGraph(gr, tb): gPad.SetLogy(1) +#draw the lumi text on the canvas +CMS_lumi.CMS_lumi(canvas, iPeriod, iPos) + +canvas.cd() +canvas.Update() +canvas.RedrawAxis() +frame = canvas.GetFrame() +frame.Draw() + if run_blind : - C.Print("BBP"+str(mH)+"_blind_v2.pdf") - C.Print("BBP"+str(mH)+"_blind_v2.png") + canvas.Print("BBP"+str(mH)+"_blind_v2.pdf") + canvas.Print("BBP"+str(mH)+"_blind_v2.png") f = TFile("narrow_BBP"+str(mH)+"_blind_v2.root","recreate") else : - C.Print("BBP"+str(mH)+"_v2.pdf") - C.Print("BBP"+str(mH)+"_v2.png") + canvas.Print("BBP"+str(mH)+"_v2.pdf") + canvas.Print("BBP"+str(mH)+"_v2.png") f = TFile("narrow_BBP"+str(mH)+"_v2.root","recreate") -C.Write() +canvas.Write() f.Close() ########### diff --git a/python/makeVBBP_300.py b/python/makeVBBP_300.py new file mode 100644 index 0000000..9ab8a66 --- /dev/null +++ b/python/makeVBBP_300.py @@ -0,0 +1,345 @@ +############### +### imports ### +############### + +import math +import os +import os.path +import numpy as numpy +from ROOT import * +from ROOT import TMath #as tmath +from ZACnC import * + +import CMS_lumi, tdrstyle + +#from ROOT.TMath import LorentzVector + +#ROOT.gSystem.Load('CMS_label_C') +#ROOT.gSystem.Load("DrawCanvas_C") + + +def drawTheoryGraph(gr, tb): + + x1 = 160 + lenght = 20 + x2 = x1+20 + y1 = gr.Eval(x1) + y2 = gr.Eval(x2) + + pente = math.atan( 300/4 *(log10(y2)-log10(y1))/(x2-x1)) + angle = pente*180/math.pi + + x3 = x1 + lenght*math.cos(pente) + y3 = gr.Eval(x3) + + print angle + print x1, y1, x3, y3 + + size = 10 + x_pl = numpy.array([x1-size, x3+size, x3+size, x1-size, x1-size], dtype=float) #[x1+10, x2+10, x2-10, x1-10, x1+10] + y_pl = numpy.array([y1+size, y3+size, y3-size, y1-size, y1+size], dtype=float) + pline = TPolyLine(5,x_pl,y_pl) + pline.SetFillColor(kWhite) + pline.Draw("f") + + l = TLatex(x1,y1,"tan #beta = "+str(tb)) + l.SetTextSize(20) + l.SetTextFont(43) + l.SetTextAlign(12) + l.SetTextAngle(angle) + l.SetTextColor(gr.GetLineColor()) + l.Draw() + SetOwnership( l, False ) + SetOwnership( pline, False ) + + + + +################### +### Definitions ### +################### + +mH=300 + +run_combine = 0 +run_blind = 0 + +TGAS_1 = TGraphAsymmErrors(0) +TGAS_2 = TGraphAsymmErrors(0) + +myTG_0 = TGraph(2) +myTG_1 = TGraph(2) +myTG_2 = TGraph(2) +myTG_3 = TGraph(2) +myTG_4 = TGraph(2) +myTG_5 = TGraph(2) + +myTG_8TeV = TGraph(2) + +myTGraph = TGraph2D(9) +#myTGraph.SetNpx(100) +#myTGraph.SetNpy(100) +myTGraph.SetName("p-value") + +#DataCards_path = "../cards/" +DataCards_path = "CARDS_combined/" +RootFiles_path = "../rootfiles_combined/" + +Signal_path = "/home/fynu/amertens/scratch/cmssw/CMSSW_7_6_3/src/cp3_llbb/ZAAnalysis/" + +options = options_() + +####################### +### Initializations ### +####################### + +n=-1 + +############################ +### Getting 8 TeV limits ### +############################ + + +f = TFile("Limit_XS_eff.root") +h_exp = f.Get("h22_lt") +h_exp.SetDirectory(0) +f.Close() + +####################### +### Running Combine ### +####################### + +for cutkey in options.cut : + + name = "HtoZAtoLLBB_mumu_13TeV" + mbb=float(options.mA_list[cutkey]) + mllbb=float(options.mH_list[cutkey]) + + DataCard = DataCards_path+str(mllbb)+"_"+str(mbb)+"/"+name+".dat" + RootFile = RootFiles_path+str(mllbb)+"_"+str(mbb)+"_"+name+".root" + + if (mbb < mllbb - 90) and (mllbb > 126.0) and mllbb == mH and os.path.exists(DataCard): + print 'entering the loop' + try : + if run_combine == 1 : + # Running combine and moving the output rootfile in the repository + + #combine_cmd = "combine -M ProfileLikelihood --signif -m "+str(int(mbb))+" "+DataCard+" --toysFreq" + #combine_cmd = "combine -M ProfileLikelihood --significance --pvalue -m "+str(int(mbb))+" "+DataCard + if run_blind == 1 : + combine_cmd = "combine -M Asymptotic -m "+str(int(mbb))+" --run=blind "+DataCard + else : + combine_cmd = "combine -M Asymptotic -m "+str(int(mbb))+" "+DataCard + os.system(str(combine_cmd)) + mv_cmd = "mv higgsCombineTest.Asymptotic.mH"+str(int(mbb))+".root "+RootFile + os.system(str(mv_cmd)) + ''' + eff_file = TFile('eff.root','READ') + eff_h = eff_file.Get('eff') + eff = eff_h.Interpolate(mbb,mllbb) + ''' + # Accessing the rootfile, get the p-value and fill a TGraph2D + fList = TFile(str(RootFile)) + mytree = fList.Get("limit") + #for entry in mytree: + mytree.GetEntry(2) + limit=mytree.limit + SignalYields = 1 #eff*2.2 + + print 'mbb : ', mbb, ' limit : ', limit #, 'eff : ', eff + + if limit > 0 and SignalYields > 0: + n+=1 + myTG_2.SetPoint(n,int(mbb),mytree.limit/SignalYields) + TGAS_1.SetPoint(n,mbb,mytree.limit/SignalYields) + TGAS_2.SetPoint(n,mbb,mytree.limit/SignalYields) + exp_limit = mytree.limit/SignalYields + mytree.GetEntry(0) + TGAS_2.SetPointEYlow(n,abs(mytree.limit/SignalYields-exp_limit)) + mytree.GetEntry(4) + TGAS_2.SetPointEYhigh(n,abs(mytree.limit/SignalYields-exp_limit)) + mytree.GetEntry(1) + TGAS_1.SetPointEYlow(n,abs(mytree.limit/SignalYields-exp_limit)) + mytree.GetEntry(3) + TGAS_1.SetPointEYhigh(n,abs(mytree.limit/SignalYields-exp_limit)) + + + mytree.GetEntry(5) + myTG_5.SetPoint(n,int(mbb),mytree.limit/SignalYields) + + limit_8TeV = h_exp.Interpolate(int(mbb),int(mllbb)) + print "8TeV : ", limit_8TeV + myTG_8TeV.SetPoint(n,int(mbb),limit_8TeV) + + #n+=1 + #myTGraph.SetPoint(n, mbb, mllbb, mytree.limit) + except: + print 'no background events' + +gStyle.SetOptStat(0) + +''' +Cname = "C_"+str(mH) +C = TCanvas(Cname,Cname,1200,1200) +C.SetLeftMargin(0.2) +C.SetBottomMargin(0.2) +C.SetTitle(" ") +''' + +#set the tdr style +tdrstyle.setTDRStyle() + +#change the CMS_lumi variables (see CMS_lumi.py) +CMS_lumi.lumi_13TeV = "2.3 fb^{-1}" +CMS_lumi.writeExtraText = 1 +CMS_lumi.extraText = "Preliminary" +CMS_lumi.lumi_sqrtS = "13 TeV" # used with iPeriod = 0, e.g. for simulation-only plots (default is an empty string) + +iPos = 11 +if( iPos==0 ): CMS_lumi.relPosX = 0.12 + +H_ref = 800; +W_ref = 800; +W = W_ref +H = H_ref + +iPeriod = 4 + +# references for T, B, L, R +T = 0.08*H_ref +B = 0.12*H_ref +L = 0.12*W_ref +R = 0.04*W_ref + +canvas = TCanvas("C_"+str(mH),"C_"+str(mH),50,50,W,H) +canvas.SetFillColor(0) +canvas.SetBorderMode(0) +canvas.SetFrameFillStyle(0) +canvas.SetFrameBorderMode(0) +canvas.SetLeftMargin( L/W ) +canvas.SetRightMargin( R/W ) +canvas.SetTopMargin( T/H ) +canvas.SetBottomMargin( B/H ) +canvas.SetTickx(0) +canvas.SetTicky(0) + +#draw the lumi text on the canvas +CMS_lumi.CMS_lumi(canvas, iPeriod, iPos) + +TGAS_2.SetFillColor(kYellow) +TGAS_2.SetFillStyle(1001) +TGAS_2.SetLineWidth(0) +TGAS_1.SetFillColor(kGreen) +TGAS_1.SetFillStyle(1001) +TGAS_1.SetLineWidth(0) + +myTG_2.SetLineColor(kBlack) +myTG_2.SetLineStyle(9) +myTG_2.SetLineWidth(2) + +myTG_8TeV.SetLineColor(kBlue) +myTG_8TeV.SetLineWidth(2) + +myTG_5.SetLineColor(kBlack) +myTG_5.SetLineWidth(2) +myTG_5.Sort() +TGAS_2.Sort() +TGAS_1.Sort() +myTG_2.Sort() +myTG_8TeV.Sort() + +#gPad.SetLogy() + +title = "M_{H} = "+str(int(mH)) +TGAS_2.SetTitle(title) +TGAS_2.GetXaxis().SetTitle("M_{A}") +#TGAS_2.GetXaxis().SetTitleSize(0.045) +#TGAS_2.GetXaxis().SetTitleOffset(0.8) +#TGAS_2.SetMinimum(0) +TGAS_2.Draw('A E3') +TGAS_1.Draw('E3') +myTG_2.Draw("L*") +myTG_2.SetMarkerStyle(20) +if run_blind == 0 : + myTG_5.Draw("L*") + myTG_5.SetMarkerStyle(20) +if run_blind == 1 : + myTG_8TeV.Draw("L") + +TGAS_2.GetXaxis().SetTitle("m_{A} [GeV]") +TGAS_2.GetYaxis().SetTitle("#sigma #times BR [fb]") +TGAS_2.GetYaxis().SetTitleOffset(1.30) +TGAS_2.SetMinimum(10) +TGAS_2.SetMaximum(100000) + + + +xAxis = TGAS_2.GetXaxis() +#xAxis.SetNdivisions(6,5,0) + +yAxis = TGAS_2.GetYaxis() +#yAxis.SetNdivisions(6,5,0) +yAxis.SetTitleOffset(1) + +## Adding theory curves + +file_theory = TFile("xsec.root") +gr_tb1 = file_theory.Get("Xsec_mH"+str(mH)+"_tb1") +gr_tb1_5 = file_theory.Get("Xsec_mH"+str(mH)+"_tb1.5") + +gr_tb1.SetLineColor(kGray+1) +gr_tb1.SetLineStyle(2) +gr_tb1.SetLineWidth(3) + +gr_tb1_5.SetLineColor(kGray+3) +gr_tb1_5.SetLineStyle(2) +gr_tb1_5.SetLineWidth(3) + +gr_tb1.Draw("C") +gr_tb1_5.Draw("C") + +drawTheoryGraph(gr_tb1, 1) +drawTheoryGraph(gr_tb1_5, 1.5) + +canvas.Update() + +#gPad().SetLogy(1) + +leg = TLegend(0.59,0.68,0.89,0.88) +leg.SetLineColor(0) +leg.SetFillColor(0) +leg.SetShadowColor(0) +leg.AddEntry(TGAS_1,"CL_{s} Expected #pm 1 #sigma","F") +leg.AddEntry(TGAS_2,"CL_{s} Expected #pm 2 #sigma","F") +leg.AddEntry(myTG_2,"CL_{s} Expected","L") +leg.AddEntry(myTG_5,"CL_{s} Observed","L") +#leg.AddEntry(myTG_8TeV,"8 TeV","L") +leg.Draw() + +gPad.SetLogy(1) + +#draw the lumi text on the canvas +CMS_lumi.CMS_lumi(canvas, iPeriod, iPos) + +canvas.cd() +canvas.Update() +canvas.RedrawAxis() +frame = canvas.GetFrame() +frame.Draw() + + +if run_blind : + canvas.Print("BBP"+str(mH)+"_blind_v2.pdf") + canvas.Print("BBP"+str(mH)+"_blind_v2.png") + f = TFile("narrow_BBP"+str(mH)+"_blind_v2.root","recreate") +else : + canvas.Print("BBP"+str(mH)+"_v2.pdf") + canvas.Print("BBP"+str(mH)+"_v2.png") + f = TFile("narrow_BBP"+str(mH)+"_v2.root","recreate") + +canvas.Write() +f.Close() + +########### +### FIN ### +########### diff --git a/python/makeVBBP_800.py b/python/makeVBBP_800.py new file mode 100644 index 0000000..58138e0 --- /dev/null +++ b/python/makeVBBP_800.py @@ -0,0 +1,345 @@ +############### +### imports ### +############### + +import math +import os +import os.path +import numpy as numpy +from ROOT import * +from ROOT import TMath #as tmath +from ZACnC import * + +import CMS_lumi, tdrstyle + +#from ROOT.TMath import LorentzVector + +#ROOT.gSystem.Load('CMS_label_C') +#ROOT.gSystem.Load("DrawCanvas_C") + + +def drawTheoryGraph(gr, tb): + + x1 = 150 + lenght = 100 + x2 = x1+50 + y1 = gr.Eval(x1) + y2 = gr.Eval(x2) + + pente = math.atan( 820/4 *(log10(y2)-log10(y1))/(x2-x1)) + angle = pente*180/math.pi + + x3 = x1 + lenght*math.cos(pente) + y3 = gr.Eval(x3) + + print angle + print x1, y1, x3, y3 + + size = 1 + x_pl = numpy.array([x1-size, x3+size, x3+size, x1-size, x1-size], dtype=float) #[x1+10, x2+10, x2-10, x1-10, x1+10] + y_pl = numpy.array([y1+size, y3+size, y3-size, y1-size, y1+size], dtype=float) + pline = TPolyLine(5,x_pl,y_pl) + pline.SetFillColor(kWhite) + pline.Draw("f") + + l = TLatex(x1,y1,"tan #beta = "+str(tb)) + l.SetTextSize(20) + l.SetTextFont(43) + l.SetTextAlign(12) + l.SetTextAngle(angle) + l.SetTextColor(gr.GetLineColor()) + l.Draw() + SetOwnership( l, False ) + SetOwnership( pline, False ) + + + + +################### +### Definitions ### +################### + +mH=800 + +run_combine = 0 +run_blind = 0 + +TGAS_1 = TGraphAsymmErrors(0) +TGAS_2 = TGraphAsymmErrors(0) + +myTG_0 = TGraph(2) +myTG_1 = TGraph(2) +myTG_2 = TGraph(2) +myTG_3 = TGraph(2) +myTG_4 = TGraph(2) +myTG_5 = TGraph(2) + +myTG_8TeV = TGraph(2) + +myTGraph = TGraph2D(9) +#myTGraph.SetNpx(100) +#myTGraph.SetNpy(100) +myTGraph.SetName("p-value") + +#DataCards_path = "../cards/" +DataCards_path = "CARDS_combined/" +RootFiles_path = "../rootfiles_combined/" + +Signal_path = "/home/fynu/amertens/scratch/cmssw/CMSSW_7_6_3/src/cp3_llbb/ZAAnalysis/" + +options = options_() + +####################### +### Initializations ### +####################### + +n=-1 + +############################ +### Getting 8 TeV limits ### +############################ + + +f = TFile("Limit_XS_eff.root") +h_exp = f.Get("h22_lt") +h_exp.SetDirectory(0) +f.Close() + +####################### +### Running Combine ### +####################### + +for cutkey in options.cut : + + name = "HtoZAtoLLBB_mumu_13TeV" + mbb=float(options.mA_list[cutkey]) + mllbb=float(options.mH_list[cutkey]) + + DataCard = DataCards_path+str(mllbb)+"_"+str(mbb)+"/"+name+".dat" + RootFile = RootFiles_path+str(mllbb)+"_"+str(mbb)+"_"+name+".root" + + if (mbb < mllbb - 90) and (mllbb > 126.0) and mllbb == mH and os.path.exists(DataCard): + print 'entering the loop' + try : + if run_combine == 1 : + # Running combine and moving the output rootfile in the repository + + #combine_cmd = "combine -M ProfileLikelihood --signif -m "+str(int(mbb))+" "+DataCard+" --toysFreq" + #combine_cmd = "combine -M ProfileLikelihood --significance --pvalue -m "+str(int(mbb))+" "+DataCard + if run_blind == 1 : + combine_cmd = "combine -M Asymptotic -m "+str(int(mbb))+" --run=blind "+DataCard + else : + combine_cmd = "combine -M Asymptotic -m "+str(int(mbb))+" "+DataCard + os.system(str(combine_cmd)) + mv_cmd = "mv higgsCombineTest.Asymptotic.mH"+str(int(mbb))+".root "+RootFile + os.system(str(mv_cmd)) + ''' + eff_file = TFile('eff.root','READ') + eff_h = eff_file.Get('eff') + eff = eff_h.Interpolate(mbb,mllbb) + ''' + # Accessing the rootfile, get the p-value and fill a TGraph2D + fList = TFile(str(RootFile)) + mytree = fList.Get("limit") + #for entry in mytree: + mytree.GetEntry(2) + limit=mytree.limit + SignalYields = 1 #eff*2.2 + + print 'mbb : ', mbb, ' limit : ', limit #, 'eff : ', eff + + if limit > 0 and SignalYields > 0: + n+=1 + myTG_2.SetPoint(n,int(mbb),mytree.limit/SignalYields) + TGAS_1.SetPoint(n,mbb,mytree.limit/SignalYields) + TGAS_2.SetPoint(n,mbb,mytree.limit/SignalYields) + exp_limit = mytree.limit/SignalYields + mytree.GetEntry(0) + TGAS_2.SetPointEYlow(n,abs(mytree.limit/SignalYields-exp_limit)) + mytree.GetEntry(4) + TGAS_2.SetPointEYhigh(n,abs(mytree.limit/SignalYields-exp_limit)) + mytree.GetEntry(1) + TGAS_1.SetPointEYlow(n,abs(mytree.limit/SignalYields-exp_limit)) + mytree.GetEntry(3) + TGAS_1.SetPointEYhigh(n,abs(mytree.limit/SignalYields-exp_limit)) + + + mytree.GetEntry(5) + myTG_5.SetPoint(n,int(mbb),mytree.limit/SignalYields) + + limit_8TeV = h_exp.Interpolate(int(mbb),int(mllbb)) + print "8TeV : ", limit_8TeV + myTG_8TeV.SetPoint(n,int(mbb),limit_8TeV) + + #n+=1 + #myTGraph.SetPoint(n, mbb, mllbb, mytree.limit) + except: + print 'no background events' + +gStyle.SetOptStat(0) + +''' +Cname = "C_"+str(mH) +C = TCanvas(Cname,Cname,1200,1200) +C.SetLeftMargin(0.2) +C.SetBottomMargin(0.2) +C.SetTitle(" ") +''' + +#set the tdr style +tdrstyle.setTDRStyle() + +#change the CMS_lumi variables (see CMS_lumi.py) +CMS_lumi.lumi_13TeV = "2.3 fb^{-1}" +CMS_lumi.writeExtraText = 1 +CMS_lumi.extraText = "Preliminary" +CMS_lumi.lumi_sqrtS = "13 TeV" # used with iPeriod = 0, e.g. for simulation-only plots (default is an empty string) + +iPos = 11 +if( iPos==0 ): CMS_lumi.relPosX = 0.12 + +H_ref = 800; +W_ref = 800; +W = W_ref +H = H_ref + +iPeriod = 4 + +# references for T, B, L, R +T = 0.08*H_ref +B = 0.12*H_ref +L = 0.12*W_ref +R = 0.04*W_ref + +canvas = TCanvas("C_"+str(mH),"C_"+str(mH),50,50,W,H) +canvas.SetFillColor(0) +canvas.SetBorderMode(0) +canvas.SetFrameFillStyle(0) +canvas.SetFrameBorderMode(0) +canvas.SetLeftMargin( L/W ) +canvas.SetRightMargin( R/W ) +canvas.SetTopMargin( T/H ) +canvas.SetBottomMargin( B/H ) +canvas.SetTickx(0) +canvas.SetTicky(0) + +#draw the lumi text on the canvas +CMS_lumi.CMS_lumi(canvas, iPeriod, iPos) + +TGAS_2.SetFillColor(kYellow) +TGAS_2.SetFillStyle(1001) +TGAS_2.SetLineWidth(0) +TGAS_1.SetFillColor(kGreen) +TGAS_1.SetFillStyle(1001) +TGAS_1.SetLineWidth(0) + +myTG_2.SetLineColor(kBlack) +myTG_2.SetLineStyle(9) +myTG_2.SetLineWidth(2) + +myTG_8TeV.SetLineColor(kBlue) +myTG_8TeV.SetLineWidth(2) + +myTG_5.SetLineColor(kBlack) +myTG_5.SetLineWidth(2) +myTG_5.Sort() +TGAS_2.Sort() +TGAS_1.Sort() +myTG_2.Sort() +myTG_8TeV.Sort() + +#gPad.SetLogy() + +title = "M_{H} = "+str(int(mH)) +TGAS_2.SetTitle(title) +TGAS_2.GetXaxis().SetTitle("M_{A}") +#TGAS_2.GetXaxis().SetTitleSize(0.045) +#TGAS_2.GetXaxis().SetTitleOffset(0.8) +#TGAS_2.SetMinimum(0) +TGAS_2.Draw('A E3') +TGAS_1.Draw('E3') +myTG_2.Draw("L*") +myTG_2.SetMarkerStyle(20) +if run_blind == 0 : + myTG_5.Draw("L*") + myTG_5.SetMarkerStyle(20) +if run_blind == 1 : + myTG_8TeV.Draw("L") + +TGAS_2.GetXaxis().SetTitle("m_{A} [GeV]") +TGAS_2.GetYaxis().SetTitle("#sigma #times BR [fb]") +TGAS_2.GetYaxis().SetTitleOffset(1.30) +TGAS_2.SetMinimum(1) +TGAS_2.SetMaximum(10000) + + + +xAxis = TGAS_2.GetXaxis() +#xAxis.SetNdivisions(6,5,0) + +yAxis = TGAS_2.GetYaxis() +#yAxis.SetNdivisions(6,5,0) +yAxis.SetTitleOffset(1) + +## Adding theory curves + +file_theory = TFile("xsec.root") +gr_tb1 = file_theory.Get("Xsec_mH"+str(mH)+"_tb1") +gr_tb1_5 = file_theory.Get("Xsec_mH"+str(mH)+"_tb1.5") + +gr_tb1.SetLineColor(kGray+1) +gr_tb1.SetLineStyle(2) +gr_tb1.SetLineWidth(3) + +gr_tb1_5.SetLineColor(kGray+3) +gr_tb1_5.SetLineStyle(2) +gr_tb1_5.SetLineWidth(3) + +gr_tb1.Draw("C") +gr_tb1_5.Draw("C") + +drawTheoryGraph(gr_tb1, 1) +drawTheoryGraph(gr_tb1_5, 1.5) + +canvas.Update() + +#gPad().SetLogy(1) + +leg = TLegend(0.59,0.68,0.89,0.88) +leg.SetLineColor(0) +leg.SetFillColor(0) +leg.SetShadowColor(0) +leg.AddEntry(TGAS_1,"CL_{s} Expected #pm 1 #sigma","F") +leg.AddEntry(TGAS_2,"CL_{s} Expected #pm 2 #sigma","F") +leg.AddEntry(myTG_2,"CL_{s} Expected","L") +leg.AddEntry(myTG_5,"CL_{s} Observed","L") +#leg.AddEntry(myTG_8TeV,"8 TeV","L") +leg.Draw() + +gPad.SetLogy(1) + +#draw the lumi text on the canvas +CMS_lumi.CMS_lumi(canvas, iPeriod, iPos) + +canvas.cd() +canvas.Update() +canvas.RedrawAxis() +frame = canvas.GetFrame() +frame.Draw() + + +if run_blind : + canvas.Print("BBP"+str(mH)+"_blind_v2.pdf") + canvas.Print("BBP"+str(mH)+"_blind_v2.png") + f = TFile("narrow_BBP"+str(mH)+"_blind_v2.root","recreate") +else : + canvas.Print("BBP"+str(mH)+"_v2.pdf") + canvas.Print("BBP"+str(mH)+"_v2.png") + f = TFile("narrow_BBP"+str(mH)+"_v2.root","recreate") + +canvas.Write() +f.Close() + +########### +### FIN ### +###########