Merge pull request cms-analysis#1 from mverzett/httbar_74

Width interpolation updates and standard flow

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

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
 
 

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])

Httbar/results/Makefile

@@ -0,0 +1,42 @@
+jobid=1D_2017Mar23
+runmode='--run blind'
+limitdir=output$(jobid)
+
+#../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
+
+$(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 $@
+
+$(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

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

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