def plot_central_and_systematics( channel, systematics, exclude = [], suffix = 'altogether' ):
    global variable, k_values, b_tag_bin, met_type

    plt.figure( figsize = ( 16, 16 ), dpi = 200, facecolor = 'white' )
    axes = plt.axes()
    axes.minorticks_on()
    
    hist_data_central = read_xsection_measurement_results( 'central', channel )[0]['unfolded_with_systematics']
    hist_data_central.markersize = 2  # points. Imagine, tangible units!
    hist_data_central.marker = 'o'
    
    
    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 )

    rplt.errorbar( hist_data_central, axes = axes, label = 'data', xerr = True )

    for systematic in sorted( systematics ):
        if systematic in exclude or systematic == 'central':
            continue

        hist_data_systematic = read_xsection_measurement_results( systematic, channel )[0]['unfolded']
        hist_data_systematic.markersize = 2
        hist_data_systematic.marker = 'o'
        colour_number = systematics.index( systematic ) + 2
        if colour_number == 10:
            colour_number = 42
        hist_data_systematic.SetMarkerColor( colour_number )
        if 'PDF' in systematic:
            rplt.errorbar( hist_data_systematic, axes = axes, label = systematic.replace( 'Weights_', ' ' ), xerr = None )
        elif met_type in systematic:
            rplt.errorbar( hist_data_systematic, axes = axes, label = met_systematics_latex[systematic.replace( met_type, '' )], xerr = None )
        else:
            rplt.errorbar( hist_data_systematic, axes = axes, label = measurements_latex[systematic], xerr = None )
            
    plt.legend( numpoints = 1, loc = 'center right', prop = {'size':25}, ncol = 2 )
    label, channel_label = get_cms_labels( channel )
    plt.title( label, CMS.title )
    # CMS text
    # note: fontweight/weight does not change anything as we use Latex text!!!
    plt.text(0.95, 0.95, r"\textbf{CMS}", transform=axes.transAxes, fontsize=42,
        verticalalignment='top',horizontalalignment='right')
    # channel text
    axes.text(0.95, 0.90, r"\emph{%s}" %channel_label, transform=axes.transAxes, fontsize=40,
        verticalalignment='top',horizontalalignment='right')
    plt.tight_layout()

    
    path = output_folder + str( measurement_config.centre_of_mass_energy ) + 'TeV/' + variable
    make_folder_if_not_exists( path )
    for output_format in output_formats:
        filename = path + '/normalised_xsection_' + channel + '_' + suffix + '_kv' + str( k_values[channel] ) + '.' + output_format
        if channel == 'combined':
            filename = filename.replace( '_kv' + str( k_values[channel] ), '' )
        plt.savefig( filename ) 

    plt.close()
    gc.collect()

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_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 make_template_plots(histograms, category, channel):
    global variable, output_folder
    
    for variable_bin in variable_bins_ROOT[variable]:
        path = output_folder + str(measurement_config.centre_of_mass) + 'TeV/' + variable + '/' + category + '/fit_templates/'
        make_folder_if_not_exists(path)
        plotname = path + channel + '_templates_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_signal = histograms[variable_bin]['signal']
        h_VJets = histograms[variable_bin]['V+Jets']
        h_QCD = histograms[variable_bin]['QCD']
        
        h_signal.linecolor = 'red'
        h_VJets.linecolor = 'green'
        h_QCD.linecolor = 'gray'
        h_VJets.linestyle = 'dashed'
        h_QCD.linestyle = 'dotted'# currently not working
        #bug report: http://trac.sagemath.org/sage_trac/ticket/13834
        
        h_signal.linewidth = 5
        h_VJets.linewidth = 5
        h_QCD.linewidth = 5
    
        plt.figure(figsize=(16, 16), dpi=200, facecolor='white')
        axes = plt.axes()
        axes.minorticks_on()
        
        plt.xlabel(r'lepton $|\eta|$', CMS.x_axis_title)
        plt.ylabel('normalised to unit area/0.2', CMS.y_axis_title)
        plt.tick_params(**CMS.axis_label_major)
        plt.tick_params(**CMS.axis_label_minor)

        rplt.hist(h_signal, axes=axes, label='signal')
        if (h_VJets.Integral() != 0):
            rplt.hist(h_VJets, axes=axes, label='V+Jets')
        else:
            print "WARNING: in %s bin %s, %s category, %s channel, V+Jets template is empty: not plotting." % (variable, variable_bin, category, channel)
        if (h_QCD.Integral() != 0):
            rplt.hist(h_QCD, axes=axes, label='QCD')
        else:
            print "WARNING: in %s bin %s, %s category, %s channel, QCD template is empty: not plotting." % (variable, variable_bin, category, channel)
        axes.set_ylim([0, 0.2])
        
        plt.legend(numpoints=1, loc='upper right', prop=CMS.legend_properties)
        plt.title(get_cms_labels(channel), CMS.title)
        plt.tight_layout()
    
        for output_format in output_formats:
            plt.savefig(plotname + '.' + output_format)
        
        plt.close()
        gc.collect()

    def submit(self):
        '''
            Submits all registered jobs to the local HTCondor scheduler using
            a job template (DailyPythonScripts/condor/job_template) description
            file and the 'condor_submit' command
        '''
        today = time.strftime("%d-%m-%Y")
        job_folder = 'jobs/{0}/'.format(today)
        make_folder_if_not_exists(job_folder)
        make_folder_if_not_exists(job_folder + 'logs')
        # construct jobs
        self._construct_jobs()
        # convert each job into a pickle file
        # construct a class ad for each job
        with open('condor/job_template', 'r') as template:
            job_template = template.read()
        condor_jobs = []

        # prepare DPS for submission
        self._dps_tar_directory_on_hdfs = '/TopQuarkGroup/condor_dps/{you}/{now}/'.format(
                                                                                            you = getpass.getuser(), 
                                                                                            now = time.strftime('%d_%m_%Y_%H_%M') 
                                                                                            )

        for i, job in enumerate(self.prepared_jobs):
            job_file = job_folder + 'job_{0}.pkl'.format(i)
            job_desc_file = job_folder + 'job_{0}.dsc'.format(i)
            job_description = job_template.replace('%pkl_file%', job_file)
            job_description = job_description.replace('%dir_of_dps_on_hdfs%',
                                                      self._dps_tar_directory_on_hdfs)
            job_description = job_description.replace('%total_memory%',
                                                      str(self.request_memory))
            job_description = job_description.replace('%n_jobs_to_run%',
                                                      str(self.n_jobs_to_run))
            job_description = job_description.replace('%n_jobs_to_split%',
                                                      str(self.n_jobs_to_split))
            input_files = []
            if hasattr(job, 'additional_input_files'):
                input_files.extend(job.additional_input_files)
            input_files_str = ','.join(input_files)
            job_description = job_description.replace('%input_files%',
                                                      input_files_str)
            job_description = job_description.replace('%today%', today)

            with open(job_file, 'w+') as jf:
                pickle.dump(job, jf)
            with open(job_desc_file, 'w+') as jdf:
                jdf.write(job_description)

            condor_jobs.append(job_desc_file)

        prepare_process = subprocess.Popen(['./condor/prepare_dps.sh',self._dps_tar_directory_on_hdfs])
        prepare_process.communicate()
        
        # # submit jobs
        for j in condor_jobs:
            p = subprocess.Popen(['condor_submit', j])
            p.communicate()  # wait until command completed

