def data_fit_staterrors(fit, prefix, verbose=0):
    outfilename  = 'bf_stat'
    
    outfile  = tools.set_file(extbase=attr.brfpath, prefix=prefix, 
                              comname=outfilename)
    fdir = outfile.replace(outfilename, '')
    com1_name = 'newfit-cat_files_data_statonly.sh'

    #bashbase = os.path.join(attr.base, 'src', attr.src, 'brf')
    com1 = tools.set_file(extbase=attr.brfpath, prefix=prefix, comname=com1_name)
    sdir = com1.replace(com1_name, '')
    com2 = create_bash_statonly(fdir, sdir, com1)

    if fit or not os.access(outfile, os.F_OK) :
        print 'Running fitter for staterrors ...'
        run_bf_fitter(com2, outfile, verbose)
    
    start_line_str   = 'Fitted parameters'
    start_column_str = '+-'
    end_column_str   = '('
    end_line_str     = 'Difference from seeds'

    data = tools.get_column_from_file(outfile,
                                      start_line_str, start_column_str,
                                      end_column_str, end_line_str,
                                      verbose = verbose)
    datafilename =  'data_fit_staterrors'

    outfile  = tools.set_file(extbase=attr.brfpath, prefix=prefix, 
                              comname=datafilename)

    fo, bakfile = tools.backup_output(outfile)
    for line in data: fo.write('%s\n' %line)
    tools.check_output(fo, outfile, bakfile)

def entries(args):
    dt_type = args[0]
    tag = args[1]
    label_A = args[2]
    label_B = args[3]

    evtpath = attr.evtpath
    rnd='.01'
    err_type=None
    label = None
    diff = 'pct'
    
    tab = DHadTable()
    tab.row_append(['Mode', label_A, label_B])
    tab.head = tab.data[0]
    for mode in attr.modes:
        for sign in [1, -1]:
            evtfile_A = tools.set_file('evt', dt_type, mode, tag, sign=sign,
                                       prefix='dir_'+label_A, extbase=evtpath)
            evtfile_B = tools.set_file('evt', dt_type, mode, tag, sign=sign,
                                       prefix='dir_'+label_B, extbase=evtpath)
            entries_A = tools.count_lines(evtfile_A)
            entries_B = tools.count_lines(evtfile_B)
            if sign == 1:
                modename = attr.modes[mode]['fname']
            else:
                modename = attr.modes[mode]['fnamebar']
            row = [modename, entries_A, entries_B]
            tab.row_append(map(str, row))

    tab.column_append_by_diff_pct('diff(%)', label_B,label_A, rnd=rnd, err_type=err_type)
    tab.output(_tabname)

def create_fig_mode(opts, tag, datatype, mode, lowmass, highmass, label):
    modekey = tools.get_modekey(mode)
    comname = '%s_%s' %(datatype.replace('/', '_'), mode)
    
    prefix = 'dir_%s/kkmass/%s_%s' % (label, lowmass, highmass)
    epsfile = tools.set_file('eps', datatype, modekey, tag,
                            prefix=prefix, extbase=attr.figpath)
    
    txtfile = tools.set_file('txt', datatype, modekey, tag,
                            prefix=prefix, extbase=attr.fitpath)

    head, figname = os.path.split(epsfile)
    figname = figname.replace('.eps', '')

    subdir = 'kkmass/%s_%s' % (lowmass, highmass)
    pdflink = './%s/%s.pdf' %(subdir, figname)
    pnglink = './%s/%s.png' %(subdir, figname)
    
    loglink = '[[../../log/%s/%s/%s.txt][log]]' %(label, subdir, comname)
    figlink = '[[%s][%s]]' %(pdflink, pnglink)
    
    modename = tools.get_orgname_from_fname(mode)
    secname = 'KK mass: %s - %s' % (lowmass, highmass)
    tablink = ''
    if os.access(txtfile, os.F_OK):
        orgtabfile = txtfile.replace('.txt', '.org')
        tab = DHadTable(txtfile)
        tab.output_org(orgtabfile)
        tablink = '#+INCLUDE: "%s"\n' % orgtabfile

    msg = '\n* %s \n  %s \n\n%s\n %s\n' % (
        secname, figlink, tablink, loglink)
        
    return msg

