def plot_fit_results(histograms, category, channel):
    global variable, b_tag_bin, output_folder
    from tools.plotting import Histogram_properties, make_data_mc_comparison_plot
    
    for variable_bin in variable_bins_ROOT[variable]:
        path = output_folder + str(measurement_config.centre_of_mass) + 'TeV/' + variable + '/' + category + '/fit_results/'
        make_folder_if_not_exists(path)
        plotname = channel + '_bin_' + variable_bin
        # check if template plots exist already
        for output_format in output_formats:
            if os.path.isfile(plotname + '.' + output_format):
                continue
            
        # plot with matplotlib
        h_data = histograms[variable_bin]['data']
        h_signal = histograms[variable_bin]['signal']
        h_background = histograms[variable_bin]['background']
        
        histogram_properties = Histogram_properties()
        histogram_properties.name = plotname
        histogram_properties.x_axis_title = channel + ' $\left|\eta\\right|$'
        histogram_properties.y_axis_title = 'events/0.2'
        histogram_properties.title = get_cms_labels(channel)
        
        make_data_mc_comparison_plot([h_data, h_background, h_signal], 
                                     ['data', 'background', 'signal'], 
                                     ['black', 'green', 'red'], histogram_properties, 
                                     save_folder = path, save_as = output_formats)    

def plot_closure(unfolded_and_truths, variable, channel, come, method):
    hp = Histogram_properties()
    hp.name = '{channel}_closure_test_for_{variable}_at_{come}TeV'.format(
        channel=channel,
        variable=variable,
        come=come,
    )
    v_latex = latex_labels.variables_latex[variable]
    unit = ''
    if variable in ['HT', 'ST', 'MET', 'WPT', 'lepton_pt']:
        unit = ' [GeV]'
    hp.x_axis_title = v_latex + unit
    # plt.ylabel( r, CMS.y_axis_title )
    hp.y_axis_title = r'$\frac{1}{\sigma}  \frac{d\sigma}{d' + v_latex + '}$' + unit
    hp.title = 'Closure tests for {variable}'.format(variable=v_latex)

    output_folder = 'plots/unfolding/closure_test/{0}/'.format(method)

    models = OrderedDict()
    measurements = OrderedDict()
    for sample in unfolded_and_truths:
        models[sample + ' truth'] = unfolded_and_truths[sample]['truth']
        measurements[sample + ' unfolded'] = unfolded_and_truths[sample]['unfolded']


    compare_measurements(
                         models = models,
                         measurements = measurements,
                         show_measurement_errors=True,
                         histogram_properties=hp,
                         save_folder=output_folder,
                         save_as=['pdf'],
                         match_models_to_measurements = True)

def plot_bias(h_unfold_model, h_data_model, unfolded_data, variable,
              channel, come, method):
    hp = Histogram_properties()
    hp.name = '{channel}_bias_test_for_{variable}_at_{come}TeV'.format(
        channel=channel,
        variable=variable,
        come=come,
    )
    v_latex = latex_labels.variables_latex[variable]
    unit = ''
    if variable in ['HT', 'ST', 'MET', 'WPT']:
        unit = ' [GeV]'
    hp.x_axis_title = v_latex + unit
    hp.y_axis_title = 'Events'
    hp.title = 'Closure tests for {variable}'.format(variable=v_latex)
    
    output_folder = 'plots/unfolding/bias_test/{0}/'.format(method)

    compare_measurements(models={'MC truth': h_data_model,
                                 'unfold model': h_unfold_model},
                         measurements={'unfolded reco': unfolded_data},
                         show_measurement_errors=True,
                         histogram_properties=hp,
                         save_folder=output_folder,
                         save_as=['png', 'pdf'])

