Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

unblind style of plotting tools update #11

Open
wants to merge 3 commits into
base: XtoYH_ggbb
Choose a base branch
from
Open

unblind style of plotting tools update #11

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b50cecd
plot bkg fit bands
YonsiG Oct 22, 2024
15bc7a9
update comments about making fit vs bkg comparison
YonsiG Oct 23, 2024
7cdd661
Observed limit plotting
YonsiG Oct 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 29 additions & 19 deletions Plots/makeSplusBModelPlot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand All @@ -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')
Expand 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
Expand Down Expand Up @@ -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")
Expand All @@ -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)
Expand Down Expand Up @@ -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:
Expand All @@ -211,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)):
Expand All @@ -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:
Expand Down Expand Up @@ -263,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)):
Expand Down Expand Up @@ -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
Expand All @@ -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:
Expand All @@ -358,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<blindingRegion[1]): continue
bval, berr = h_data.GetBinContent(ibin), h_data.GetBinError(ibin)
bkgval = h_bpdf['nBins'].GetBinContent(ibin)
h_data_ratio.SetBinContent(ibin,bval-bkgval)
h_data_ratio.SetBinError(ibin,berr)
bkgerr = h_bpdf['nBins'].GetBinError(ibin)
if berr!=0:
h_data_ratio.SetBinContent(ibin,(bval-bkgval)/berr)
h_data_ratio.SetBinError(ibin,berr)
if opt.doCatWeights:
h_wbpdf_ratio = h_wbpdf['pdfNBins']-h_wbpdf['pdfNBins']
h_wspdf_ratio = h_wspdf['pdfNBins'].Clone()
Expand Down Expand Up @@ -420,7 +430,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
Expand Down
43 changes: 21 additions & 22 deletions Plots/makeToysHH.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand All @@ -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
Expand All @@ -43,24 +42,24 @@ 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()

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")
Expand All @@ -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")
Expand All @@ -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"
Loading