def fit_single_mode(datatype, mode, label, lowmass, highmass, test):
    evtpath = os.path.join(attr.datpath, 'sel', label, 'kkmass')
    modekey = tools.get_modekey(mode)
    tag = 'single'
    evtfile = tools.set_file('evt', datatype, modekey, tag,  prefix='',
                             forceCombine=1, extbase=evtpath)
    
    load_roofit_lib(datatype, label)
    cuts = '%f<kkmass && kkmass<%f' % (lowmass, highmass)
    #err_type = 'ASYM'
    err_type = 'SYMM'
    setGamma, setMres, setR, mc = get_common_parameters(datatype, label)
    Sigma, Fa, Fb, Sa, Sb = get_resolution_paras(modekey, label)

    N1, N2, Nbkgd1, Nbkgd2, md, p, sigmap1, xi = init_paras_single(
        label, datatype, 's', modekey)

    title1 = '%s : %s < KK mass < %s ' %(attr.modes[modekey]['uname'],
                                         lowmass, highmass)
    title2 = '%s : %s < KK mass < %s ' %(attr.modes[modekey]['unamebar'],
                                         lowmass, highmass)
    prefix='dir_%s/kkmass2/%s_%s' % (label, lowmass, highmass)
    epsfile = tools.set_file('eps', datatype, modekey, tag,
                             prefix=prefix, extbase=attr.figpath)
    txtfile = tools.set_file('txt', datatype, modekey, tag,
                             prefix=prefix, extbase=attr.fitpath)

    #mbc_gau_che(evtfile, mc, setMres, setGamma, setR, sigmap1, Sa, Sb, Fa,
    mbc_dline_che(evtfile, mc, setMres, setGamma, setR, sigmap1, Sa, Sb, Fa,
                  Fb, md, p, xi, N1, N2, Nbkgd1, Nbkgd2, title1, title2,
                  epsfile, txtfile, cuts=cuts, err_type=err_type, test=test)
    
    if not test:
        tools.eps2png(epsfile)
        tools.eps2pdf(epsfile)

def signal_single_efficiencies(opts, prefix, prefix_signal, verbose=0):
    effmatrix = []
    errmatrix = []
    mode_sign_list  = [(i,j) for i in modes for j in (1, -1)]
    dt_type = 's'
    tag     = 's'

    label = prefix.replace('dir_', '')
    
    prefix_yld = prefix
    if ('/nofsr' in prefix) or ('/widede' in prefix) or (
        '/gamma/' in prefix) or ('/p/0.5' in prefix):
        prefix_yld = prefix_yld.split('/')[0]

    if '281ipbv0' in label:
        yldfile = os.path.join(attr.datpath, 'brf', '281ipbv0',
                               'yields_and_efficiencies')
    else:
        yldfile = yields_and_efficiencies(opts, [label])

    for i in mode_sign_list:
        mode = i[0]
        sign = i[1]
        effline = []
        errline = []
        tot_gen  = get_generated_numbers_single(i, prefix_signal)
        for j in  mode_sign_list:
            if i == j: 
                eff  = get_diag_eff(prefix_signal, dt_type,
                                    tag, mode, sign, tot_gen)
                err  = get_diag_err(prefix_signal, dt_type,
                                    tag, mode, sign, tot_gen, eff)
            else:      
                eff  = get_off_diag_eff(i, j, yldfile, verbose)
                err  = get_off_diag_err(i, j, yldfile, verbose)

            effline.append(eff)
            errline.append(err)

        effmatrix.append(effline)
        errmatrix.append(errline)

    eff_file_name = 'signal_single_efficiencies_for_werner'
    err_file_name = 'signal_single_efficiencies_errors_for_werner'

    eff_file  = tools.set_file(
        extbase=attr.brfpath, prefix=prefix, 
        comname=eff_file_name)

    err_file  = tools.set_file(
        extbase=attr.brfpath, prefix=prefix, 
        comname=err_file_name)

    output_matrix(effmatrix, eff_file)
    output_matrix(errmatrix, err_file)