def compare_vjets_templates( variable = 'MET', met_type = 'patType1CorrectedPFMet',
                             title = 'Untitled', channel = 'electron' ):
    ''' Compares the V+jets templates in different bins
     of the current variable'''
    global fit_variable_properties, b_tag_bin, save_as
    variable_bins = variable_bins_ROOT[variable]
    histogram_template = get_histogram_template( variable )
    
    for fit_variable in electron_fit_variables:
        all_hists = {}
        inclusive_hist = None
        save_path = 'plots/%dTeV/fit_variables/%s/%s/' % ( measurement_config.centre_of_mass_energy, variable, fit_variable )
        make_folder_if_not_exists( save_path + '/vjets/' )
        
        max_bins = len( variable_bins )
        for bin_range in variable_bins[0:max_bins]:
            
            params = {'met_type': met_type, 'bin_range':bin_range, 'fit_variable':fit_variable, 'b_tag_bin':b_tag_bin, 'variable':variable}
            fit_variable_distribution = histogram_template % params
            # format: histograms['data'][qcd_fit_variable_distribution]
            histograms = get_histograms_from_files( [fit_variable_distribution], histogram_files )
            prepare_histograms( histograms, rebin = fit_variable_properties[fit_variable]['rebin'], scale_factor = measurement_config.luminosity_scale )
            all_hists[bin_range] = histograms['V+Jets'][fit_variable_distribution]
    
        # create the inclusive distributions
        inclusive_hist = deepcopy( all_hists[variable_bins[0]] )
        for bin_range in variable_bins[1:max_bins]:
            inclusive_hist += all_hists[bin_range]
        for bin_range in variable_bins[0:max_bins]:
            if not all_hists[bin_range].Integral() == 0:
                all_hists[bin_range].Scale( 1 / all_hists[bin_range].Integral() )
        # normalise all histograms
        inclusive_hist.Scale( 1 / inclusive_hist.Integral() )
        # now compare inclusive to all bins
        histogram_properties = Histogram_properties()
        histogram_properties.x_axis_title = fit_variable_properties[fit_variable]['x-title']
        histogram_properties.y_axis_title = fit_variable_properties[fit_variable]['y-title']
        histogram_properties.y_axis_title = histogram_properties.y_axis_title.replace( 'Events', 'a.u.' )
        histogram_properties.x_limits = [fit_variable_properties[fit_variable]['min'], fit_variable_properties[fit_variable]['max']]
        histogram_properties.title = title
        histogram_properties.additional_text = channel_latex[channel] + ', ' + b_tag_bins_latex[b_tag_bin]
        histogram_properties.name = variable + '_' + fit_variable + '_' + b_tag_bin + '_VJets_template_comparison'
        histogram_properties.y_max_scale = 1.5
        measurements = {bin_range + ' GeV': histogram for bin_range, histogram in all_hists.iteritems()}
        measurements = OrderedDict( sorted( measurements.items() ) )
        fit_var = fit_variable.replace( 'electron_', '' )
        fit_var = fit_var.replace( 'muon_', '' )
        graphs = spread_x( measurements.values(), fit_variable_bin_edges[fit_var] )
        for key, graph in zip( sorted( measurements.keys() ), graphs ):
            measurements[key] = graph
        compare_measurements( models = {'inclusive' : inclusive_hist},
                             measurements = measurements,
                             show_measurement_errors = True,
                             histogram_properties = histogram_properties,
                             save_folder = save_path + '/vjets/',
                             save_as = save_as )

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 plot( self, plot_options = {} ):
     '''
         Plots the stored histograms based on the plot options
     '''
     file_name = self.name
     if plot_options.has_key('output_file'):
         file_name = plot_options['output_file']
     output_format = plot_options['output_format'] 
     output_folder = plot_options['output_folder']
     
     plot_type = plot_options['plot_type']
     
     histogram_properties = Histogram_properties(plot_options)
     histogram_properties.name = file_name
     if plot_options.has_key('rebin') and plot_options['rebin'] > 1:
         for hist in self.histograms:
             hist.Rebin(plot_options['rebin'])
     
     colours = ['green', 'yellow', 'magenta', 'red', 'black']
     if plot_options.has_key('colours'):
         colours = plot_options['colours']
     
     if plot_type == 'shape_comparison':
         make_shape_comparison_plot( shapes = self.histograms,
                                    names = self.labels,
                                    colours = colours,
                                    histogram_properties = histogram_properties,
                                    fill_area = True,
                                    make_ratio = True,
                                    alpha = 0.5,
                                    save_folder = output_folder,
                                    save_as = output_format )
     elif plot_type == 'data_mc_comparison':
         make_data_mc_comparison_plot( histograms = self.histograms,
                                      histogram_lables = self.labels,
                                      histogram_colors = colours,
                                      histogram_properties = histogram_properties,
                                      data_index = 0,
                                      save_folder = output_folder,
                                      save_as = output_format,
                                      normalise = False,
                                      show_ratio = True,
                                      show_stat_errors_on_mc = True,
                                      draw_vertical_line = 0 )
     else:
         print 'Plot type "%s" not known, exiting.' % plot_type
         sys.exit()

