def read_fit_templates_and_results_as_histograms( category, channel ):
    global path_to_JSON, variable, met_type, phase_space
    templates = read_data_from_JSON( path_to_JSON + '/fit_results/' + category + '/templates_' + channel + '_' + met_type + '.txt' )

    data_values = read_data_from_JSON( path_to_JSON + '/fit_results/' + category + '/initial_values_' + channel + '_' + met_type + '.txt' )['data']
    fit_results = read_data_from_JSON( path_to_JSON + '/fit_results/' + category + '/fit_results_' + channel + '_' + met_type + '.txt' )
    fit_variables = templates.keys()
    template_histograms = {fit_variable: {} for fit_variable in fit_variables}
    fit_results_histograms = {fit_variable: {} for fit_variable in fit_variables}

    variableBins = None
    if phase_space == 'VisiblePS':
        variableBins = variable_bins_visiblePS_ROOT
    elif phase_space == 'FullPS':
        variableBins = variable_bins_ROOT

    for bin_i, variable_bin in enumerate( variableBins[variable] ):
        for fit_variable in fit_variables:
            h_template_data = value_tuplelist_to_hist( templates[fit_variable]['data'][bin_i], fit_variable_bin_edges[fit_variable] )
            h_template_ttjet =  value_tuplelist_to_hist( templates[fit_variable]['TTJet'][bin_i], fit_variable_bin_edges[fit_variable] )
            h_template_singletop =  value_tuplelist_to_hist( templates[fit_variable]['SingleTop'][bin_i], fit_variable_bin_edges[fit_variable] )
            h_template_VJets = value_tuplelist_to_hist( templates[fit_variable]['V+Jets'][bin_i], fit_variable_bin_edges[fit_variable] )
            h_template_QCD = value_tuplelist_to_hist( templates[fit_variable]['QCD'][bin_i], fit_variable_bin_edges[fit_variable] )
            template_histograms[fit_variable][variable_bin] = {
                                        'TTJet' : h_template_ttjet,
                                        'SingleTop' : h_template_singletop,
                                        'V+Jets':h_template_VJets,
                                        'QCD':h_template_QCD
                                        }
            h_data = h_template_data.Clone()
            h_ttjet = h_template_ttjet.Clone()
            h_singletop = h_template_singletop.Clone()
            h_VJets = h_template_VJets.Clone()
            h_QCD = h_template_QCD.Clone()

            data_normalisation = data_values[bin_i][0]
            n_ttjet = fit_results['TTJet'][bin_i][0]
            n_singletop = fit_results['SingleTop'][bin_i][0]
            VJets_normalisation = fit_results['V+Jets'][bin_i][0]
            QCD_normalisation = fit_results['QCD'][bin_i][0]

            h_data.Scale( data_normalisation )
            h_ttjet.Scale( n_ttjet )
            h_singletop.Scale( n_singletop )
            h_VJets.Scale( VJets_normalisation )
            h_QCD.Scale( QCD_normalisation )
            h_background = h_VJets + h_QCD + h_singletop

            for bin_i_data in range( len( h_data ) ):
                h_data.SetBinError( bin_i_data + 1, sqrt( h_data.GetBinContent( bin_i_data + 1 ) ) )

            fit_results_histograms[fit_variable][variable_bin] = {
                                                    'data' : h_data,
                                                    'signal' : h_ttjet,
                                                    'background' : h_background
                                                    }

    return template_histograms, fit_results_histograms