def make_plots_matplotlib(histograms, category, output_folder, histname):
    global variable, variables_latex_matplotlib, measurements_latex_matplotlib, k_value
    
    channel = 'electron'
    if 'electron' in histname:
        channel = 'electron'
    elif 'muon' in histname:
        channel = 'muon'
    else:
        channel = 'combined'
        
    # plot with matplotlib
    hist_data = histograms['unfolded']
    hist_measured = histograms['measured']
    
    hist_data.markersize = 2
    hist_measured.markersize = 2
    hist_data.marker = 'o'
    hist_measured.marker = 'o'
    hist_measured.color = 'red'

    plt.figure(figsize=(14, 10), dpi=200, facecolor='white')
    axes = plt.axes()
    axes.minorticks_on()
    
    plt.xlabel('$%s$ [GeV]' % variables_latex_matplotlib[variable], CMS.x_axis_title)
    plt.ylabel(r'$\frac{1}{\sigma} \times \frac{d\sigma}{d' + variables_latex_matplotlib[variable] + '} \left[\mathrm{GeV}^{-1}\\right]$', CMS.y_axis_title)
    plt.tick_params(**CMS.axis_label_major)
    plt.tick_params(**CMS.axis_label_minor)

    rplt.errorbar(hist_data, axes=axes, label='unfolded', xerr=False)
    rplt.errorbar(hist_measured, axes=axes, label='measured', xerr=False)
    
    for key, hist in histograms.iteritems():
        if not 'unfolded' in key and not 'measured' in key:
            hist.linestyle = 'dashed'
            hist.linewidth = 2