def bf_stat_sys(fit , prefix):
    outfilename = 'bf_stat_sys'
    outfile  = tools.set_file(extbase=attr.brfpath, prefix=prefix, 
                              comname=outfilename)
    fdir = outfile.replace(outfilename, '')
    com1_name = 'newfit-cat_files_data.sh'
    com1 = tools.set_file(extbase=attr.brfpath, prefix=prefix, comname=com1_name)
    sdir = com1.replace(com1_name, '')
    com2  = create_bash_stat_syst(fdir, sdir, com1)

    if fit or not os.access(outfile, os.F_OK) :
        print 'Running fitter for stat and syst ...'
        run_bf_fitter(com2, outfile)

def signal_double_efficiencies(prefix, prefix_signal, verbose=0):
    effmatrix = []
    errmatrix = []
    dt_type = 's'
    tag     = 'd'

    mode_pair_list = attr.PossibleDoubleTags
    for i in mode_pair_list:
        txtfile = tools.set_file('txt', dt_type, i, tag, 
                                 prefix=prefix_signal,
                                 extbase=attr.fitpath)

        N      = float(DHadTable(txtfile).cell_get(1, 'Value'))
        N_err  = float(DHadTable(txtfile).cell_get(1, 'Error'))

        tot_gen  = get_generated_numbers_double(i, prefix_signal)

        effline = []
        errline = []

        for j in  mode_pair_list:
            if i == j:
                if '281ipbv0.1' in prefix and i == (1, 1):
                    tot_gen = tot_gen*(0.257402179684/0.140119760479)
                
                eff  = N/tot_gen
                err  = N_err/tot_gen*math.sqrt(1-eff)
            else:
                eff = 0
                err = 0

            effline.append(eff)
            errline.append(err)

        effmatrix.append(effline)
        errmatrix.append(errline)
                
    eff_file_name = 'signal_double_efficiencies_for_werner'
    err_file_name = 'signal_double_efficiencies_errors_for_werner'

    eff_file  = tools.set_file(
        extbase   = attr.brfpath, prefix  = prefix, 
        comname   = eff_file_name)

    err_file  = tools.set_file(
        extbase   = attr.brfpath, prefix  = prefix, 
        comname   = err_file_name)

    output_matrix(effmatrix, eff_file)
    output_matrix(errmatrix, err_file)

def create_script_logfile_jobname(var, datatype, tag, mode, label, test):
    logname = '%s_%s.log' %(datatype, mode)
    logpath = os.path.join(attr.yldlogpath, label, 'var', var)
    logfile = tools.set_file(extbase=logpath, comname=logname)

    content =  '''#!/usr/bin/env python

from yld import var

var.single_tag_mode("%s", "%s", "%s", "%s", test=%s)

'''% (var, datatype, mode, label, test)

    mode, sign = tools.get_modekey_sign(mode)

    if tag == 'single':
        if sign == 1:
            sign = 'p'
        else:
            sign = 'm'
        filename = 'var-%s-%s-m%s-%s-%s.py' % (datatype, tag, mode, sign, label)
        qjobname = 'var%s,%s' % (mode, sign)
    else:
        filename = 'var-%s-%s-m%s-%s-%s.py' % (datatype, tag, mode[0], mode[1], label)
        qjobname = 'var%s,%s' % (mode[0], mode[1])

    file_ = os.path.join(attr.base, 'src', attr.src, 'yld', filename)

    f = UserFile()
    f.data.append(content) 
    f.output(file_)
    os.chmod(file_, 0755)

    return filename, logfile, qjobname

