Wed 17 Apr 09:28:30 CDT 2019

6 files changed, 231 insertions, 183 deletions

diff --git a/plot_sens.py b/plot_sens.py
index 3fb585e..42b711f 100755
--- a/plot_sens.py
+++ b/plot_sens.py
@@ -27,7 +27,7 @@ from utils.misc import gen_identifier, SortingHelpFormatter
 from utils.param import Param, ParamSet
 
 
-MASK_X = (0.3, 0.8)
+MASK_X = (0.3, 0.89)
 
 
 def process_args(args):
@@ -227,13 +227,18 @@ def main():
         *map(parse_enum, [args.stat_method, args.data])
     ) + r'{0}'.format(prefix)
 
+    argsc = deepcopy(args)
     if args.plot_x:
-        plot_utils.plot_x(
-            data      = data,
-            outfile   = baseoutfile + '/hese_x',
-            outformat = ['png', 'pdf'],
-            args      = args,
-        )
+        for idim, dim in enumerate(args.dimensions):
+            print '|||| DIM = {0}'.format(dim)
+            argsc.dimension = dim
+            plot_utils.plot_x(
+                data      = data[idim],
+                outfile   = baseoutfile + '/hese_x_DIM{0}'.format(dim),
+                outformat = ['png', 'pdf'],
+                args      = argsc,
+                normalise = True
+            )
 
     if args.plot_table:
         plot_utils.plot_table_sens(
diff --git a/sens.py b/sens.py
index 7213381..cafbf10 100755
--- a/sens.py
+++ b/sens.py
@@ -205,7 +205,8 @@ def main():
     ) + misc_utils.gen_identifier(args)
 
     # Setup Golemfit.
-    gf_utils.setup_fitter(args, asimov_paramset)
+    if args.run_mn:
+        gf_utils.setup_fitter(args, asimov_paramset)
 
     # Initialise data structure.
     stat_arr = np.full((eval_dim, 2), np.nan)
@@ -254,13 +255,15 @@ def main():
         # Cleanup.
         if reset_range is not None:
             scale_prm.ranges = reset_range
-        try:
-            for f in glob.glob(prefix + '*'):
-                print 'cleaning file {0}'.format(f)
-                os.remove(f)
-        except:
-            print 'got error trying to cleanup, continuing'
-            pass
+
+        if args.run_mn and not args.debug:
+            try:
+                for f in glob.glob(prefix + '*'):
+                    print 'cleaning file {0}'.format(f)
+                    os.remove(f)
+            except:
+                print 'got error trying to cleanup, continuing'
+                pass
 
     misc_utils.make_dir(outfile)
     print 'Saving to {0}'.format(outfile+'.npy')
diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py
index a4ec65c..601985d 100644
--- a/submitter/sens_dag.py
+++ b/submitter/sens_dag.py
@@ -9,26 +9,27 @@ import numpy as np
 #     (0, 1, 0),
 # ]
 
-MASK_X = (0.3, 0.8)
-x_segments = 50
+XLIMS = (0., 0.3)
+x_segments = 100
 x_array = np.linspace(0, 1, x_segments)
 sources = []
 for x in x_array:
-    if x > MASK_X[0] and x < MASK_X[1]: continue
-    sources.append([x, 1-x, 0])
+    if x >= XLIMS[0] and x <= XLIMS[1]:
+        sources.append([x, 1-x, 0])
 
 dims = [
     6
 ]
 
 textures = [
-    'OEU', 'OET', 'OUT'
+    # 'OEU', 'OET', 'OUT'
+    'OET'
 ]
 
 datadir = '/data/user/smandalia/flavour_ratio/data/sensitivity'
 
-prefix = ''
-# prefix = '_noprior'
+# prefix = ''
+prefix = '_OET'
 
 golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
 condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub'