#            hist.SetLineStyle(7)
#            hist.SetLineWidth(2)
            #setting colours
            if 'POWHEG' in key or 'matchingdown' in key:
                hist.SetLineColor(kBlue)
            elif 'MADGRAPH' in key or 'matchingup' in key:
                hist.SetLineColor(kRed + 1)
            elif 'MCATNLO'  in key or 'scaleup' in key:
                hist.SetLineColor(kMagenta + 3)
            elif 'scaledown' in key:
                hist.SetLineColor(kGreen)
            rplt.hist(hist, axes=axes, label=measurements_latex_matplotlib[key])
    
    plt.legend(numpoints=1, loc='upper right', prop=CMS.legend_properties)
    plt.title(get_cms_labels_matplotlib(channel), CMS.title)
    plt.tight_layout()

    path = output_folder + str(measurement_config.centre_of_mass) + 'TeV/' + variable + '/' + category
    make_folder_if_not_exists(path)
    for output_format in output_formats:
        plt.savefig(path + '/' + histname + '_kv' + str(k_value) + '.' + output_format)

def main():
    '''
        Main function for this script
    '''
    set_root_defaults(msg_ignore_level=3001)

    parser = OptionParser()
    parser.add_option("-o", "--output",
                      dest="output_folder", default='data/pull_data/',
                      help="output folder for pull data files")
    parser.add_option("-n", "--n_input_mc", type=int,
                      dest="n_input_mc", default=100,
                      help="number of toy MC used for the tests")
    parser.add_option("-k", "--k_value", type=int,
                      dest="k_value", default=3,
                      help="k-value for SVD unfolding")
    parser.add_option("--tau", type='float',
                      dest="tau_value", default=-1.,
                      help="tau-value for SVD unfolding")
    parser.add_option("-m", "--method", type='string',
                      dest="method", default='RooUnfoldSvd',
                      help="unfolding method")
    parser.add_option("-f", "--file", type='string',
                      dest="file", default='data/toy_mc/unfolding_toy_mc.root',
                      help="file with toy MC")
    parser.add_option("-v", "--variable", dest="variable", default='MET',
                      help="set the variable to analyse (MET, HT, ST, MT, WPT)")
    parser.add_option("-s", "--centre-of-mass-energy", dest="CoM", default=13,
                      help='''set the centre of mass energy for analysis.
                      Default = 8 [TeV]''', type=int)
    parser.add_option("-c", "--channel", type='string',
                      dest="channel", default='combined',
                      help="channel to be analysed: electron|muon|combined")

    parser.add_option("--offset_toy_mc", type=int,
                      dest="offset_toy_mc", default=0,
                      help="offset of the toy MC used to response matrix")
    parser.add_option("--offset_toy_data", type=int,
                      dest="offset_toy_data", default=0,
                      help="offset of the toy MC used as data for unfolding")
    (options, _) = parser.parse_args()

    centre_of_mass = options.CoM
    make_folder_if_not_exists(options.output_folder)

    # set the number of toy MC for error calculation
    k_value = options.k_value
    tau_value = options.tau_value
    use_n_toy = options.n_input_mc
    offset_toy_mc = options.offset_toy_mc
    offset_toy_data = options.offset_toy_data
    method = options.method
    variable = options.variable

    create_unfolding_pull_data(options.file, method, options.channel,
                               centre_of_mass, variable, use_n_toy, use_n_toy,
                               options.output_folder, offset_toy_mc,
                               offset_toy_data, k_value, tau_value)

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 save(self):
     make_folder_if_not_exists(self._path)
     for f in self.__properties.formats:
         file_name = '{path}{name}.{format}'
         file_name = file_name.format(
                         path = self._path,
                         name = self.__properties.name,
                         format = f,
                         )
         plt.savefig(file_name)