def sqrt_fig_mode(tag, dt_type, mode, label):
    if tag == 'single':
        modename = attr.modes[mode]['orgname']
    else:
        if mode == 'double_all_d0s':
            modename = 'Double all D0s'
        else:
            mode1 = mode[0]
            mode2 = mode[1]
            modename = '%s, %s' % (attr.modes[mode1]['orgname'],
                                   attr.modes[mode2]['orgnamebar'])
    prefix = 'dir_' + label
    rootfile = tools.set_file('root', dt_type, mode, tag,
                              prefix=prefix, extbase=attr.figpath)
    print rootfile
    #load_roofit_lib(dt_type)
    f = ROOT.TFile(rootfile)
    print f
    canvas = f.Get('canvas')
    print canvas

    canvas_1 = canvas.GetListOfPrimitives().FindObject('canvas_1')
    canvas_1.SetLogy(0)

    pdffile = rootfile.replace('.root', '.pdf')
    
    canvas.Print(pdffile)

    
    sys.exit()

def get_generated_numbers_single(i, prefix=''):
    if '281ipb' in prefix:
        mult_fact = 10
    elif '537ipb' in prefix:
        mult_fact = 20
    elif '818ipb' in prefix or 'v13' in prefix:
        mult_fact = 30
    else:
        raise NameError(prefix)
    
    mode = i[0]
    sign = i[1]
    tag_number = modes[mode]['tag_num_s']
    count = tag_number*mult_fact
    if mode == 205: 
        if sign == 1:
            count = 29972 # 30000
        if sign == -1: 
            count = 35969 # 36000

    if 'noFSR' in prefix:
        mode = i[0]
        sign = i[1]
        logfile =  tools.set_file('log', 's',  mode, 's', sign, prefix,
                                      extbase=evtlogbase)
        tab = DHadTable(logfile)
        FSR = tab.cell_get('FSR', 'Value')
        count = count - int(FSR)
        
    return count

def get_data_statonly_crosssectionsdef(prefix):
    filename = 'data_statonly_crosssectionsdef'
    file_ = tools.set_file(extbase=attr.brfpath, prefix=prefix, 
                           comname=filename)
    if '281ipb' in prefix:
        lumi = '281.50'
    elif '537ipb' in prefix:
        lumi = '536.57'
    elif '818ipb' in prefix or 'v13' in prefix:
        lumi = '818.10'
    else:
        raise NameError(prefix)

    content = '''0
4
%s
0.
0
n
0
n
n
n
''' % lumi

    f = UserFile()
    f.append(content)
    f.output(file_, verbose=1)

def get_external_bkg_files(prefix):
    label = prefix.replace('dir_', '')
    bkg_file_list = ['data_external_bkg_bfs_for_werner',
                     'data_statonly_external_bkg_bfs_for_werner']

    for bkg_file in bkg_file_list:
        source_file = os.path.join(attr.srcbrfpath, bkg_file)
        if label in ['818ipbv10', '818ipbv7/generic', '818ipbv12.1',
                     '818ipbv12.2']:
            
            source_file += '_'+label
            if label == '818ipbv7/generic':
                source_file = source_file.replace('/generic', '_generic')
            if label == '818ipbv12.2':
                source_file = source_file.replace('v12.2', 'v12.1')

        dest_file   = tools.set_file(
            extbase   = attr.brfpath, prefix  = prefix, 
            comname   =  bkg_file)

        if '281ipb' in label or label in ['818ipbv12.0', '818ipbv12.0/generic',
                                          '818ipbv12.1', '818ipbv12.2',
                                          '818ipbv12.3', '818ipbv12.4', 
                                          '818ipbv12.8', 'v13.1', 'v13.2']:
            shutil.copy2(source_file, dest_file) 

        elif '537ipb' in label:
            factor = 1.91
            scale_absolute_bkgs(source_file, dest_file, factor)

        elif label in ['818ipbv7', '818ipbv8', '818ipbv9']:
            factor = 2.91
            scale_absolute_bkgs(source_file, dest_file, factor)
        else:
            raise NameError(prefix)

