From b50cecd713db2379cff5df522f65c4ddae843994 Mon Sep 17 00:00:00 2001 From: YonsiG Date: Tue, 22 Oct 2024 06:46:36 -0700 Subject: [PATCH 1/3] plot bkg fit bands --- Plots/makeSplusBModelPlot.py | 32 ++++++++++------- Plots/makeToysHH.py | 43 ++++++++++++----------- Plots/make_plots.sh | 66 ++++++++++++++++++++---------------- Plots/plottingToolsHH.py | 34 +++++++++---------- 4 files changed, 94 insertions(+), 81 deletions(-) diff --git a/Plots/makeSplusBModelPlot.py b/Plots/makeSplusBModelPlot.py index 4bda2c88..a99f5168 100644 --- a/Plots/makeSplusBModelPlot.py +++ b/Plots/makeSplusBModelPlot.py @@ -26,10 +26,11 @@ def leave(): def get_options(): parser = OptionParser() parser.add_option("--inputWSFile", dest="inputWSFile", default=None, help="Input RooWorkspace file. If loading snapshot then use a post-fit workspace where the option --saveWorkspace was set") + parser.add_option("--inputSpecialFile", dest="inputSpecialFile", default=None, help="Input RooWorkspace file. If loading snapshot then use a post-fit workspace where the option --saveWorkspace was set") parser.add_option("--loadSnapshot", dest="loadSnapshot", default=None, help="Load best-fit snapshot name") parser.add_option("--cats", dest="cats", default=None, help="Analysis categories. all = loop over cats and plot sum") - parser.add_option("--unblind", dest="unblind", default=False, action="store_true", help="Unblind signal region") - parser.add_option("--blindingRegion", dest="blindingRegion", default="116,134", help="Region in xvar to blind") + parser.add_option("--unblind", dest="unblind", default=False, action="store_true", help="Unblind signal region") + parser.add_option("--blindingRegion", dest="blindingRegion", default="120,130", help="Region in xvar to blind") # plot with --unblind to avoid using parser.add_option("--dataScaler", dest="dataScaler", default=1., type='float', help="Scaling term for data histogram") parser.add_option("--doBkgRenormalization", dest="doBkgRenormalization", default=False, action="store_true", help="Do Bkg renormalization") parser.add_option("--doBands", dest="doBands", default=False, action="store_true", help="Do +-1/2sigma bands for bkg model") @@ -46,8 +47,8 @@ def get_options(): parser.add_option("--ext", dest="ext", default='', help="Extension for saving") parser.add_option("--mass", dest="mass", default=None, help="Higgs mass") parser.add_option("--xvar", dest="xvar", default="CMS_hgg_mass,m_{#gamma#gamma},GeV", help="X-variable: name,title,units") - parser.add_option("--nBins", dest="nBins", default=80, type='int', help="Number of bins") - parser.add_option("--pdfNBins", dest="pdfNBins", default=3200, type='int', help="Number of bins") + parser.add_option("--nBins", dest="nBins", default=160, type='int', help="Number of bins") + parser.add_option("--pdfNBins", dest="pdfNBins", default=160, type='int', help="Number of bins") parser.add_option("--translateCats", dest="translateCats", default=None, help="JSON to store cat translations") parser.add_option("--translatePOIs", dest="translatePOIs", default=None, help="JSON to store poi translations") parser.add_option("--problematicCats", dest="problematicCats", default='', help='Problematic analysis categories to skip when processing all') @@ -58,8 +59,10 @@ def get_options(): # Open WS if opt.inputWSFile is not None: print " --> Opening workspace: %s"%opt.inputWSFile - f = ROOT.TFile(opt.inputWSFile) + f = ROOT.TFile(opt.inputSpecialFile) + f2 = ROOT.TFile(opt.inputWSFile) w = f.Get("w") + w2 = f2.Get("w") # If required loadSnapshot if opt.loadSnapshot is not None: print " * Loading snapshot: %s"%opt.loadSnapshot @@ -106,6 +109,7 @@ def get_options(): # Extract the total SB/B models sb_model, b_model = w.pdf("model_s"), w.pdf("model_b") +nonres_model = w2.pdf("model_b") # Extract dataset for opt.cats d_obs = w.data("data_obs") @@ -114,7 +118,7 @@ def get_options(): for cidx in range(chan.numTypes()): chan.setIndex(cidx) c = chan.getLabel() - if(opt.cats!='all')&(c not in opt.cats.split(",")): continue + if opt.cats != 'all' and (opt.cats not in c or "cr" in c): continue if( opt.doHHMjjFix )&( c in catsfix ): _xvar_argset, _wxvar_argset = xvarfix_argset, wxvarfix_argset else: _xvar_argset, _wxvar_argset = xvar_argset, wxvar_argset data_cats[c] = ROOT.RooDataSet("d_%s"%c,"d_%s"%c,_xvar_argset) @@ -187,7 +191,7 @@ def get_options(): # Loop over bins and add entry for each "weight" to cat datasets for i in range(d_obs.numEntries()): p = d_obs.get(i) - if(opt.cats!='all')&(p.getCatLabel("CMS_channel") not in opt.cats.split(",")): continue + if(opt.cats!='all')&(opt.cats not in p.getCatLabel("CMS_channel") or "cr" in p.getCatLabel("CMS_channel")): continue nent = int(d_obs.weight()) for ient in range(nent): data_cats[p.getCatLabel("CMS_channel")].add(p) if opt.doCatWeights: @@ -235,7 +239,7 @@ def get_options(): for cidx in range(chan.numTypes()): chan.setIndex(cidx) c = chan.getLabel() - if( opt.cats == 'all' )|( c in opt.cats.split(",") ): + if( opt.cats == 'all' )|(opt.cats in c and "cr" not in c): if( opt.doHHMjjFix )&( c in catsfix ): _xvar, _xvar_arglist = xvarfix, xvarfix_arglist else: @@ -311,12 +315,16 @@ def get_options(): # Extract pdfs for category and create histograms print " * creating pdf histograms: S+B, B" sbpdf, bpdf = sb_model.getPdf(c), b_model.getPdf(c) + nonrespdf = nonres_model.getPdf(c) h_sbpdf = {'pdfNBins':sbpdf.createHistogram("h_sb_pdfNBins_%s"%c,_xvar,ROOT.RooFit.Binning(opt.pdfNBins,xvar.getMin(),xvar.getMax())), 'nBins':sbpdf.createHistogram("h_sb_nBins_%s"%c,_xvar,ROOT.RooFit.Binning(opt.nBins,xvar.getMin(),xvar.getMax())) } h_bpdf = {'pdfNBins':bpdf.createHistogram("h_b_pdfNBins_%s"%c,_xvar,ROOT.RooFit.Binning(opt.pdfNBins,xvar.getMin(),xvar.getMax())), 'nBins':bpdf.createHistogram("h_b_nBins_%s"%c,_xvar,ROOT.RooFit.Binning(opt.nBins,xvar.getMin(),xvar.getMax())) } + h_nonrespdf = {'pdfNBins':nonrespdf.createHistogram("h_b_pdfNBins_%s"%c,_xvar,ROOT.RooFit.Binning(opt.pdfNBins,xvar.getMin(),xvar.getMax())), + 'nBins':nonrespdf.createHistogram("h_b_nBins_%s"%c,_xvar,ROOT.RooFit.Binning(opt.nBins,xvar.getMin(),xvar.getMax())) + } # Calculate yields SB, B = sbpdf.expectedEvents(_xvar_argset), bpdf.expectedEvents(_xvar_argset) S = SB-B @@ -332,11 +340,11 @@ def get_options(): else: print " * Yield for category: S = %.2f, B=%.2f"%(S,B) # Extract signal pdf + signalScale = 1 print " * creating pdf histogram: S" - h_spdf = {'pdfNBins':h_sbpdf['pdfNBins']-h_bpdf['pdfNBins'], - 'nBins':h_sbpdf['nBins']-h_bpdf['nBins'] + h_spdf = {'pdfNBins':(h_sbpdf['pdfNBins']-h_bpdf['pdfNBins'])*signalScale, + 'nBins':(h_sbpdf['nBins']-h_bpdf['nBins'])*signalScale } - # Scale pdf histograms to match binning used xvar_range = int(xvar.getBinning().highBound()-xvar.getBinning().lowBound()) if opt.nBins != xvar_range: @@ -420,7 +428,7 @@ def get_options(): if not opt.skipIndividualCatPlots: print " * making plot" if not os.path.isdir("./SplusBModels%s"%(opt.ext)): os.system("mkdir ./SplusBModels%s"%(opt.ext)) - if opt.doBands: makeSplusBPlot(w,h_data,h_sbpdf,h_bpdf,h_spdf,h_data_ratio,h_bpdf_ratio,h_spdf_ratio,c,opt,df_bands,_reduceRange) + if opt.doBands: makeSplusBPlot(w,h_data,h_sbpdf,h_bpdf,h_spdf,h_nonrespdf,h_data_ratio,h_bpdf_ratio,h_spdf_ratio,c,opt,df_bands,_reduceRange) else: makeSplusBPlot(w,h_data,h_sbpdf,h_bpdf,h_spdf,h_data_ratio,h_bpdf_ratio,h_spdf_ratio,c,opt,None,_reduceRange) # Delete histograms diff --git a/Plots/makeToysHH.py b/Plots/makeToysHH.py index 732069fb..13e9fdd5 100644 --- a/Plots/makeToysHH.py +++ b/Plots/makeToysHH.py @@ -11,7 +11,6 @@ def get_options(): parser.add_option("--inputWSFile", dest="inputWSFile", default=None, help="Input RooWorkspace file. If loading snapshot then use a post-fit workspace where the option --saveWorkspace was set") parser.add_option("--loadSnapshot", dest="loadSnapshot", default=None, help="Load best-fit snapshot name") parser.add_option("--MX", dest="MX", default=None, help="Resonance mass") - parser.add_option("--MY", dest="MY", default=None, help="Resonance mass") parser.add_option("--ext", dest="ext", default='', help="Extension for saving") parser.add_option("--dropResonantBkg", dest="dropResonantBkg", default=False, action="store_true") parser.add_option("--nToys", dest="nToys", default=500, type='int', help="Number of toys") @@ -23,12 +22,12 @@ def get_options(): (opt,args) = get_options() # Create jobs directory -if not os.path.isdir("./SplusBModels_%s"%(opt.ext)): os.system("mkdir ./SplusBModels_%s"%(opt.ext)) -if not os.path.isdir("./SplusBModels_%s/toys"%(opt.ext)): os.system("mkdir ./SplusBModels_%s/toys"%(opt.ext)) -if not os.path.isdir("./SplusBModels_%s/toys/jobs"%(opt.ext)): os.system("mkdir ./SplusBModels_%s/toys/jobs"%(opt.ext)) +if not os.path.isdir("./SplusBModels%s"%(opt.ext)): os.system("mkdir ./SplusBModels%s"%(opt.ext)) +if not os.path.isdir("./SplusBModels%s/toys"%(opt.ext)): os.system("mkdir ./SplusBModels%s/toys"%(opt.ext)) +if not os.path.isdir("./SplusBModels%s/toys/jobs"%(opt.ext)): os.system("mkdir ./SplusBModels%s/toys/jobs"%(opt.ext)) # Delete all old jobs -for job in glob.glob("./SplusBModels_%s/toys/jobs/sub*.sh"%opt.ext): os.system("rm %s"%job) +for job in glob.glob("./SplusBModels%s/toys/jobs/sub*.sh"%opt.ext): os.system("rm %s"%job) # Open workspace and extract best fit mass and signal strength inputWSFile = opt.inputWSFile @@ -43,14 +42,14 @@ def get_options(): setParamStr += "%s=%.3f,"%(p,v) setParam0Str += "%s=0,"%p if opt.dropResonantBkg: - setParamStr += "res_bkg_scaler=0," - setParam0Str += "res_bkg_scaler=0," + setParamStr += "singleH=0," + setParamStr += "doubleH=0," + setParam0Str += "singleH=0," + setParam0Str += "doubleH=0," + if opt.MX is not None: setParamStr += "MX=%s,"%opt.MX setParam0Str += "MX=%s,"%opt.MX -if opt.MY is not None: - setParamStr += "MY=%s,"%opt.MY - setParam0Str += "MY=%s,"%opt.MY setParamStr = setParamStr[:-1] setParam0Str = setParam0Str[:-1] mh_bf = w.var("MH").getVal() @@ -58,9 +57,9 @@ def get_options(): if opt.batch == 'IC': # Create submission file for itoy in range(0,opt.nToys): - fsub = open("./SplusBModels_%s/toys/jobs/sub_toy_%g.sh"%(opt.ext,itoy),'w') + fsub = open("./SplusBModels%s/toys/jobs/sub_toy_%g.sh"%(opt.ext,itoy),'w') fsub.write("#!/bin/bash\n\n") - fsub.write("cd %s/src/flashggFinalFit/Plots/SplusBModels_%s/toys\n\n"%(os.environ['CMSSW_BASE'],opt.ext)) + fsub.write("cd %s/src/flashggFinalFit/Plots/SplusBModels%s/toys\n\n"%(os.environ['CMSSW_BASE'],opt.ext)) fsub.write("eval `scramv1 runtime -sh`\n\n") # Generate command fsub.write("#Generate command\n") @@ -85,20 +84,20 @@ def get_options(): fsub.close() # Change permission of all files and set running on batch - os.system("chmod 775 ./SplusBModels_%s/toys/jobs/sub*.sh"%opt.ext) + os.system("chmod 775 ./SplusBModels%s/toys/jobs/sub*.sh"%opt.ext) if not opt.dryRun: - subs = glob.glob("./SplusBModels_%s/toys/jobs/sub*"%opt.ext) + subs = glob.glob("./SplusBModels%s/toys/jobs/sub*"%opt.ext) for fsub in subs: os.system("qsub -q hep.q -l h_rt=0:20:0 -l h_vmem=12G %s"%fsub) else: print " --> [DRY-RUN] jobs have not been submitted" elif opt.batch == 'condor': # Create submission file - fsub = open("./SplusBModels_%s/toys/jobs/sub_toys.sh"%opt.ext,'w') + fsub = open("./SplusBModels%s/toys/jobs/sub_toys.sh"%opt.ext,'w') fsub.write("#!/bin/bash\n") fsub.write("ulimit -s unlimited\n") fsub.write("set -e\n") - fsub.write("cd %s/src/flashggFinalFit/Plots/SplusBModels_%s/toys\n"%(os.environ['CMSSW_BASE'],opt.ext)) + fsub.write("cd %s/src/flashggFinalFit/Plots/SplusBModels%s/toys\n"%(os.environ['CMSSW_BASE'],opt.ext)) fsub.write("source /cvmfs/cms.cern.ch/cmsset_default.sh\n") fsub.write("eval `scramv1 runtime -sh`\n\n") fsub.write("itoy=$1\n\n") @@ -123,18 +122,18 @@ def get_options(): fsub.close() # Write condor submission file - fcondor = open("./SplusBModels_%s/toys/jobs/sub_toys.sub"%opt.ext,'w') + fcondor = open("./SplusBModels%s/toys/jobs/sub_toys.sub"%opt.ext,'w') fcondor.write("executable = sub_toys.sh\n") fcondor.write("arguments = $(ProcId)\n") - fcondor.write("output = %s/src/flashggFinalFit/Plots/SplusBModels_%s/toys/jobs/sub_toys.$(ClusterId).$(ProcId).out\n"%(os.environ['CMSSW_BASE'],opt.ext)) - fcondor.write("error = %s/src/flashggFinalFit/Plots/SplusBModels_%s/toys/jobs/sub_toys.$(ClusterId).$(ProcId).err\n"%(os.environ['CMSSW_BASE'],opt.ext)) - fcondor.write("log = %s/src/flashggFinalFit/Plots/SplusBModels_%s/toys/jobs/sub_toys.$(ClusterId).log\n\n"%(os.environ['CMSSW_BASE'],opt.ext)) + fcondor.write("output = %s/src/flashggFinalFit/Plots/SplusBModels%s/toys/jobs/sub_toys.$(ClusterId).$(ProcId).out\n"%(os.environ['CMSSW_BASE'],opt.ext)) + fcondor.write("error = %s/src/flashggFinalFit/Plots/SplusBModels%s/toys/jobs/sub_toys.$(ClusterId).$(ProcId).err\n"%(os.environ['CMSSW_BASE'],opt.ext)) + fcondor.write("log = %s/src/flashggFinalFit/Plots/SplusBModels%s/toys/jobs/sub_toys.$(ClusterId).log\n\n"%(os.environ['CMSSW_BASE'],opt.ext)) fcondor.write("on_exit_hold = (ExitBySignal == True) || (ExitCode != 0)\n") fcondor.write("periodic_release = (NumJobStarts < 3) && ((CurrentTime - EnteredCurrentStatus) > 600)\n\n") fcondor.write("+JobFlavour = \"%s\"\n"%opt.queue) fcondor.write("queue %s\n"%opt.nToys) # Submission - os.system("chmod 775 ./SplusBModels_%s/toys/jobs/sub_toys.sh"%opt.ext) - if not opt.dryRun: os.system("cd ./SplusBModels_%s/toys/jobs; source /cvmfs/cms.cern.ch/cmsset_default.sh; eval `scramv1 runtime -sh`; condor_submit sub_toys.sub; cd ../../.."%opt.ext) + os.system("chmod 775 ./SplusBModels%s/toys/jobs/sub_toys.sh"%opt.ext) + if not opt.dryRun: os.system("cd ./SplusBModels%s/toys/jobs; source /cvmfs/cms.cern.ch/cmsset_default.sh; eval `scramv1 runtime -sh`; condor_submit sub_toys.sub; cd ../../.."%opt.ext) else: print " --> [DRY-RUN] jobs have not been submitted" diff --git a/Plots/make_plots.sh b/Plots/make_plots.sh index f302fd8b..72a9d02e 100644 --- a/Plots/make_plots.sh +++ b/Plots/make_plots.sh @@ -1,46 +1,52 @@ #!/usr/bin/env bash -set -x -#source /cvmfs/cms.cern.ch/cmsset_default.sh -#source /vols/grid/cms/setup.sh +tag=mx1000my125 -tag=XtoYH_mx360my90 - -#cmsenv -#source setup.sh - -nToys=10 -#nToys=500 -#Need run toys after generating a Datacard.root file with only the signals used +nToys=5000 make_toys(){ pushd Plots - rm -rf SplusBModels_$tag - python makeToysHH.py --inputWSFile ../Combine/Datacard_ggbbres_mx360my90_ggbbres.root --ext $tag --dryRun --nToys $nToys --dropResonantBkg + rm -rf SplusBModels$tag + #python makeToysHH.py --inputWSFile /home/users/yagu/XYH/XtoYH_sys/plots_2023-12-19_CategorizationAndPreselSFCorrected/CMSSW_10_2_13/src/flashggFinalFit/higgsCombine_MultiDimFit_r_ggbbres_mx1000my125mh125.MultiDimFit.mH125.root --ext $tag --dryRun --nToys $nToys + python makeToysHH.py --inputWSFile /home/users/yagu/XYH/XtoYH_sys/plots_2023-12-19_CategorizationAndPreselSFCorrected/CMSSW_10_2_13/src/flashggFinalFit/higgsCombine_MultiDimFit_r_ggbbres_mx1000my125mh125.MultiDimFit.mH125.root --ext $tag --dryRun --nToys $nToys #--loadSnapshot MultiDimFit #--dropResonantBkg + #python makeToysHH.py --inputWSFile /home/users/iareed/CMSSW_10_2_13/src/flashggFinalFit/Combine/Datacard_ttHHggXX.root --ext $tag --dryRun --nToys $nToys #--dropResonantBkg iter=0 while [ $iter -lt $nToys ] - do - ./SplusBModels_${tag}/toys/jobs/sub_toy_$(($iter+1)).sh & - ./SplusBModels_${tag}/toys/jobs/sub_toy_$(($iter+2)).sh & - ./SplusBModels_${tag}/toys/jobs/sub_toy_$(($iter+3)).sh & - ./SplusBModels_${tag}/toys/jobs/sub_toy_$(($iter+4)).sh & - ./SplusBModels_${tag}/toys/jobs/sub_toy_$(($iter+5)).sh & - ./SplusBModels_${tag}/toys/jobs/sub_toy_$(($iter+6)).sh & - ./SplusBModels_${tag}/toys/jobs/sub_toy_$(($iter+7)).sh & - ./SplusBModels_${tag}/toys/jobs/sub_toy_$(($iter+8)).sh & - ./SplusBModels_${tag}/toys/jobs/sub_toy_$(($iter+9)).sh & - ./SplusBModels_${tag}/toys/jobs/sub_toy_$iter.sh + do + ./SplusBModels${tag}/toys/jobs/sub_toy_$(($iter+1)).sh & + ./SplusBModels${tag}/toys/jobs/sub_toy_$(($iter+2)).sh & + ./SplusBModels${tag}/toys/jobs/sub_toy_$(($iter+3)).sh & + ./SplusBModels${tag}/toys/jobs/sub_toy_$(($iter+4)).sh & + ./SplusBModels${tag}/toys/jobs/sub_toy_$(($iter+5)).sh & + ./SplusBModels${tag}/toys/jobs/sub_toy_$(($iter+6)).sh & + ./SplusBModels${tag}/toys/jobs/sub_toy_$(($iter+7)).sh & + ./SplusBModels${tag}/toys/jobs/sub_toy_$(($iter+8)).sh & + ./SplusBModels${tag}/toys/jobs/sub_toy_$(($iter+9)).sh & + ./SplusBModels${tag}/toys/jobs/sub_toy_$iter.sh iter=$(($iter+10)) done - popd + popd } -#Should be run with a Datacard.root file with everything make_SpB(){ - pushd Plots - python makeSplusBModelPlot.py --inputWSFile ../Combine/Datacard_ggbbres_mx360my90_ggbbres.root --cat "all" --doBands --ext - popd + pushd Plots + #python makeSplusBModelPlot.py --inputWSFile /home/users/yagu/XYH/XtoYH_sys/plots_2023-12-19_CategorizationAndPreselSFCorrected/CMSSW_10_2_13/src/flashggFinalFit/higgsCombine_MultiDimFit_r_ggbbres_mx1000my125mh125.MultiDimFit.mH125.root --inputSpecialFile /home/users/yagu/XYH/XtoYH_sys/plots_2023-12-19_CategorizationAndPreselSFCorrected/CMSSW_10_2_13/src/flashggFinalFit/higgsCombine_MultiDimFit_r_ggbbres_mx1000my125mh125.MultiDimFit.mH125.root --cat "cat0" --doBands --ext $tag --unblind # --parameterMap "singleH:1,doubleH:1" + python makeSplusBModelPlot.py --inputWSFile /home/users/yagu/XYH/XtoYH_sys/plots_2023-12-19_CategorizationAndPreselSFCorrected/CMSSW_10_2_13/src/flashggFinalFit/higgsCombine_MultiDimFit_r_ggbbres_mx1000my125mh125.MultiDimFit.mH125.root --inputSpecialFile /home/users/yagu/XYH/XtoYH_sys/plots_2023-12-19_CategorizationAndPreselSFCorrected/CMSSW_10_2_13/src/flashggFinalFit/higgsCombine_MultiDimFit_r_ggbbres_mx1000my125mh125.MultiDimFit.mH125.root --cat "cat0" --doBands --ext $tag --unblind #--loadSnapshot MultiDimFit # --parameterMap "singleH:1,doubleH:1" + ##python makeSplusBModelPlot.py --inputWSFile /home/users/iareed/CMSSW_10_2_13/src/flashggFinalFit/Combine/Datacard_ttHHggXX.root --cat "SR1" --doBands --ext $tag # --parameterMap "singleH:1,doubleH:1" + #mv SplusBModels${tag}/SR1_CMS_hgg_mass.pdf SplusBModels${tag}/SR1_CMS_hgg_mass_blind.pdf + #mv SplusBModels${tag}/SR1_CMS_hgg_mass.png SplusBModels${tag}/SR1_CMS_hgg_mass_blind.png + #mv SplusBModels${tag}/SR2_CMS_hgg_mass.pdf SplusBModels${tag}/SR2_CMS_hgg_mass_blind.pdf + #mv SplusBModels${tag}/SR2_CMS_hgg_mass.png SplusBModels${tag}/SR2_CMS_hgg_mass_blind.png + + ##python makeSplusBModelPlot.py --inputWSFile /home/users/iareed/CMSSW_10_2_13/src/flashggFinalFit/Combine/Datacard_ttHHggXX.root --cat "SR1" --doBands --ext $tag --unblind # --parameterMap "singleH:1,doubleH:1" + ##python makeSplusBModelPlot.py --inputWSFile /home/users/iareed/CMSSW_10_2_13/src/flashggFinalFit/Combine/Datacard_Tprime_M500.root --cat "SR1" --doBands --ext $tag --unblind # --parameterMap "singleH:1,doubleH:1" + #mv SplusBModels${tag}/SR1_CMS_hgg_mass.pdf SplusBModels${tag}/SR1_CMS_hgg_mass_unblind.pdf + #mv SplusBModels${tag}/SR1_CMS_hgg_mass.png SplusBModels${tag}/SR1_CMS_hgg_mass_unblind.png + #mv SplusBModels${tag}/SR2_CMS_hgg_mass.pdf SplusBModels${tag}/SR2_CMS_hgg_mass_unblind.pdf + #mv SplusBModels${tag}/SR2_CMS_hgg_mass.png SplusBModels${tag}/SR2_CMS_hgg_mass_unblind.png + #cp SplusBModels${tag}/* /home/users/iareed/public_html/ttHH/flashggFinalFit/${tag}/Bands/ + popd } make_toys -#make_SpB +make_SpB diff --git a/Plots/plottingToolsHH.py b/Plots/plottingToolsHH.py index ea1ed1af..61bf8d0e 100644 --- a/Plots/plottingToolsHH.py +++ b/Plots/plottingToolsHH.py @@ -77,9 +77,9 @@ def extractBandProperties(data,category,bidx): props['down2sigma'] = np.percentile(data['%s_%g'%(c,bidx)].values,50*(1+math.erf(-2./math.sqrt(2)))) return props -def makeSplusBPlot(workspace,hD,hSB,hB,hS,hDr,hBr,hSr,cat,options,dB=None,reduceRange=None,limit=None): - translateCats = {} if options.translateCats is None else LoadTranslations(options.translateCats) - translatePOIs = {} if options.translatePOIs is None else LoadTranslations(options.translatePOIs) +def makeSplusBPlot(workspace,hD,hSB,hB,hS,hNonRes,hDr,hBr,hSr,cat,options,dB=None,reduceRange=None,limit=None): + translateCats = {} #if options.translateCats is None else LoadTranslations(options.translateCats) + translatePOIs = {} #if options.translatePOIs is None else LoadTranslations(options.translatePOIs) blindingRegion = [float(options.blindingRegion.split(",")[0]),float(options.blindingRegion.split(",")[1])] if reduceRange is not None: for h in [hD,hDr,hBr,hSr]: h.GetXaxis().SetRangeUser(reduceRange[0],reduceRange[1]) @@ -166,7 +166,7 @@ def makeSplusBPlot(workspace,hD,hSB,hB,hS,hDr,hBr,hSr,cat,options,dB=None,reduce leg.SetTextSize(0.045) leg.AddEntry(hD,"Data","ep") if options.unblind: - leg.AddEntry(hSB['pdfNBins'],"S+B fit","l") +# leg.AddEntry(hSB['pdfNBins'],"S+B fit","l") leg.AddEntry(hB['pdfNBins'],"B component","l") else: leg.AddEntry(hB['pdfNBins'],"B fit","l") @@ -177,12 +177,12 @@ def makeSplusBPlot(workspace,hD,hSB,hB,hS,hDr,hBr,hSr,cat,options,dB=None,reduce leg.Draw("Same") # Set pdf style if options.unblind: - hSB['pdfNBins'].SetLineWidth(3) - hSB['pdfNBins'].SetLineColor(2) - hSB['pdfNBins'].Draw("Hist same c") +# hSB['pdfNBins'].SetLineWidth(3) +# hSB['pdfNBins'].SetLineColor(2) +# hSB['pdfNBins'].Draw("Hist same c") +# hSB['pdfNBins'].SetLineStyle(2) hB['pdfNBins'].SetLineWidth(3) hB['pdfNBins'].SetLineColor(2) - hB['pdfNBins'].SetLineStyle(2) hB['pdfNBins'].Draw("Hist same c") else: hS['pdfNBins'].SetLineWidth(3) @@ -211,9 +211,9 @@ def makeSplusBPlot(workspace,hD,hSB,hB,hS,hDr,hBr,hSr,cat,options,dB=None,reduce lat0.DrawLatex(0.6,0.8,"#scale[0.6]{%s}"%Translate(cat,translateCats)) #lat0.DrawLatex(0.15,0.83,"#scale[0.75]{H#rightarrow#gamma#gamma}") #lat0.DrawLatex(0.15,0.83,"#scale[0.75]{Non-resonant HH #rightarrow #gamma#gamma#tau#tau}") - if(options.loadSnapshot is not None): +# if(options.loadSnapshot is not None): #lat0.DrawLatex(0.15,0.77,"#scale[0.6]{#vec{#alpha} = STXS stage 1.2 minimal}") - lat0.DrawLatex(0.15,0.77,"#scale[0.6]{#vec{#alpha} = (#mu_{ggH}, #mu_{VBF}, #mu_{VH}, #mu_{top})}") + #lat0.DrawLatex(0.15,0.77,"#scale[0.6]{#vec{#alpha} = (#mu_{ggH}, #mu_{VBF}, #mu_{VH}, #mu_{top})}") #lat0.DrawLatex(0.15,0.77,"#scale[0.5]{(#hat{#mu}_{ggH},#hat{#mu}_{VBF},#hat{#mu}_{VH},#hat{#mu}_{top}) = (1.07,1.04,1.34,1.35)}") #lat0.DrawLatex(0.15,0.77,"#scale[0.75]{#hat{#mu} = 1.03}") #muhat_ggh, muhat_vbf, muhat_vh, muhat_top, mhhat = workspace.var("r_ggH").getVal(), workspace.var("r_VBF").getVal(), workspace.var("r_VH").getVal(), workspace.var("r_top").getVal(), workspace.var("MH").getVal() @@ -227,14 +227,14 @@ def makeSplusBPlot(workspace,hD,hSB,hB,hS,hDr,hBr,hSr,cat,options,dB=None,reduce # poiStr += ' %s = %.1f,'%(Translate(k,translatePOIs),float(v)) # poiStr = poiStr[:-1] # lat0.DrawLatex(0.13,0.77,"#scale[0.75]{%s}"%poiStr) - else: - if limit is not None: lat0.DrawLatex(0.15,0.77,"#scale[0.5]{m_{X} = %s GeV, #sigma #bf{#it{#Beta}} = %.3f fb}"%(options.MX,limit)) - else: lat0.DrawLatex(0.15,0.77,"#scale[0.75]{m_{H} = 125.38 GeV}") +# else: +# if limit is not None: lat0.DrawLatex(0.15,0.77,"#scale[0.5]{m_{X} = %s GeV, #sigma #bf{#it{#Beta}} = %.3f fb}"%(options.MX,limit)) # Ratio plot pad2.cd() h_axes_ratio = hDr.Clone() h_axes_ratio.Reset() - h_axes_ratio.SetMaximum(max((hDr.GetMaximum()+hDr.GetBinError(hDr.GetMaximumBin()))*1.5,hSr.GetMaximum()*1.2)) +### h_axes_ratio.SetMaximum(max((hDr.GetMaximum()+hDr.GetBinError(hDr.GetMaximumBin()))*1.5,hSr.GetMaximum()*1.2)) + h_axes_ratio.SetMaximum((hDr.GetMaximum()+hDr.GetBinError(hDr.GetMaximumBin()))*1.5) h_axes_ratio.SetMinimum((hDr.GetMinimum()-hDr.GetBinError(hDr.GetMinimumBin()))*1.3) h_axes_ratio.SetTitle("") h_axes_ratio.GetXaxis().SetTitleSize(0.05*padSizeRatio) @@ -251,9 +251,9 @@ def makeSplusBPlot(workspace,hD,hSB,hB,hS,hDr,hBr,hSr,cat,options,dB=None,reduce gr_1sig_r.Draw("LE3SAME") # Set pdf style if options.unblind: - hSr.SetLineWidth(3) - hSr.SetLineColor(2) - hSr.Draw("Hist same c") + ### hSr.SetLineWidth(3) + ### hSr.SetLineColor(2) + ### hSr.Draw("Hist same c") hBr.SetLineWidth(3) hBr.SetLineStyle(2) hBr.SetLineColor(2) From 15bc7a99f99a283236f54b0743d50cac73f0e044 Mon Sep 17 00:00:00 2001 From: YonsiG Date: Wed, 23 Oct 2024 08:57:34 -0700 Subject: [PATCH 2/3] update comments about making fit vs bkg comparison --- Plots/makeSplusBModelPlot.py | 16 +++++++++------- Plots/plottingToolsHH.py | 12 +++++++++--- 2 files changed, 18 insertions(+), 10 deletions(-) diff --git a/Plots/makeSplusBModelPlot.py b/Plots/makeSplusBModelPlot.py index a99f5168..cd435d55 100644 --- a/Plots/makeSplusBModelPlot.py +++ b/Plots/makeSplusBModelPlot.py @@ -215,9 +215,9 @@ def get_options(): # Create dataframe df_bands = pd.DataFrame(columns=_columns) # Loop over toys file and add row for each toy dataset - toyFiles = glob.glob("./SplusBModels%s/toys/toy_*.root"%opt.ext) + toyFiles = glob.glob("./postfit/SplusBModels%s/toys/toy_*.root"%opt.ext) if len(toyFiles) == 0: - print " * [ERROR] No toys files of form ./SplusBModels%s/toys/toy_*.root. Skipping bands"%opt.ext + print " * [ERROR] No toys files of form ./postfit/SplusBModels%s/toys/toy_*.root. Skipping bands"%opt.ext opt.doBands = False else: for tidx in range(len(toyFiles)): @@ -267,8 +267,8 @@ def get_options(): else: print " --> Toy veto: zero entries in first bin" # Savin toy yields dataframe to pickle file if opt.saveToyYields: - print " * Saving toy yields to: SplusBModels%s/toyYields_%s.pkl"%(opt.ext,opt.xvar.split(",")[0]) - with open("SplusBModels%s/toyYields_%s.pkl"%(opt.ext,opt.xvar.split(",")[0]),"w") as fD: pickle.dump(df_bands,fD) + print " * Saving toy yields to: postfit/SplusBModels%s/toyYields_%s.pkl"%(opt.ext,opt.xvar.split(",")[0]) + with open("postfit/SplusBModels%s/toyYields_%s.pkl"%(opt.ext,opt.xvar.split(",")[0]),"w") as fD: pickle.dump(df_bands,fD) # Process each category separately for cidx in range(len(cats)): @@ -366,14 +366,16 @@ def get_options(): h_bpdf_ratio = h_bpdf['pdfNBins']-h_bpdf['pdfNBins'] h_spdf_ratio = h_spdf['pdfNBins'].Clone() h_data_ratio = h_data.Clone() - h_data_ratio.Reset() + #h_data_ratio.Reset() for ibin in range(1,h_data.GetNbinsX()+1): bcenter = h_data.GetBinCenter(ibin) if(not opt.unblind)&(bcenter>blindingRegion[0])&(bcenter Date: Wed, 23 Oct 2024 08:57:55 -0700 Subject: [PATCH 3/3] Observed limit plotting --- plot_limits.py | 179 +++++++++++++++++++++++++++++-------------------- 1 file changed, 105 insertions(+), 74 deletions(-) diff --git a/plot_limits.py b/plot_limits.py index ee9da220..4f72ff93 100644 --- a/plot_limits.py +++ b/plot_limits.py @@ -10,6 +10,7 @@ import tabulate import pandas as pd import os +from optparse import OptionParser import scipy.interpolate as spi @@ -26,6 +27,13 @@ # "MY": [70, 250, 550, 800, 70, 100, 70, 70, 125, 300, 90, 400, 70, 100, 300, 250, 500, 70, 190, 450, 550, 70, 500, 600, 70, 170, 400, 650, 70, 80, 190, 650, 700, 70, 500], # "limit":[0.0559,0.1305,0.1922,0.4891,0.7906,0.8094,0.3141,0.1508,0.3344,0.7500,0.1289,0.6969,0.0871,0.1262,0.3438,0.2656,0.6969,0.0699,0.1953,0.3734,0.7312,0.0645,0.3750,0.7719,0.0633,0.1414,0.2227,0.7469,0.1063,0.0621,0.1266,0.3562,0.4969,0.0535,0.1914]}) +def get_options(): + parser = OptionParser() + parser.add_option("--inputFile", dest="inputFile", default=None, help="Input limit result file") + parser.add_option("--outputFile", dest="outputFile", default=None, help="Output folder to save results") + parser.add_option("--doObserved", dest="doObserved", default=False, action="store_true", help="do observed limit") + return parser.parse_args() +(opt,args) = get_options() def getLimits(results_path): with open(results_path, "r") as f: @@ -45,6 +53,7 @@ def getLimits(results_path): limits_no_sys = np.zeros((5, len(masses))) limits_no_res_bkg = np.zeros((5, len(masses))) limits_no_dy_bkg = np.zeros((5, len(masses))) + limits_observed = np.zeros(len(masses)) for line in results: m = line.split(".")[0].split("_")[-1] @@ -66,8 +75,10 @@ def getLimits(results_path): limit = float(line.split("r < ")[1]) - if "no_" not in line: + if "no_" not in line and "Expected" in line: limits[idx2][idx1] = limit + elif "no_" not in line and "Observed" in line: + limits_observed[idx1] = limit else: if "no_sys" in line: limits_no_sys[idx2][idx1] = limit @@ -76,9 +87,6 @@ def getLimits(results_path): elif "no_dy_bkg" in line: limits_no_dy_bkg[idx2][idx1] = limit - #print(limits[2]) - #print(limits_no_sys[2]) - masses = np.array(masses) #sort out scan over mh (mgg) if len(np.unique(np.array(masses)[:,2])) != 1: #if more than 125 in mh @@ -100,12 +108,7 @@ def getLimits(results_path): limits_no_res_bkg = np.delete(limits_no_res_bkg, to_delete, axis=1) limits_no_dy_bkg = np.delete(limits_no_dy_bkg, to_delete, axis=1) - # masses[:,1] = masses[:,2] #set my to be mh - - #masses = masses[:,:2] - #print(masses) - - return masses, limits, limits_no_sys, limits_no_res_bkg, limits_no_dy_bkg + return masses, limits, limits_no_sys, limits_no_res_bkg, limits_no_dy_bkg, limits_observed def plotLimits(mX, limits, ylabel, nominal_masses, savename=None, xlabel=r"$m_X$"): plt.scatter(mX, limits[2], zorder=3, facecolors="none", edgecolors="blue") @@ -129,17 +132,22 @@ def plotLimits(mX, limits, ylabel, nominal_masses, savename=None, xlabel=r"$m_X$ plt.savefig(savename+"_log.pdf") plt.clf() -def plotLimitsStackMX(masses, limits, ylabel, nominal_mx, nominal_my, savename): +def plotLimitsStackMX(masses, limits, ylabel, nominal_mx, nominal_my, savename, doObserved = False, limits_observed = None): label1 = "Nominal masses" label2 = "Expected 95% CL limit" label3 = r"$\pm$ $1\sigma$" label4 = r"$\pm$ $2\sigma$" + label5 = "Observed 95% CL limits" for i, mx in enumerate(np.sort(np.unique(masses[:,0]))): my = masses[masses[:,0]==mx,1] limits_slice = limits[:,masses[:,0]==mx] - limits_slice = limits_slice[:,np.argsort(my)] + + if doObserved: + limits_observed_slice = limits_observed[masses[:,0]==mx] + limits_observed_slice = limits_observed_slice[np.argsort(my)] + my = my[np.argsort(my)] limits_slice *= 10**i @@ -153,7 +161,26 @@ def plotLimitsStackMX(masses, limits, ylabel, nominal_mx, nominal_my, savename): label1 = label2 = label3 = label4 = None plt.text(my[-1]+10, limits_slice[2][-1], r"$m_X=%d$ GeV $(\times 10^{%d})$"%(mx, i), fontsize=12, verticalalignment="center") - + + if doObserved: + limits_observed_slice *= 10**i + if mx == 600: + # Apply the mask only when mx is 600 to exclude my=400 + mask = my != 400 + else: + # No mask applied, include all points + mask = np.ones_like(my, dtype=bool) + + plt.scatter(my[mask], limits_observed_slice[mask], zorder=3, facecolors="none", edgecolors="black") + #plt.scatter(my, limits_observed_slice, zorder=3, facecolors="none", edgecolors="black") + if mx in nominal_mx: + plt.scatter(my[np.isin(my, nominal_my) & mask], + limits_observed_slice[np.isin(my, nominal_my) & mask], + zorder=5, facecolors="black", edgecolors="black") + #plt.scatter(my[np.isin(my, nominal_my)], limits_observed_slice[np.isin(my, nominal_my)], zorder=5, facecolors="black", edgecolors="black") + #plt.plot(my, limits_observed_slice, 'k-', zorder=5, label=label5) + plt.plot(my[mask], limits_observed_slice[mask], 'k-', zorder=5, label=label5) + label5 = None plt.xlabel(r"$m_Y$") plt.ylabel(ylabel) @@ -259,9 +286,6 @@ def plotLimits2D(masses, limits, ylabel, savename): plt.savefig(savename+"_exclude.pdf") #s = limits[2] < max_allowed_values # s = limits[2] -# print(limits) -# print(s) -# print(masses[s,0]) # plt.scatter(masses[s,0], masses[s,1], marker='x', color="r", label="Limit below maximally allowed in NMSSM") # plt.savefig(savename+"_exclude_points.png") @@ -375,28 +399,30 @@ def tabulateLimitsAll(masses, limits, limits_no_sys, limits_no_res_bkg, path): f.write(df.to_latex(float_format="%.4f", index=False)) df.to_csv(os.path.join(path, "limits.csv"), float_format="%.4f") -masses, limits, limits_no_sys, limits_no_res_bkg, limits_no_dy_bkg = getLimits(sys.argv[1]) -print(masses) -os.makedirs(os.path.join(sys.argv[2], "Limits_xs_br"), exist_ok=True) -os.makedirs(os.path.join(sys.argv[2], "Limits_xs"), exist_ok=True) -os.makedirs(os.path.join(sys.argv[2], "Limits_xs_br_no_sys"), exist_ok=True) -os.makedirs(os.path.join(sys.argv[2], "Limits_xs_no_sys"), exist_ok=True) -os.makedirs(os.path.join(sys.argv[2], "Limits_xs_br_no_res_bkg"), exist_ok=True) -os.makedirs(os.path.join(sys.argv[2], "Limits_xs_no_res_bkg"), exist_ok=True) -os.makedirs(os.path.join(sys.argv[2], "Limits_systematics_comparison"), exist_ok=True) -os.makedirs(os.path.join(sys.argv[2], "Limits_res_bkg_comparison"), exist_ok=True) -os.makedirs(os.path.join(sys.argv[2], "Limits_dy_bkg_comparison"), exist_ok=True) +masses, limits, limits_no_sys, limits_no_res_bkg, limits_no_dy_bkg, limits_observed = getLimits(opt.inputFile) +os.makedirs(os.path.join(opt.outputFile, "Limits_xs_br"), exist_ok=True) +os.makedirs(os.path.join(opt.outputFile, "Limits_xs"), exist_ok=True) +if not opt.doObserved: + os.makedirs(os.path.join(opt.outputFile, "Limits_xs_br_no_sys"), exist_ok=True) + os.makedirs(os.path.join(opt.outputFile, "Limits_xs_no_sys"), exist_ok=True) + os.makedirs(os.path.join(opt.outputFile, "Limits_xs_br_no_res_bkg"), exist_ok=True) + os.makedirs(os.path.join(opt.outputFile, "Limits_xs_no_res_bkg"), exist_ok=True) + os.makedirs(os.path.join(opt.outputFile, "Limits_systematics_comparison"), exist_ok=True) + os.makedirs(os.path.join(opt.outputFile, "Limits_res_bkg_comparison"), exist_ok=True) + os.makedirs(os.path.join(opt.outputFile, "Limits_dy_bkg_comparison"), exist_ok=True) -tabulateLimits(masses, limits, os.path.join(sys.argv[2], "Limits_xs_br")) -tabulateLimitsAll(masses, limits, limits_no_sys, limits_no_res_bkg, os.path.join(sys.argv[2], "Limits_xs_br")) -tabulateLimits(masses, limits / BR_HH_GGBB, os.path.join(sys.argv[2], "Limits_xs")) +tabulateLimits(masses, limits, os.path.join(opt.outputFile, "Limits_xs_br")) +if not opt.doObserved: + tabulateLimitsAll(masses, limits, limits_no_sys, limits_no_res_bkg, os.path.join(opt.outputFile, "Limits_xs_br")) +tabulateLimits(masses, limits / BR_HH_GGBB, os.path.join(opt.outputFile, "Limits_xs")) -tabulateLimits(masses, limits_no_sys, os.path.join(sys.argv[2], "Limits_xs_br_no_sys")) -tabulateLimits(masses, limits_no_sys / BR_HH_GGBB, os.path.join(sys.argv[2], "Limits_xs_no_sys")) +if not opt.doObserved: + tabulateLimits(masses, limits_no_sys, os.path.join(opt.outputFile, "Limits_xs_br_no_sys")) + tabulateLimits(masses, limits_no_sys / BR_HH_GGBB, os.path.join(opt.outputFile, "Limits_xs_no_sys")) -tabulateLimits(masses, limits_no_res_bkg, os.path.join(sys.argv[2], "Limits_xs_br_no_res_bkg")) -tabulateLimits(masses, limits_no_res_bkg / BR_HH_GGBB, os.path.join(sys.argv[2], "Limits_xs_no_res_bkg")) + tabulateLimits(masses, limits_no_res_bkg, os.path.join(opt.outputFile, "Limits_xs_br_no_res_bkg")) + tabulateLimits(masses, limits_no_res_bkg / BR_HH_GGBB, os.path.join(opt.outputFile, "Limits_xs_no_res_bkg")) if len(np.unique(masses[:,1])) == 1: #if 1D (graviton or radion) mx = masses[:,0] @@ -404,31 +430,32 @@ def tabulateLimitsAll(masses, limits, limits_no_sys, limits_no_res_bkg, path): limits_no_sys = limits_no_sys[:,np.argsort(mx)] limits_no_res_bkg = limits_no_res_bkg[:,np.argsort(mx)] limits_no_dy_bkg = limits_no_dy_bkg[:,np.argsort(mx)] - + limits_observed = limits_observed[np.argsort(mx)] mx = mx[np.argsort(mx)] nominal_masses = [260,270,280,290,300,320,350,400,450,500,550,600,650,700,750,800,900,1000] ylabel = r"$\sigma(pp \rightarrow X) B(X \rightarrow HH \rightarrow \gamma\gamma bb)$ [fb]" - plotLimits(mx, limits, ylabel, nominal_masses, os.path.join(sys.argv[2], "Limits_xs_br", "limits")) - plotLimits(mx, limits_no_sys, ylabel, nominal_masses, os.path.join(sys.argv[2], "Limits_xs_br_no_sys", "limits_no_sys")) - plotLimits(mx, limits_no_res_bkg, ylabel, nominal_masses, os.path.join(sys.argv[2], "Limits_xs_br_no_res_bkg", "limits_no_res_bkg")) + if opt.doObserved: + plotLimits(mx, limits, limits_observed, ylabel, nominal_masses, os.path.join(opt.outputFile, "Limits_xs_br", "limits")) + else: + plotLimits(mx, limits, ylabel, nominal_masses, os.path.join(opt.outputFile, "Limits_xs_br", "limits")) + plotLimits(mx, limits_no_sys, ylabel, nominal_masses, os.path.join(opt.outputFile, "Limits_xs_br_no_sys", "limits_no_sys")) + plotLimits(mx, limits_no_res_bkg, ylabel, nominal_masses, os.path.join(opt.outputFile, "Limits_xs_br_no_res_bkg", "limits_no_res_bkg")) ylabel = r"$\sigma(pp \rightarrow X) B(X \rightarrow HH)$ [fb]" - plotLimits(mx, limits / BR_HH_GGBB, ylabel, nominal_masses, os.path.join(sys.argv[2], "Limits_xs", "limits")) - plotLimits(mx, limits_no_sys / BR_HH_GGBB, ylabel, nominal_masses, os.path.join(sys.argv[2], "Limits_xs_no_sys", "limits_no_sys")) - plotLimits(mx, limits_no_res_bkg / BR_HH_GGBB, ylabel, nominal_masses, os.path.join(sys.argv[2], "Limits_xs_no_res_bkg", "limits_no_res_bkg")) + if opt.doObserved: + plotLimits(mx, limits / BR_HH_GGBB, limits_observed / BR_HH_GGBB, ylabel, nominal_masses, os.path.join(opt.outputFile, "Limits_xs", "limits")) + else: + plotLimits(mx, limits_no_sys / BR_HH_GGBB, ylabel, nominal_masses, os.path.join(opt.outputFile, "Limits_xs_no_sys", "limits_no_sys")) + plotLimits(mx, limits_no_res_bkg / BR_HH_GGBB, ylabel, nominal_masses, os.path.join(opt.outputFile, "Limits_xs_no_res_bkg", "limits_no_res_bkg")) - plotSystematicComparison(mx, limits, limits_no_sys, nominal_masses, os.path.join(sys.argv[2], "Limits_systematics_comparison", "125")) - ylabel = r"$\sigma(pp \rightarrow X) B(X \rightarrow HH \rightarrow \gamma\gamma bb)$ [fb]" - plotSystematicComparison2(mx, limits, limits_no_sys, ylabel, nominal_masses, os.path.join(sys.argv[2], "Limits_systematics_comparison", "125_2")) - plotResBkgComparison2(mx, limits, limits_no_res_bkg, ylabel, nominal_masses, os.path.join(sys.argv[2], "Limits_res_bkg_comparison", "125_2")) + plotSystematicComparison(mx, limits, limits_no_sys, nominal_masses, os.path.join(opt.outputFile, "Limits_systematics_comparison", "125")) + ylabel = r"$\sigma(pp \rightarrow X) B(X \rightarrow HH \rightarrow \gamma\gamma bb)$ [fb]" + plotSystematicComparison2(mx, limits, limits_no_sys, ylabel, nominal_masses, os.path.join(opt.outputFile, "Limits_systematics_comparison", "125_2")) + plotResBkgComparison2(mx, limits, limits_no_res_bkg, ylabel, nominal_masses, os.path.join(opt.outputFile, "Limits_res_bkg_comparison", "125_2")) else: - #nominal_mx = [300,400,500,600,700,800,900,1000] - #nominal_my = [70,80,90,100,125] - #nominal_my = [70,80,90,100,125,150,200,250,300,400,500,600,700,800] - #nominal_my = [125,150,200,250,300,400,500,600,700,800] nominal_mx = [240,280,300,320,360,400,450,500,550,600,650,700,750,800,850,900,950,1000] nominal_my = [70,80,90,100,125,150,170,190,250,300,350,400,450,500,550,600,650,700,800] @@ -438,15 +465,19 @@ def tabulateLimitsAll(masses, limits, limits_no_sys, limits_no_res_bkg, path): # limits_no_sys = limits_no_sys[:, s] # masses = masses[s] ylabel = r"$\sigma(pp \rightarrow X) B(X \rightarrow YH \rightarrow \gamma\gamma bb)$ [fb]" - plotLimitsStackMX(masses, limits, ylabel, nominal_mx, nominal_my, os.path.join(sys.argv[2], "Limits_xs_br", "limits_stack_mx")) - plotLimitsStackMX(masses, limits_no_sys, ylabel, nominal_mx, nominal_my, os.path.join(sys.argv[2], "Limits_xs_br_no_sys", "limits_stack_mx_no_sys")) - plotLimitsStackMX(masses, limits_no_res_bkg, ylabel, nominal_mx, nominal_my, os.path.join(sys.argv[2], "Limits_xs_br_no_res_bkg", "limits_stack_mx_no_res_bkg")) -# plotLimitsStackMY(masses, limits, ylabel, nominal_mx, nominal_my, os.path.join(sys.argv[2], "Limits_xs_br", "limits_stack_my")) -# plotLimitsStackMY(masses, limits_no_sys, ylabel, nominal_mx, nominal_my, os.path.join(sys.argv[2], "Limits_xs_br_no_sys", "limits_stack_my_no_sys")) -# plotLimitsStackMY(masses, limits_no_res_bkg, ylabel, nominal_mx, nominal_my, os.path.join(sys.argv[2], "Limits_xs_br_no_res_bkg", "limits_stack_my_no_res_bkg")) - plotLimits2D(masses, limits, ylabel, os.path.join(sys.argv[2], "Limits_xs_br", "limits_2d")) - plotLimits2D(masses, limits_no_sys, ylabel, os.path.join(sys.argv[2], "Limits_xs_br_no_sys", "limits_2d_no_sys")) - plotLimits2D(masses, limits_no_res_bkg, ylabel, os.path.join(sys.argv[2], "Limits_xs_br_no_res_bkg", "limits_2d_no_res_bkg")) + if opt.doObserved: + plotLimitsStackMX(masses, limits, ylabel, nominal_mx, nominal_my, os.path.join(opt.outputFile, "Limits_xs_br", "limits_stack_mx"), opt.doObserved, limits_observed) + else: + plotLimitsStackMX(masses, limits, ylabel, nominal_mx, nominal_my, os.path.join(opt.outputFile, "Limits_xs_br", "limits_stack_mx"), opt.doObserved) + plotLimitsStackMX(masses, limits_no_sys, ylabel, nominal_mx, nominal_my, os.path.join(opt.outputFile, "Limits_xs_br_no_sys", "limits_stack_mx_no_sys"), opt.doObserved) + plotLimitsStackMX(masses, limits_no_res_bkg, ylabel, nominal_mx, nominal_my, os.path.join(opt.outputFile, "Limits_xs_br_no_res_bkg", "limits_stack_mx_no_res_bkg"), opt.doObserved) +# plotLimitsStackMY(masses, limits, ylabel, nominal_mx, nominal_my, os.path.join(opt.outputFile, "Limits_xs_br", "limits_stack_my")) +# plotLimitsStackMY(masses, limits_no_sys, ylabel, nominal_mx, nominal_my, os.path.join(opt.outputFile, "Limits_xs_br_no_sys", "limits_stack_my_no_sys")) +# plotLimitsStackMY(masses, limits_no_res_bkg, ylabel, nominal_mx, nominal_my, os.path.join(opt.outputFile, "Limits_xs_br_no_res_bkg", "limits_stack_my_no_res_bkg")) + plotLimits2D(masses, limits, ylabel, os.path.join(opt.outputFile, "Limits_xs_br", "limits_2d")) + if not opt.doObserved: + plotLimits2D(masses, limits_no_sys, ylabel, os.path.join(opt.outputFile, "Limits_xs_br_no_sys", "limits_2d_no_sys")) + plotLimits2D(masses, limits_no_res_bkg, ylabel, os.path.join(opt.outputFile, "Limits_xs_br_no_res_bkg", "limits_2d_no_res_bkg")) for mx in np.unique(masses[:,0]): my = masses[masses[:,0]==mx,1] @@ -466,17 +497,16 @@ def tabulateLimitsAll(masses, limits, limits_no_sys, limits_no_res_bkg, path): else: nm = [] - print(mx) - print(my, limits_slice) ylabel = r"$\sigma(pp \rightarrow X(%d)) B(X \rightarrow YH \rightarrow \gamma\gamma bb)$ [fb]"%mx - plotLimits(my, limits_slice, ylabel, nm, os.path.join(sys.argv[2], "Limits_xs_br", "limits_mx%d"%mx), xlabel=r"$m_Y$") - plotLimits(my, limits_no_sys_slice, ylabel, nm, os.path.join(sys.argv[2], "Limits_xs_br_no_sys", "limits_mx%d_no_sys"%mx), xlabel=r"$m_Y$") - plotLimits(my, limits_no_res_bkg_slice, ylabel, nm, os.path.join(sys.argv[2], "Limits_xs_br_no_res_bkg", "limits_mx%d_no_res_bkg"%mx), xlabel=r"$m_Y$") - plotSystematicComparison(my, limits_slice, limits_no_sys_slice, nm, os.path.join(sys.argv[2], "Limits_systematics_comparison", "mx%d"%mx), xlabel=r"$m_Y$") - plotSystematicComparison2(my, limits_slice, limits_no_sys_slice, ylabel, nm, os.path.join(sys.argv[2], "Limits_systematics_comparison", "mx%d_2"%mx), xlabel=r"$m_Y$") - plotResBkgComparison2(my, limits_slice, limits_no_res_bkg_slice, ylabel, nm, os.path.join(sys.argv[2], "Limits_res_bkg_comparison", "mx%d_2"%mx), xlabel=r"$m_Y$") - plotDYBkgComparison2(my, limits_slice, limits_no_dy_bkg_slice, ylabel, nm, os.path.join(sys.argv[2], "Limits_dy_bkg_comparison", "mx%d_2"%mx), xlabel=r"$m_Y$") + plotLimits(my, limits_slice, ylabel, nm, os.path.join(opt.outputFile, "Limits_xs_br", "limits_mx%d"%mx), xlabel=r"$m_Y$") + if not opt.doObserved: + plotLimits(my, limits_no_sys_slice, ylabel, nm, os.path.join(opt.outputFile, "Limits_xs_br_no_sys", "limits_mx%d_no_sys"%mx), xlabel=r"$m_Y$") + plotLimits(my, limits_no_res_bkg_slice, ylabel, nm, os.path.join(opt.outputFile, "Limits_xs_br_no_res_bkg", "limits_mx%d_no_res_bkg"%mx), xlabel=r"$m_Y$") + plotSystematicComparison(my, limits_slice, limits_no_sys_slice, nm, os.path.join(opt.outputFile, "Limits_systematics_comparison", "mx%d"%mx), xlabel=r"$m_Y$") + plotSystematicComparison2(my, limits_slice, limits_no_sys_slice, ylabel, nm, os.path.join(opt.outputFile, "Limits_systematics_comparison", "mx%d_2"%mx), xlabel=r"$m_Y$") + plotResBkgComparison2(my, limits_slice, limits_no_res_bkg_slice, ylabel, nm, os.path.join(opt.outputFile, "Limits_res_bkg_comparison", "mx%d_2"%mx), xlabel=r"$m_Y$") + plotDYBkgComparison2(my, limits_slice, limits_no_dy_bkg_slice, ylabel, nm, os.path.join(opt.outputFile, "Limits_dy_bkg_comparison", "mx%d_2"%mx), xlabel=r"$m_Y$") for my in np.unique(masses[:,1]): mx = masses[masses[:,1]==my,0] @@ -498,12 +528,13 @@ def tabulateLimitsAll(masses, limits, limits_no_sys, limits_no_res_bkg, path): nm = [] ylabel = r"$\sigma(pp \rightarrow X) B(X \rightarrow Y(%d)H \rightarrow \gamma\gamma bb)$ [fb]"%my - plotLimits(mx, limits_slice, ylabel, nm, os.path.join(sys.argv[2], "Limits_xs_br", "limits_my%d"%my)) - plotLimits(mx, limits_no_sys_slice, ylabel, nm, os.path.join(sys.argv[2], "Limits_xs_br_no_sys", "limits_my%d_no_sys"%my)) - plotLimits(mx, limits_no_res_bkg_slice, ylabel, nm, os.path.join(sys.argv[2], "Limits_xs_br_no_res_bkg", "limits_my%d_no_res_bkg"%my)) - plotSystematicComparison(mx, limits_slice, limits_no_sys_slice, nm, os.path.join(sys.argv[2], "Limits_systematics_comparison", "my%d"%my)) - plotSystematicComparison2(mx, limits_slice, limits_no_sys_slice, ylabel, nm, os.path.join(sys.argv[2], "Limits_systematics_comparison", "my%d_2"%my)) - plotResBkgComparison2(mx, limits_slice, limits_no_res_bkg_slice, ylabel, nm, os.path.join(sys.argv[2], "Limits_res_bkg_comparison", "my%d_2"%my)) - plotDYBkgComparison2(mx, limits_slice, limits_no_dy_bkg_slice, ylabel, nm, os.path.join(sys.argv[2], "Limits_dy_bkg_comparison", "my%d_2"%my)) + plotLimits(mx, limits_slice, ylabel, nm, os.path.join(opt.outputFile, "Limits_xs_br", "limits_my%d"%my)) + if not opt.doObserved: + plotLimits(mx, limits_no_sys_slice, ylabel, nm, os.path.join(opt.outputFile, "Limits_xs_br_no_sys", "limits_my%d_no_sys"%my)) + plotLimits(mx, limits_no_res_bkg_slice, ylabel, nm, os.path.join(opt.outputFile, "Limits_xs_br_no_res_bkg", "limits_my%d_no_res_bkg"%my)) + plotSystematicComparison(mx, limits_slice, limits_no_sys_slice, nm, os.path.join(opt.outputFile, "Limits_systematics_comparison", "my%d"%my)) + plotSystematicComparison2(mx, limits_slice, limits_no_sys_slice, ylabel, nm, os.path.join(opt.outputFile, "Limits_systematics_comparison", "my%d_2"%my)) + plotResBkgComparison2(mx, limits_slice, limits_no_res_bkg_slice, ylabel, nm, os.path.join(opt.outputFile, "Limits_res_bkg_comparison", "my%d_2"%my)) + plotDYBkgComparison2(mx, limits_slice, limits_no_dy_bkg_slice, ylabel, nm, os.path.join(opt.outputFile, "Limits_dy_bkg_comparison", "my%d_2"%my))