def plot_central_and_systematics(channel, systematics, exclude=[], suffix='altogether'):
    global variable, variables_latex, k_value, b_tag_bin, maximum, ttbar_generator_systematics

    canvas = Canvas(width=700, height=500)
    canvas.SetLeftMargin(0.15)
    canvas.SetBottomMargin(0.15)
    canvas.SetTopMargin(0.05)
    canvas.SetRightMargin(0.05)
    legend = plotting.create_legend(x0=0.6, y1=0.5)
    
    hist_data_central = read_xsection_measurement_results('central', channel)[0]['unfolded']
    
    hist_data_central.GetXaxis().SetTitle(variables_latex[variable] + ' [GeV]')
    hist_data_central.GetYaxis().SetTitle('#frac{1}{#sigma} #frac{d#sigma}{d' + variables_latex[variable] + '} [GeV^{-1}]')
    hist_data_central.GetXaxis().SetTitleSize(0.05)
    hist_data_central.GetYaxis().SetTitleSize(0.05)
    hist_data_central.SetMinimum(0)
    hist_data_central.SetMaximum(maximum[variable])
    hist_data_central.SetMarkerSize(1)
    hist_data_central.SetMarkerStyle(20)

    gStyle.SetEndErrorSize(20)
    hist_data_central.Draw('P')
    legend.AddEntry(hist_data_central, 'measured (unfolded)', 'P')
    
    for systematic in systematics:
        if systematic in exclude or systematic == 'central':
            continue

        hist_data_systematic = read_xsection_measurement_results(systematic, channel)[0]['unfolded']
        hist_data_systematic.SetMarkerSize(0.5)
        hist_data_systematic.SetMarkerStyle(20)
        colour_number = systematics.index(systematic)+1
        if colour_number == 10:
            colour_number = 42
        hist_data_systematic.SetMarkerColor(colour_number)
        hist_data_systematic.Draw('same P')
        legend.AddEntry(hist_data_systematic, systematic, 'P')
    
    legend.Draw()
    
    cms_label, channel_label = get_cms_labels(channel)
    cms_label.Draw()
    
    channel_label.Draw()
    
    canvas.Modified()
    canvas.Update()
    
    path = output_folder + str(measurement_config.centre_of_mass) + 'TeV/' + variable
    make_folder_if_not_exists(path)
    
    for output_format in output_formats:
        canvas.SaveAs(path + '/normalised_xsection_' + channel + '_' + suffix + '_kv' + str(k_value) + '.' + output_format)

def check_save_folder(save_folder):
    '''
        Checks and fixes (if necessary) the save folder
    '''
    # save_folder should end with an '/'
    if not save_folder.endswith('/'):
        save_folder += '/'
    # save_folder should exist
    make_folder_if_not_exists(save_folder)
    
    return save_folder