def process_generic_single_mode(datatype, mode, label, test):
    prefix = 'dir_'+label
    sys.stdout.write('Processing backgrounds for %s ...\n' %mode)
    sys.stdout.flush()
    
    files = {}

    inputlabel = label
    if 'noxfeed' in label:
        inputlabel = label.split('/')[0]

    rootfile = tools.get_rootfile(datatype, mode, inputlabel)

    if test:
        sys.stdout.write('\ntest: rootfile %s\n' % rootfile)
    pt = tools.add_rootfile(rootfile)

    bkgname = '%s_%s.evt' %(datatype, mode)
    bkgname = bkgname.replace('/', '_')
    
    bkgfile = tools.set_file(
        extbase=attr.bkgpath(), prefix=prefix, comname=bkgname)

    t_iterstart = time.time()
    nfaked, npassed, ntotal = output_mbc_ebeam(
        pt, mode, label, bkgfile, test)
    dur = tools.duration_human(time.time()-t_iterstart)

    sys.stdout.write('faked %s , passed %s , total %s. done in %s.\n'
                     %(nfaked, npassed, ntotal, dur))
    sys.stdout.flush()

def process_peak_generic_single_mode(datatype, mode, label, test):
    prefix = 'dir_'+label
    sys.stdout.write('Processing backgrounds peak for %s ...\n' %mode)
    sys.stdout.flush()
    
    files = {}

    rootname = '*.root'
    rootfile = os.path.join(attr.datpath, datatype, label, rootname)

    if test:
        sys.stdout.write('\ntest: rootfile %s\n' % rootfile)
    pt = tools.add_rootfile(rootfile)

    bkgname = 'peak_%s_%s.evt' %(datatype, mode)
    bkgname = bkgname.replace('/', '_')
    
    bkgfile = tools.set_file(
        extbase=attr.bkgpath(), prefix=prefix, comname=bkgname)

    t_iterstart = time.time()
    npeak, nfaked, npassed, ntotal = output_peak_mbc_ebeam(
        pt, mode, bkgfile, test)
    dur = tools.duration_human(time.time()-t_iterstart)

    sys.stdout.write('peak %s, faked %s , passed %s , total %s. done in %s.\n'
                     %(npeak, nfaked, npassed, ntotal, dur))
    sys.stdout.flush()

def get_other_files_generic(prefix):
    file_list = ['generic_seeds',
                 'newfit-modedef',
                 'generic_statonly_brratiodef',
                 'generic_statonly_crosssectionsdef'] 

    for _file in file_list:
        source_file = os.path.join(attr.srcbrfpath, _file)
        if '281ipbv0' in prefix:
            if _file == 'generic_seeds':
                source_file += '_281ipbv0'
        elif '281ipbv12' in prefix:
            if _file == 'generic_seeds':
                source_file += '_281ipbv12'
            if 'crosssectionsdef' in _file: 
                source_file += '_281ipb'
        elif '818ipbv12' in prefix:
            if _file == 'generic_seeds':
                source_file += '_818ipbv12'
        else:
            raise NameError(prefix)
        
        dest_file   = tools.set_file(
            extbase   = attr.brfpath, prefix  = prefix, 
            comname   =  _file)

        tools.check_and_copy(source_file, dest_file, verbose=1) 

def get_other_files(prefix):

    file_list = ['newfit-modedef',
                 'zero_single_efficiencies_errors_for_werner',
                 'zero_double_efficiencies_errors_for_werner',
                 'data_seeds',
                 'data_brratiodef', 
                 'data_statonly_brratiodef'] 

    for _file in file_list:
        source_file = os.path.join(attr.srcbrfpath, _file)
        if '/generic' in prefix and _file == 'data_seeds':
            source_file += '_generic'
            if '281ipbv0' in prefix:
                source_file += '_281ipbv0'

        dest_file   = tools.set_file(
            extbase   = attr.brfpath, prefix  = prefix, 
            comname   =  _file)

        if os.access(dest_file, os.F_OK) :
            if filecmp.cmp(source_file, dest_file):
                print 'up-to-date: %s' % _file
                continue
            else:
                print 'Updating %s ...' %_file
                shutil.copy2(dest_file, dest_file+'.bak') 
                
        else:
            print 'Copying %s ...' %_file
        shutil.copy2(source_file, dest_file) 

