From 37c53dfdd4eaacc1cd18b70366148488082966a1 Mon Sep 17 00:00:00 2001 From: mauro verzetti Date: Thu, 30 Mar 2017 03:42:31 -0500 Subject: [PATCH 1/7] added options to freeze nuisances and to run expected only limits --- CombineTools/python/combine/EnhancedCombine.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/CombineTools/python/combine/EnhancedCombine.py b/CombineTools/python/combine/EnhancedCombine.py index 80c85b27c40..cfd43795a03 100755 --- a/CombineTools/python/combine/EnhancedCombine.py +++ b/CombineTools/python/combine/EnhancedCombine.py @@ -38,6 +38,10 @@ def attach_intercept_args(self, group): help='Name used to label the combine output file, can be modified by other options') group.add_argument( '--setPhysicsModelParameterRanges', help='Some other options will modify or add to the list of parameter ranges') + group.add_argument( + '--run', choices=['both', 'observed', 'expected', 'blind', ''], default='', help='what and how to run') + group.add_argument( + '--freezeNuisances', default='', help='coma-separated list of nuisances to freeze') def attach_args(self, group): @@ -116,6 +120,12 @@ def run_method(self): current_ranges = self.args.setPhysicsModelParameterRanges put_back_ranges = current_ranges is not None + if self.args.run: + self.passthru.extend(['--run', self.args.run]) + + if self.args.freezeNuisances: + self.passthru.extend(['--freezeNuisances=%s' % self.args.freezeNuisances]) + if self.args.boundlist is not None: # We definitely don't need to put the parameter ranges back # into the args because they're going in via the boundlist From f06fcb6d6d5f83308d04d9bf6491809ab5833f31 Mon Sep 17 00:00:00 2001 From: mauro verzetti Date: Thu, 30 Mar 2017 03:44:18 -0500 Subject: [PATCH 2/7] allow remapping of the limit value to different interpretations --- CombineTools/python/plotting.py | 26 +++++++++++++------------- CombineTools/scripts/plotLimits.py | 15 ++++++++++++--- 2 files changed, 25 insertions(+), 16 deletions(-) diff --git a/CombineTools/python/plotting.py b/CombineTools/python/plotting.py index 72e2e81227a..96abf1b3a7a 100644 --- a/CombineTools/python/plotting.py +++ b/CombineTools/python/plotting.py @@ -677,21 +677,21 @@ def MakeErrorBand(LowerGraph, UpperGraph): return errorBand -def LimitTGraphFromJSON(js, label): +def LimitTGraphFromJSON(js, label, mapping=lambda x: x): xvals = [] yvals = [] for key in js: xvals.append(float(key)) - yvals.append(js[key][label]) + yvals.append(mapping(js[key][label])) graph = R.TGraph(len(xvals), array('d', xvals), array('d', yvals)) graph.Sort() return graph -def LimitTGraphFromJSONFile(jsfile, label): +def LimitTGraphFromJSONFile(jsfile, label, mapping=lambda x: x): with open(jsfile) as jsonfile: js = json.load(jsonfile) - return LimitTGraphFromJSON(js, label) + return LimitTGraphFromJSON(js, label, mapping) def ToyTGraphFromJSON(js, label): xvals = [] @@ -720,35 +720,35 @@ def ToyTGraphFromJSONFile(jsfile, label): js = json.load(jsonfile) return ToyTGraphFromJSON(js, label) -def LimitBandTGraphFromJSON(js, central, lo, hi): +def LimitBandTGraphFromJSON(js, central, lo, hi, mapping=lambda x: x): xvals = [] yvals = [] yvals_lo = [] yvals_hi = [] for key in js: xvals.append(float(key)) - yvals.append(js[key][central]) - yvals_lo.append(js[key][central] - js[key][lo]) - yvals_hi.append(js[key][hi] - js[key][central]) + yvals.append(mapping(js[key][central])) + yvals_lo.append(mapping(js[key][central]) - mapping(js[key][lo])) + yvals_hi.append(mapping(js[key][hi]) - mapping(js[key][central])) graph = R.TGraphAsymmErrors(len(xvals), array('d', xvals), array('d', yvals), array( 'd', [0]), array('d', [0]), array('d', yvals_lo), array('d', yvals_hi)) graph.Sort() return graph -def StandardLimitsFromJSONFile(json_file, draw=['obs', 'exp0', 'exp1', 'exp2']): +def StandardLimitsFromJSONFile(json_file, draw=['obs', 'exp0', 'exp1', 'exp2'], mapping=lambda x: x): graphs = {} data = {} with open(json_file) as jsonfile: data = json.load(jsonfile) if 'obs' in draw: - graphs['obs'] = LimitTGraphFromJSON(data, 'obs') + graphs['obs'] = LimitTGraphFromJSON(data, 'obs', mapping) if 'exp0' in draw or 'exp' in draw: - graphs['exp0'] = LimitTGraphFromJSON(data, 'exp0') + graphs['exp0'] = LimitTGraphFromJSON(data, 'exp0', mapping) if 'exp1' in draw or 'exp' in draw: - graphs['exp1'] = LimitBandTGraphFromJSON(data, 'exp0', 'exp-1', 'exp+1') + graphs['exp1'] = LimitBandTGraphFromJSON(data, 'exp0', 'exp-1', 'exp+1', mapping) if 'exp2' in draw or 'exp' in draw: - graphs['exp2'] = LimitBandTGraphFromJSON(data, 'exp0', 'exp-2', 'exp+2') + graphs['exp2'] = LimitBandTGraphFromJSON(data, 'exp0', 'exp-2', 'exp+2', mapping) return graphs diff --git a/CombineTools/scripts/plotLimits.py b/CombineTools/scripts/plotLimits.py index 94f39cfed51..f647336a3f6 100755 --- a/CombineTools/scripts/plotLimits.py +++ b/CombineTools/scripts/plotLimits.py @@ -2,6 +2,8 @@ import ROOT import CombineHarvester.CombineTools.plotting as plot import argparse +from math import * +from pdb import set_trace # import CombineHarvester.CombineTools.maketable as maketable parser = argparse.ArgumentParser() @@ -38,11 +40,15 @@ '--ratio-to', default=None) parser.add_argument( '--pad-style', default=None, help="""Extra style options for the pad, e.g. Grid=(1,1)""") +parser.add_argument( + '--grid', action='store_true', help="""Extra style options for the pad, e.g. Grid=(1,1)""") +parser.add_argument( + '--mapping', default='lambda x: x', help="""mapping for the obtained limit values""") parser.add_argument( '--auto-style', nargs='?', const='', default=None, help="""Take line colors and styles from a pre-defined list""") parser.add_argument('--table_vals', help='Amount of values to be written in a table for different masses', default=10) args = parser.parse_args() - +mapping = eval(args.mapping) def DrawAxisHists(pads, axis_hists, def_pad=None): for i, pad in enumerate(pads): @@ -69,7 +75,10 @@ def DrawAxisHists(pads, axis_hists, def_pad=None): # Set the style options of the pads for padx in pads: # Use tick marks on oppsite axis edges - plot.Set(padx, Tickx=1, Ticky=1, Logx=args.logx) + plot.Set(padx, Tickx=1, Ticky=1, Logx=args.logx) + if args.grid: + padx.SetGridx() + padx.SetGridy() if args.pad_style is not None: settings = {x.split('=')[0]: eval(x.split('=')[1]) for x in args.pad_style.split(',')} print 'Applying style options to the TPad(s):' @@ -105,7 +114,7 @@ def DrawAxisHists(pads, axis_hists, def_pad=None): file = splitsrc[0] # limit.json => Draw as full obs + exp limit band if len(splitsrc) == 1: - graph_sets.append(plot.StandardLimitsFromJSONFile(file, args.show.split(','))) + graph_sets.append(plot.StandardLimitsFromJSONFile(file, args.show.split(','), mapping)) if axis is None: axis = plot.CreateAxisHists(len(pads), graph_sets[-1].values()[0], True) DrawAxisHists(pads, axis, pads[0]) From 49b6f29cb2c1cbfedba075b5c075af60834b079e Mon Sep 17 00:00:00 2001 From: mauro verzetti Date: Thu, 30 Mar 2017 03:44:43 -0500 Subject: [PATCH 3/7] refactored code --- CombineTools/scripts/setup_andrey.py | 235 ++++++++++++++------------- 1 file changed, 120 insertions(+), 115 deletions(-) mode change 100644 => 100755 CombineTools/scripts/setup_andrey.py diff --git a/CombineTools/scripts/setup_andrey.py b/CombineTools/scripts/setup_andrey.py old mode 100644 new mode 100755 index 44493d51cd4..c7d7fdb7ac1 --- a/CombineTools/scripts/setup_andrey.py +++ b/CombineTools/scripts/setup_andrey.py @@ -2,6 +2,12 @@ import CombineHarvester.CombineTools.ch as ch import os +from argparse import ArgumentParser + +parser = ArgumentParser() +parser.add_argument('inputfile') +parser.add_argument('limit', choices=['electrons', 'muons', 'cmb'], help='choose leptonic decay type') +args = parser.parse_args() cb = ch.CombineHarvester() @@ -11,123 +17,122 @@ addBBB = False masses = ['400', '500', '600', '750'] +widths = ['5', '10', '25', '50'] # for mode in ['scalar']:#, 'pseudoscalar']: -for mode in ['pseudoscalar']:#, 'pseudoscalar']: - - patterns = ['HToTT-sgn-{mode}-m', 'HToTT-int-{mode}-m', 'HToTT-int-{mode}_neg-m'] - procs = { - # 'sig' : ['S0_{mode}_M'.format(mode=mode), 'S0_neg_{mode}_M'.format(mode=mode)], - 'sig' : [pattern.format(mode=mode) for pattern in patterns], - # 'sim' : ['WZ', 'ZZ', 'ttW', 'ttZ', 'ttH'], - 'bkg' : ['Wjets', 'single-top', 'VV', 'DY', 'ttbar', 'QCD'] - } - - cats = [(0, '1b'), (1, '2b')] - - - - cb.AddObservations(['*'], ['andrey'], ["8TeV"], ['htt'], cats) - cb.AddProcesses(['*'], ['andrey'], ["8TeV"], ['htt'], procs['bkg'], cats, False) - cb.AddProcesses(masses, ['andrey'], ["8TeV"], ['htt'], procs['sig'], cats, True) - # cb.AddProcesses(['*'], ['andrey'], ["8TeV"], ['htt'], procs['sig'], cats, True) - - print '>> Adding systematic uncertainties...' - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'lumi', 'lnN', ch.SystMap()(1.026)) - - cb.cp().process(['VV']).AddSyst( - cb, 'diboson_rate', 'lnN', ch.SystMap()(1.3)) - - cb.cp().process(['ttbar']).AddSyst( - cb, 'ttbar_rate', 'lnN', ch.SystMap()(1.05)) - - cb.cp().process(['single-top']).AddSyst( - cb, 'st_rate', 'lnN', ch.SystMap()(1.05)) - - cb.cp().process(['Wjets']).AddSyst( - cb, 'w_rate', 'lnN', ch.SystMap()(1.03)) - - cb.cp().process(['DY']).AddSyst( - cb, 'dy_rate', 'lnN', ch.SystMap()(1.03)) - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'trigger', 'lnN', ch.SystMap()(1.0115)) - - cb.cp().process(['ttbar']).AddSyst( - cb, 'scaleTT', 'shape', ch.SystMap()(1.)) - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'bTag', 'shape', ch.SystMap()(1.)) - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'jec', 'shape', ch.SystMap()(1.)) - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'jer', 'shape', ch.SystMap()(1.)) - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'pu', 'shape', ch.SystMap()(1.)) - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'muID', 'shape', ch.SystMap()(1.)) - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'eID', 'shape', ch.SystMap()(1.)) - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'metUncl', 'shape', ch.SystMap()(1.)) - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'bMistag', 'shape', ch.SystMap()(1.)) - - cb.cp().process(['Wjets']).AddSyst( - cb, 'scaleW', 'shape', ch.SystMap()(1.)) - - cb.cp().process(['DY']).AddSyst( - cb, 'scaleDY', 'shape', ch.SystMap()(1.)) - - - print '>> Extracting histograms from input root files...' - in_file = aux_shapes + 'andrey_tth.inputs-bsm-8TeV.root' - cb.cp().backgrounds().ExtractShapes( - in_file, '$BIN/$PROCESS', '$BIN/$PROCESS__$SYSTEMATIC') - cb.cp().signals().ExtractShapes( - in_file, '$BIN/$PROCESS$MASS', '$BIN/$PROCESS$MASS__$SYSTEMATIC') - # in_file, '$BIN/$PROCESS', '$BIN/$PROCESS__$SYSTEMATIC') - - - # print '>> Generating bbb uncertainties...' - # bbb = ch.BinByBinFactory() - # bbb.SetAddThreshold(0.1).SetFixNorm(True) - # bbb.AddBinByBin(cb.cp().process(['reducible']), cb) - - - # for mass in masses: - # norm_initial = norms[mode][int(mass)] - - # cb.cp().process(['S0_neg_{mode}_M{mass}'.format(mode=mode, mass=mass), 'S0_{mode}_M{mass}'.format(mode=mode, mass=mass)]).ForEachProc(lambda p: p.set_rate(p.rate()/norm_initial)) - - if addBBB: - bbb = ch.BinByBinFactory().SetAddThreshold(0.).SetFixNorm(False) - bbb.MergeBinErrors(cb.cp().backgrounds()) - bbb.AddBinByBin(cb.cp().backgrounds(), cb) - - print '>> Setting standardised bin names...' - ch.SetStandardBinNames(cb) - cb.PrintAll() - - writer = ch.CardWriter('$TAG/$MASS/$ANALYSIS_$CHANNEL_$BINID.txt', - # writer = ch.CardWriter('$TAG/$ANALYSIS_$CHANNEL_$BINID_$ERA.txt', - '$TAG/$ANALYSIS_$CHANNEL.input.root') - # writer.SetVerbosity(100) - writer.WriteCards('output/{mode}'.format(mode=mode), cb) - print 'Try writing cards...' - # import ROOT - # f_out = ROOT.TFile('andrey_out.root', 'RECREATE') - # cb.WriteDatacard("andrey_out.txt", 'andrey_out.root') - # writer.WriteCards('output/andrey_cards/', cb) +for parity in ['pseudoscalar']:#, 'pseudoscalar']: + for width in widths: + mode='{mode}_{width}pc'.format(mode=parity, width=width) + patterns = [i % width for i in ['ggA_pos-int-%spc-M', 'ggA_neg-int-%spc-M', 'ggA_pos-sgn-%spc-M']] + procs = { + # 'sig' : ['S0_{mode}_M'.format(mode=mode), 'S0_neg_{mode}_M'.format(mode=mode)], + 'sig' : patterns, + # 'sim' : ['WZ', 'ZZ', 'ttW', 'ttZ', 'ttH'], + 'bkg' : ['TT', 'VV', 'TTV', 'tChannel', 'tWChannel', 'WJets', 'ZJets'] + } + if args.limit == 'electrons' or args.limit == 'cmb': + procs['bkg'].append('QCDejets') + if args.limit == 'muons' or args.limit == 'cmb': + procs['bkg'].append('QCDmujets') + + cats = [] + if args.limit == 'muons' or args.limit == 'cmb': + cats.append((len(cats), 'mujets')) + if args.limit == 'electrons' or args.limit == 'cmb': + cats.append((len(cats), 'ejets')) + + cb.AddObservations(['*'], ['mauro'], ["13TeV"], ['htt'], cats) + cb.AddProcesses(['*'], ['mauro'], ["13TeV"], ['htt'], procs['bkg'], cats, False) + cb.AddProcesses(masses, ['mauro'], ["13TeV"], ['htt'], procs['sig'], cats, True) + # cb.AddProcesses(['*'], ['mauro'], ["13TeV"], ['htt'], procs['sig'], cats, True) + + print '>> Adding systematic uncertainties...' + # + # Normalization uncertainty + # + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( + cb, 'lumi_13TeV', 'lnN', ch.SystMap()(1.026)) + + cb.cp().process(['VV']).AddSyst( + cb, 'CMS_httbar_vvXsec_13TeV', 'lnN', ch.SystMap()(1.3)) + + cb.cp().process(['TT']).AddSyst( + cb, 'ttbar_rate', 'lnN', ch.SystMap()(1.05)) + + cb.cp().process(['tChannel']).AddSyst( + cb, 'st_rate', 'lnN', ch.SystMap()(1.05)) + + cb.cp().process(['tWChannel']).AddSyst( + cb, 'st_rate', 'lnN', ch.SystMap()(1.05)) + + cb.cp().process(['WJets']).AddSyst( + cb, 'w_rate', 'lnN', ch.SystMap()(1.03)) + cb.cp().process(['ZJets']).AddSyst( + cb, 'w_rate', 'lnN', ch.SystMap()(1.03)) + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( + cb, 'trigger', 'lnN', ch.SystMap()(1.0115)) + + # + # Shape uncetainties + # + #on everything + for name in [ + 'CMS_scale_j_13TeV', + 'CMS_eff_b_13TeV', + 'CMS_fake_b_13TeV', + 'CMS_pileup', + 'CMS_METunclustered_13TeV' + ]: + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( + cb, name, 'shape', ch.SystMap()(1.)) + + #on TT only + for name in [ + 'QCDscaleMERenorm_TT', + 'QCDscaleMEFactor_TT', + 'QCDscaleMERenormFactor_TT', + 'QCDscalePS_TT', + 'TMass', + 'pdf' + ]: + cb.cp().process(['TT']).AddSyst( + cb, name, 'shape', ch.SystMap()(1.)) + + if args.limit == 'muons' or args.limit == 'cmb': + cb.cp().process(['QCDmujets']).AddSyst( + cb, 'QCD_mu_norm', 'lnN', ch.SystMap()(2.)) + if args.limit == 'electrons' or args.limit == 'cmb': + cb.cp().process(['QCDejets']).AddSyst( + cb, 'QCD_el_norm', 'lnN', ch.SystMap()(2.)) + + + print '>> Extracting histograms from input root files...' + in_file = args.inputfile + cb.cp().backgrounds().ExtractShapes( + in_file, '$BIN/$PROCESS', '$BIN/$PROCESS_$SYSTEMATIC') + cb.cp().signals().ExtractShapes( + in_file, '$BIN/$PROCESS$MASS', '$BIN/$PROCESS$MASS_$SYSTEMATIC') + # in_file, '$BIN/$PROCESS', '$BIN/$PROCESS__$SYSTEMATIC') + + if addBBB: + bbb = ch.BinByBinFactory().SetAddThreshold(0.).SetFixNorm(False) + bbb.MergeBinErrors(cb.cp().backgrounds()) + bbb.AddBinByBin(cb.cp().backgrounds(), cb) + + print '>> Setting standardised bin names...' + ch.SetStandardBinNames(cb) + cb.PrintAll() + + writer = ch.CardWriter( + '$TAG/$MASS/$ANALYSIS_$CHANNEL_$BINID.txt', + '$TAG/$ANALYSIS_$CHANNEL.input.root') + # writer.SetVerbosity(100) + writer.WriteCards('output/{mode}'.format(mode=mode), cb) + print 'Try writing cards...' + # import ROOT + # f_out = ROOT.TFile('andrey_out.root', 'RECREATE') + # cb.WriteDatacard("andrey_out.txt", 'andrey_out.root') + # writer.WriteCards('output/andrey_cards/', cb) print '>> Done!' From 79731e6da3849fd340096356ec48a65737ae904f Mon Sep 17 00:00:00 2001 From: mauro verzetti Date: Thu, 30 Mar 2017 09:27:13 -0500 Subject: [PATCH 4/7] Make script more flexible and update recipes --- CombineTools/scripts/setup_andrey.py | 153 ------------ Httbar/scripts/setup_andrey.py | 337 +++++++++++++++------------ 2 files changed, 186 insertions(+), 304 deletions(-) delete mode 100755 CombineTools/scripts/setup_andrey.py diff --git a/CombineTools/scripts/setup_andrey.py b/CombineTools/scripts/setup_andrey.py deleted file mode 100755 index c7d7fdb7ac1..00000000000 --- a/CombineTools/scripts/setup_andrey.py +++ /dev/null @@ -1,153 +0,0 @@ -#!/usr/bin/env python - -import CombineHarvester.CombineTools.ch as ch -import os -from argparse import ArgumentParser - -parser = ArgumentParser() -parser.add_argument('inputfile') -parser.add_argument('limit', choices=['electrons', 'muons', 'cmb'], help='choose leptonic decay type') -args = parser.parse_args() - -cb = ch.CombineHarvester() - -auxiliaries = os.environ['CMSSW_BASE'] + '/src/CombineHarvester/CombineTools/scripts/' -aux_shapes = auxiliaries - -addBBB = False - -masses = ['400', '500', '600', '750'] -widths = ['5', '10', '25', '50'] - -# for mode in ['scalar']:#, 'pseudoscalar']: -for parity in ['pseudoscalar']:#, 'pseudoscalar']: - for width in widths: - mode='{mode}_{width}pc'.format(mode=parity, width=width) - patterns = [i % width for i in ['ggA_pos-int-%spc-M', 'ggA_neg-int-%spc-M', 'ggA_pos-sgn-%spc-M']] - procs = { - # 'sig' : ['S0_{mode}_M'.format(mode=mode), 'S0_neg_{mode}_M'.format(mode=mode)], - 'sig' : patterns, - # 'sim' : ['WZ', 'ZZ', 'ttW', 'ttZ', 'ttH'], - 'bkg' : ['TT', 'VV', 'TTV', 'tChannel', 'tWChannel', 'WJets', 'ZJets'] - } - if args.limit == 'electrons' or args.limit == 'cmb': - procs['bkg'].append('QCDejets') - if args.limit == 'muons' or args.limit == 'cmb': - procs['bkg'].append('QCDmujets') - - cats = [] - if args.limit == 'muons' or args.limit == 'cmb': - cats.append((len(cats), 'mujets')) - if args.limit == 'electrons' or args.limit == 'cmb': - cats.append((len(cats), 'ejets')) - - cb.AddObservations(['*'], ['mauro'], ["13TeV"], ['htt'], cats) - cb.AddProcesses(['*'], ['mauro'], ["13TeV"], ['htt'], procs['bkg'], cats, False) - cb.AddProcesses(masses, ['mauro'], ["13TeV"], ['htt'], procs['sig'], cats, True) - # cb.AddProcesses(['*'], ['mauro'], ["13TeV"], ['htt'], procs['sig'], cats, True) - - print '>> Adding systematic uncertainties...' - # - # Normalization uncertainty - # - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'lumi_13TeV', 'lnN', ch.SystMap()(1.026)) - - cb.cp().process(['VV']).AddSyst( - cb, 'CMS_httbar_vvXsec_13TeV', 'lnN', ch.SystMap()(1.3)) - - cb.cp().process(['TT']).AddSyst( - cb, 'ttbar_rate', 'lnN', ch.SystMap()(1.05)) - - cb.cp().process(['tChannel']).AddSyst( - cb, 'st_rate', 'lnN', ch.SystMap()(1.05)) - - cb.cp().process(['tWChannel']).AddSyst( - cb, 'st_rate', 'lnN', ch.SystMap()(1.05)) - - cb.cp().process(['WJets']).AddSyst( - cb, 'w_rate', 'lnN', ch.SystMap()(1.03)) - cb.cp().process(['ZJets']).AddSyst( - cb, 'w_rate', 'lnN', ch.SystMap()(1.03)) - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'trigger', 'lnN', ch.SystMap()(1.0115)) - - # - # Shape uncetainties - # - #on everything - for name in [ - 'CMS_scale_j_13TeV', - 'CMS_eff_b_13TeV', - 'CMS_fake_b_13TeV', - 'CMS_pileup', - 'CMS_METunclustered_13TeV' - ]: - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, name, 'shape', ch.SystMap()(1.)) - - #on TT only - for name in [ - 'QCDscaleMERenorm_TT', - 'QCDscaleMEFactor_TT', - 'QCDscaleMERenormFactor_TT', - 'QCDscalePS_TT', - 'TMass', - 'pdf' - ]: - cb.cp().process(['TT']).AddSyst( - cb, name, 'shape', ch.SystMap()(1.)) - - if args.limit == 'muons' or args.limit == 'cmb': - cb.cp().process(['QCDmujets']).AddSyst( - cb, 'QCD_mu_norm', 'lnN', ch.SystMap()(2.)) - if args.limit == 'electrons' or args.limit == 'cmb': - cb.cp().process(['QCDejets']).AddSyst( - cb, 'QCD_el_norm', 'lnN', ch.SystMap()(2.)) - - - print '>> Extracting histograms from input root files...' - in_file = args.inputfile - cb.cp().backgrounds().ExtractShapes( - in_file, '$BIN/$PROCESS', '$BIN/$PROCESS_$SYSTEMATIC') - cb.cp().signals().ExtractShapes( - in_file, '$BIN/$PROCESS$MASS', '$BIN/$PROCESS$MASS_$SYSTEMATIC') - # in_file, '$BIN/$PROCESS', '$BIN/$PROCESS__$SYSTEMATIC') - - if addBBB: - bbb = ch.BinByBinFactory().SetAddThreshold(0.).SetFixNorm(False) - bbb.MergeBinErrors(cb.cp().backgrounds()) - bbb.AddBinByBin(cb.cp().backgrounds(), cb) - - print '>> Setting standardised bin names...' - ch.SetStandardBinNames(cb) - cb.PrintAll() - - writer = ch.CardWriter( - '$TAG/$MASS/$ANALYSIS_$CHANNEL_$BINID.txt', - '$TAG/$ANALYSIS_$CHANNEL.input.root') - # writer.SetVerbosity(100) - writer.WriteCards('output/{mode}'.format(mode=mode), cb) - print 'Try writing cards...' - # import ROOT - # f_out = ROOT.TFile('andrey_out.root', 'RECREATE') - # cb.WriteDatacard("andrey_out.txt", 'andrey_out.root') - # writer.WriteCards('output/andrey_cards/', cb) - -print '>> Done!' - -# Post instructions: -''' -combineTool.py -M T2W -i {scalar,pseudoscalar}/* -o workspace.root -P CombineHarvester.CombineTools.InterferenceModel:interferenceModel -combineTool.py -M Asymptotic -d */*/workspace.root --there -n .limit --parallel 4 -combineTool.py -M CollectLimits */*/*.limit.* --use-dirs -o limits.json -plotLimits.py --y-title="Coupling modifier" --x-title="M_{A} (GeV)" limits_default.json - -combineTool.py -M Impacts -d workspace.root -m 600 --doInitialFit --robustFit 1 -combineTool.py -M Impacts -d workspace.root -m 600 --robustFit 1 --doFits -# combineTool.py -M ImpactsFromScans -d workspace.root -m 600 --robustFit 1 --doFits --robustFit on -combineTool.py -M Impacts -d workspace.root -m 600 -o impacts.json -plotImpacts.py -i impacts.json -o impacts -''' - - diff --git a/Httbar/scripts/setup_andrey.py b/Httbar/scripts/setup_andrey.py index 490849322e5..645a1f76da4 100644 --- a/Httbar/scripts/setup_andrey.py +++ b/Httbar/scripts/setup_andrey.py @@ -1,4 +1,17 @@ #!/usr/bin/env python +from argparse import ArgumentParser +from pdb import set_trace + +parser = ArgumentParser() +parser.add_argument('jobid') +parser.add_argument('--limit' , choices=['electrons', 'muons', 'cmb'], help='choose leptonic decay type', default='cmb') +parser.add_argument('--masses', default='400,500,600,750', help='coma separated list of masses') +parser.add_argument('--widths', default='5,10,25,50', help='coma separated list of widths') +parser.add_argument('--parity', default='A,H', help='coma separated list of parity (A,H only)') +parser.add_argument('--noBBB', action='store_true') +parser.add_argument('--noMorph', action='store_true') +args = parser.parse_args() + import os from ROOT import RooWorkspace, TFile, RooRealVar @@ -6,160 +19,182 @@ import CombineHarvester.CombineTools.ch as ch from CombineHarvester.CombinePdfs.morphing import BuildRooMorphing -cb = ch.CombineHarvester() - auxiliaries = os.environ['CMSSW_BASE'] + '/src/CombineHarvester/Httbar/data/' aux_shapes = auxiliaries -addBBB = True - -masses = ['400', '500', '600', '750'] - -width = '5' # in percent - -doMorph = True - -for mode in ['A']: - - patterns = ['gg{mode}_pos-sgn-{width}pc-M', 'gg{mode}_pos-int-{width}pc-M', 'gg{mode}_neg-int-{width}pc-M'] - - procs = { - 'sig': [pattern.format(mode=mode, width=width) for pattern in patterns], - 'bkg': ['WJets', 'tWChannel', 'tChannel', 'VV', 'ZJets', 'TT','TTV'], - # 'bkg_mu':['QCDmujets'], # JAN: Ignore QCD for now because of extreme bbb uncertainties - 'bkg_mu':[], - 'bkg_e':[] - } - - cats = [(0, 'mujets'), (1, 'ejets')] - cat_to_id = {a:b for b, a in cats} - - cb.AddObservations(['*'], ['httbar'], ["8TeV"], [''], cats) - cb.AddProcesses(['*'], ['httbar'], ["8TeV"], [''], procs['bkg'], cats, False) - cb.AddProcesses(masses, ['httbar'], ["8TeV"], [''], procs['sig'], cats, True) - - print '>> Adding systematic uncertainties...' - - ### RATE UNCERTAINTIES - - # THEORY - cb.cp().process(['VV']).AddSyst( - cb, 'CMS_httbar_vvNorm_13TeV', 'lnN', ch.SystMap()(1.1)) - - cb.cp().process(['TT']).AddSyst( - cb, 'QCDscale_ttbar', 'lnN', ch.SystMap()(1.059)) - - cb.cp().process(['tWChannel']).AddSyst( - cb, 'CMS_httbar_tWChannelNorm_13TeV', 'lnN', ch.SystMap()(1.15)) - - cb.cp().process(['tChannel']).AddSyst( - cb, 'CMS_httbar_tChannelNorm_13TeV', 'lnN', ch.SystMap()(1.20)) - - cb.cp().process(['WJets']).AddSyst( - cb, 'CMS_httbar_WNorm_13TeV', 'lnN', ch.SystMap()(1.5)) - - cb.cp().process(['ZJets']).AddSyst( - cb, 'CMS_httbar_ZNorm_13TeV', 'lnN', ch.SystMap()(1.5)) - - cb.cp().process(['TTV']).AddSyst( - cb, 'CMS_httbar_TTVNorm_13TeV', 'lnN', ch.SystMap()(1.3)) - - # EXPERIMENT - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'lumi', 'lnN', ch.SystMap()(1.058)) - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'CMS_eff_trigger_m', 'lnN', ch.SystMap('bin_id')([cat_to_id['mujets']], 1.02)) - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'CMS_eff_trigger_e', 'lnN', ch.SystMap('bin_id')([cat_to_id['ejets']], 1.02)) - - - # GENERIC SHAPE UNCERTAINTIES - shape_uncertainties = ['CMS_pileup', 'CMS_eff_b_13TeV', 'CMS_fake_b_13TeV', 'CMS_scale_j_13TeV', 'CMS_res_j_13TeV', 'CMS_METunclustered_13TeV'] - - for shape_uncertainty in shape_uncertainties: - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst(cb, shape_uncertainty, 'shape', ch.SystMap()(1.)) - - shape_uncertainties_mu = ['CMS_eff_m'] - for shape_uncertainty in shape_uncertainties_mu: - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst(cb, shape_uncertainty, 'shape', ch.SystMap('bin_id')([cat_to_id['mujets']], 1.)) - - # SPECIFIC SHAPE UNCERTAINTIES - - shape_uncertainties_tt = ['QCDscaleMERenorm_TT', 'QCDscaleMEFactor_TT', 'QCDscaleMERenormFactor_TT', - 'Hdamp_TT', 'QCDscalePS_TT', 'TMass', 'pdf'] - - for shape_uncertainty in shape_uncertainties_tt: - cb.cp().process(['TT']).AddSyst( - cb, shape_uncertainty, 'shape', ch.SystMap()(1.)) - - print '>> Extracting histograms from input root files...' - in_file = aux_shapes + 'templates1D_161110.root' - cb.cp().backgrounds().ExtractShapes( - in_file, '$BIN/$PROCESS', '$BIN/$PROCESS_$SYSTEMATIC') - cb.cp().signals().ExtractShapes( - in_file, '$BIN/$PROCESS$MASS', '$BIN/$PROCESS$MASS_$SYSTEMATIC') - # in_file, '$BIN/$PROCESS', '$BIN/$PROCESS__$SYSTEMATIC') - - # print '>> Generating bbb uncertainties...' - # bbb = ch.BinByBinFactory() - # bbb.SetAddThreshold(0.1).SetFixNorm(True) - # bbb.AddBinByBin(cb.cp().process(['reducible']), cb) - - # for mass in masses: - # norm_initial = norms[mode][int(mass)] - - # cb.cp().process(['S0_neg_{mode}_M{mass}'.format(mode=mode, mass=mass), 'S0_{mode}_M{mass}'.format(mode=mode, mass=mass)]).ForEachProc(lambda p: p.set_rate(p.rate()/norm_initial)) - - if addBBB: - bbb = ch.BinByBinFactory().SetAddThreshold(0.).SetFixNorm(False).SetMergeThreshold(0.5) - bbb.MergeAndAdd(cb.cp().backgrounds(), cb) - - if doMorph: - mA = RooRealVar('MH', 'MH', 400., 750.); - mA.setConstant(True) - - mass_debug = False - if mass_debug: - f_debug = TFile('morph_debug.root', 'RECREATE') - print 'Try to morph between masses' - ws = RooWorkspace('httbar', 'httbar') - bins = cb.bin_set() - for bin in bins: - for proc in procs['sig']: - BuildRooMorphing(ws, cb, bin, proc, mA, "norm", True, True, False, f_debug if mass_debug else None) - - if mass_debug: - f_debug.Close() - cb.AddWorkspace(ws, False) - cb.cp().process(procs['sig']).ExtractPdfs(cb, "httbar", "$BIN_$PROCESS_morph", "") - - #BuildRooMorphing(ws, cb, bin, process, mass_var, norm_postfix='norm', allow_morph=True, verbose=False, force_template_limit=False, file=None) - # void BuildRooMorphing(RooWorkspace& ws, CombineHarvester& cb, - # std::string const& bin, std::string const& process, - # RooAbsReal& mass_var, std::string norm_postfix, - # bool allow_morph, bool verbose, bool force_template_limit, TFile * file) - - print '>> Setting standardised bin names...' - ch.SetStandardBinNames(cb) - # cb.PrintAll() - - if not doMorph: - writer = ch.CardWriter('$TAG/$MASS/$ANALYSIS_$CHANNEL_$BINID.txt', - # writer = ch.CardWriter('$TAG/$ANALYSIS_$CHANNEL_$BINID_$ERA.txt', - '$TAG/$ANALYSIS_$CHANNEL.input.root') - else: - writer = ch.CardWriter('$TAG/MORPH/$ANALYSIS_$CHANNEL_$BINID.txt', - '$TAG/$ANALYSIS_$CHANNEL.input.root') - writer.SetWildcardMasses([]) - writer.SetVerbosity(100) - writer.WriteCards('output/{mode}_{width}'.format(mode=mode, width=width), cb) - print 'Writing cards...' - # import ROOT - # f_out = ROOT.TFile('andrey_out.root', 'RECREATE') - # cb.WriteDatacard("andrey_out.txt", 'andrey_out.root') - # writer.WriteCards('output/andrey_cards/', cb) +addBBB = not args.noBBB + +masses = args.masses.split(',') +widths = args.widths.split(',')#['5', '10', '25', '50'] +modes = args.parity.split(',') + +doMorph = not args.noMorph + +for mode in modes: + for width in widths: + cb = ch.CombineHarvester() + patterns = ['gg{mode}_pos-sgn-{width}pc-M', 'gg{mode}_pos-int-{width}pc-M', 'gg{mode}_neg-int-{width}pc-M'] + + procs = { + 'sig': [pattern.format(mode=mode, width=width) for pattern in patterns], + 'bkg': ['WJets', 'tWChannel', 'tChannel', 'VV', 'ZJets', 'TT','TTV'], + 'bkg_mu':['QCDmujets'], + 'bkg_e' :['QCDejets' ] + } + + cats = [] + if args.limit == 'muons' or args.limit == 'cmb': + cats.append((len(cats), 'mujets')) + if args.limit == 'electrons' or args.limit == 'cmb': + cats.append((len(cats), 'ejets')) + cat_to_id = {a:b for b, a in cats} + + cb.AddObservations(['*'], ['httbar'], ["8TeV"], [''], cats) + cb.AddProcesses(['*'], ['httbar'], ["8TeV"], [''], procs['bkg'], cats, False) + cb.AddProcesses(masses, ['httbar'], ["8TeV"], [''], procs['sig'], cats, True) + if args.limit == 'muons' or args.limit == 'cmb': + cb.AddProcesses( + ['*'], ['httbar'], ["8TeV"], + [''], procs['bkg_mu'], + [(cat_to_id['mujets'], 'mujets')], False) + if args.limit == 'electrons' or args.limit == 'cmb': + cb.AddProcesses( + ['*'], ['httbar'], ["8TeV"], + [''], procs['bkg_e'], + [(cat_to_id['ejets'], 'ejets')], False) + + print '>> Adding systematic uncertainties...' + + ### RATE UNCERTAINTIES + + # THEORY + cb.cp().process(['VV']).AddSyst( + cb, 'CMS_httbar_vvNorm_13TeV', 'lnN', ch.SystMap()(1.5)) + + cb.cp().process(['TT']).AddSyst( + cb, 'CMS_httbar_ttbarNorm_13TeV', 'lnN', ch.SystMap()(1.059)) + + cb.cp().process(['tWChannel']).AddSyst( + cb, 'CMS_httbar_tWChannelNorm_13TeV', 'lnN', ch.SystMap()(1.15)) + + cb.cp().process(['tChannel']).AddSyst( + cb, 'CMS_httbar_tChannelNorm_13TeV', 'lnN', ch.SystMap()(1.20)) + + cb.cp().process(['WJets']).AddSyst( + cb, 'CMS_httbar_WNorm_13TeV', 'lnN', ch.SystMap()(1.5)) + + cb.cp().process(['ZJets']).AddSyst( + cb, 'CMS_httbar_ZNorm_13TeV', 'lnN', ch.SystMap()(1.5)) + + cb.cp().process(['TTV']).AddSyst( + cb, 'CMS_httbar_TTVNorm_13TeV', 'lnN', ch.SystMap()(1.2)) + + # Experiment + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( + cb, 'lumi', 'lnN', ch.SystMap()(1.058)) + + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( + cb, 'CMS_eff_trigger_m', 'lnN', ch.SystMap('bin_id')([cat_to_id['mujets']], 1.02)) + + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( + cb, 'CMS_eff_trigger_e', 'lnN', ch.SystMap('bin_id')([cat_to_id['ejets']], 1.02)) + + if args.limit == 'muons' or args.limit == 'cmb': + cb.cp().process(['QCDmujets']).AddSyst( + cb, 'CMS_httbar_mujets_QCDNorm', 'lnN', ch.SystMap('bin_id')([cat_to_id['mujets']], 2.0)) + if args.limit == 'electrons' or args.limit == 'cmb': + cb.cp().process(['QCDejets']).AddSyst( + cb, 'CMS_httbar_ejets_QCDNorm', 'lnN', ch.SystMap('bin_id')([cat_to_id['ejets']], 2.0)) + + # GENERIC SHAPE UNCERTAINTIES + shape_uncertainties = [ + 'CMS_pileup', 'CMS_eff_b_13TeV', 'CMS_fake_b_13TeV', + 'CMS_scale_j_13TeV', 'CMS_res_j_13TeV', 'CMS_METunclustered_13TeV'] + + for shape_uncertainty in shape_uncertainties: + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst(cb, shape_uncertainty, 'shape', ch.SystMap()(1.)) + + shape_uncertainties_mu = ['CMS_eff_m'] + for shape_uncertainty in shape_uncertainties_mu: + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst(cb, shape_uncertainty, 'shape', ch.SystMap('bin_id')([cat_to_id['mujets']], 1.)) + + # SPECIFIC SHAPE UNCERTAINTIES + + shape_uncertainties_tt = [ + 'QCDscaleMERenorm_TT', 'QCDscaleMEFactor_TT', + 'QCDscaleISR_TT', 'QCDscaleFSR_TT', + 'Hdamp_TT', 'TMass', 'pdf'] + + for shape_uncertainty in shape_uncertainties_tt: + cb.cp().process(['TT']).AddSyst( + cb, shape_uncertainty, 'shape', ch.SystMap()(1.)) + + print '>> Extracting histograms from input root files...' + in_file = aux_shapes + 'templates_%s.root' % args.jobid #1D_161110 + cb.cp().backgrounds().ExtractShapes( + in_file, '$BIN/$PROCESS', '$BIN/$PROCESS_$SYSTEMATIC') + cb.cp().signals().ExtractShapes( + in_file, '$BIN/$PROCESS$MASS', '$BIN/$PROCESS$MASS_$SYSTEMATIC') + # in_file, '$BIN/$PROCESS', '$BIN/$PROCESS__$SYSTEMATIC') + + # print '>> Generating bbb uncertainties...' + # bbb = ch.BinByBinFactory() + # bbb.SetAddThreshold(0.1).SetFixNorm(True) + # bbb.AddBinByBin(cb.cp().process(['reducible']), cb) + + # for mass in masses: + # norm_initial = norms[mode][int(mass)] + + # cb.cp().process(['S0_neg_{mode}_M{mass}'.format(mode=mode, mass=mass), 'S0_{mode}_M{mass}'.format(mode=mode, mass=mass)]).ForEachProc(lambda p: p.set_rate(p.rate()/norm_initial)) + + if addBBB: + bbb = ch.BinByBinFactory().SetAddThreshold(0.).SetFixNorm(False).SetMergeThreshold(0.5) + bbb.MergeAndAdd(cb.cp().process(['TT']), cb) + + if doMorph: + mA = RooRealVar('MH', 'MH', 400., 750.); + mA.setConstant(True) + + mass_debug = False + if mass_debug: + f_debug = TFile('morph_debug.root', 'RECREATE') + print 'Try to morph between masses' + ws = RooWorkspace('httbar', 'httbar') + bins = cb.bin_set() + for bin in bins: + for proc in procs['sig']: + BuildRooMorphing(ws, cb, bin, proc, mA, "norm", True, True, False, f_debug if mass_debug else None) + + if mass_debug: + f_debug.Close() + cb.AddWorkspace(ws, False) + cb.cp().process(procs['sig']).ExtractPdfs(cb, "httbar", "$BIN_$PROCESS_morph", "") + + #BuildRooMorphing(ws, cb, bin, process, mass_var, norm_postfix='norm', allow_morph=True, verbose=False, force_template_limit=False, file=None) + # void BuildRooMorphing(RooWorkspace& ws, CombineHarvester& cb, + # std::string const& bin, std::string const& process, + # RooAbsReal& mass_var, std::string norm_postfix, + # bool allow_morph, bool verbose, bool force_template_limit, TFile * file) + + print '>> Setting standardised bin names...' + ch.SetStandardBinNames(cb) + #cb.PrintAll() + + if not doMorph: + writer = ch.CardWriter('$TAG/$MASS/$ANALYSIS_$CHANNEL_$BINID.txt', + # writer = ch.CardWriter('$TAG/$ANALYSIS_$CHANNEL_$BINID_$ERA.txt', + '$TAG/$ANALYSIS_$CHANNEL.input.root') + else: + writer = ch.CardWriter('$TAG/MORPH/$ANALYSIS_$CHANNEL_$BINID.txt', + '$TAG/$ANALYSIS_$CHANNEL.input.root') + #writer.SetWildcardMasses([]) + writer.SetVerbosity(100) + writer.WriteCards('output{jobid}/{mode}_{width}'.format(jobid=args.jobid, mode=mode, width=width), cb) + print 'Writing cards...' + # import ROOT + # f_out = ROOT.TFile('andrey_out.root', 'RECREATE') + # cb.WriteDatacard("andrey_out.txt", 'andrey_out.root') + # writer.WriteCards('output/andrey_cards/', cb) print '>> Done!' From ed4953986231ef006a0d72a8991d3b74fa72ef27 Mon Sep 17 00:00:00 2001 From: mauro verzetti Date: Fri, 31 Mar 2017 04:58:10 -0500 Subject: [PATCH 5/7] wip --- Httbar/scripts/setup_andrey.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) mode change 100644 => 100755 Httbar/scripts/setup_andrey.py diff --git a/Httbar/scripts/setup_andrey.py b/Httbar/scripts/setup_andrey.py old mode 100644 new mode 100755 index 645a1f76da4..1a788765f41 --- a/Httbar/scripts/setup_andrey.py +++ b/Httbar/scripts/setup_andrey.py @@ -131,6 +131,7 @@ print '>> Extracting histograms from input root files...' in_file = aux_shapes + 'templates_%s.root' % args.jobid #1D_161110 + cb.cp().backgrounds().ExtractShapes( in_file, '$BIN/$PROCESS', '$BIN/$PROCESS_$SYSTEMATIC') cb.cp().signals().ExtractShapes( @@ -182,19 +183,24 @@ if not doMorph: writer = ch.CardWriter('$TAG/$MASS/$ANALYSIS_$CHANNEL_$BINID.txt', - # writer = ch.CardWriter('$TAG/$ANALYSIS_$CHANNEL_$BINID_$ERA.txt', '$TAG/$ANALYSIS_$CHANNEL.input.root') else: writer = ch.CardWriter('$TAG/MORPH/$ANALYSIS_$CHANNEL_$BINID.txt', '$TAG/$ANALYSIS_$CHANNEL.input.root') #writer.SetWildcardMasses([]) + set_trace() writer.SetVerbosity(100) - writer.WriteCards('output{jobid}/{mode}_{width}'.format(jobid=args.jobid, mode=mode, width=width), cb) print 'Writing cards...' - # import ROOT - # f_out = ROOT.TFile('andrey_out.root', 'RECREATE') - # cb.WriteDatacard("andrey_out.txt", 'andrey_out.root') - # writer.WriteCards('output/andrey_cards/', cb) + writer.WriteCards('output{jobid}/{mode}_{width}'.format(jobid=args.jobid, mode=mode, width=width), cb) + #writer.WriteCards('output/{mode}_{width}'.format(mode=mode, width=width), cb) + ## import ROOT + ## dirname = 'output{jobid}/{mode}_{width}'.format(jobid=args.jobid, mode=mode, width=width) + ## if os.path.exists(dirname): + ## os.makedirs(dirname) + ## f_out = ROOT.TFile('%s/httbar.input.root', 'RECREATE') + ## split = cb.cp().mass(["*"]) + ## split.WriteDatacard(dirname + "/httbar.txt", f_out) + ## f_out.Close() print '>> Done!' From c0ba260c52d045a557c359202ee7ef5d150a86d5 Mon Sep 17 00:00:00 2001 From: mauro verzetti Date: Fri, 31 Mar 2017 08:04:50 -0500 Subject: [PATCH 6/7] make system mode flexible and automated plot production --- Httbar/BuildFile.xml | 0 Httbar/results/Makefile | 29 +++++++++++++++++++++++++ Httbar/scripts/Httbar_plots.sh | 12 +++++++++++ Httbar/scripts/Httbar_workspace.sh | 8 +++++++ Httbar/scripts/setup_andrey.py | 34 ++++-------------------------- 5 files changed, 53 insertions(+), 30 deletions(-) create mode 100644 Httbar/BuildFile.xml create mode 100644 Httbar/results/Makefile create mode 100755 Httbar/scripts/Httbar_plots.sh create mode 100755 Httbar/scripts/Httbar_workspace.sh diff --git a/Httbar/BuildFile.xml b/Httbar/BuildFile.xml new file mode 100644 index 00000000000..e69de29bb2d diff --git a/Httbar/results/Makefile b/Httbar/results/Makefile new file mode 100644 index 00000000000..a9325386c4e --- /dev/null +++ b/Httbar/results/Makefile @@ -0,0 +1,29 @@ +jobid=1D_2017Mar23 +runmode='--run blind' +limitdir=output$(jobid) + +$(limitdir)/.cards: ../data/templates_$(jobid).root ../scripts/setup_andrey.py + python ../scripts/setup_andrey.py 1D_2017Mar23 --parity=A && touch $@ + +cards: $(limitdir)/.cards + +$(limitdir)/.workspace: $(limitdir)/.cards ../scripts/Httbar_workspace.sh + cd .. && scram b && cd - && cd $(limitdir) && Httbar_workspace.sh && cd - && touch $@ + +workspace: $(limitdir)/.workspace + +$(limitdir)/.limit: $(limitdir)/.workspace + cd $(limitdir) && combineTool.py -M Asymptotic -m 400:750:50 --freezeNuisances MH -d */*/workspace.root --there -n .limit --parallel 8 $(runmode) && cd - && touch $@ + +limits: $(limitdir)/.limit + +$(limitdir)/.jsons: $(limitdir)/.limit + cd $(limitdir) && combineTool.py -M CollectLimits */*/*.limit.* --use-dirs && cd - && touch $@ + +jsons: $(limitdir)/.jsons + +$(limitdir)/.plot: $(limitdir)/.jsons ../scripts/Httbar_plots.sh + cd .. && scram b && cd - && cd $(limitdir) && mkdir -p plots && Httbar_plots.sh && cd - && touch $@ + +plots: $(limitdir)/.plot + diff --git a/Httbar/scripts/Httbar_plots.sh b/Httbar/scripts/Httbar_plots.sh new file mode 100755 index 00000000000..6977165147c --- /dev/null +++ b/Httbar/scripts/Httbar_plots.sh @@ -0,0 +1,12 @@ +#! /bin/bash + +set -o nounset +set -o errexit + +for json in $(ls -d limits_*.json); do + #remove "limits" + mode_width=`echo "${json#*_}" | sed 's|.json||g'` + ah="${mode_width%_*}" + width="${mode_width#*_}" + plotLimits.py --y-title='Coupling modifier' --x-title='M_{'$ah'} (GeV)' $json -o plots/limit_$ah'_'$width'_pc' --show=exp --grid --mapping='lambda x: sqrt(x)' +done diff --git a/Httbar/scripts/Httbar_workspace.sh b/Httbar/scripts/Httbar_workspace.sh new file mode 100755 index 00000000000..5d2d4e8941c --- /dev/null +++ b/Httbar/scripts/Httbar_workspace.sh @@ -0,0 +1,8 @@ +#! /bin/bash + +set -o nounset +set -o errexit + +for wdir in $(ls -d [AH]_[0-9]*); do + combineTool.py -M T2W -i $wdir/* -o workspace.root -P CombineHarvester.CombineTools.InterferenceModel:interferenceModel +done diff --git a/Httbar/scripts/setup_andrey.py b/Httbar/scripts/setup_andrey.py index 1a788765f41..5de1e93603a 100755 --- a/Httbar/scripts/setup_andrey.py +++ b/Httbar/scripts/setup_andrey.py @@ -47,6 +47,7 @@ cats.append((len(cats), 'mujets')) if args.limit == 'electrons' or args.limit == 'cmb': cats.append((len(cats), 'ejets')) + cats = [(0, 'mujets'), (1, 'ejets')] cat_to_id = {a:b for b, a in cats} cb.AddObservations(['*'], ['httbar'], ["8TeV"], [''], cats) @@ -136,17 +137,6 @@ in_file, '$BIN/$PROCESS', '$BIN/$PROCESS_$SYSTEMATIC') cb.cp().signals().ExtractShapes( in_file, '$BIN/$PROCESS$MASS', '$BIN/$PROCESS$MASS_$SYSTEMATIC') - # in_file, '$BIN/$PROCESS', '$BIN/$PROCESS__$SYSTEMATIC') - - # print '>> Generating bbb uncertainties...' - # bbb = ch.BinByBinFactory() - # bbb.SetAddThreshold(0.1).SetFixNorm(True) - # bbb.AddBinByBin(cb.cp().process(['reducible']), cb) - - # for mass in masses: - # norm_initial = norms[mode][int(mass)] - - # cb.cp().process(['S0_neg_{mode}_M{mass}'.format(mode=mode, mass=mass), 'S0_{mode}_M{mass}'.format(mode=mode, mass=mass)]).ForEachProc(lambda p: p.set_rate(p.rate()/norm_initial)) if addBBB: bbb = ch.BinByBinFactory().SetAddThreshold(0.).SetFixNorm(False).SetMergeThreshold(0.5) @@ -170,13 +160,7 @@ f_debug.Close() cb.AddWorkspace(ws, False) cb.cp().process(procs['sig']).ExtractPdfs(cb, "httbar", "$BIN_$PROCESS_morph", "") - - #BuildRooMorphing(ws, cb, bin, process, mass_var, norm_postfix='norm', allow_morph=True, verbose=False, force_template_limit=False, file=None) - # void BuildRooMorphing(RooWorkspace& ws, CombineHarvester& cb, - # std::string const& bin, std::string const& process, - # RooAbsReal& mass_var, std::string norm_postfix, - # bool allow_morph, bool verbose, bool force_template_limit, TFile * file) - + print '>> Setting standardised bin names...' ch.SetStandardBinNames(cb) #cb.PrintAll() @@ -185,22 +169,12 @@ writer = ch.CardWriter('$TAG/$MASS/$ANALYSIS_$CHANNEL_$BINID.txt', '$TAG/$ANALYSIS_$CHANNEL.input.root') else: - writer = ch.CardWriter('$TAG/MORPH/$ANALYSIS_$CHANNEL_$BINID.txt', + writer = ch.CardWriter('$TAG/000/$ANALYSIS_$CHANNEL_$BINID.txt', '$TAG/$ANALYSIS_$CHANNEL.input.root') - #writer.SetWildcardMasses([]) - set_trace() + writer.SetWildcardMasses([]) writer.SetVerbosity(100) print 'Writing cards...' writer.WriteCards('output{jobid}/{mode}_{width}'.format(jobid=args.jobid, mode=mode, width=width), cb) - #writer.WriteCards('output/{mode}_{width}'.format(mode=mode, width=width), cb) - ## import ROOT - ## dirname = 'output{jobid}/{mode}_{width}'.format(jobid=args.jobid, mode=mode, width=width) - ## if os.path.exists(dirname): - ## os.makedirs(dirname) - ## f_out = ROOT.TFile('%s/httbar.input.root', 'RECREATE') - ## split = cb.cp().mass(["*"]) - ## split.WriteDatacard(dirname + "/httbar.txt", f_out) - ## f_out.Close() print '>> Done!' From 37738cb3d1d54d887bcdecfccb53df53ee7311a1 Mon Sep 17 00:00:00 2001 From: mauro verzetti Date: Wed, 10 May 2017 03:44:33 -0500 Subject: [PATCH 7/7] WIP: width interpolation and standard flow --- Httbar/results/Makefile | 15 +++- Httbar/scripts/compare_limits.py | 66 +++++++++++++++++ Httbar/scripts/morph_widths.py | 97 +++++++++++++++++++++++++ Httbar/scripts/setup_andrey.py | 39 +++++----- Httbar/scripts/systematics_breakdown.py | 83 +++++++++++++++++++++ 5 files changed, 281 insertions(+), 19 deletions(-) create mode 100755 Httbar/scripts/compare_limits.py create mode 100755 Httbar/scripts/morph_widths.py create mode 100755 Httbar/scripts/systematics_breakdown.py diff --git a/Httbar/results/Makefile b/Httbar/results/Makefile index a9325386c4e..f3998e3572c 100644 --- a/Httbar/results/Makefile +++ b/Httbar/results/Makefile @@ -2,7 +2,10 @@ jobid=1D_2017Mar23 runmode='--run blind' limitdir=output$(jobid) -$(limitdir)/.cards: ../data/templates_$(jobid).root ../scripts/setup_andrey.py +#../data/templates_$(jobid)_morphed.root: ../data/templates_$(jobid).root ../scripts/morph_widths.py +# ../scripts/morph_widths.py ../data/templates_$(jobid).root + +$(limitdir)/.cards: ../data/templates_$(jobid).root ../scripts/setup_andrey.py #../data/templates_$(jobid)_morphed.root ../scripts/setup_andrey.py python ../scripts/setup_andrey.py 1D_2017Mar23 --parity=A && touch $@ cards: $(limitdir)/.cards @@ -25,5 +28,15 @@ jsons: $(limitdir)/.jsons $(limitdir)/.plot: $(limitdir)/.jsons ../scripts/Httbar_plots.sh cd .. && scram b && cd - && cd $(limitdir) && mkdir -p plots && Httbar_plots.sh && cd - && touch $@ +$(limitdir)/.jsons: $(limitdir)/.limit + cd $(limitdir) && combineTool.py -M CollectLimits */*/*.limit.* --use-dirs && cd - && touch $@ + plots: $(limitdir)/.plot +$(limitdir)/.breakdown: $(limitdir)/.workspace + cd $(limitdir) && systematics_breakdown.py $(runmode) && cd - && touch $@ + +breakdown: $(limitdir)/.breakdown + + +all: $(limitdir)/.plot diff --git a/Httbar/scripts/compare_limits.py b/Httbar/scripts/compare_limits.py new file mode 100755 index 00000000000..690611332e9 --- /dev/null +++ b/Httbar/scripts/compare_limits.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python +from argparse import ArgumentParser + +parser = ArgumentParser() +parser.add_argument('jsons', nargs='+', help='json_file:label') +#FIXME add auto labels +parser.add_argument('-x', help='x label') +parser.add_argument('-y', help='y label') +parser.add_argument('-t', help='title') +parser.add_argument('-o', default='out.png', help='output') +parser.add_argument('--autolabels', action='store_true', help='output') +parser.add_argument( + '--use', choices=["exp+1", "exp+2", "exp-1", "exp-2", "exp0", 'obs'], + default='exp0', help='to use') +args = parser.parse_args() + +from ROOT import gROOT, TCanvas, TGraph, TLegend +import json +import uuid +import os + +def json2graph(jfile, topick): + jinfo = json.loads(open(jfile).read()) + points = [(float(i), j[topick]) for i, j in jinfo.iteritems()] + points.sort() + gr = TGraph(len(points)) + for i, xy in enumerate(points): + gr.SetPoint(i, xy[0], xy[1]) + gr.SetMarkerStyle(20) + return gr, min(i[1] for i in points), max(i[1] for i in points) + +uid = str(uuid.uuid1()) +canvas = TCanvas(uid, uid, 800, 800) +legend = TLegend(0.1,0.7,0.48,0.9); +keep = [] +ms = [] +Ms = [] +colors = [2,4,6,8,28,46,14,31,1] +first=True +for jfilelabel, color in zip(args.jsons, colors): + if not args.autolabels: + jfile, label = tuple(jfilelabel.split(':')) + else: + jfile = jfilelabel + label = os.path.basename(jfile).split('_')[0] + graph, m, M = json2graph(jfile, args.use) + ms.append(m) + Ms.append(M) + graph.SetMarkerColor(color) + graph.SetLineColor(color) + legend.AddEntry(graph, label, 'lp') + graph.SetTitle(args.t) + keep.append(graph) + +for graph in keep: + if first: + graph.Draw('ALP') + graph.GetYaxis().SetRangeUser(min(ms)*0.8, max(Ms)*1.2) + graph.GetXaxis().SetTitle(args.x) + graph.GetYaxis().SetTitle(args.y) + first = False + else: + graph.Draw('LP same') + +legend.Draw() +canvas.SaveAs(args.o) diff --git a/Httbar/scripts/morph_widths.py b/Httbar/scripts/morph_widths.py new file mode 100755 index 00000000000..a977811a239 --- /dev/null +++ b/Httbar/scripts/morph_widths.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python +from argparse import ArgumentParser + +parser = ArgumentParser() +parser.add_argument('inputfile') +parser.add_argument('--forchecks', action='store_true') +parser.add_argument('--hyperbolic', action='store_true', help='use hyperbolic interpolation') +args = parser.parse_args() + +import ROOT +import shutil +import re +from pdb import set_trace + +mapping = { + '2p5pc' : 2.5, + '5pc' : 5.0, + '10pc' : 10., + '25pc' : 25., + '50pc' : 50., + } + +checks = { + 5.0 : 'checks_5pc' , + 10. : 'checks_10pc', + 25. : 'checks_25pc', +} + +new_points = { + 7.0 : '7pc' , + 15. : '15pc', + 20. : '20pc', + 30. : '30pc', + 40. : '40pc', +} + +if args.forchecks: + new_points.update(checks) + +available = set(mapping.values()) +to_make = set(new_points.keys()) +intersection = available.intersection(to_make) +if len(intersection) and not args.forchecks: + raise ValueError('the points: %s are already computed, I will not do them again!' % ','.join(list(intersection))) + +interpolation = {} +for point in to_make: + above = sorted([i for i in available if point < i]) + below = sorted([i for i in available if point > i]) + if not above or not below: + raise ValueError('I cannot interpolate %s as it is outside the provided range!' % point) + interpolation[point] = (below[-1], above[0]) + +def get_info(name): + try: + process, sigint, width, mass_sys = tuple(name.split('-')) + except: + set_trace() + key = (process, sigint, mass_sys) + return key, mapping[width] + +def make_name(key, width): + process, sigint, mass_sys = key + return '-'.join([process, sigint, new_points[width], mass_sys]) + +outname = args.inputfile.replace('.root', '_morphed.root') +shutil.copyfile(args.inputfile, outname) +infile = ROOT.TFile(outname, 'UPDATE') + +for category in [i.GetName() for i in infile.GetListOfKeys()]: + tdir = infile.Get(category) + tdir.cd() + shapes = [i.GetName() for i in tdir.GetListOfKeys()] + shapes = [i for i in shapes if i.startswith('ggA_') or i.startswith('ggH_')] + #there are still a lot of spurious shapes, remove them + shapes = [i for i in shapes if '_pos-sgn' in i or '-int-' in i] + shapes_map = {} + for shape in shapes: + k, w = get_info(shape) + if k not in shapes_map: + shapes_map[k] = {} + shapes_map[k][w] = tdir.Get(shape) + + #compute new histograms + for key in shapes_map: + for width in new_points: + b, a = interpolation[width] + factor = (width-b)/(a-b) + if args.hyperbolic: + factor = (1./width-1./b)/(1./a-1./b) + below = shapes_map[key][b] + above = shapes_map[key][a] + new_name = make_name(key, width) + new_hist = below.Clone(new_name) + new_hist.Scale(1-factor) + new_hist.Add(above, factor) + new_hist.Write() diff --git a/Httbar/scripts/setup_andrey.py b/Httbar/scripts/setup_andrey.py index 5de1e93603a..cafbb08133b 100755 --- a/Httbar/scripts/setup_andrey.py +++ b/Httbar/scripts/setup_andrey.py @@ -6,7 +6,7 @@ parser.add_argument('jobid') parser.add_argument('--limit' , choices=['electrons', 'muons', 'cmb'], help='choose leptonic decay type', default='cmb') parser.add_argument('--masses', default='400,500,600,750', help='coma separated list of masses') -parser.add_argument('--widths', default='5,10,25,50', help='coma separated list of widths') +parser.add_argument('--widths', default='2p5,5,10,25,50', help='coma separated list of widths') parser.add_argument('--parity', default='A,H', help='coma separated list of parity (A,H only)') parser.add_argument('--noBBB', action='store_true') parser.add_argument('--noMorph', action='store_true') @@ -93,31 +93,34 @@ # Experiment cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( cb, 'lumi', 'lnN', ch.SystMap()(1.058)) + + # GENERIC SHAPE UNCERTAINTIES + shape_uncertainties = [ + 'CMS_pileup', 'CMS_eff_b_13TeV', 'CMS_fake_b_13TeV', + 'CMS_scale_j_13TeV', 'CMS_res_j_13TeV', 'CMS_METunclustered_13TeV'] - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'CMS_eff_trigger_m', 'lnN', ch.SystMap('bin_id')([cat_to_id['mujets']], 1.02)) - - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( - cb, 'CMS_eff_trigger_e', 'lnN', ch.SystMap('bin_id')([cat_to_id['ejets']], 1.02)) - + for shape_uncertainty in shape_uncertainties: + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst(cb, shape_uncertainty, 'shape', ch.SystMap()(1.)) + + #E/Mu uncertainties if args.limit == 'muons' or args.limit == 'cmb': cb.cp().process(['QCDmujets']).AddSyst( cb, 'CMS_httbar_mujets_QCDNorm', 'lnN', ch.SystMap('bin_id')([cat_to_id['mujets']], 2.0)) + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( + cb, 'CMS_eff_trigger_m', 'lnN', ch.SystMap('bin_id')([cat_to_id['mujets']], 1.02)) + shape_uncertainties_mu = ['CMS_eff_m'] + for shape_uncertainty in shape_uncertainties_mu: + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst(cb, shape_uncertainty, 'shape', ch.SystMap('bin_id')([cat_to_id['mujets']], 1.)) if args.limit == 'electrons' or args.limit == 'cmb': cb.cp().process(['QCDejets']).AddSyst( cb, 'CMS_httbar_ejets_QCDNorm', 'lnN', ch.SystMap('bin_id')([cat_to_id['ejets']], 2.0)) + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst( + cb, 'CMS_eff_trigger_e', 'lnN', ch.SystMap('bin_id')([cat_to_id['ejets']], 1.02)) + shape_uncertainties_mu = ['CMS_eff_e'] + for shape_uncertainty in shape_uncertainties_mu: + cb.cp().process(procs['sig'] + procs['bkg']).AddSyst(cb, shape_uncertainty, 'shape', ch.SystMap('bin_id')([cat_to_id['ejets']], 1.)) - # GENERIC SHAPE UNCERTAINTIES - shape_uncertainties = [ - 'CMS_pileup', 'CMS_eff_b_13TeV', 'CMS_fake_b_13TeV', - 'CMS_scale_j_13TeV', 'CMS_res_j_13TeV', 'CMS_METunclustered_13TeV'] - - for shape_uncertainty in shape_uncertainties: - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst(cb, shape_uncertainty, 'shape', ch.SystMap()(1.)) - shape_uncertainties_mu = ['CMS_eff_m'] - for shape_uncertainty in shape_uncertainties_mu: - cb.cp().process(procs['sig'] + procs['bkg']).AddSyst(cb, shape_uncertainty, 'shape', ch.SystMap('bin_id')([cat_to_id['mujets']], 1.)) # SPECIFIC SHAPE UNCERTAINTIES @@ -131,7 +134,7 @@ cb, shape_uncertainty, 'shape', ch.SystMap()(1.)) print '>> Extracting histograms from input root files...' - in_file = aux_shapes + 'templates_%s.root' % args.jobid #1D_161110 + in_file = aux_shapes + 'templates_%s_morphed.root' % args.jobid #1D_161110 cb.cp().backgrounds().ExtractShapes( in_file, '$BIN/$PROCESS', '$BIN/$PROCESS_$SYSTEMATIC') diff --git a/Httbar/scripts/systematics_breakdown.py b/Httbar/scripts/systematics_breakdown.py new file mode 100755 index 00000000000..9a14fb34d3e --- /dev/null +++ b/Httbar/scripts/systematics_breakdown.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python +from argparse import ArgumentParser +from pdb import set_trace +from fnmatch import fnmatch +from glob import glob +import shutil +import os +from subprocess import Popen, PIPE, STDOUT +from pdb import set_trace + +def call(cmd): + proc = Popen(cmd, stdout=PIPE, stderr=STDOUT, shell=True) + stdout, nothing = proc.communicate() + return proc.wait() + +def backup(src): + ret = '%s.bak' % src + shutil.copyfile(src, ret) + return ret + +def unbackup(src): + if not src.endswith('.bak'): + raise ValueError('%s is not a backup file!') + shutil.copyfile(src, src[:-4]) + return src[:-4] + +def tag(src, value): + dname, bname = os.path.dirname(src), os.path.basename(src) + ret = '%s/%s_%s' % (dname, value, bname) + shutil.copyfile(src, ret) + return ret + +parser = ArgumentParser() +parser.add_argument('--run', choices=['blind']) +args = parser.parse_args() + +txts = glob('[AH]_[0-9]*/000/combined.txt.cmb')# % args.jobid) +if not txts: + raise RuntimeError('I could not find any datacard available') +txt = txts[0] + +#get all nuisances +nuisances = set() +with open(txt) as card: + for line in card: + #the space in ' param ' is important to not select spurious things + if any((i in line for i in [' shape ', ' lnN ', ' param '])): + nuisances.add(line.split()[0]) + +groups = { + 'BinByBin' : ['*_TT_bin_[0-9]*'], + 'BTag' : ['CMS_*_b_13TeV'], + 'Leptons' : ['CMS_eff_[em]', 'CMS_eff_trigger_[em]'], + 'JetMET' : ['CMS_*_j_13TeV', 'CMS_METunclustered_13TeV'], + 'XSections' : ['CMS_httbar_*_13TeV'], + 'QCDNorm' : ['CMS_httbar_*_QCDNorm'], + 'PU' : ['CMS_pileup'], + 'Theory' : ['*_TT', 'TMass', 'pdf'], +} + +#backup files +backed_up = [backup(i) for i in glob('[AH]_[0-9]*/000/higgsCombine.limit.Asymptotic.*.root')] + +#run limits +for group, names in groups.iteritems(): + print 'running for group: ', group + tofreeze = [i for i in nuisances if any((fnmatch(i, j) for j in names))] + retval = call( + ('combineTool.py -M Asymptotic -m 400:750:50 --freezeNuisances MH,%s' + ' -d */*/workspace.root --there -n .limit --parallel 8 %s') % \ + (','.join(tofreeze), '--run blind' if args.run == 'blind' else '') + ) + if retval != 0: + raise RuntimeError('combineTool.py crashed with exit code %d' % retval) + #tag outputs + [tag(i, group) for i in glob('[AH]_[0-9]*/000/higgsCombine.limit.Asymptotic.*.root')] + retval = call( + 'combineTool.py -M CollectLimits */*/{0}_*.limit.* --use-dirs -o {0}.json'.format(group) + ) + if retval != 0: + raise RuntimeError('limit harvesting crashed with exit code %d' % retval) + +[unbackup(i) for i in backed_up]