def unfold_results( results, category, channel, k_value, h_truth, h_measured, h_response, h_fakes, method ):
    global variable, path_to_JSON, options
    h_data = value_error_tuplelist_to_hist( results, bin_edges[variable] )
    unfolding = Unfolding( h_truth, h_measured, h_response, h_fakes, method = method, k_value = k_value )
    
    # turning off the unfolding errors for systematic samples
    if not category == 'central':
        unfoldCfg.Hreco = 0
    else:
        unfoldCfg.Hreco = options.Hreco
        
    h_unfolded_data = unfolding.unfold( h_data )
    
    if options.write_unfolding_objects:
        # export the D and SV distributions
        SVD_path = path_to_JSON + '/unfolding_objects/' + channel + '/kv_' + str( k_value ) + '/'
        make_folder_if_not_exists( SVD_path )
        if method == 'TSVDUnfold':
            SVDdist = File( SVD_path + method + '_SVDdistributions_' + category + '.root', 'recreate' )
            directory = SVDdist.mkdir( 'SVDdist' )
            directory.cd()
            unfolding.unfoldObject.GetD().Write()
            unfolding.unfoldObject.GetSV().Write()
            #    unfolding.unfoldObject.GetUnfoldCovMatrix(data_covariance_matrix(h_data), unfoldCfg.SVD_n_toy).Write()
            SVDdist.Close()
        else:
            SVDdist = File( SVD_path + method + '_SVDdistributions_Hreco' + str( unfoldCfg.Hreco ) + '_' + category + '.root', 'recreate' )
            directory = SVDdist.mkdir( 'SVDdist' )
            directory.cd()
            unfolding.unfoldObject.Impl().GetD().Write()
            unfolding.unfoldObject.Impl().GetSV().Write()
            h_truth.Write()
            h_measured.Write()
            h_response.Write()
            #    unfolding.unfoldObject.Impl().GetUnfoldCovMatrix(data_covariance_matrix(h_data), unfoldCfg.SVD_n_toy).Write()
            SVDdist.Close()
    
        # export the whole unfolding object if it doesn't exist
        if method == 'TSVDUnfold':
            unfolding_object_file_name = SVD_path + method + '_unfoldingObject_' + category + '.root'
        else:
            unfolding_object_file_name = SVD_path + method + '_unfoldingObject_Hreco' + str( unfoldCfg.Hreco ) + '_' + category + '.root'
        if not os.path.isfile( unfolding_object_file_name ):
            unfoldingObjectFile = File( unfolding_object_file_name, 'recreate' )
            directory = unfoldingObjectFile.mkdir( 'unfoldingObject' )
            directory.cd()
            if method == 'TSVDUnfold':
                unfolding.unfoldObject.Write()
            else:
                unfolding.unfoldObject.Impl().Write()
            unfoldingObjectFile.Close()
    
    del unfolding
    return hist_to_value_error_tuplelist( h_unfolded_data )

def main():
    '''
        Main function for this script
    '''
    set_root_defaults(msg_ignore_level=3001)

    parser = OptionParser()
    parser.add_option("-o", "--output",
                      dest="output_folder", default='data/pull_data/',
                      help="output folder for pull data files")
    parser.add_option("-n", "--n_input_mc", type=int,
                      dest="n_input_mc", default=100,
                      help="number of toy MC used for the tests")
    parser.add_option("--tau", type='float',
                      dest="tau_value", default=-1.,
                      help="tau-value for SVD unfolding")
    parser.add_option("-m", "--method", type='string',
                      dest="method", default='TUnfold',
                      help="unfolding method")
    parser.add_option("-f", "--file", type='string',
                      dest="file", default='data/toy_mc/unfolding_toy_mc.root',
                      help="file with toy MC")
    parser.add_option("-v", "--variable", dest="variable", default='MET',
                      help="set the variable to analyse (defined in config/variable_binning.py)")
    parser.add_option("--com", "--centre-of-mass-energy", dest="CoM", default=13,
                      help='''set the centre of mass energy for analysis.
                      Default = 8 [TeV]''', type=int)
    parser.add_option("-c", "--channel", type='string',
                      dest="channel", default='combined',
                      help="channel to be analysed: electron|muon|combined")
    parser.add_option("-s", type='string',
                      dest="sample", default='madgraph',
                      help="channel to be analysed: electron|muon|combined")

    (options, _) = parser.parse_args()

    centre_of_mass = options.CoM
    measurement_config = XSectionConfig(centre_of_mass)
    make_folder_if_not_exists(options.output_folder)

    use_n_toy = options.n_input_mc
    method = options.method
    variable = options.variable
    sample = options.sample
    tau_value = options.tau_value

    create_unfolding_pull_data(options.file, method, options.channel,
                               centre_of_mass, variable,
                               sample,
                               measurement_config.unfolding_central,
                               use_n_toy,
                               options.output_folder,
                               tau_value)