def plot_fit_results(fit_results, initial_values, channel):
    global variable, output_folder

    title = electron_histogram_title if channel == "electron" else muon_histogram_title

    histogram_properties = Histogram_properties()
    histogram_properties.title = title

    histogram_properties.x_axis_title = variable + " [GeV]"
    histogram_properties.mc_error = 0.0
    histogram_properties.legend_location = "upper right"
    # we will need 4 histograms: TTJet, SingleTop, QCD, V+Jets
    for sample in ["TTJet", "SingleTop", "QCD", "V+Jets"]:
        histograms = {}
        # absolute eta measurement as baseline
        h_absolute_eta = None
        h_before = None
        histogram_properties.y_axis_title = "Fitted number of events for " + samples_latex[sample]

        for fit_var_input in fit_results.keys():
            latex_string = create_latex_string(fit_var_input)
            fit_data = fit_results[fit_var_input][sample]
            h = value_error_tuplelist_to_hist(fit_data, bin_edges[variable])
            if fit_var_input == "absolute_eta":
                h_absolute_eta = h
            elif fit_var_input == "before":
                h_before = h
            else:
                histograms[latex_string] = h
        graphs = spread_x(histograms.values(), bin_edges[variable])
        for key, graph in zip(histograms.keys(), graphs):
            histograms[key] = graph
        filename = sample.replace("+", "_") + "_fit_var_comparison_" + channel
        histogram_properties.name = filename
        histogram_properties.y_limits = 0, limit_range_y(h_absolute_eta)[1] * 1.3
        histogram_properties.x_limits = bin_edges[variable][0], bin_edges[variable][-1]

        h_initial_values = value_error_tuplelist_to_hist(initial_values[sample], bin_edges[variable])
        h_initial_values.Scale(closure_tests["simple"][sample])

        compare_measurements(
            models={
                fit_variables_latex["absolute_eta"]: h_absolute_eta,
                "initial values": h_initial_values,
                "before": h_before,
            },
            measurements=histograms,
            show_measurement_errors=True,
            histogram_properties=histogram_properties,
            save_folder=output_folder,
            save_as=["png", "pdf"],
        )

def plot_fit_results( fit_results, initial_values, channel ):
    global variable, output_folder
    
    title = electron_histogram_title if channel == 'electron' else muon_histogram_title
    
    
    histogram_properties = Histogram_properties()
    histogram_properties.title = title
    
    histogram_properties.x_axis_title = variable + ' [GeV]'
    histogram_properties.mc_error = 0.0
    histogram_properties.legend_location = 'upper right'
    # we will need 4 histograms: TTJet, SingleTop, QCD, V+Jets
    for sample in ['TTJet', 'SingleTop', 'QCD', 'V+Jets']:
        histograms = {}
        # absolute eta measurement as baseline
        h_absolute_eta = None
        h_before = None
        histogram_properties.y_axis_title = 'Fitted number of events for ' + samples_latex[sample]
        
        for fit_var_input in fit_results.keys():
            latex_string = create_latex_string( fit_var_input )
            fit_data = fit_results[fit_var_input][sample]
            h = value_error_tuplelist_to_hist( fit_data,
                                              bin_edges[variable] )
            if fit_var_input == 'absolute_eta':
                h_absolute_eta = h
            elif fit_var_input == 'before':
                h_before = h
            else:
                histograms[latex_string] = h
        graphs = spread_x( histograms.values(), bin_edges[variable] )
        for key, graph in zip( histograms.keys(), graphs ):
            histograms[key] = graph
        filename = sample.replace( '+', '_' ) + '_fit_var_comparison_' + channel
        histogram_properties.name = filename
        histogram_properties.y_limits = 0, limit_range_y( h_absolute_eta )[1] * 1.3
        histogram_properties.x_limits = bin_edges[variable][0], bin_edges[variable][-1]
        
        h_initial_values = value_error_tuplelist_to_hist( initial_values[sample],
                                                         bin_edges[variable] )
        h_initial_values.Scale(closure_tests['simple'][sample])
        
        compare_measurements( models = {fit_variables_latex['absolute_eta']:h_absolute_eta,
                                        'initial values' : h_initial_values,
                                        'before': h_before},
                             measurements = histograms,
                             show_measurement_errors = True,
                             histogram_properties = histogram_properties,
                             save_folder = output_folder,
                             save_as = ['png', 'pdf'] )

