Merge pull request #90 from senkin/ratio_plots_etc

Added the ratio plots, updated k-values, fixed some cosmetics

bin/convert_unfolding

@@ -5,11 +5,13 @@
 # in read speed between fine-binned and asymmetric is a factor of 200!
 # TODO: create the combined histograms as well.
 from src.cross_section_measurement.lib import convert_unfolding_histograms
-import config.cross_section_measurement_7TeV as config_7TeV
-import config.cross_section_measurement_8TeV as config_8TeV
+from config import XSectionConfig
 from tools.Timer import Timer
 from multiprocessing import Pool
 
+config_7TeV = XSectionConfig(7)
+config_8TeV = XSectionConfig(8)
+
 # it takes ~ 20 min to do all files
 files_to_load = [config_7TeV.unfolding_madgraph_raw,
                  config_7TeV.unfolding_matching_down_raw,

bin/x_02_all_vars

@@ -1,12 +1,12 @@
 #!/bin/bash
-python src/cross_section_measurement/02_unfold_and_measure.py --combine-before-unfolding &> logs/MET_unfold_8TeV.log &
-python src/cross_section_measurement/02_unfold_and_measure.py --combine-before-unfolding -v HT &> logs/HT_unfold_8TeV.log &
-python src/cross_section_measurement/02_unfold_and_measure.py --combine-before-unfolding -v ST &> logs/ST_unfold_8TeV.log &
-python src/cross_section_measurement/02_unfold_and_measure.py --combine-before-unfolding -v MT &> logs/MT_unfold_8TeV.log &
-python src/cross_section_measurement/02_unfold_and_measure.py --combine-before-unfolding -v WPT &> logs/WPT_unfold_8TeV.log &
+python src/cross_section_measurement/02_unfold_and_measure.py &> logs/MET_unfold_8TeV.log &
+python src/cross_section_measurement/02_unfold_and_measure.py -v HT &> logs/HT_unfold_8TeV.log &
+python src/cross_section_measurement/02_unfold_and_measure.py -v ST &> logs/ST_unfold_8TeV.log &
+python src/cross_section_measurement/02_unfold_and_measure.py -v MT &> logs/MT_unfold_8TeV.log &
+python src/cross_section_measurement/02_unfold_and_measure.py -v WPT &> logs/WPT_unfold_8TeV.log &
 # 7 TeV
-python src/cross_section_measurement/02_unfold_and_measure.py --combine-before-unfolding -c 7 &> logs/MET_unfold_7TeV.log &
-python src/cross_section_measurement/02_unfold_and_measure.py --combine-before-unfolding -c 7 -v HT &> logs/HT_unfold_7TeV.log &
-python src/cross_section_measurement/02_unfold_and_measure.py --combine-before-unfolding -c 7 -v ST &> logs/ST_unfold_7TeV.log &
-python src/cross_section_measurement/02_unfold_and_measure.py --combine-before-unfolding -c 7 -v MT &> logs/MT_unfold_7TeV.log &
-python src/cross_section_measurement/02_unfold_and_measure.py --combine-before-unfolding -c 7 -v WPT &> logs/WPT_unfold_7TeV.log &
+python src/cross_section_measurement/02_unfold_and_measure.py -c 7 &> logs/MET_unfold_7TeV.log &
+python src/cross_section_measurement/02_unfold_and_measure.py -c 7 -v HT &> logs/HT_unfold_7TeV.log &
+python src/cross_section_measurement/02_unfold_and_measure.py -c 7 -v ST &> logs/ST_unfold_7TeV.log &
+python src/cross_section_measurement/02_unfold_and_measure.py -c 7 -v MT &> logs/MT_unfold_7TeV.log &
+python src/cross_section_measurement/02_unfold_and_measure.py -c 7 -v WPT &> logs/WPT_unfold_7TeV.log &

bin/x_05_all_vars

@@ -1,5 +1,4 @@
 #!/bin/bash
-mkdir -p plots
 python src/cross_section_measurement/05_make_tables.py &> logs/MET_table_8TeV.log &
 python src/cross_section_measurement/05_make_tables.py -v HT &> logs/HT_table_8TeV.log &
 python src/cross_section_measurement/05_make_tables.py -v ST &> logs/ST_table_8TeV.log &

config/CMS.py

@@ -21,6 +21,13 @@
               'horizontalalignment': 'right'
               }
 