def init_paras_single(label, dt_type, tag, mode):

    if '281ipb' in label:
        init_para_prefix = 'dir_281ipbv12'
    elif '537ipb' in label:
        init_para_prefix = 'dir_537ipbv12'
    elif '818ipb' in label or 'v13' in label:
        init_para_prefix = 'dir_818ipbv12'
    else:
        raise ValueError(label)

    init_para_base = attr.fitpath

    if dt_type == 'generic':
        if label in ['281ipbv7', '537ipbv7', '818ipbv7']:
            dt_type = 'signal'
    
    tabfile =  tools.set_file('txt', dt_type, mode, tag,
                              prefix=init_para_prefix,
                              extbase=init_para_base)

    sys.stdout.write('Loading init paras from %s ... \n' %tabfile)

    tab     = DHadTable(tabfile)
    
    N1      = float(tab.cell_get('N1', 'Value'))
    N2      = float(tab.cell_get('N2', 'Value'))
    Nbkgd1  = float(tab.cell_get('Nbkgd1', 'Value'))
    Nbkgd2  = float(tab.cell_get('Nbkgd2', 'Value'))
    md      = float(tab.cell_get('md', 'Value'))
    p       = float(tab.cell_get('p', 'Value'))
    sigmap1 = float(tab.cell_get('sigmap1', 'Value'))
    xi      = float(tab.cell_get('xi', 'Value'))

    return N1, N2, Nbkgd1, Nbkgd2, md, p, sigmap1, xi

def create_fig_mode(tag, dt_type, mode, label):
    if tag == 'single':
        modename = attr.modes[mode]['orgname']
    else:
        if mode == 'double_all_d0s':
            modename = 'Double all D0s'
        elif mode == 'double_all_dps':
            modename = 'Double all Dps'
        else:
            mode1 = mode[0]
            mode2 = mode[1]
            modename = '%s, %s' % (attr.modes[mode1]['orgname'],
                                   attr.modes[mode2]['orgnamebar'])

    prefix = 'dir_' + label
    epsfile = tools.set_file('eps', dt_type, mode, tag,
                             prefix=prefix, extbase=attr.figpath)

    head, tail = os.path.split(epsfile)
    comname = tail.replace('.eps', '')
    figlabel = label.split('/')[0]
    relpath = label.replace(figlabel, '')
    
    pdflink = '.%s/%s.pdf' %(relpath, comname)
    pnglink = '.%s/%s.png' %(relpath, comname)
    loglink = '[[../../log/%s/%s.txt][log]]' %(label, comname)
    figlink = '[[%s][%s]]' %(pdflink, pnglink)
    fitpath = attr.fitpath
    tabfile = tools.set_file('txt', dt_type, mode, tag,
                             prefix=prefix, extbase=fitpath)

    if os.access(tabfile, os.F_OK):
        orgtabfile =  tools.set_file('org', dt_type, mode, tag,
                                     prefix=prefix, extbase=fitpath)
        tab = DHadTable(tabfile)
        tab.output_org(orgtabfile)
        abspath = os.path.join(attr.base)#, attr.analysis)
        orgtabfile = orgtabfile.replace(abspath, '../..')
        tablink = '#+INCLUDE: "%s"\n' % orgtabfile
        msg = '\n* %s \n  %s \n\n%s\n %s\n' % (
            modename, figlink, tablink, loglink)
    else:
        sys.stdout.write('File does not exist: %s \n' %tabfile)
        sys.stdout.write('Skipping mode %s ... \n' % str(mode))
        msg = '\n* %s \n Skipped.\n'  % modename
        
    return msg