def plot_fit_results( histograms, category, channel ):
    global variable, b_tag_bin, output_folder
    from tools.plotting import Histogram_properties, make_data_mc_comparison_plot
    fit_variables = histograms.keys()
    for variable_bin in variable_bins_ROOT[variable]:
        path = output_folder + str( measurement_config.centre_of_mass_energy ) + 'TeV/' + variable + '/' + category + '/fit_results/'
        make_folder_if_not_exists( path )
        for fit_variable in fit_variables:
            plotname = channel + '_' + fit_variable + '_bin_' + variable_bin
            # check if template plots exist already
            for output_format in output_formats:
                if os.path.isfile( plotname + '.' + output_format ):
                    continue
                
            # plot with matplotlib
            h_data = histograms[fit_variable][variable_bin]['data']
            h_signal = histograms[fit_variable][variable_bin]['signal']
            h_background = histograms[fit_variable][variable_bin]['background']
            
            histogram_properties = Histogram_properties()
            histogram_properties.name = plotname
            histogram_properties.x_axis_title = fit_variables_latex[fit_variable]
            histogram_properties.y_axis_title = 'Events/(%s)' % get_unit_string(fit_variable)
            label, _ = get_cms_labels( channel )
            histogram_properties.title = label
            histogram_properties.x_limits = measurement_config.fit_boundaries[fit_variable]
            
            make_data_mc_comparison_plot( [h_data, h_background, h_signal],
                                         ['data', 'background', 'signal'],
                                         ['black', 'green', 'red'], histogram_properties,
                                         save_folder = path, save_as = output_formats )    

def plot_bias(unfolded_and_truths, variable, channel, come, method):
    hp = Histogram_properties()
    hp.name = 'Bias_{channel}_{variable}_at_{come}TeV'.format(
        channel=channel,
        variable=variable,
        come=come,
    )
    v_latex = latex_labels.variables_latex[variable]
    unit = ''
    if variable in ['HT', 'ST', 'MET', 'WPT', 'lepton_pt']:
        unit = ' [GeV]'
    hp.x_axis_title = v_latex + unit
    # plt.ylabel( r, CMS.y_axis_title )
    hp.y_axis_title = 'Unfolded / Truth'
    hp.y_limits = [0.92, 1.08]
    hp.title = 'Bias for {variable}'.format(variable=v_latex)

    output_folder = 'plots/unfolding/bias_test/'

    measurements = { 'Central' : unfolded_and_truths['Central']['bias'] }

    models = {}
    for sample in unfolded_and_truths:
        if sample == 'Central' : continue
        models[sample] = unfolded_and_truths[sample]['bias']


    compare_measurements(
                         models = models,
                         measurements = measurements,
                         show_measurement_errors=True,
                         histogram_properties=hp,
                         save_folder=output_folder,
                         save_as=['pdf'],
                         match_models_to_measurements = True)