def make_template_plots(histograms, category, channel):
    global variable, output_folder
    
    for variable_bin in variable_bins_ROOT[variable]:
        path = output_folder + str(measurement_config.centre_of_mass) + 'TeV/' + variable + '/' + category + '/fit_templates/'
        make_folder_if_not_exists(path)
        plotname = path + channel + '_templates_bin_' + variable_bin
        #check if template plots exist already
        for output_format in output_formats:
            if os.path.isfile(plotname + '.' + output_format):
                continue
        canvas = Canvas(width=700, height=500)
        canvas.SetLeftMargin(0.15)
        canvas.SetBottomMargin(0.15)
        canvas.SetTopMargin(0.05)
        canvas.SetRightMargin(0.05)
        legend = plotting.create_legend(x0=0.7, y1=0.8)
        h_signal = histograms[variable_bin]['signal']
        h_VJets = histograms[variable_bin]['V+Jets']
        h_QCD = histograms[variable_bin]['QCD']
        
        h_signal.GetXaxis().SetTitle('Lepton #eta')
        h_signal.GetYaxis().SetTitle('Normalised Events')
        h_signal.GetXaxis().SetTitleSize(0.05)
        h_signal.GetYaxis().SetTitleSize(0.05)
        h_signal.SetMinimum(0)
        h_signal.SetMaximum(0.2)
        h_signal.SetLineWidth(2)
        h_VJets.SetLineWidth(2)
        h_QCD.SetLineWidth(2)
        h_signal.SetLineColor(kRed + 1)
        h_VJets.SetLineColor(kBlue)
        h_QCD.SetLineColor(kYellow)
        h_signal.Draw('hist')
        h_VJets.Draw('hist same')
        h_QCD.Draw('hist same')
        legend.AddEntry(h_signal, 'signal', 'l')
        legend.AddEntry(h_VJets, 'V+Jets', 'l')
        legend.AddEntry(h_QCD, 'QCD', 'l')
        legend.Draw()
        
        cms_label, channel_label = get_cms_labels(channel)
        cms_label.Draw()
        channel_label.Draw()
        
        canvas.Modified()
        canvas.Update()
        for output_format in output_formats:
            canvas.SaveAs(plotname + '.' + output_format)

def plot_central_and_systematics(channel, systematics, exclude=[], suffix='altogether'):
    global variable, k_value, b_tag_bin, met_type

    plt.figure(figsize=(16, 16), dpi=200, facecolor='white')
    axes = plt.axes()
    axes.minorticks_on()
    
    hist_data_central = read_xsection_measurement_results('central', channel)[0]['unfolded_with_systematics']
    hist_data_central.markersize = 2  # points. Imagine, tangible units!
    hist_data_central.marker = 'o'
    
    
    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)

    rplt.errorbar(hist_data_central, axes=axes, label='data', xerr=True)

    for systematic in sorted(systematics):
        if systematic in exclude or systematic == 'central':
            continue

        hist_data_systematic = read_xsection_measurement_results(systematic, channel)[0]['unfolded']
        hist_data_systematic.markersize = 2
        hist_data_systematic.marker = 'o'
        colour_number = systematics.index(systematic) + 2
        if colour_number == 10:
            colour_number = 42
        hist_data_systematic.SetMarkerColor(colour_number)
        if 'PDF' in systematic:
            rplt.errorbar(hist_data_systematic, axes=axes, label=systematic.replace('Weights_', ' '), xerr=False)
        elif met_type in systematic:
            rplt.errorbar(hist_data_systematic, axes=axes, label=met_systematics_latex[systematic.replace(met_type, '')], xerr=False)
        else:
            rplt.errorbar(hist_data_systematic, axes=axes, label=measurements_latex[systematic], xerr=False)
            
    plt.legend(numpoints=1, loc='upper right', prop={'size':25}, ncol=2)
    plt.title(get_cms_labels(channel), CMS.title)
    plt.tight_layout()

    
    path = output_folder + str(measurement_config.centre_of_mass) + 'TeV/' + variable
    make_folder_if_not_exists(path)
    for output_format in output_formats:
        plt.savefig(path + '/normalised_xsection_' + channel + '_' + suffix + '_kv' + str(k_value) + '.' + output_format) 

    plt.close()
    gc.collect()

