diff options
| -rwxr-xr-x | plot_sens.py | 41 | ||||
| -rw-r--r-- | submitter/mcmc_dag.py | 2 | ||||
| -rw-r--r-- | submitter/sens_dag.py | 61 | ||||
| -rw-r--r-- | utils/fr.py | 2 | ||||
| -rw-r--r-- | utils/gf.py | 2 | ||||
| -rw-r--r-- | utils/param.py | 4 | ||||
| -rw-r--r-- | utils/plot.py | 63 |
7 files changed, 107 insertions, 68 deletions
diff --git a/plot_sens.py b/plot_sens.py index 9aa5089..de7b348 100755 --- a/plot_sens.py +++ b/plot_sens.py @@ -19,6 +19,7 @@ import numpy as np import numpy.ma as ma from utils import fr as fr_utils +from utils import gf as gf_utils from utils import likelihood as llh_utils from utils import misc as misc_utils from utils import plot as plot_utils @@ -28,6 +29,18 @@ from utils.param import Param, ParamSet, get_paramsets from utils import multinest as mn_utils +import matplotlib as mpl +print mpl.rcParams.keys() +mpl.rcParams['text.usetex'] = True +# mpl.rcParams['text.latex.unicode'] = True +# mpl.rcParams['text.latex.preamble'] = r'\usepackage{cmbright}' +mpl.rcParams['font.family'] = 'sans-serif' +# mpl.rcParams['font.sans-serif'] = 'DejaVu Sans' +# mpl.rcParams['mathtext.fontset'] = 'stixsans' +# mpl.rcParams['mathtext.rm'] = 'DejaVu Sans' +# mpl.rcParams['mathtext.it'] = 'DejaVu Sans:italic' +# mpl.rcParams['mathtext.bf'] = 'DejaVu Sans:bold' + def define_nuisance(): """Define the nuisance parameters.""" @@ -149,6 +162,7 @@ def parse_args(args=None): help='Plot MultiNest evidence or LLH value' ) fr_utils.fr_argparse(parser) + gf_utils.gf_argparse(parser) llh_utils.likelihood_argparse(parser) mn_utils.mn_argparse(parser) nuisance_argparse(parser) @@ -229,8 +243,8 @@ def main(): infile += '/gaussian/' if args.likelihood is Likelihood.GAUSSIAN: infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_')) - infile += '/{0}/{1}/{2}/fr_stat'.format( - *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data]) + infile += '/DIM{0}/fix_ifr/{1}/{2}/{3}/fr_stat'.format( + dim, *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data]) ) + misc_utils.gen_identifier(argsc) print '== {0:<25} = {1}'.format('infile', infile) @@ -264,28 +278,27 @@ def main(): print 'Plotting statistic' argsc = deepcopy(args) - base_infile = args.infile - if args.likelihood is Likelihood.GOLEMFIT: - base_infile += '/golemfit/' - elif args.likelihood is Likelihood.GAUSSIAN: - base_infile += '/gaussian/' for idim, dim in enumerate(args.dimensions): argsc.dimension = dim _, scale_region = fr_utils.estimate_scale(argsc) argsc.scale_region = scale_region + base_infile = args.infile + if args.likelihood is Likelihood.GOLEMFIT: + base_infile += '/golemfit/' + elif args.likelihood is Likelihood.GAUSSIAN: + base_infile += '/gaussian/' if args.likelihood is Likelihood.GAUSSIAN: - infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_')) + base_infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_')) + base_infile += '/DIM{0}/fix_ifr'.format(dim) for isrc, src in enumerate(args.source_ratios): argsc.source_ratio = src - finfile = infile +'/{0}/{1}/{2}/fr_stat'.format( + infile = base_infile +'/{0}/{1}/{2}/fr_stat'.format( *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data]) ) + misc_utils.gen_identifier(argsc) - basename = os.path.dirname(finfile) + '/statrun/' + basename = os.path.dirname(infile) baseoutfile = basename[:5]+basename[5:].replace('data', 'plots') - if args.data: - baseoutfile += '/data/' - baseoutfile += os.path.basename(finfile) + baseoutfile += '/' + os.path.basename(infile) if args.run_method is SensitivityCateg.FULL: outfile = baseoutfile plot_utils.plot_statistic( @@ -334,7 +347,7 @@ def main(): print 'Plotting sensitivities' basename = args.infile[:5]+args.infile[5:].replace('data', 'plots') - baseoutfile = basename + '/{1}/{2}/{3}'.format( + baseoutfile = basename + '/{0}/{1}/{2}/'.format( *map(misc_utils.parse_enum, [args.likelihood, args.stat_method, args.data]) ) diff --git a/submitter/mcmc_dag.py b/submitter/mcmc_dag.py index e40c043..bad73ac 100644 --- a/submitter/mcmc_dag.py +++ b/submitter/mcmc_dag.py @@ -57,7 +57,7 @@ GLOBAL_PARAMS.update(dict( # GolemFit GLOBAL_PARAMS.update(dict( ast = 'p2_0', - data = 'asimov' + data = 'real' )) # Plot diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py index 5e00dfd..0652705 100644 --- a/submitter/sens_dag.py +++ b/submitter/sens_dag.py @@ -27,9 +27,9 @@ GLOBAL_PARAMS = {} sens_eval_bin = 'true' # set to 'all' to run normally GLOBAL_PARAMS.update(dict( sens_run = 'True', - run_method = 'fixed_angle', # full, fixed_angle, corr_angle + run_method = 'corr_angle', # full, fixed_angle, corr_angle stat_method = 'bayesian', - sens_bins = 20, + sens_bins = 15, seed = 'None' )) @@ -46,7 +46,7 @@ dimension = [3] GLOBAL_PARAMS.update(dict( threads = 1, binning = '6e4 1e7 20', - # binning = '1e5 1e7 5', + # binning = '1e5 1e7 20', no_bsm = 'False', scale_region = "1E10", energy_dependance = 'spectral', @@ -65,7 +65,7 @@ GLOBAL_PARAMS.update(dict( # GolemFit GLOBAL_PARAMS.update(dict( ast = 'p2_0', - data = 'asimov' + data = 'real' )) # Plot @@ -73,10 +73,12 @@ GLOBAL_PARAMS.update(dict( plot_statistic = 'True' )) -outfile = 'dagman_FR_SENS_{0}_{1}_{2}_{3}.submit'.format( +outfile = 'dagman_FR_SENS_{0}_{1}_{2}_{3}'.format( GLOBAL_PARAMS['stat_method'], GLOBAL_PARAMS['run_method'], GLOBAL_PARAMS['likelihood'], GLOBAL_PARAMS['data'] ) +# outfile += '_100TeV' +outfile += '.submit' golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit' condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub' @@ -99,6 +101,7 @@ with open(outfile, 'w') as f: output = outchain_head + '/fix_ifr/' if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian': output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_')) + # output += '100TeV/' for r in xrange(sens_runs): print 'run', r f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script)) @@ -119,29 +122,29 @@ with open(outfile, 'w') as f: f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, output)) job_number += 1 - for frs in full_scan_mfr: - print 'frs', frs - output = outchain_head + '/full/' - if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian': - output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_')) - for r in xrange(sens_runs): - print 'run', r - f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script)) - f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim)) - f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0])) - f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1])) - f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2])) - f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, False)) - f.write('VARS\tjob{0}\tsr0="{1}"\n'.format(job_number, 0)) - f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, 0)) - f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, 0)) - if sens_eval_bin.lower() != 'all': - f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, r)) - else: - f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, 'all')) - for key in GLOBAL_PARAMS.iterkeys(): - f.write('VARS\tjob{0}\t{1}="{2}"\n'.format(job_number, key, GLOBAL_PARAMS[key])) - f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, output)) - job_number += 1 + # for frs in full_scan_mfr: + # print 'frs', frs + # output = outchain_head + '/full/' + # if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian': + # output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_')) + # for r in xrange(sens_runs): + # print 'run', r + # f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script)) + # f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim)) + # f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0])) + # f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1])) + # f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2])) + # f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, False)) + # f.write('VARS\tjob{0}\tsr0="{1}"\n'.format(job_number, 0)) + # f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, 0)) + # f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, 0)) + # if sens_eval_bin.lower() != 'all': + # f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, r)) + # else: + # f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, 'all')) + # for key in GLOBAL_PARAMS.iterkeys(): + # f.write('VARS\tjob{0}\t{1}="{2}"\n'.format(job_number, key, GLOBAL_PARAMS[key])) + # f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, output)) + # job_number += 1 print 'dag file = {0}'.format(outfile) diff --git a/utils/fr.py b/utils/fr.py index 639927f..f98c730 100644 --- a/utils/fr.py +++ b/utils/fr.py @@ -248,7 +248,7 @@ def fr_argparse(parser): help='Fold in the spectral index when using GolemFit' ) parser.add_argument( - '--binning', default=[6e4, 1e7, 5], type=float, nargs=3, + '--binning', default=[6e4, 1e7, 20], type=float, nargs=3, help='Binning for spectral energy dependance' ) parser.add_argument( diff --git a/utils/gf.py b/utils/gf.py index cc093e1..c316ed2 100644 --- a/utils/gf.py +++ b/utils/gf.py @@ -65,7 +65,7 @@ def steering_params(args): # For Tianlu # params.years = [999] - # params.minFitEnergy = 1.0e5 # GeV + params.minFitEnergy = args.binning[0] # GeV return params diff --git a/utils/param.py b/utils/param.py index 25558c0..5684e91 100644 --- a/utils/param.py +++ b/utils/param.py @@ -255,12 +255,12 @@ def get_paramsets(args, nuisance_paramset): if hasattr(args, 'dimension'): llh_paramset.append( Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, - tex=r'{\rm log}_{10}\Lambda^{-1}'+get_units(args.dimension), tag=tag) + tex=r'{\rm log}_{10}\left (\Lambda^{-1}'+get_units(args.dimension)+r'\right )', tag=tag) ) elif hasattr(args, 'dimensions'): llh_paramset.append( Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, - tex=r'{\rm log}_{10}\Lambda^{-1} / GeV^{-d+4}', tag=tag) + tex=r'{\rm log}_{10}\left (\Lambda^{-1} / GeV^{-d+4}\right )', tag=tag) ) if not args.fix_source_ratio: tag = ParamTag.SRCANGLES diff --git a/utils/plot.py b/utils/plot.py index a51b0f1..65635bc 100644 --- a/utils/plot.py +++ b/utils/plot.py @@ -22,7 +22,8 @@ from matplotlib.offsetbox import AnchoredText from getdist import plots, mcsamples from utils import misc as misc_utils -from utils.enums import EnergyDependance, Likelihood, ParamTag, StatCateg +from utils.enums import DataType, EnergyDependance +from utils.enums import Likelihood, ParamTag, StatCateg from utils.fr import angles_to_u, angles_to_fr rc('text', usetex=False) @@ -34,12 +35,12 @@ def centers(x): def get_units(dimension): - if dimension == 3: return r' / GeV' + if dimension == 3: return r' / \:GeV' if dimension == 4: return r'' - if dimension == 5: return r' / GeV^{-1}' - if dimension == 6: return r' / GeV^{-2}' - if dimension == 7: return r' / GeV^{-3}' - if dimension == 8: return r' / GeV^{-4}' + if dimension == 5: return r' / \:GeV^{-1}' + if dimension == 6: return r' / \:GeV^{-2}' + if dimension == 7: return r' / \:GeV^{-3}' + if dimension == 8: return r' / \:GeV^{-4}' def calc_nbins(x): @@ -230,7 +231,7 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset): g.export(outfile+'_elements.'+of) -def myround(x, base=5, up=False, down=False): +def myround(x, base=1, up=False, down=False): if up == down and up is True: assert 0 if up: return int(base * np.round(float(x)/base-0.5)) elif down: return int(base * np.round(float(x)/base+0.5)) @@ -271,7 +272,7 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None): ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.3, linewidth=1) at = AnchoredText( - '\n'+fig_text, prop=dict(size=7), frameon=True, loc=2 + fig_text, prop=dict(size=10), frameon=True, loc=2 ) at.patch.set_boxstyle("round,pad=0.1,rounding_size=0.5") ax.add_artist(at) @@ -296,7 +297,7 @@ def plot_sens_full(data, outfile, outformat, args): ax.set_xlim(args.dimensions[0]-1, args.dimensions[-1]+1) ax.set_xticklabels([''] + xticks + ['']) ax.set_xlabel(r'BSM operator dimension ' + r'$d$') - ax.set_ylabel(r'${\rm log}_{10} \Lambda^{-1} / GeV^{-d+4}$') + ax.set_ylabel(r'${\rm log}_{10} \left (\Lambda^{-1} / GeV^{-d+4} \right )$') argsc = deepcopy(args) for idim in xrange(len(data)): @@ -360,7 +361,6 @@ def plot_sens_fixed_angle(data, outfile, outformat, args): print 'Making FIXED_ANGLE sensitivity plot' colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'} - # xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$'] xticks = [r'$\mathcal{O}_{e\mu}$', r'$\mathcal{O}_{e\tau}$', r'$\mathcal{O}_{\mu\tau}$'] argsc = deepcopy(args) for idim in xrange(len(data)): @@ -376,7 +376,7 @@ def plot_sens_fixed_angle(data, outfile, outformat, args): ax.set_xlim(0, len(xticks)+1) ax.set_xticklabels([''] + xticks + ['']) ax.set_xlabel(r'BSM operator angle') - ax.set_ylabel(r'${\rm log}_{10} \Lambda^{-1} / GeV^{-d+4}$') + ax.set_ylabel(r'${\rm log}_{10} \left (\Lambda^{-1}' + get_units(dim) +r'\right )$') for isrc in xrange(len(data[idim])): src = args.source_ratios[isrc] @@ -391,7 +391,6 @@ def plot_sens_fixed_angle(data, outfile, outformat, args): if args.stat_method is StatCateg.BAYESIAN: reduced_ev = -(statistic - null) al = scales[reduced_ev > np.log(10**(3/2.))] # Strong degree of belief - # al = scales[reduced_ev > np.log(10**(1/2.))] elif args.stat_method is StatCateg.FREQUENTIST: reduced_ev = -2*(statistic - null) al = scales[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks @@ -401,22 +400,24 @@ def plot_sens_fixed_angle(data, outfile, outformat, args): dim, src, reduced_ev ) continue + arr_len = dim-2 lim = al[0] print 'limit = {0}'.format(lim) - label = '[{0}, {1}, {2}]'.format(*misc_utils.solve_ratio(src)) + label = '{0} : {1} : {2}'.format(*misc_utils.solve_ratio(src)) if lim < yranges[0]: yranges[0] = lim - if lim > yranges[1]: yranges[1] = lim+5 + if lim > yranges[1]: yranges[1] = lim+arr_len+1 # if lim > yranges[1]: yranges[1] = lim + xoff = 0.15 line = plt.Line2D( - (ian+1-0.1, ian+1+0.1), (lim, lim), lw=3, color=colour[isrc], label=label + (ian+1-xoff, ian+1+xoff), (lim, lim), lw=2., color=colour[isrc], label=label ) ax.add_line(line) if len(legend_handles) < isrc+1: legend_handles.append(line) - x_offset = isrc*0.05 - 0.05 + x_offset = isrc*xoff/2. - xoff/2. ax.annotate( - s='', xy=(ian+1+x_offset, lim), xytext=(ian+1+x_offset, lim+3), - arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[isrc]} + s='', xy=(ian+1+x_offset, lim-0.01), xytext=(ian+1+x_offset, lim+arr_len), + arrowprops={'arrowstyle': '<-', 'lw': 2., 'color':colour[isrc]} ) try: @@ -424,7 +425,29 @@ def plot_sens_fixed_angle(data, outfile, outformat, args): ax.set_ylim(yranges) except: pass - ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right') + legend = ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right', + title=r'$\nu_e:\nu_\mu:\nu_\tau{\rm\:\:at\:\:source}$', + framealpha=1., edgecolor='black') + plt.setp(legend.get_title(), fontsize='8') + legend.get_frame().set_linestyle('-') + + an_text = 'Dimension {0}'.format(dim) + at = AnchoredText( + an_text, prop=dict(size=10), frameon=True, loc=2 + ) + at.patch.set_boxstyle("round,pad=0.1,rounding_size=0.5") + ax.add_artist(at) + + fig.text(0.42, 0.8, 'Excluded', color='red', fontsize=20, ha='center', + va='center', fontweight='bold') + + if args.data is DataType.REAL: + fig.text(0.805, 0.14, 'IceCube Preliminary', color='red', fontsize=11, + ha='center', va='center') + elif args.data is DataType.ASIMOV: + fig.text(0.805, 0.14, 'IceCube Simulation', color='red', fontsize=11, + ha='center', va='center') + for ymaj in ax.yaxis.get_majorticklocs(): ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1) for xmaj in ax.xaxis.get_majorticklocs(): @@ -469,7 +492,7 @@ def plot_sens_corr_angle(data, outfile, outformat, args): ax = fig.add_subplot(111) ax.set_ylim(0, 1) ax.set_ylabel(labels[ian]) - ax.set_xlabel(r'${\rm log}_{10} \Lambda^{-1}'+get_units(dim)+r'$') + ax.set_xlabel(r'${\rm log}_{10} \left (\Lambda^{-1}'+get_units(dim)+r'\right )$') xranges = [np.inf, -np.inf] legend_handles = [] |