def compare( central_mc, expected_result = None, measured_result = None, results = {}, variable = 'MET',
             channel = 'electron', bin_edges = [] ):
    global input_file, plot_location, ttbar_xsection, luminosity, centre_of_mass, method, test, log_plots

    channel_label = ''
    if channel == 'electron':
        channel_label = 'e+jets, $\geq$4 jets'
    elif channel == 'muon':
        channel_label = '$\mu$+jets, $\geq$4 jets'
    else:
        channel_label = '$e, \mu$ + jets combined, $\geq$4 jets'

    if test == 'data':
        title_template = 'CMS Preliminary, $\mathcal{L} = %.1f$ fb$^{-1}$  at $\sqrt{s}$ = %d TeV \n %s'
        title = title_template % ( luminosity / 1000., centre_of_mass, channel_label )
    else:
        title_template = 'CMS Simulation at $\sqrt{s}$ = %d TeV \n %s'
        title = title_template % ( centre_of_mass, channel_label )

    models = {latex_labels.measurements_latex['MADGRAPH'] : central_mc}
    if expected_result and test == 'data':
        models.update({'fitted data' : expected_result})
        # scale central MC to lumi
        nEvents = input_file.EventFilter.EventCounter.GetBinContent( 1 )  # number of processed events 
        lumiweight = ttbar_xsection * luminosity / nEvents
        central_mc.Scale( lumiweight )
    elif expected_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():
        measurements['k = ' + str( key )] = value
    
    # get some spread in x    
    graphs = spread_x( measurements.values(), bin_edges )
    for key, graph in zip( measurements.keys(), graphs ):
        measurements[key] = graph

    histogram_properties = Histogram_properties()
    histogram_properties.name = channel + '_' + variable + '_' + method + '_' + test
    histogram_properties.title = title + ', ' + latex_labels.b_tag_bins_latex['2orMoreBtags']
    histogram_properties.x_axis_title = '$' + latex_labels.variables_latex[variable] + '$'
    histogram_properties.y_axis_title = r'Events'
#     histogram_properties.y_limits = [0, 0.03]
    histogram_properties.x_limits = [bin_edges[0], bin_edges[-1]]

    if log_plots:
        histogram_properties.set_log_y = True
        histogram_properties.name += '_log'

    compare_measurements( models, measurements, show_measurement_errors = True,
                          histogram_properties = histogram_properties,
                          save_folder = plot_location, save_as = ['pdf'] )

def plot_results ( results ):
    '''
    Takes results fo the form:
        {centre-of-mass-energy: {
            channel : {
                variable : {
                    fit_variable : {
                        test : { sample : []},
                        }
                    }
                }
            }
        }
    '''
    global options
    output_base = 'plots/fit_checks/chi2'
    for COMEnergy in results.keys():
        tmp_result_1 = results[COMEnergy]
        for channel in tmp_result_1.keys():
            tmp_result_2 = tmp_result_1[channel]
            for variable in tmp_result_2.keys():
                tmp_result_3 = tmp_result_2[variable]
                for fit_variable in tmp_result_3.keys():
                    tmp_result_4 = tmp_result_3[fit_variable]
                    # histograms should be {sample: {test : histogram}}
                    histograms = {}
                    for test, chi2 in tmp_result_4.iteritems():
                        for sample in chi2.keys():
                            if not histograms.has_key(sample):
                                histograms[sample] = {}
                            # reverse order of test and sample
                            histograms[sample][test] = value_tuplelist_to_hist(chi2[sample], bin_edges[variable])
                    for sample in histograms.keys():
                        hist_properties = Histogram_properties()
                        hist_properties.name = sample.replace('+', '') + '_chi2'
                        hist_properties.title = '$\\chi^2$ distribution for fit output (' + sample + ')'
                        hist_properties.x_axis_title = '$' + latex_labels.variables_latex[variable] + '$ [TeV]'
                        hist_properties.y_axis_title = '$\chi^2 = \\left({N_{fit}} - N_{{exp}}\\right)^2$'
                        hist_properties.set_log_y = True
                        hist_properties.y_limits = (1e-20, 1e20)
                        path = output_base + '/' + COMEnergy + 'TeV/' + channel + '/' + variable + '/' + fit_variable + '/'
                        if options.test:
                            path = output_base + '/test/'
                        
                        measurements = {}
                        for test, histogram in histograms[sample].iteritems():
                            measurements[test.replace('_',' ')] = histogram
                        compare_measurements({}, 
                                             measurements, 
                                             show_measurement_errors = False, 
                                             histogram_properties = hist_properties, 
                                             save_folder = path, 
                                             save_as = ['pdf'])