+y_axis_title_small = {
+              'fontsize':40,
+              'position' : (0, 1.),
+              'verticalalignment': 'bottom',
+              'horizontalalignment': 'right'
+              }
+
 axis_label_major = {
               'which':'major',
               'labelsize':40,
@@ -36,7 +43,7 @@
               'pad': 12
               }
 
-legend_properties = {'size':40}
+legend_properties = {'size':35}
 
 figsize = (16,16)
 dpi = 400

config/__init__.py

@@ -182,18 +182,18 @@ def __fill_defaults__( self ):
         # optimal regularisation parameters
         self.k_values_electron = {
                    'MET' : 3,
-                   'HT' : 7,
-                   'ST' : 6,
+                   'HT' : 4,
+                   'ST' : 3,
                    'MT' : 2,
                    'WPT' : 3
                    }
 
         self.k_values_muon = {
                    'MET' : 3,
-                   'HT' : 6,
-                   'ST' : 6,
+                   'HT' : 4,
+                   'ST' : 3,
                    'MT' : 2,
-                   'WPT' : 6
+                   'WPT' : 3
                    }
 
         self.k_values_combined = {
@@ -242,7 +242,7 @@ def __fill_defaults_8TeV__( self ):
         middle = self.middle
         path_to_files = self.path_to_files
 
-        self.new_luminosity = self.luminosity  # pb^-1
+        self.new_luminosity = 19712  # pb^-1
         self.ttbar_xsection = 245.8  # pb
 
         self.data_file_electron = path_to_files + 'central/SingleElectron' + middle + '.root'

config/latex_labels.py

@@ -18,9 +18,9 @@
 
 measurements_latex = {'unfolded': 'unfolded',
                         'measured': 'measured',
-                        'MADGRAPH': '$t\\bar{t}$ (MadGraph)',
-                        'MCATNLO': '$t\\bar{t}$ (MC@NLO)',
-                        'POWHEG': '$t\\bar{t}$ (POWHEG)',
+                        'MADGRAPH': '$t\\bar{t}$ (MadGraph+Pythia)',
+                        'MCATNLO': '$t\\bar{t}$ (MC@NLO+Herwig)',
+                        'POWHEG': '$t\\bar{t}$ (POWHEG+Pythia)',
                         'matchingdown': '$t\\bar{t}$ (matching down)',
                         'matchingup': '$t\\bar{t}$ (matching up)',
                         'scaledown': '$t\\bar{t}$ ($Q^{2}$ down)',

src/cross_section_measurement/04_make_plots_matplotlib.py

@@ -9,7 +9,7 @@
 from config import XSectionConfig
 from tools.file_utilities import read_data_from_JSON, make_folder_if_not_exists
 from tools.hist_utilities import value_error_tuplelist_to_hist, \
-value_tuplelist_to_hist, value_errors_tuplelist_to_graph
+value_tuplelist_to_hist, value_errors_tuplelist_to_graph, graph_to_value_errors_tuplelist
 from math import sqrt
 # rootpy & matplotlib
 from ROOT import kRed, kGreen, kMagenta, kBlue
@@ -18,6 +18,8 @@
 mpl.use( 'agg' )
 import rootpy.plotting.root2matplotlib as rplt
 import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+from matplotlib.ticker import MultipleLocator
 from config import CMS
 from matplotlib import rc
 rc( 'font', **CMS.font )
@@ -258,7 +260,7 @@ def get_cms_labels( channel ):
     return label
 
 
-def make_plots( histograms, category, output_folder, histname, show_before_unfolding = False ):
+def make_plots( histograms, category, output_folder, histname, show_ratio = True, show_before_unfolding = False ):
     global variable, k_values
 
     channel = 'electron'
@@ -286,11 +288,16 @@ def make_plots( histograms, category, output_folder, histname, show_before_unfol
     hist_measured.marker = 'o'
     hist_measured.color = 'red'
 
-    plt.figure( figsize = ( 16, 16 ), dpi = 200, facecolor = 'white' )
-    axes = plt.axes()
+    plt.figure( figsize = CMS.figsize, dpi = CMS.dpi, facecolor = CMS.facecolor )
+    if show_ratio:
+        gs = gridspec.GridSpec( 2, 1, height_ratios = [5, 1] ) 
+        axes = plt.subplot( gs[0] )
+    else:
+        axes = plt.axes()
+        plt.xlabel( '$%s$ [GeV]' % variables_latex[variable], CMS.x_axis_title )
+
     axes.minorticks_on()
 
-    plt.xlabel( '$%s$ [GeV]' % variables_latex[variable], CMS.x_axis_title )
     plt.ylabel( r'$\frac{1}{\sigma}  \frac{d\sigma}{d' + variables_latex[variable] + '} \left[\mathrm{GeV}^{-1}\\right]$', CMS.y_axis_title )
     plt.tick_params( **CMS.axis_label_major )
     plt.tick_params( **CMS.axis_label_minor )
@@ -341,9 +348,64 @@ def make_plots( histograms, category, output_folder, histname, show_before_unfol
             new_handles.append( handle )
             new_labels.append( label )
 
-    plt.legend( new_handles, new_labels, numpoints = 1, loc = 'upper right', prop = CMS.legend_properties )
+    legend_location = 'upper right'
+    if variable == 'MT':
+        legend_location = 'upper left'
+    plt.legend( new_handles, new_labels, numpoints = 1, loc = legend_location, prop = CMS.legend_properties )
     plt.title( get_cms_labels( channel ), CMS.title )
-    plt.tight_layout()
+
+    if show_ratio:
+        plt.setp( axes.get_xticklabels(), visible = False )
+        ax1 = plt.subplot( gs[1] )
+        ax1.minorticks_on()
+        #ax1.grid( True, 'major', linewidth = 1 )
+        # setting the x_limits identical to the main plot
+        x_limits = axes.get_xlim()
+        ax1.set_xlim(x_limits)
+        ax1.yaxis.set_major_locator( MultipleLocator( 0.5 ) )
+        ax1.yaxis.set_minor_locator( MultipleLocator( 0.1 ) )
+
+        plt.xlabel( '$%s$ [GeV]' % variables_latex[variable], CMS.x_axis_title )
+        plt.tick_params( **CMS.axis_label_major )
+        plt.tick_params( **CMS.axis_label_minor ) 
+        plt.ylabel( '$\\frac{\\textrm{theory}}{\\textrm{data}}$', CMS.y_axis_title_small )
+        ax1.yaxis.set_label_coords(-0.115, 0.8)
+        #draw a horizontal line at y=1 for data
+        plt.axhline(y = 1, color = 'black', linewidth = 1)
+
+        for key, hist in sorted( histograms.iteritems() ):
+            if not 'unfolded' in key and not 'measured' in key:
+                ratio = hist.Clone()
+                ratio.Divide( hist_data ) #divide by data
+                rplt.hist( ratio, axes = ax1, label = 'do_not_show' )
+
+        stat_lower = hist_data.Clone()
+        stat_upper = hist_data.Clone()
+        syst_lower = hist_data.Clone()
+        syst_upper = hist_data.Clone()
+
+        # plot error bands on data in the ratio plot
+        for bin_i in range( 1, hist_data.GetNbinsX() + 1 ):
+            stat_errors = graph_to_value_errors_tuplelist(hist_data)
+            stat_lower.SetBinContent( bin_i, 1 - stat_errors[bin_i-1][1]/stat_errors[bin_i-1][0] )
+            stat_upper.SetBinContent( bin_i, 1 + stat_errors[bin_i-1][2]/stat_errors[bin_i-1][0] )
+            if category == 'central':
+                syst_errors = graph_to_value_errors_tuplelist(hist_data_with_systematics)
+                syst_lower.SetBinContent( bin_i, 1 - syst_errors[bin_i-1][1]/syst_errors[bin_i-1][0] )
+                syst_upper.SetBinContent( bin_i, 1 + syst_errors[bin_i-1][2]/syst_errors[bin_i-1][0] )
+
+        if category == 'central':
+            rplt.fill_between( syst_lower, syst_upper, ax1, facecolor = 'yellow', alpha = 0.5 )
+
+        rplt.fill_between( stat_upper, stat_lower, ax1, facecolor = '0.75', alpha = 0.5 )
+
+        # p1 = plt.Rectangle((0, 0), 1, 1, fc="yellow")
+        # p2 = plt.Rectangle((0, 0), 1, 1, fc="0.75") 
+        # plt.legend([p1, p2], ['Stat. $\\oplus$ Syst.', 'Stat.'], loc = 'upper left', prop = {'size':20})
+        ax1.set_ylim( ymin = 0.5, ymax = 1.5 )
+
+    if CMS.tight_layout:
+        plt.tight_layout()
 
     path = output_folder + str( measurement_config.centre_of_mass_energy ) + 'TeV/' + variable + '/' + category
     make_folder_if_not_exists( path )

src/unfolding_tests/compare_unfolding_parameters.py

@@ -68,7 +68,7 @@ def run_test( h_truth, h_measured, h_response, h_data, h_fakes = None, variable
         result = calculate_normalised_xsection( 
                         hist_to_value_error_tuplelist( unfolded_data ),
                         bin_widths[variable],
-                        normalise_to_one = True )
+                        normalise_to_one = False )
         h_result = value_error_tuplelist_to_hist( result, bin_edges[variable] )
         k_value_results[k_value] = deepcopy( h_result )
 
@@ -87,7 +87,7 @@ def run_test( h_truth, h_measured, h_response, h_data, h_fakes = None, variable
             result = calculate_normalised_xsection( 
                             hist_to_value_error_tuplelist( unfolded_data ),
                             bin_widths[variable],
-                            normalise_to_one = True )
+                            normalise_to_one = False)
             h_result = value_error_tuplelist_to_hist( result, bin_edges[variable] )
             tau_value_results[tau_value] = deepcopy( h_result )
 
@@ -112,11 +112,11 @@ def compare( central_mc, expected_result = None, measured_result = None, results
         title_template = 'CMS Simulation at $\sqrt{s}$ = %d TeV \n %s'
         title = title_template % ( centre_of_mass, channel_label )
 
-    models = {'central' : central_mc}
+    models = {latex_labels.measurements_latex['MADGRAPH'] : central_mc}
     if expected_result:
-        models['expected'] = expected_result
-    if measured_result:
-        models['measured'] = measured_result
+        models.update({'expected' : expected_result})
+    if measured_result and test != 'data':
+        models.update({'measured' : measured_result})
 
     measurements = collections.OrderedDict()
     for key, value in results['k_value_results'].iteritems():
@@ -213,6 +213,7 @@ def compare( central_mc, expected_result = None, measured_result = None, results
 
             if test == 'data':
                 h_data = get_data_histogram( channel, variable, met_type )
+                h_expected = deepcopy( h_data )
             elif test == 'bias':
                 h_truth_bias, h_measured_bias, _, h_fakes = get_unfold_histogram_tuple( 
                                 inputfile = input_file_bias,
@@ -233,19 +234,19 @@ def compare( central_mc, expected_result = None, measured_result = None, results
             central_mc_result = calculate_normalised_xsection( 
                             hist_to_value_error_tuplelist( h_truth ),
                             bin_widths[variable],
-                            normalise_to_one = True )
+                            normalise_to_one = False )
             central_mc = value_error_tuplelist_to_hist( central_mc_result, bin_edges[variable] )
             if h_expected:
                 expected_result = calculate_normalised_xsection( 
                             hist_to_value_error_tuplelist( h_expected ),
                             bin_widths[variable],
-                            normalise_to_one = True )
+                            normalise_to_one = False )
                 h_expected = value_error_tuplelist_to_hist( expected_result, bin_edges[variable] )
 
             data_result = calculate_normalised_xsection( 
                             hist_to_value_error_tuplelist( h_data ),
                             bin_widths[variable],
-                            normalise_to_one = True )
+                            normalise_to_one = False )
             h_data = value_error_tuplelist_to_hist( data_result, bin_edges[variable] )
 
             compare( central_mc = central_mc, expected_result = h_expected, measured_result = h_data,

src/unfolding_tests/tau_value_determination.py

@@ -163,7 +163,7 @@ def draw_l_shape( l_shape, best_tau, x_value, y_value, channel, variable ):
         plt.savefig( 'plots/tau_from_L_shape_%s_channel_%s_MC.png' % ( channel, variable ) )
 
 def get_data_histogram( channel, variable, met_type ):
-    fit_result_input = '../cross_section_measurement/data/8TeV/%(variable)s/fit_results/central/fit_results_%(channel)s_%(met_type)s.txt'
+    fit_result_input = '../../data/8TeV/%(variable)s/fit_results/central/fit_results_%(channel)s_%(met_type)s.txt'
     fit_results = read_data_from_JSON( fit_result_input % {'channel': channel, 'variable': variable, 'met_type':met_type} )
     fit_data = fit_results['TTJet']
     h_data = value_error_tuplelist_to_hist( fit_data, bin_edges[variable] )