def read_xsection_measurement_results(category, channel):
    global path_to_JSON, variable, k_value, met_type
    
    normalised_xsection_unfolded = None
    if category in met_uncertainties and variable == 'HT':
        normalised_xsection_unfolded = read_data_from_JSON(path_to_JSON + '/xsection_measurement_results' + '/kv' + str(k_value) + '/' 
                                                       + 'central' + '/normalised_xsection_' + channel + '_' + met_type + '.txt')
    else:
        normalised_xsection_unfolded = read_data_from_JSON(path_to_JSON + '/xsection_measurement_results' + '/kv' + str(k_value) + '/' 
                                                       + category + '/normalised_xsection_' + channel + '_' + met_type + '.txt')
        
    h_normalised_xsection = value_error_tuplelist_to_hist(normalised_xsection_unfolded['TTJet_measured'], bin_edges[variable])
    h_normalised_xsection_unfolded = value_error_tuplelist_to_hist(normalised_xsection_unfolded['TTJet_unfolded'], bin_edges[variable])
    
    
    histograms_normalised_xsection_different_generators = {'measured':h_normalised_xsection,
                                                           'unfolded':h_normalised_xsection_unfolded}
    
    histograms_normalised_xsection_systematics_shifts = {'measured':h_normalised_xsection,
                                                         'unfolded':h_normalised_xsection_unfolded}
    
    if category == 'central':
        # true distributions
        h_normalised_xsection_MADGRAPH = value_error_tuplelist_to_hist(normalised_xsection_unfolded['MADGRAPH'], bin_edges[variable])
        h_normalised_xsection_POWHEG = value_error_tuplelist_to_hist(normalised_xsection_unfolded['POWHEG'], bin_edges[variable])
        h_normalised_xsection_MCATNLO = value_error_tuplelist_to_hist(normalised_xsection_unfolded['MCATNLO'], bin_edges[variable])
        h_normalised_xsection_mathchingup = value_error_tuplelist_to_hist(normalised_xsection_unfolded['matchingup'], bin_edges[variable])
        h_normalised_xsection_mathchingdown = value_error_tuplelist_to_hist(normalised_xsection_unfolded['matchingdown'], bin_edges[variable])
        h_normalised_xsection_scaleup = value_error_tuplelist_to_hist(normalised_xsection_unfolded['scaleup'], bin_edges[variable])
        h_normalised_xsection_scaledown = value_error_tuplelist_to_hist(normalised_xsection_unfolded['scaledown'], bin_edges[variable])
        
        histograms_normalised_xsection_different_generators.update({'MADGRAPH':h_normalised_xsection_MADGRAPH,
                                                                    'POWHEG':h_normalised_xsection_POWHEG,
                                                                    'MCATNLO':h_normalised_xsection_MCATNLO})
        
        histograms_normalised_xsection_systematics_shifts.update({'MADGRAPH':h_normalised_xsection_MADGRAPH,
                                                                  'matchingdown': h_normalised_xsection_mathchingdown,
                                                                  'matchingup': h_normalised_xsection_mathchingup,
                                                                  'scaledown': h_normalised_xsection_scaledown,
                                                                  'scaleup': h_normalised_xsection_scaleup})
        
        normalised_xsection_unfolded_with_errors = read_data_from_JSON(path_to_JSON + '/xsection_measurement_results' + '/kv' + 
                                                                   str(k_value) + '/' + category + '/normalised_xsection_' + 
                                                                   channel + '_' + met_type + '_with_errors.txt')
        # a rootpy.Graph with asymmetric errors!
        h_normalised_xsection_with_systematics = value_errors_tuplelist_to_graph(normalised_xsection_unfolded_with_errors['TTJet_measured'], bin_edges[variable])
        h_normalised_xsection_with_systematics_unfolded = value_errors_tuplelist_to_graph(normalised_xsection_unfolded_with_errors['TTJet_unfolded'], bin_edges[variable])
        
        histograms_normalised_xsection_different_generators['measured_with_systematics'] = h_normalised_xsection_with_systematics
        histograms_normalised_xsection_different_generators['unfolded_with_systematics'] = h_normalised_xsection_with_systematics_unfolded
        
        histograms_normalised_xsection_systematics_shifts['measured_with_systematics'] = h_normalised_xsection_with_systematics
        histograms_normalised_xsection_systematics_shifts['unfolded_with_systematics'] = h_normalised_xsection_with_systematics_unfolded
    
    return histograms_normalised_xsection_different_generators, histograms_normalised_xsection_systematics_shifts

def compare_combine_before_after_unfolding(measurement='normalised_xsection',
                              add_before_unfolding=False):
    file_template = 'data/normalisation/background_subtraction/13TeV/'
    file_template += '{variable}/VisiblePS/central/'
    file_template += '{measurement}_{channel}_RooUnfold{method}.txt'

    variables = ['MET', 'HT', 'ST', 'NJets',
                 'lepton_pt', 'abs_lepton_eta', 'WPT']
    for variable in variables:
        combineBefore = file_template.format(
            variable=variable,
            method='Svd',
            channel='combinedBeforeUnfolding',
            measurement=measurement)
        combineAfter = file_template.format(
            variable=variable,
            method='Svd',
            channel='combined',
            measurement=measurement)
        data = read_data_from_JSON(combineBefore)
        before_unfolding = data['TTJet_measured']
        combineBefore_data = data['TTJet_unfolded']
        combineAfter_data = read_data_from_JSON(combineAfter)['TTJet_unfolded']
        h_combineBefore = value_error_tuplelist_to_hist(
            combineBefore_data, bin_edges_vis[variable])
        h_combineAfter = value_error_tuplelist_to_hist(
            combineAfter_data, bin_edges_vis[variable])
        h_before_unfolding = value_error_tuplelist_to_hist(
            before_unfolding, bin_edges_vis[variable])

        properties = Histogram_properties()
        properties.name = '{0}_compare_combine_before_after_unfolding_{1}'.format(
            measurement, variable)
        properties.title = 'Comparison of combining before/after unfolding'
        properties.path = 'plots'
        properties.has_ratio = True
        properties.xerr = True
        properties.x_limits = (
            bin_edges_vis[variable][0], bin_edges_vis[variable][-1])
        properties.x_axis_title = variables_latex[variable]
        if 'xsection' in measurement:
            properties.y_axis_title = r'$\frac{1}{\sigma}  \frac{d\sigma}{d' + \
                variables_latex[variable] + '}$'
        else:
            properties.y_axis_title = r'$t\bar{t}$ normalisation'

        histograms = {'Combine before unfolding': h_combineBefore, 'Combine after unfolding': h_combineAfter}
        if add_before_unfolding:
            histograms['before unfolding'] = h_before_unfolding
            properties.name += '_ext'
            properties.has_ratio = False
        plot = Plot(histograms, properties)
        plot.draw_method = 'errorbar'
        compare_histograms(plot)

def compare_unfolding_methods(measurement='normalised_xsection',
                              add_before_unfolding=False, channel='combined'):
    file_template = '/hdfs/TopQuarkGroup/run2/dpsData/'
    file_template += 'data/normalisation/background_subtraction/13TeV/'
    file_template += '{variable}/VisiblePS/central/'
    file_template += '{measurement}_{channel}_RooUnfold{method}.txt'

    variables = ['MET', 'HT', 'ST', 'NJets',
                 'lepton_pt', 'abs_lepton_eta', 'WPT']
    for variable in variables:
        svd = file_template.format(
            variable=variable,
            method='Svd',
            channel=channel,
            measurement=measurement)
        bayes = file_template.format(
            variable=variable,
            method='Bayes', channel=channel,
            measurement=measurement)
        data = read_data_from_JSON(svd)
        before_unfolding = data['TTJet_measured_withoutFakes']
        svd_data = data['TTJet_unfolded']
        bayes_data = read_data_from_JSON(bayes)['TTJet_unfolded']
        h_svd = value_error_tuplelist_to_hist(
            svd_data, bin_edges_vis[variable])
        h_bayes = value_error_tuplelist_to_hist(
            bayes_data, bin_edges_vis[variable])
        h_before_unfolding = value_error_tuplelist_to_hist(
            before_unfolding, bin_edges_vis[variable])

        properties = Histogram_properties()
        properties.name = '{0}_compare_unfolding_methods_{1}_{2}'.format(
            measurement, variable, channel)
        properties.title = 'Comparison of unfolding methods'
        properties.path = 'plots'
        properties.has_ratio = True
        properties.xerr = True
        properties.x_limits = (
            bin_edges_vis[variable][0], bin_edges_vis[variable][-1])
        properties.x_axis_title = variables_latex[variable]
        if 'xsection' in measurement:
            properties.y_axis_title = r'$\frac{1}{\sigma}  \frac{d\sigma}{d' + \
                variables_latex[variable] + '}$'
        else:
            properties.y_axis_title = r'$t\bar{t}$ normalisation'

        histograms = {'SVD': h_svd, 'Bayes': h_bayes}
        if add_before_unfolding:
            histograms['before unfolding'] = h_before_unfolding
            properties.name += '_ext'
            properties.has_ratio = False
        plot = Plot(histograms, properties)
        plot.draw_method = 'errorbar'
        compare_histograms(plot)

def read_fit_templates_and_results_as_histograms(category, channel):
    global path_to_JSON, variable, met_type
    templates = read_data_from_JSON(
        path_to_JSON + "/" + variable + "/fit_results/" + category + "/templates_" + channel + "_" + met_type + ".txt"
    )
    data_values = read_data_from_JSON(
        path_to_JSON
        + "/"
        + variable
        + "/fit_results/"
        + category
        + "/initial_values_"
        + channel
        + "_"
        + met_type
        + ".txt"
    )["data"]
    fit_results = read_data_from_JSON(
        path_to_JSON + "/" + variable + "/fit_results/" + category + "/fit_results_" + channel + "_" + met_type + ".txt"
    )
    template_histograms = {}
    fit_results_histograms = {}
    for bin_i, variable_bin in enumerate(variable_bins_ROOT[variable]):
        h_template_data = value_tuplelist_to_hist(templates["data"][bin_i], eta_bin_edges)
        h_template_signal = value_tuplelist_to_hist(templates["signal"][bin_i], eta_bin_edges)
        h_template_VJets = value_tuplelist_to_hist(templates["V+Jets"][bin_i], eta_bin_edges)
        h_template_QCD = value_tuplelist_to_hist(templates["QCD"][bin_i], eta_bin_edges)
        template_histograms[variable_bin] = {
            "signal": h_template_signal,
            "V+Jets": h_template_VJets,
            "QCD": h_template_QCD,
        }
        h_data = h_template_data.Clone()
        h_signal = h_template_signal.Clone()
        h_VJets = h_template_VJets.Clone()
        h_QCD = h_template_QCD.Clone()

        data_normalisation = data_values[bin_i]
        signal_normalisation = fit_results["signal"][bin_i][0]
        VJets_normalisation = fit_results["V+Jets"][bin_i][0]
        QCD_normalisation = fit_results["QCD"][bin_i][0]

        h_data.Scale(data_normalisation)
        h_signal.Scale(signal_normalisation)
        h_VJets.Scale(VJets_normalisation)
        h_QCD.Scale(QCD_normalisation)
        h_background = h_VJets + h_QCD

        for bin_i in range(len(h_data)):
            h_data.SetBinError(bin_i + 1, sqrt(h_data.GetBinContent(bin_i + 1)))

        fit_results_histograms[variable_bin] = {"data": h_data, "signal": h_signal, "background": h_background}

    return template_histograms, fit_results_histograms

def read_xsection_measurement_results_with_errors(channel):
    global path_to_JSON, variable, k_values, met_type
    category = 'central'

    file_template = path_to_JSON + '/' + variable +  '/xsection_measurement_results/' + channel + '/kv' + str(k_values[channel]) + '/' + category + '/normalised_xsection_' + met_type + '.txt' 
    if channel == 'combined':
        file_template = file_template.replace('kv' + str(k_values[channel]), '')

    file_name = file_template
    normalised_xsection_unfolded = read_data_from_JSON( file_name )
    
    normalised_xsection_measured_unfolded = {'measured':normalised_xsection_unfolded['TTJet_measured'],
                                            'unfolded':normalised_xsection_unfolded['TTJet_unfolded']}
    
    file_name = file_template.replace('.txt', '_with_errors.txt')
    normalised_xsection_unfolded_with_errors = read_data_from_JSON( file_name )

    file_name = file_template.replace('.txt', '_ttbar_generator_errors.txt')
    normalised_xsection_ttbar_generator_errors = read_data_from_JSON( file_name )

    file_name = file_template.replace('.txt', '_MET_errors.txt')
    normalised_xsection_MET_errors = read_data_from_JSON( file_name )

    file_name = file_template.replace('.txt', '_topMass_errors.txt')
    normalised_xsection_topMass_errors = read_data_from_JSON( file_name )

    file_name = file_template.replace('.txt', '_kValue_errors.txt')
    normalised_xsection_kValue_errors = read_data_from_JSON( file_name )

    file_name = file_template.replace('.txt', '_PDF_errors.txt')
    normalised_xsection_PDF_errors = read_data_from_JSON( file_name )

    file_name = file_template.replace('.txt', '_other_errors.txt')
    normalised_xsection_other_errors = read_data_from_JSON( file_name )

    file_name = file_template.replace('.txt', '_new_errors.txt')
    normalised_xsection_new_errors = read_data_from_JSON( file_name )
    
    normalised_xsection_measured_unfolded.update({'measured_with_systematics':normalised_xsection_unfolded_with_errors['TTJet_measured'],
                                                'unfolded_with_systematics':normalised_xsection_unfolded_with_errors['TTJet_unfolded']})
    
    normalised_xsection_measured_errors = normalised_xsection_ttbar_generator_errors['TTJet_measured']
    normalised_xsection_measured_errors.update(normalised_xsection_PDF_errors['TTJet_measured'])
    normalised_xsection_measured_errors.update(normalised_xsection_MET_errors['TTJet_measured'])
    normalised_xsection_measured_errors.update(normalised_xsection_topMass_errors['TTJet_measured'])
    normalised_xsection_measured_errors.update(normalised_xsection_kValue_errors['TTJet_measured'])
    normalised_xsection_measured_errors.update(normalised_xsection_other_errors['TTJet_measured'])
    normalised_xsection_measured_errors.update(normalised_xsection_new_errors['TTJet_measured'])

    normalised_xsection_unfolded_errors = normalised_xsection_ttbar_generator_errors['TTJet_unfolded']
    normalised_xsection_unfolded_errors.update(normalised_xsection_PDF_errors['TTJet_unfolded'])
    normalised_xsection_unfolded_errors.update(normalised_xsection_MET_errors['TTJet_unfolded'])
    normalised_xsection_unfolded_errors.update(normalised_xsection_topMass_errors['TTJet_unfolded'])
    normalised_xsection_unfolded_errors.update(normalised_xsection_kValue_errors['TTJet_unfolded'])
    normalised_xsection_unfolded_errors.update(normalised_xsection_other_errors['TTJet_unfolded'])
    normalised_xsection_unfolded_errors.update(normalised_xsection_new_errors['TTJet_unfolded'])

    return normalised_xsection_measured_unfolded, normalised_xsection_measured_errors, normalised_xsection_unfolded_errors

def read_normalised_xsection_measurement(category, channel):
    global path_to_JSON, met_type, met_uncertainties
    normalised_xsection = None
    
    if category in met_uncertainties and variable == 'HT':
        normalised_xsection = read_data_from_JSON(path_to_JSON + 'central' + '/normalised_xsection_' + channel + '_' + met_type + '.txt')
    else:
        normalised_xsection = read_data_from_JSON(path_to_JSON + category + '/normalised_xsection_' + channel + '_' + met_type + '.txt')
    
    measurement = normalised_xsection['TTJet_measured']
    measurement_unfolded = normalised_xsection['TTJet_unfolded']
    
    return measurement, measurement_unfolded

def compare_combine_before_after_unfolding_uncertainties():
    file_template = 'data/normalisation/background_subtraction/13TeV/'
    file_template += '{variable}/VisiblePS/central/'
    file_template += 'unfolded_normalisation_{channel}_RooUnfoldSvd.txt'

    variables = ['MET', 'HT', 'ST', 'NJets',
                 'lepton_pt', 'abs_lepton_eta', 'WPT']
#     variables = ['ST']
    for variable in variables:
        beforeUnfolding = file_template.format(
            variable=variable, channel='combinedBeforeUnfolding')
        afterUnfolding = file_template.format(
            variable=variable, channel='combined')
        data = read_data_from_JSON(beforeUnfolding)
        before_unfolding = data['TTJet_measured']
        beforeUnfolding_data = data['TTJet_unfolded']
        afterUnfolding_data = read_data_from_JSON(afterUnfolding)['TTJet_unfolded']

        before_unfolding = [e / v * 100 for v, e in before_unfolding]
        beforeUnfolding_data = [e / v * 100 for v, e in beforeUnfolding_data]
        afterUnfolding_data = [e / v * 100 for v, e in afterUnfolding_data]

        h_beforeUnfolding = value_tuplelist_to_hist(
            beforeUnfolding_data, bin_edges_vis[variable])
        h_afterUnfolding = value_tuplelist_to_hist(
            afterUnfolding_data, bin_edges_vis[variable])
        h_before_unfolding = value_tuplelist_to_hist(
            before_unfolding, bin_edges_vis[variable])

        properties = Histogram_properties()
        properties.name = 'compare_combine_before_after_unfolding_uncertainties_{0}'.format(
            variable)
        properties.title = 'Comparison of unfolding uncertainties'
        properties.path = 'plots'
        properties.has_ratio = False
        properties.xerr = True
        properties.x_limits = (
            bin_edges_vis[variable][0], bin_edges_vis[variable][-1])
        properties.x_axis_title = variables_latex[variable]
        properties.y_axis_title = 'relative uncertainty (\\%)'
        properties.legend_location = (0.98, 0.95)

        histograms = {'Combine before unfolding': h_beforeUnfolding, 'Combine after unfolding': h_afterUnfolding,
                      # 'before unfolding': h_before_unfolding
                      }
        plot = Plot(histograms, properties)
        plot.draw_method = 'errorbar'
        compare_histograms(plot)

def compare_unfolding_uncertainties():
    file_template = '/hdfs/TopQuarkGroup/run2/dpsData/'
    file_template += 'data/normalisation/background_subtraction/13TeV/'
    file_template += '{variable}/VisiblePS/central/'
    file_template += 'unfolded_normalisation_combined_RooUnfold{method}.txt'

    variables = ['MET', 'HT', 'ST', 'NJets',
                 'lepton_pt', 'abs_lepton_eta', 'WPT']
#     variables = ['ST']
    for variable in variables:
        svd = file_template.format(
            variable=variable, method='Svd')
        bayes = file_template.format(
            variable=variable, method='Bayes')
        data = read_data_from_JSON(svd)
        before_unfolding = data['TTJet_measured_withoutFakes']
        svd_data = data['TTJet_unfolded']
        bayes_data = read_data_from_JSON(bayes)['TTJet_unfolded']

        before_unfolding = [e / v * 100 for v, e in before_unfolding]
        svd_data = [e / v * 100 for v, e in svd_data]
        bayes_data = [e / v * 100 for v, e in bayes_data]

        h_svd = value_tuplelist_to_hist(
            svd_data, bin_edges_vis[variable])
        h_bayes = value_tuplelist_to_hist(
            bayes_data, bin_edges_vis[variable])
        h_before_unfolding = value_tuplelist_to_hist(
            before_unfolding, bin_edges_vis[variable])

        properties = Histogram_properties()
        properties.name = 'compare_unfolding_uncertainties_{0}'.format(
            variable)
        properties.title = 'Comparison of unfolding uncertainties'
        properties.path = 'plots'
        properties.has_ratio = False
        properties.xerr = True
        properties.x_limits = (
            bin_edges_vis[variable][0], bin_edges_vis[variable][-1])
        properties.x_axis_title = variables_latex[variable]
        properties.y_axis_title = 'relative uncertainty (\\%)'
        properties.legend_location = (0.98, 0.95)

        histograms = {'SVD': h_svd, 'Bayes': h_bayes,
                      'before unfolding': h_before_unfolding}
        plot = Plot(histograms, properties)
        plot.draw_method = 'errorbar'
        compare_histograms(plot)