def compare_vjets_btag_regions( variable = 'MET', met_type = 'patType1CorrectedPFMet',
                                title = 'Untitled', channel = 'electron' ):
    ''' Compares the V+Jets template in different b-tag bins'''
    global fit_variable_properties, b_tag_bin, save_as, b_tag_bin_ctl
    b_tag_bin_ctl = '0orMoreBtag'
    variable_bins = variable_bins_ROOT[variable]
    histogram_template = get_histogram_template( variable )
    
    for fit_variable in electron_fit_variables:
        if '_bl' in fit_variable:
                b_tag_bin_ctl = '1orMoreBtag'
        else:
            b_tag_bin_ctl = '0orMoreBtag'
        save_path = 'plots/%dTeV/fit_variables/%s/%s/' % ( measurement_config.centre_of_mass_energy, variable, fit_variable )
        make_folder_if_not_exists( save_path + '/vjets/' )
        histogram_properties = Histogram_properties()
        histogram_properties.x_axis_title = fit_variable_properties[fit_variable]['x-title']
        histogram_properties.y_axis_title = fit_variable_properties[fit_variable]['y-title']
        histogram_properties.y_axis_title = histogram_properties.y_axis_title.replace( 'Events', 'a.u.' )
        histogram_properties.x_limits = [fit_variable_properties[fit_variable]['min'], fit_variable_properties[fit_variable]['max']]
        histogram_properties.title = title
        histogram_properties.additional_text = channel_latex[channel] + ', ' + b_tag_bins_latex[b_tag_bin_ctl]
        histogram_properties.y_max_scale = 1.5
        for bin_range in variable_bins:
            params = {'met_type': met_type, 'bin_range':bin_range, 'fit_variable':fit_variable, 'b_tag_bin':b_tag_bin, 'variable':variable}
            fit_variable_distribution = histogram_template % params
            fit_variable_distribution_ctl = fit_variable_distribution.replace( b_tag_bin, b_tag_bin_ctl )
            # format: histograms['data'][qcd_fit_variable_distribution]
            histograms = get_histograms_from_files( [fit_variable_distribution, fit_variable_distribution_ctl], {'V+Jets' : histogram_files['V+Jets']} )
            prepare_histograms( histograms, rebin = fit_variable_properties[fit_variable]['rebin'], scale_factor = measurement_config.luminosity_scale )
            histogram_properties.name = variable + '_' + bin_range + '_' + fit_variable + '_' + b_tag_bin_ctl + '_VJets_template_comparison'
            histograms['V+Jets'][fit_variable_distribution].Scale( 1 / histograms['V+Jets'][fit_variable_distribution].Integral() )
            histograms['V+Jets'][fit_variable_distribution_ctl].Scale( 1 / histograms['V+Jets'][fit_variable_distribution_ctl].Integral() )
            compare_measurements( models = {'no b-tag' : histograms['V+Jets'][fit_variable_distribution_ctl]},
                             measurements = {'$>=$ 2 b-tags': histograms['V+Jets'][fit_variable_distribution]},
                             show_measurement_errors = True,
                             histogram_properties = histogram_properties,
                             save_folder = save_path + '/vjets/',
                             save_as = save_as )

