diff --git a/bin/convert_unfolding b/bin/convert_unfolding index a6c1e46d..14f0bc12 100755 --- a/bin/convert_unfolding +++ b/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, diff --git a/bin/x_02_all_vars b/bin/x_02_all_vars index e2827ad6..80396155 100755 --- a/bin/x_02_all_vars +++ b/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 & diff --git a/bin/x_05_all_vars b/bin/x_05_all_vars index 5d8d1250..75534b37 100755 --- a/bin/x_05_all_vars +++ b/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 & diff --git a/config/CMS.py b/config/CMS.py index eb13a2ae..e8594947 100644 --- a/config/CMS.py +++ b/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 diff --git a/config/__init__.py b/config/__init__.py index 54a028ec..7a973c87 100644 --- a/config/__init__.py +++ b/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' diff --git a/config/latex_labels.py b/config/latex_labels.py index b4af5396..4114b580 100644 --- a/config/latex_labels.py +++ b/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)', diff --git a/src/cross_section_measurement/04_make_plots_matplotlib.py b/src/cross_section_measurement/04_make_plots_matplotlib.py index 0c871ded..eb1488af 100644 --- a/src/cross_section_measurement/04_make_plots_matplotlib.py +++ b/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 ) diff --git a/src/unfolding_tests/compare_unfolding_parameters.py b/src/unfolding_tests/compare_unfolding_parameters.py index 98b5906f..4e0ec61c 100644 --- a/src/unfolding_tests/compare_unfolding_parameters.py +++ b/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, diff --git a/src/unfolding_tests/tau_value_determination.py b/src/unfolding_tests/tau_value_determination.py index 55eb0114..42bdeeff 100644 --- a/src/unfolding_tests/tau_value_determination.py +++ b/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] )