@@ -44,17 +45,18 @@ GLOBAL_PARAMS.update(dict(
 
 # MultiNest
 GLOBAL_PARAMS.update(dict(
-    mn_live_points = 1000,
+    # mn_live_points = 1000,
     # mn_live_points = 600,
-    # mn_live_points = 200,
+    mn_live_points = 200,
     # mn_tolerance   = 0.1,
     mn_tolerance   = 0.3,
-    mn_output      = './mnrun'
+    mn_output      = './mnrun',
+    run_mn         = 'True'
 ))
 
 # FR
 GLOBAL_PARAMS.update(dict(
-    threads = 2,
+    threads = 1,
     binning = '6e4 1e7 20',
     no_bsm  = 'False'
 ))
diff --git a/submitter/sens_submit.sub b/submitter/sens_submit.sub
index 7d21c96..315ef39 100644
--- a/submitter/sens_submit.sub
+++ b/submitter/sens_submit.sub
@@ -1,5 +1,5 @@
 Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/sens.py
-Arguments = "--ast $(ast) --data $(data) --dimension $(dimension) --no-bsm $(no_bsm) --datadir $(datadir) --seed $(seed) --source-ratio $(sr0) $(sr1) $(sr2) --threads $(threads) --binning $(binning) --texture $(texture) --segments $(segments) --eval-segment $(eval_segment) --stat-method $(stat_method) --mn-live-points $(mn_live_points) --mn-tolerance $(mn_tolerance) --mn-output $(mn_output)"
+Arguments = "--ast $(ast) --data $(data) --dimension $(dimension) --no-bsm $(no_bsm) --datadir $(datadir) --seed $(seed) --source-ratio $(sr0) $(sr1) $(sr2) --threads $(threads) --binning $(binning) --texture $(texture) --segments $(segments) --eval-segment $(eval_segment) --stat-method $(stat_method) --mn-live-points $(mn_live_points) --mn-tolerance $(mn_tolerance) --mn-output $(mn_output) --run-mn $(run_mn)"
 
 # All logs will go to a single file 
 log    = /scratch/smandalia/flavour_ratio/submitter/logs/job_$(Cluster).log
@@ -10,13 +10,14 @@ getenv = True
 # environment = "X509_USER_PROXY=x509up_u14830" 
 
 request_memory = 3GB
-request_cpus = 2
+request_cpus = 1
 
 initialdir = /home/smandalia/condor
 
 Universe        = vanilla 
 Notification    = never 
 
+# +AccountingGroup="quicktest.$ENV(USER)"
 # +AccountingGroup="sanctioned.$ENV(USER)"
 # run on both SL5 and 6
 # +WantRHEL6 = True
diff --git a/utils/mn.py b/utils/mn.py
index 335df96..ac42858 100644
--- a/utils/mn.py
+++ b/utils/mn.py
@@ -16,7 +16,7 @@ import numpy as np
 from pymultinest import analyse, run
 
 from utils import llh as llh_utils
-from utils.misc import gen_identifier, make_dir, solve_ratio
+from utils.misc import gen_identifier, make_dir, parse_bool, solve_ratio
 
 
 def CubePrior(cube, ndim, n_params):
@@ -58,6 +58,10 @@ def mn_argparse(parser):
         '--mn-output', type=str, default='./mnrun/',
         help='Folder to store MultiNest evaluations'
     )
+    parser.add_argument(
+        '--run-mn', type=parse_bool, default='True',
+        help='Run MultiNest'
+    )
 
 
 def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, prefix='mn'):
@@ -75,19 +79,21 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, prefix='mn'):
         args            = args,
     )
 
-    make_dir(prefix)
-    print 'Running evidence calculation for {0}'.format(prefix)
-    run(
-        LogLikelihood              = lnProbEval,
-        Prior                      = CubePrior,
-        n_dims                     = n_params,
-        n_live_points              = args.mn_live_points,
-        evidence_tolerance         = args.mn_tolerance,
-        outputfiles_basename       = prefix,
-        importance_nested_sampling = True,
-        resume                     = False,
-        verbose                    = True
-    )
+    if args.run_mn:
+        make_dir(prefix)
+        print 'Running evidence calculation for {0}'.format(prefix)
+        run(
+            LogLikelihood              = lnProbEval,
+            Prior                      = CubePrior,
+            n_dims                     = n_params,
+            n_live_points              = args.mn_live_points,
+            evidence_tolerance         = args.mn_tolerance,
+            outputfiles_basename       = prefix,
+            importance_nested_sampling = True,
+            # resume                     = False,
+            resume                     = True,
+            verbose                    = True
+        )
 
     analyser = analyse.Analyzer(
         outputfiles_basename=prefix, n_params=n_params
diff --git a/utils/plot.py b/utils/plot.py
index abd4548..db65dda 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -28,6 +28,7 @@ from matplotlib import pyplot as plt
 from matplotlib.offsetbox import AnchoredText
 from matplotlib.lines import Line2D
 from matplotlib.patches import Patch
+from matplotlib.patches import Arrow
 
 from getdist import plots, mcsamples
 
@@ -111,6 +112,7 @@ def cmap_discretize(cmap, N):
 
 def get_limit(scales, statistic, args, mask_initial=False):
     max_st = np.max(statistic)
+    print 'scales, stat', zip(scales, statistic)
     if args.stat_method is StatCateg.BAYESIAN:
         if (statistic[0] - max_st) > np.log(10**(BAYES_K)):
             raise AssertionError('Discovered LV!')
@@ -118,8 +120,10 @@ def get_limit(scales, statistic, args, mask_initial=False):
     try:
         tck, u = splprep([scales, statistic], s=0)
     except:
-        return None
-    sc, st = splev(np.linspace(0, 1, 10000), tck)
+        print 'Failed to spline'
+        # return None
+        raise
+    sc, st = splev(np.linspace(0, 1, 1000), tck)
 
     if mask_initial:
         scales_rm = sc[sc >= scales[1]]
@@ -132,7 +136,7 @@ def get_limit(scales, statistic, args, mask_initial=False):
     null = statistic[min_idx]
     if args.stat_method is StatCateg.BAYESIAN:
         reduced_ev = -(statistic_rm - null)
-        print 'reduced_ev', reduced_ev
+        print '[reduced_ev > np.log(10**(BAYES_K))]', np.sum([reduced_ev > np.log(10**(BAYES_K))])
         al = scales_rm[reduced_ev > np.log(10**(BAYES_K))]
     else:
         assert 0
@@ -142,11 +146,11 @@ def get_limit(scales, statistic, args, mask_initial=False):
         )
         return None
     if reduced_ev[-1] < np.log(10**(BAYES_K)) - 0.1:
-        print 'Peaked contour does not exclude large scales! For ' \
+        print 'Warning, peaked contour does not exclude large scales! For ' \
             'DIM {0} [{1}, {2}, {3}]!'.format(
                 args.dimension, *args.source_ratio
             )
-        return None
+        # return None
     lim = al[0]
     print 'limit = {0}'.format(lim)
     return lim
@@ -623,13 +627,11 @@ def plot_table_sens(data, outfile, outformat, args):
         fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150)
 
 
-def plot_x(data, outfile, outformat, args):
+def plot_x(data, outfile, outformat, args, normalise=False):
     """Limit plot as a function of the source flavour ratio for each operator
     texture."""
     print 'Making X sensitivity plot'
-    argsc = deepcopy(args)
-
-    dims = args.dimensions
+    dim = args.dimension
     srcs = args.source_ratios
     x_arr = np.array([i[0] for i in srcs])
     if args.texture is Texture.NONE:
@@ -639,137 +641,167 @@ def plot_x(data, outfile, outformat, args):
 
     # Rearrange data structure
     r_data = np.full((
-        data.shape[0], data.shape[2], data.shape[1], data.shape[3], data.shape[4]
+        data.shape[1], data.shape[0], data.shape[2], data.shape[3]
     ), np.nan)
 
-    for idim in xrange(data.shape[0]):
-        for isrc in xrange(data.shape[1]):
-            for itex in xrange(data.shape[2]):
-                r_data[idim][itex][isrc] = data[idim][isrc][itex]
+    for isrc in xrange(data.shape[0]):
+        for itex in xrange(data.shape[1]):
+            r_data[itex][isrc] = data[isrc][itex]
     r_data = ma.masked_invalid(r_data)
     print r_data.shape, 'r_data.shape'
 
-    fig = plt.figure(figsize=(8, 6))
-    gs = gridspec.GridSpec(len(dims), 1)
-    gs.update(hspace=0.15)
+    if normalise:
+        fig = plt.figure(figsize=(7, 6))
+    else:
+        fig = plt.figure(figsize=(8, 6))
+    ax = fig.add_subplot(111)
 
-    xlims = (-60, -20)
-    ylims = (0, 1)
+    if normalise:
+        ylims = (-10, 8)
+    else:
+        ylims = (-60, -20)
+    xlims = (0, 1)
 
-    colour = {0:'red', 1:'blue', 2:'green'}
-    rgb_co = {0:[1,0,0], 1:[0,0,1], 2:[0.0, 0.5019607843137255, 0.0]}
+    colour = {0:'red', 1:'green', 2:'blue'}
+    rgb_co = {0:[1,0,0], 1:[0.0, 0.5019607843137255, 0.0], 2:[0,0,1]}
 
-    first_ax = None
     legend_log = []
     legend_elements = []