def make_template_plots_matplotlib(histograms, category, channel):
    global variable, output_folder
    from matplotlib import rc
    rc('text', usetex=True)
    
    for variable_bin in variable_bins_ROOT[variable]:
        path = output_folder + str(measurement_config.centre_of_mass) + 'TeV/' + variable + '/' + category + '/fit_templates/'
        make_folder_if_not_exists(path)
        plotname = path + channel + '_templates_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_signal = histograms[variable_bin]['signal']
        h_VJets = histograms[variable_bin]['V+Jets']
        h_QCD = histograms[variable_bin]['QCD']
        
        h_signal.linecolor = 'red'
        h_VJets.linecolor = 'green'
        h_QCD.linecolor = 'yellow'
        
        h_signal.linewidth = 5
        h_VJets.linewidth = 5
        h_QCD.linewidth = 5
    
        plt.figure(figsize=(14, 10), dpi=200, facecolor='white')
        axes = plt.axes()
        axes.minorticks_on()
        
        plt.xlabel(r'lepton $|\eta|$', CMS.x_axis_title)
        plt.ylabel('normalised to unit area/0.2', CMS.y_axis_title)
        plt.tick_params(**CMS.axis_label_major)
        plt.tick_params(**CMS.axis_label_minor)

        rplt.hist(h_signal, axes=axes, label='signal')
        rplt.hist(h_VJets, axes=axes, label='V+Jets')
        rplt.hist(h_QCD, axes=axes, label='QCD')
        axes.set_ylim([0,0.2])
        
        plt.legend(numpoints=1, loc='upper right', prop=CMS.legend_properties)
        plt.title(get_cms_labels_matplotlib(channel), CMS.title)
        plt.tight_layout()
    
        for output_format in output_formats:
            plt.savefig(plotname + '.' + output_format) 

def create_toy_mc(input_file, sample, output_folder, n_toy, centre_of_mass, ttbar_xsection):
    from tools.file_utilities import make_folder_if_not_exists
    from tools.toy_mc import generate_toy_MC_from_distribution, generate_toy_MC_from_2Ddistribution
    from tools.Unfolding import get_unfold_histogram_tuple
    make_folder_if_not_exists(output_folder)
    input_file_hists = File(input_file)
    output_file_name = get_output_file_name(output_folder, sample, n_toy, centre_of_mass)
    variable_bins = bin_edges_vis.copy()
    with root_open(output_file_name, 'recreate') as f_out:
        for channel in ['combined']:
            for variable in variable_bins:
                output_dir = f_out.mkdir(channel + '/' + variable, recurse=True)
                cd = output_dir.cd
                mkdir = output_dir.mkdir
                h_truth, h_measured, h_response, _ = get_unfold_histogram_tuple(input_file_hists,
                                                                        variable,
                                                                        channel,
                                                                        centre_of_mass = centre_of_mass,
                                                                        ttbar_xsection = ttbar_xsection,
                                                                        visiblePS = True,
                                                                        load_fakes=False)

                cd()

                mkdir('Original')
                cd ('Original')
                h_truth.Write('truth')
                h_measured.Write('measured')
                h_response.Write('response')

                for i in range(1, n_toy+1):
                    toy_id = 'toy_{0}'.format(i)
                    mkdir(toy_id)
                    cd(toy_id)
                    # create histograms
                    # add tuples (truth, measured, response) of histograms
                    truth = generate_toy_MC_from_distribution(h_truth)
                    measured = generate_toy_MC_from_distribution(h_measured)
                    response = generate_toy_MC_from_2Ddistribution(h_response)

                    truth.SetName('truth')
                    measured.SetName('measured')
                    response.SetName('response')

                    truth.Write()
                    measured.Write()
                    response.Write()

def print_xsections(xsections, channel, toFile = True, print_before_unfolding = False):
    global output_folder, variable, k_value, met_type, b_tag_bin
    printout = '=' * 60
    printout += '\n'
    printout += 'Results for %s variable, %s channel, k-value %s, met type %s, %s b-tag region\n' % (variable, channel, k_value, met_type, b_tag_bin)
    if print_before_unfolding:
        printout += 'BEFORE UNFOLDING\n'
    printout += '=' * 60
    printout += '\n'
    printout += '$%s$ bin & $\sigma_{meas}$' % variables_latex[variable]
    printout += '\\\\ \n\hline\n'
    scale = 100
    
    bins = variable_bins_ROOT[variable]
    assert(len(bins) == len(xsections['unfolded_with_systematics']))
    
    for bin_i, variable_bin in enumerate(bins):
        if print_before_unfolding:
            value, error_up, error_down = xsections['measured_with_systematics'][bin_i]
        else:
            value, error_up, error_down = xsections['unfolded_with_systematics'][bin_i]
        r