def plot_fit_results(fit_results, centre_of_mass, channel, variable, k_value,
                     tau_value, output_folder, output_formats, bin_edges):
    h_mean = Hist(bin_edges, type='D')
    h_sigma = Hist(bin_edges, type='D')
    n_bins = h_mean.nbins()
    assert len(fit_results) == n_bins
    for i, fr in enumerate(fit_results):
        h_mean.SetBinContent(i + 1, fr.mean)
        h_mean.SetBinError(i + 1, fr.meanError)
        h_sigma.SetBinContent(i + 1, fr.sigma)
        h_sigma.SetBinError(i + 1, fr.sigmaError)
    histogram_properties = Histogram_properties()
    name_mpt = 'pull_distribution_mean_and_sigma_{0}_{1}_{2}TeV'
    histogram_properties.name = name_mpt.format(
        variable,
        channel,
        centre_of_mass
    )
    histogram_properties.y_axis_title = r'$\mu_{\text{pull}}$ ($\sigma_{\text{pull}}$)'
    histogram_properties.x_axis_title = latex_labels.variables_latex[variable]
    value = get_value_title(k_value, tau_value)
    title = 'pull distribution mean \& sigma for {0}'.format(value)
    histogram_properties.title = title
    histogram_properties.y_limits = [-0.5, 2]
    histogram_properties.xerr = True

    compare_measurements(
        models={
            'ideal $\mu$': make_line_hist(bin_edges, 0),
            'ideal $\sigma$': make_line_hist(bin_edges, 1)
        },
        measurements={
            r'$\mu_{\text{pull}}$': h_mean,
            r'$\sigma_{\text{pull}}$': h_sigma
        },
        show_measurement_errors=True,
        histogram_properties=histogram_properties,
        save_folder=output_folder,
        save_as=output_formats)

def draw_regularisation_histograms( h_truth, h_measured, h_response, h_fakes = None, h_data = None ):
    global method, variable, output_folder, output_formats, test
    k_max = h_measured.nbins()
    unfolding = Unfolding( h_truth,
                           h_measured,
                           h_response,
                           h_fakes,
                           method = method,
                           k_value = k_max,
                           error_treatment = 4,
                           verbose = 1 )
    
    RMSerror, MeanResiduals, RMSresiduals, Chi2 = unfolding.test_regularisation ( h_data, k_max )

    histogram_properties = Histogram_properties()
    histogram_properties.name = 'chi2_%s_channel_%s' % ( channel, variable )
    histogram_properties.title = '$\chi^2$ for $%s$ in %s channel, %s test' % ( variables_latex[variable], channel, test )
    histogram_properties.x_axis_title = '$i$'
    histogram_properties.y_axis_title = '$\chi^2$'
    histogram_properties.set_log_y = True
    make_plot(Chi2, 'chi2', histogram_properties, output_folder, output_formats, draw_errorbar = True, draw_legend = False)

    histogram_properties = Histogram_properties()
    histogram_properties.name = 'RMS_error_%s_channel_%s' % ( channel, variable )
    histogram_properties.title = 'Mean error for $%s$ in %s channel, %s test' % ( variables_latex[variable], channel, test )
    histogram_properties.x_axis_title = '$i$'
    histogram_properties.y_axis_title = 'Mean error'
    make_plot(RMSerror, 'RMS', histogram_properties, output_folder, output_formats, draw_errorbar = True, draw_legend = False)

    histogram_properties = Histogram_properties()
    histogram_properties.name = 'RMS_residuals_%s_channel_%s' % ( channel, variable )
    histogram_properties.title = 'RMS of residuals for $%s$ in %s channel, %s test' % ( variables_latex[variable], channel, test )
    histogram_properties.x_axis_title = '$i$'
    histogram_properties.y_axis_title = 'RMS of residuals'
    if test == 'closure':
        histogram_properties.set_log_y = True
    make_plot(RMSresiduals, 'RMSresiduals', histogram_properties, output_folder, output_formats, draw_errorbar = True, draw_legend = False)

    histogram_properties = Histogram_properties()
    histogram_properties.name = 'mean_residuals_%s_channel_%s' % ( channel, variable )
    histogram_properties.title = 'Mean of residuals for $%s$ in %s channel, %s test' % ( variables_latex[variable], channel, test )
    histogram_properties.x_axis_title = '$i$'
    histogram_properties.y_axis_title = 'Mean of residuals'
    make_plot(MeanResiduals, 'MeanRes', histogram_properties, output_folder, output_formats, draw_errorbar = True, draw_legend = False)