-
-    for idim, dim in enumerate(dims):
-        print '|||| DIM = {0}'.format(dim)
-        argsc.dimension = dim
-        ax = fig.add_subplot(gs[idim])
-
-        if first_ax is None:
-            first_ax = ax
-
-        ax.set_xlim(xlims)
-        ax.set_ylim(ylims)
-        yticks = []
-        ylabels = []
-        if idim == 0:
-            yticks += [0]
-            ylabels += [r'$0$']
-        yticks += [1/3., 0.5, 2/3.]
-        ylabels += [r'$\frac{1}{3}$', r'$\frac{1}{2}$', r'$\frac{2}{3}$']
-        if idim == len(dims)-1:
-            yticks += [1]
-            ylabels += [r'$1$']
-        ax.set_yticks([], minor=True)
-        ax.set_yticks(yticks, minor=False)
-        ax.set_yticklabels(ylabels, fontsize=13)
-        for xmaj in ax.xaxis.get_majorticklocs():
+    labelsize = 13
+    largesize = 17
+
+    ax.set_xlim(xlims)
+    ax.set_ylim(ylims)
+    xticks = [0, 1/3., 0.5, 2/3., 1]
+    # xlabels = [r'$0$', r'$\frac{1}{3}$', r'$\frac{1}{2}$', r'$\frac{2}{3}$', r'$1$']
+    xlabels = [r'$0$', r'$1 / 3$', r'$1/2$', r'$2/3$', r'$1$']
+    ax.set_xticks([], minor=True)
+    ax.set_xticks(xticks, minor=False)
+    ax.set_xticklabels(xlabels, fontsize=largesize)
+    for ymaj in ax.yaxis.get_majorticklocs():
+        ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.2, linewidth=1)
+    for xmaj in xticks:
+        if xmaj == 1/3.:
+            ax.axvline(x=xmaj, ls='-', color='gray', alpha=0.5, linewidth=1)
+        else:
             ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.2, linewidth=1)
-        for ymaj in yticks:
-            if ymaj == 1/3.:
-                ax.axhline(y=ymaj, ls='-', color='gray', alpha=0.5, linewidth=1)
-            else:
-                ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.2, linewidth=1)
-        ax.get_xaxis().set_visible(False)
-
-        ax.text(
-            1.01, 0.5, r'$d = {0}$'.format(dim), fontsize=16,
-            rotation='90', verticalalignment='center', transform=ax.transAxes
-        )
-
-        ax.text(
-            0.01, (1/3.)-0.05, r'$f_{\alpha}^S=(1:2:0)$', fontsize=13,
-            transform=ax.transAxes
-        )
 
-        for itex, tex in enumerate(textures):
-            print '|||| TEX = {0}'.format(tex)
-            lims = np.full(len(srcs), np.nan)
-
-            for isrc, src in enumerate(srcs):
-                x = src[0]
-                print '|||| X = {0}'.format(x)
-                argsc.source_ratio = src
-                try:
-                    scales, statistic = ma.compress_rows(r_data[idim][itex][isrc]).T
-                except:
-                    continue
-                lim = get_limit(scales, statistic, argsc, mask_initial=True)
-                if lim is None: continue
-                lims[isrc] = lim
-
-            lims = ma.masked_invalid(lims)
-            size = np.sum(~lims.mask)
-            if size == 0: continue
-
-            tck, u = splprep([x_arr[~lims.mask], lims[~lims.mask]], s=0, k=1)
-            y, x = splev(np.linspace(0, 1, 100), tck)
-            ax.scatter(lims, x_arr, marker='o', s=10, alpha=1, zorder=5,
-                       color='red')
-            ax.fill_betweenx(y, x, 0, color=colour[itex], alpha=1.)
-
-            if itex not in legend_log:
-                legend_log.append(itex)
-                label = texture_label(tex)[:-1] + r'\:{\rm\:texture}$'
-                legend_elements.append(
-                    Patch(facecolor=rgb_co[itex]+[0.3],
-                          edgecolor=rgb_co[itex]+[1], label=label)
-                )
+    ax.text(
+        (1/3.)+0.01, 0.01,  r'$f_{\alpha}^S=(1:2:0)$', fontsize=labelsize,
+        transform=ax.transAxes, rotation='vertical', va='bottom'
+    )
+    ax.text(
+        0.96, 0.01, r'$f_{\alpha}^S=(1:0:0)$', fontsize=labelsize,
+        transform=ax.transAxes, rotation='vertical', va='bottom', ha='left'
+    )
+    ax.text(
+        0.01, 0.01, r'$f_{\alpha}^S=(0:1:0)$', fontsize=labelsize,
+        transform=ax.transAxes, rotation='vertical', va='bottom'
+    )
+    ax.text(
+        0.07, 0.46, r'${\rm \bf Excluded}$', fontsize=largesize,
+        transform=ax.transAxes, color = 'g', rotation='vertical', zorder=10
+    )
+    ax.text(
+        0.95, 0.46, r'${\rm \bf Excluded}$', fontsize=largesize,
+        transform=ax.transAxes, color = 'b', rotation='vertical', zorder=10
+    )
 