def read_unfolded_xsections(channel):
    global path_to_JSON, variable, k_value, met_type, b_tag_bin
    TTJet_xsection_unfolded = {}
    for category in categories:
        normalised_xsections = read_data_from_JSON(path_to_JSON  + '/xsection_measurement_results' + '/kv' + str(k_value) + '/' + category + '/normalised_xsection_' + channel + '_' + met_type + '.txt')
        TTJet_xsection_unfolded[category] = normalised_xsections['TTJet_unfolded']
    return TTJet_xsection_unfolded

def get_fitted_normalisation(variable, channel):
    global path_to_JSON, category, met_type
    fit_results = read_data_from_JSON(path_to_JSON + variable + '/fit_results/' + category + '/fit_results_' + channel + '_' + met_type + '.txt')

    N_fit_ttjet = [0, 0]
    N_fit_singletop = [0, 0]
    N_fit_vjets = [0, 0]
    N_fit_qcd = [0, 0]

    bins = variable_bins_ROOT[variable]
    for bin_i, _ in enumerate(bins):
        #central values
        N_fit_ttjet[0] += fit_results['TTJet'][bin_i][0]
        N_fit_singletop[0] += fit_results['SingleTop'][bin_i][0]
        N_fit_vjets[0] += fit_results['V+Jets'][bin_i][0]
        N_fit_qcd[0] += fit_results['QCD'][bin_i][0]

        #errors
        N_fit_ttjet[1] += fit_results['TTJet'][bin_i][1]
        N_fit_singletop[1] += fit_results['SingleTop'][bin_i][1]
        N_fit_vjets[1] += fit_results['V+Jets'][bin_i][1]
        N_fit_qcd[1] += fit_results['QCD'][bin_i][1]

    fitted_normalisation = {
                'TTJet': N_fit_ttjet,
                'SingleTop': N_fit_singletop,
                'V+Jets': N_fit_vjets,
                'QCD': N_fit_qcd
                }
    return fitted_normalisation

def get_fitted_normalisation( variable, channel, path_to_JSON, category, met_type ):
    '''
    This function now gets the error on the fit correctly,
    so that it can be applied if the --normalise_to_fit option is used
    '''
    import config.variable_binning
    variable_bins_ROOT = config.variable_binning.variable_bins_ROOT 
    fit_results = read_data_from_JSON( path_to_JSON + variable + '/fit_results/' + category + '/fit_results_' + channel + '_' + met_type + '.txt' )

    N_fit_ttjet = [0, 0]
    N_fit_singletop = [0, 0]
    N_fit_vjets = [0, 0]
    N_fit_qcd = [0, 0]

    bins = variable_bins_ROOT[variable]
    for bin_i, _ in enumerate( bins ):
        # central values
        N_fit_ttjet[0] += fit_results['TTJet'][bin_i][0]
        N_fit_singletop[0] += fit_results['SingleTop'][bin_i][0]
        N_fit_vjets[0] += fit_results['V+Jets'][bin_i][0]
        N_fit_qcd[0] += fit_results['QCD'][bin_i][0]

        # errors
        N_fit_ttjet[1] += fit_results['TTJet'][bin_i][1]
        N_fit_singletop[1] += fit_results['SingleTop'][bin_i][1]
        N_fit_vjets[1] += fit_results['V+Jets'][bin_i][1]
        N_fit_qcd[1] += fit_results['QCD'][bin_i][1]

    fitted_normalisation = {
                'TTJet': N_fit_ttjet,
                'SingleTop': N_fit_singletop,
                'V+Jets': N_fit_vjets,
                'QCD': N_fit_qcd
                }
    return fitted_normalisation

def get_data_histogram( channel, variable, met_type ):
    fit_result_input = 'data/M3_angle_bl/13TeV/%(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']
    print fit_data
    print bin_edges[variable]
    h_data = value_error_tuplelist_to_hist( fit_data, bin_edges[variable] )
    return h_data

def main(options, args):
    config = XSectionConfig(options.CoM)
    variables = ['MET', 'HT', 'ST', 'WPT']
    channels = ['electron', 'muon', 'combined']
    m_file = 'normalised_xsection_patType1CorrectedPFMet.txt'
    m_with_errors_file = 'normalised_xsection_patType1CorrectedPFMet_with_errors.txt'
    path_template = args[0]
    output_file = 'measurement_{0}TeV.root'.format(options.CoM)
    f = File(output_file, 'recreate')
    for channel in channels:
        d = f.mkdir(channel)
        d.cd()
        for variable in variables:
            dv = d.mkdir(variable)
            dv.cd()
            if channel == 'combined':
                path = path_template.format(variable=variable,
                                            channel=channel,
                                            centre_of_mass_energy=options.CoM)
            else:
                kv = channel + \
                    '/kv{0}/'.format(config.k_values[channel][variable])
                path = path_template.format(variable=variable,
                                            channel=kv,
                                            centre_of_mass_energy=options.CoM)

            m = read_data_from_JSON(path + '/' + m_file)
            m_with_errors = read_data_from_JSON(
                path + '/' + m_with_errors_file)

            for name, result in m.items():
                h = make_histogram(result, bin_edges[variable])
                h.SetName(name)
                h.write()

            for name, result in m_with_errors.items():
                if not 'TTJet' in name:
                    continue
                h = make_histogram(result, bin_edges[variable])
                h.SetName(name + '_with_syst')
                h.write()
            dv.write()
            d.cd()
        d.write()
    f.write()
    f.close()

def main():
    global config, options
    parser = OptionParser()
    parser.add_option( "-p", "--path", dest = "path", default = 'data/fit_checks/no_merging',
                  help = "set path to JSON files" )
    parser.add_option( '--create_fit_data', dest = "create_fit_data", action = "store_true",
                      help = "create the fit data for testing." )
    parser.add_option( '--refit', dest = "refit", action = "store_true",
                      help = "Fit again even if the output already exists" )
    parser.add_option( '--test', dest = "test", action = "store_true",
                      help = "Test only: run just one selected sample" )
    variables = config.histogram_path_templates.keys()
    fit_variables = fit_var_inputs
    mc_samples = ['TTJet', 'SingleTop', 'QCD', 'V+Jets']
    tests = closure_tests
    channels = ['electron', 'muon']
    COMEnergies = ['7', '8']
    
    ( options, _ ) = parser.parse_args()
    print 'Running from path', options.path 
    
    if ( options.create_fit_data ):
        create_fit_data( options.path, variables, fit_variables, mc_samples,
                        COMEnergies = COMEnergies, channels = channels )
     
    output_file = options.path + '/fit_test_output.txt'
    if options.test:
        output_file = options.path + '/fit_test_output_test.txt'
    if options.refit or not os.path.isfile( output_file ) or options.test:
        if os.path.isfile( options.path + '/fit_check_data.txt' ):
            fit_data = read_data_from_JSON( options.path + '/fit_check_data.txt' )
            results = run_tests( fit_data,
                                COMEnergies,
                                variables,
                                fit_variables,
                                mc_samples,
                                channels,
                                tests )
            write_data_to_JSON(results, output_file )
        else:
            print 'Please run bin/prepare_data_for_fit_checks first'
            print 'Then run this script with the option --create_fit_data.'
     
    results = read_data_from_JSON(output_file)
    plot_results( results )

def read_fit_templates_and_results_as_histograms(category, channel):
    global path_to_JSON, variable, met_type
    templates = read_data_from_JSON(path_to_JSON  + '/fit_results/' + category + '/templates_' + channel + '_' + met_type + '.txt')
    data_values = read_data_from_JSON(path_to_JSON  + '/fit_results/' + category + '/initial_values_' + channel + '_' + met_type + '.txt')['data']
    fit_results = read_data_from_JSON(path_to_JSON  + '/fit_results/' + category + '/fit_results_' + channel + '_' + met_type + '.txt')
    template_histograms = {}
    fit_results_histograms = {}
    for bin_i, variable_bin in enumerate(variable_bins_ROOT[variable]):
        h_template_data = value_tuplelist_to_hist(templates['data'][bin_i], eta_bin_edges)
        h_template_signal = value_tuplelist_to_hist(templates['signal'][bin_i], eta_bin_edges)
        h_template_VJets = value_tuplelist_to_hist(templates['V+Jets'][bin_i], eta_bin_edges)
        h_template_QCD = value_tuplelist_to_hist(templates['QCD'][bin_i], eta_bin_edges)
        template_histograms[variable_bin] = {
                                    'signal':h_template_signal,
                                    'V+Jets':h_template_VJets,
                                    'QCD':h_template_QCD
                                    }
        h_data = h_template_data.Clone()
        h_signal = h_template_signal.Clone()
        h_VJets = h_template_VJets.Clone()
        h_QCD = h_template_QCD.Clone()
        
        data_normalisation = data_values[bin_i]
        signal_normalisation = fit_results['signal'][bin_i][0]
        VJets_normalisation = fit_results['V+Jets'][bin_i][0]
        QCD_normalisation = fit_results['QCD'][bin_i][0]
        
        h_data.Scale(data_normalisation)
        h_signal.Scale(signal_normalisation)
        h_VJets.Scale(VJets_normalisation)
        h_QCD.Scale(QCD_normalisation)
        h_background = h_VJets + h_QCD
        
        for bin_i in range(len(h_data)):
            h_data.SetBinError(bin_i+1, sqrt(h_data.GetBinContent(bin_i+1)))
        
        fit_results_histograms[variable_bin] = {
                                                'data':h_data,
                                                'signal':h_signal,
                                                'background':h_background
                                                }
        
    return template_histograms, fit_results_histograms

def read_xsection_measurement_results(category, channel):
    global path_to_JSON, variable, k_value, met_type
    
    normalised_xsection_unfolded = read_data_from_JSON(path_to_JSON  + '/xsection_measurement_results' + '/kv' + str(k_value) + '/' 
                                                       + category + '/normalised_xsection_' + channel + '_' + met_type + '.txt')
    h_normalised_xsection = value_error_tuplelist_to_hist(normalised_xsection_unfolded['TTJet_measured'], bin_edges[variable])
    h_normalised_xsection_unfolded = value_error_tuplelist_to_hist(normalised_xsection_unfolded['TTJet_unfolded'], bin_edges[variable])
    
    if category == 'central':    
        normalised_xsection_unfolded_with_errors = read_data_from_JSON(path_to_JSON  + '/xsection_measurement_results' + '/kv' + 
                                                                   str(k_value) + '/' + category + '/normalised_xsection_' + 
                                                                   channel + '_' + met_type + '_with_errors.txt')
        h_normalised_xsection = value_errors_tuplelist_to_graph(normalised_xsection_unfolded_with_errors['TTJet_measured'], bin_edges[variable])
        h_normalised_xsection_unfolded = value_errors_tuplelist_to_graph(normalised_xsection_unfolded_with_errors['TTJet_unfolded'], bin_edges[variable])
    
    
    #true distributions
    h_normalised_xsection_MADGRAPH = value_error_tuplelist_to_hist(normalised_xsection_unfolded['MADGRAPH'], bin_edges[variable])
    h_normalised_xsection_POWHEG = value_error_tuplelist_to_hist(normalised_xsection_unfolded['POWHEG'], bin_edges[variable])
    h_normalised_xsection_MCATNLO = value_error_tuplelist_to_hist(normalised_xsection_unfolded['MCATNLO'], bin_edges[variable])
    h_normalised_xsection_mathchingup = value_error_tuplelist_to_hist(normalised_xsection_unfolded['matchingup'], bin_edges[variable])
    h_normalised_xsection_mathchingdown = value_error_tuplelist_to_hist(normalised_xsection_unfolded['matchingdown'], bin_edges[variable])
    h_normalised_xsection_scaleup = value_error_tuplelist_to_hist(normalised_xsection_unfolded['scaleup'], bin_edges[variable])
    h_normalised_xsection_scaledown = value_error_tuplelist_to_hist(normalised_xsection_unfolded['scaledown'], bin_edges[variable])
    
    histograms_normalised_xsection_different_generators = {
                  'measured':h_normalised_xsection,
                  'unfolded':h_normalised_xsection_unfolded,
                  'MADGRAPH':h_normalised_xsection_MADGRAPH,
                  'POWHEG':h_normalised_xsection_POWHEG,
                  'MCATNLO':h_normalised_xsection_MCATNLO
                  }
    
    histograms_normalised_xsection_systematics_shifts = {
                  'measured':h_normalised_xsection,
                  'unfolded':h_normalised_xsection_unfolded,
                  'matchingdown': h_normalised_xsection_mathchingdown,
                  'matchingup': h_normalised_xsection_mathchingup,
                  'scaledown': h_normalised_xsection_scaledown,
                  'scaleup': h_normalised_xsection_scaleup
                  }
    
    return histograms_normalised_xsection_different_generators, histograms_normalised_xsection_systematics_shifts

def get_fit_results_histogram( data_path = 'data/M3_angle_bl',
                               centre_of_mass = 8,
                               channel = 'electron',
                               variable = 'MET',
                               met_type = 'patType1CorrectedPFMet',
                               bin_edges = [] ):
    fit_result_input = data_path + '/%(CoM)dTeV/%(variable)s/fit_results/central/fit_results_%(channel)s_%(met_type)s.txt'
    fit_results = read_data_from_JSON( fit_result_input % {'CoM': centre_of_mass, 'channel': channel, 'variable': variable, 'met_type':met_type} )
    fit_data = fit_results['TTJet']
    h_data = value_error_tuplelist_to_hist( fit_data, bin_edges )
    return h_data