def get_data_crosssectionsdef(prefix):
    filename = 'data_crosssectionsdef'
    file_ = tools.set_file(extbase=attr.brfpath, prefix=prefix, 
                           comname=filename)
    if '281ipb' in prefix:
        lumi = '281.50'
        lumi_err = '2.8150'
    elif '537ipb' in prefix:
        lumi = '536.57'
        lumi_err = '5.3657'
    elif '818ipb' in prefix or 'v13' in prefix:
        lumi = '818.10'
        lumi_err = '8.1810'
    else:
        raise NameError(prefix)

    bffilename = 'bf_stat'
    bffile = tools.set_file(extbase=attr.brfpath, prefix=prefix,
                            comname=bffilename)
    brf = BrfFile(bffile)

    ne_xsec_err = brf.parsed['sigma(D0D0bar)'].split('+-')[-1].strip()
    ch_xsec_err = brf.parsed['sigma(D+D-)'].split('+-')[-1].strip()
    xsec_err =  brf.parsed['sigma(DDbar)'].split('+-')[-1].strip()
    xsec_ratio_err = brf.parsed['chg/neu'].split('+-')[-1].strip()

    content = '''0
4
%s
%s
0
y
%s
0
y
%s
y
%s
y
%s
''' % (lumi, lumi_err, ne_xsec_err, ch_xsec_err, xsec_err, xsec_ratio_err)

    f = UserFile()
    f.append(content)
    f.output(file_, verbose=1)

def draw_momenta_kkpi(figpath, datatype, mode, label, test):
    modekey = tools.get_modekey(mode)
    sname = attr.modes[modekey]['sname'].lower()
    figname = '%s_momenta' % sname

    epsfile = set_file(extbase=figpath, comname=figname, ext='eps')

    #tab = DHadTable()
    #tab.row_append(['Name', 'Data/MC'])
                   
    # --------------------------------------------------
    h_pkm = {}; h_pkp = {}; h_ppi = {}
    h_pkm_c = {}; h_pkp_c = {}; h_ppi_c = {}

    for datatype in [datatype]:
        selfile = get_selfile(datatype, mode, label, test=test)
        selfile = selfile.replace('/trkmtm/', '/trkmtm2/')

        f = TFile(selfile)
        if datatype == 'signal': # or datatype == 'generic':
            tp = 'mc'
        #if datatype == 'data':
        #    tp = 'data'

        h_pkm[tp]= f.Get('h_pk'+tp)
        h_pkp[tp]= f.Get('h_ppi1'+tp)
        h_ppi[tp]= f.Get('h_ppiz'+tp)

        h_pkm_c[tp]= f.Get('h_pk_c'+tp)
        h_pkp_c[tp]= f.Get('h_ppi1_c'+tp)
        h_ppi_c[tp]= f.Get('h_ppiz_c'+tp)

        f.Clear()

    #ratio = h_pkm['data'].Integral()/h_pkm['mc'].Integral()

    c1 = TCanvas('c1', 'canvas', 900, 900)
    c1.Divide(2,2)
    c1.cd(1)
    h_pkm['mc'].Draw()
    #h_pkm['data'].Draw('PE')
    #h_pkm['mc'].Scale(ratio)
    #h_pkm['mc'].Draw('SAME')
    c1.cd(2)
    h_pkp['mc'].Draw()
    #h_pkp['data'].Draw('PE')
    #h_pkp['mc'].Scale(ratio)
    #h_pkp['mc'].Draw('SAME')
    c1.cd(3)
    c1.cd(4)
    h_ppi['mc'].Draw()
    #h_ppi['data'].Draw('PE')
    #h_ppi['mc'].Scale(ratio)
    #h_ppi['mc'].Draw('SAME')

    c1.Print(epsfile)
    tools.eps2pdf(epsfile)