-        LV_lim = np.log10(LV_ATMO_90PC_LIMITS[dim])
-        ax.add_patch(patches.Rectangle(
-            (LV_lim[1], ylims[0]), LV_lim[0]-LV_lim[1], np.diff(ylims),
-            fill=False, hatch='\\\\'
-        ))
-
-        if dim in PLANCK_SCALE.iterkeys():
-            ps = np.log10(PLANCK_SCALE[dim])
-            if ps < xlims[0]:
-                ax.annotate(
-                    s='', xy=(xlims[0], 1), xytext=(xlims[0]+1, 1),
-                    arrowprops={'arrowstyle': '->, head_length=0.2',
-                                'lw': 1, 'color':'purple'}
-                )
-            elif ps > xlims[1]:
-                ax.annotate(
-                    s='', xy=(xlims[1]-1, 1), xytext=(xlims[1], 1),
-                    arrowprops={'arrowstyle': '<-, head_length=0.2',
-                                'lw': 1, 'color':'purple'}
-                )
-            else:
-                ax.axvline(x=ps, color='purple', alpha=1., linewidth=1.5)
+    for itex, tex in enumerate(textures):
+        print '|||| TEX = {0}'.format(tex)
+        lims = np.full(len(srcs), np.nan)
+
+        for isrc, src in enumerate(srcs):
+            x = src[0]
+            print '|||| X = {0}'.format(x)
+            args.source_ratio = src
+            d = r_data[itex][isrc]
+            if np.sum(d.mask) > 0:
+                continue
+            scales, statistic = ma.compress_rows(d).T
+            lim = get_limit(scales, statistic, args, mask_initial=True)
+            if lim is None: continue
+            if normalise:
+                lim -= np.log10(PLANCK_SCALE[dim])
+            lims[isrc] = lim
+
+        lims = ma.masked_invalid(lims)
+        size = np.sum(~lims.mask)
+        if size == 0: continue
+
+        print 'x_arr, lims', zip(x_arr, lims)
+        if normalise:
+            zeropoint = 100
+        else:
+            zeropoint = 0
+        lims[lims.mask] = zeropoint
+        tck, u = splprep([x_arr, lims], s=0, k=1)
+        x, y = splev(np.linspace(0, 1, 1000), tck)
+        y = gaussian_filter(y, sigma=4)
+        ax.fill_between(x, y, zeropoint, color=colour[itex], alpha=0.3)
+        # ax.scatter(x, y, color='black', s=1)
+        # ax.scatter(x_arr, lims, color=colour[itex], s=8)
+
+        if itex not in legend_log:
+            legend_log.append(itex)
+            label = texture_label(tex)[:-1] + r'\:{\rm\:texture}$'
+            legend_elements.append(
+                Patch(facecolor=rgb_co[itex]+[0.3],
+                      edgecolor=rgb_co[itex]+[1], label=label)
+            )
 