def do_shape_check(channel, control_region_1, control_region_2, variable, normalisation, title, x_title, y_title, x_limits, y_limits,
                   name_region_1='conversions' , name_region_2='non-isolated electrons', name_region_3='fit results', rebin=1):
    global b_tag_bin
    # QCD shape comparison
    if channel == 'electron':
        histograms = get_histograms_from_files([control_region_1, control_region_2], histogram_files)
        
        region_1 = histograms[channel][control_region_1].Clone() - histograms['TTJet'][control_region_1].Clone() - histograms['V+Jets'][control_region_1].Clone() - histograms['SingleTop'][control_region_1].Clone()
        region_2 = histograms[channel][control_region_2].Clone() - histograms['TTJet'][control_region_2].Clone() - histograms['V+Jets'][control_region_2].Clone() - histograms['SingleTop'][control_region_2].Clone()
        
        region_1.Rebin(rebin)
        region_2.Rebin(rebin)
        
        histogram_properties = Histogram_properties()
        histogram_properties.name = 'QCD_control_region_comparison_' + channel + '_' + variable + '_' + b_tag_bin
        histogram_properties.title = title + ', ' + b_tag_bins_latex[b_tag_bin]
        histogram_properties.x_axis_title = x_title
        histogram_properties.y_axis_title = 'arbitrary units/(0.1)'
        histogram_properties.x_limits = x_limits
        histogram_properties.y_limits = y_limits[0]
        histogram_properties.mc_error = 0.0
        histogram_properties.legend_location = 'upper right'
        make_control_region_comparison(region_1, region_2,
                                       name_region_1=name_region_1, name_region_2=name_region_2,
                                       histogram_properties=histogram_properties, save_folder=output_folder)
        
        # QCD shape comparison to fit results
        histograms = get_histograms_from_files([control_region_1], histogram_files)
        
        region_1_tmp = histograms[channel][control_region_1].Clone() - histograms['TTJet'][control_region_1].Clone() - histograms['V+Jets'][control_region_1].Clone() - histograms['SingleTop'][control_region_1].Clone()
        region_1 = rebin_asymmetric(region_1_tmp, bin_edges[variable])
        
        fit_results_QCD = normalisation[variable]['QCD']
        region_2 = value_error_tuplelist_to_hist(fit_results_QCD, bin_edges_vis[variable])
        
        histogram_properties = Histogram_properties()
        histogram_properties.name = 'QCD_control_region_comparison_' + channel + '_' + variable + '_fits_with_conversions_' + b_tag_bin
        histogram_properties.title = title + ', ' + b_tag_bins_latex[b_tag_bin]
        histogram_properties.x_axis_title = x_title
        histogram_properties.y_axis_title = 'arbitrary units/(0.1)'
        histogram_properties.x_limits = x_limits
        histogram_properties.y_limits = y_limits[1]
        histogram_properties.mc_error = 0.0
        histogram_properties.legend_location = 'upper right'
        make_control_region_comparison(region_1, region_2,
                                       name_region_1=name_region_1, name_region_2=name_region_3,
                                       histogram_properties=histogram_properties, save_folder=output_folder)
    
    histograms = get_histograms_from_files([control_region_2], histogram_files)
    
    region_1_tmp = histograms[channel][control_region_2].Clone() - histograms['TTJet'][control_region_2].Clone() - histograms['V+Jets'][control_region_2].Clone() - histograms['SingleTop'][control_region_2].Clone()
    region_1 = rebin_asymmetric(region_1_tmp, bin_edges_vis[variable])    
    
    fit_results_QCD = normalisation[variable]['QCD']
    region_2 = value_error_tuplelist_to_hist(fit_results_QCD, bin_edges[variable])
    
    histogram_properties = Histogram_properties()
    histogram_properties.name = 'QCD_control_region_comparison_' + channel + '_' + variable + '_fits_with_noniso_' + b_tag_bin
    histogram_properties.title = title + ', ' + b_tag_bins_latex[b_tag_bin]
    histogram_properties.x_axis_title = x_title
    histogram_properties.y_axis_title = 'arbitrary units/(0.1)'
    histogram_properties.x_limits = x_limits
    histogram_properties.y_limits = y_limits[1]
    histogram_properties.mc_error = 0.0
    histogram_properties.legend_location = 'upper right'
    make_control_region_comparison(region_1, region_2,
                                   name_region_1=name_region_2, name_region_2=name_region_3,
                                   histogram_properties=histogram_properties, save_folder=output_folder)