-    fig.text(0.06, 0.5,
-             r'$x\: : \: {\rm Source\:Ratio}\:[\:f_{\alpha}^S=\left (\:x,\left (1-x\right ),0\:\right )\:]$',
-             va='center', rotation='vertical', fontsize=19)
+    LV_lim = np.log10(LV_ATMO_90PC_LIMITS[dim])
+    if normalise:
+        LV_lim -= np.log10(PLANCK_SCALE[dim])
+    ax.add_patch(patches.Rectangle(
+        (xlims[0], LV_lim[1]), np.diff(xlims), LV_lim[0]-LV_lim[1],
+        fill=False, hatch='\\\\'
+    ))
+
+    if dim in PLANCK_SCALE:
+        ps = np.log10(PLANCK_SCALE[dim])
+        if normalise:
+            ps -= np.log10(PLANCK_SCALE[dim])
+            ax.add_patch(Arrow(
+                0.27, -0.009, 0, -5, width=0.12, capstyle='butt',
+                facecolor='purple', fill=True, alpha=0.8,
+                edgecolor='darkmagenta'
+            ))
+            ax.add_patch(Arrow(
+                0.82, -0.009, 0, -5, width=0.12, capstyle='butt',
+                facecolor='purple', fill=True, alpha=0.8,
+                edgecolor='darkmagenta'
+            ))
+
+            ax.text(
+                0.3, 0.4, r'${\rm \bf Quantum\:Gravity\:Frontier}$',
+                fontsize=largesize-2, transform=ax.transAxes, va='top',
+                ha='left', color='purple'
+            )
+        ax.axhline(y=ps, color='purple', alpha=1., linewidth=1.5)
+
+    if normalise:
+        fig.text(
+            0.02, 0.5,
+            r'${\rm New\:Physics\:Scale}\:[\:{\rm log}_{10} \left (\Lambda_{' +
+            r'\:d={0}'.format(args.dimension)+r'}\:/\:{\rm M}_{\:\rm Planck}^{\:'+
+            r'{0}'.format(args.dimension-4)+ r'}\right )\: ]$', ha='left',
+            va='center', rotation='vertical', fontsize=largesize
+        )
+    else:
+        fig.text(
+            0.02, 0.5,
+            r'${\rm New\:Physics\:Scale}\:[\:{\rm log}_{10} \left (\Lambda_{' +
+            r'\:d={0}'.format(args.dimension)+r'}^{-1}\:' + get_units(args.dimension) +
+            r'\right )\: ]$', ha='left',
+            va='center', rotation='vertical', fontsize=largesize
+        )
 
-    ax.get_xaxis().set_visible(True)
-    ax.set_xlabel(r'${\rm New\:Physics\:Scale}\:[\:{\rm log}_{10} (\Lambda^{-1}\:/\:{\rm GeV}^{-d+4})\: ]$',
-                  fontsize=19)
-    ax.tick_params(axis='x', labelsize=16)
+    ax.set_xlabel(
+        r'${\rm Source\:Flavour\:Ratio}\:[\:f_{\alpha}^S=\left (\:x:1-x:0\:\right )\:]$',
+        fontsize=largesize
+    )
+    ax.tick_params(axis='x', labelsize=largesize-1)
 
     purple = [0.5019607843137255, 0.0, 0.5019607843137255]
     legend_elements.append(
@@ -778,25 +810,24 @@ def plot_x(data, outfile, outformat, args):
     legend_elements.append(
         Patch(facecolor='none', hatch='\\\\', edgecolor='k', label='IceCube, Nature.Phy.14,961(2018)')
     )
-    legend = first_ax.legend(
-        handles=legend_elements, prop=dict(size=11), loc='center left',#loc='upper left',
-        title='Excluded regions', framealpha=1., edgecolor='black',
-        frameon=True
+    legend = ax.legend(
+        handles=legend_elements, prop=dict(size=labelsize-2),
+        loc='upper center', title='Excluded regions', framealpha=1.,
+        edgecolor='black', frameon=True, bbox_to_anchor=(0.55, 1)
     )
-    first_ax.set_zorder(10)
-    plt.setp(legend.get_title(), fontsize='11')
+    plt.setp(legend.get_title(), fontsize=labelsize)
     legend.get_frame().set_linestyle('-')
 
-    ybound = 0.61
-    if args.data is DataType.REAL:
-        fig.text(0.278, ybound, r'\bf IceCube Preliminary', color='red', fontsize=13,
-                 ha='center', va='center', zorder=11)
-    elif args.data is DataType.REALISATION:
-        fig.text(0.278, ybound-0.05, r'\bf IceCube Simulation', color='red', fontsize=13,
-                 ha='center', va='center', zorder=11)
-    else:
-        fig.text(0.278, ybound, r'\bf IceCube Simulation', color='red', fontsize=13,
-                 ha='center', va='center', zorder=11)
+    # ybound = 0.61
+    # if args.data is DataType.REAL:
+    #     fig.text(0.278, ybound, r'\bf IceCube Preliminary', color='red', fontsize=13,
+    #              ha='center', va='center', zorder=11)
+    # elif args.data is DataType.REALISATION:
+    #     fig.text(0.278, ybound-0.05, r'\bf IceCube Simulation', color='red', fontsize=13,
+    #              ha='center', va='center', zorder=11)
+    # else:
+    #     fig.text(0.278, ybound, r'\bf IceCube Simulation', color='red', fontsize=13,
+    #              ha='center', va='center', zorder=11)
 
     make_dir(outfile)
     for of in outformat: