Wed May 23 16:23:12 CDT 2018

11 files changed, 273 insertions, 916 deletions

diff --git a/bout/plot.py b/bout/plot.py
deleted file mode 100644
index 717cf81..0000000
--- a/bout/plot.py
+++ /dev/null
@@ -1,129 +0,0 @@
-import os
-
-import numpy as np
-import numpy.ma as ma
-
-import matplotlib as mpl
-mpl.use('Agg')
-from matplotlib import pyplot as plt
-from matplotlib.offsetbox import AnchoredText
-from matplotlib import rc
-
-rc('text', usetex=False)
-rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
-
-fix_sfr_mfr = [
-    (1, 1, 1, 1, 2, 0),
-    # (1, 1, 1, 1, 0, 0),
-    (1, 1, 1, 0, 1, 0),
-]
-
-# FR
-# dimension         = [3, 6]
-dimension         = [3, 6]
-sigma_ratio       = ['0.01']
-energy_dependance = 'spectral'
-spectral_index    = -2
-binning           = [1e4, 1e7, 5]
-fix_mixing        = 'False'
-fix_mixing_almost = 'False'
-scale_region      = "1E10"
-
-# Likelihood
-likelihood = 'golemfit'
-confidence = 2.71 # 90% for 1DOF
-outformat = ['png']
-
-
-def gen_identifier(measured_ratio, source_ratio, dimension, sigma_ratio=0.01):
-    mr = np.array(measured_ratio) / float(np.sum(measured_ratio))
-    sr = np.array(source_ratio) / float(np.sum(source_ratio))
-    si = sigma_ratio
-    out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_single_scale'.format(
-        int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000),
-        int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), dimension
-    )
-    return out
-
-
-def get_units(dimension):
-    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}'
-
-
-def myround(x, base=5, 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))
-    else: int(base * np.round(float(x)/base))
-
-
-colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
-
-for i_dim, dim in enumerate(dimension):
-    fig = plt.figure(figsize=(7, 5))
-    ax = fig.add_subplot(111)
-    yranges = [np.inf, -np.inf]
-    legend_handles = []
-    xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$']
-    ax.set_xlim(0, len(xticks)+1)
-    ax.set_xticklabels([''] + xticks + [''])
-    ax.set_xlabel(r'BSM operator angle')
-    ylabel = r'${\rm log}_{10} \Lambda' + get_units(dim) + r'$'
-    ax.set_ylabel(ylabel)
-    for i_frs, frs in enumerate(fix_sfr_mfr):
-        print '== DIM{0}'.format(dim)
-        print '== FRS = {0}'.format(frs)
-        outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}/fix_ifr/0_01/'.format(likelihood, dim, spectral_index)
-        infile = outchain_head + '/angles_limit/fr_anfr_evidence'+ gen_identifier(frs[:3], frs[-3:], dim) + '.npy'
-        try:
-            array = np.load(infile)
-        except IOError:
-            print 'failed to open {0}'.format(infile)
-            continue
-        print 'array', array
-        print 'array', array.shape
-        for i_th in xrange(len(xticks)):
-            scale, llhs = array[i_th].T
-            min_llh = np.min(llhs)
-            delta_llh = 2*(llhs - min_llh)
-            print 'scale', scale
-            print 'delta_llh', delta_llh
-            al = scale[delta_llh < confidence]
-            if len(al) > 0:
-                label = '[{0}, {1}, {2}]'.format(frs[3], frs[4], frs[5])
-                lim = al[0]
-                print 'frs, dim, lim = ', frs, dim, lim
-                if lim < yranges[0]: yranges[0] = lim
-                if lim > yranges[1]: yranges[1] = lim+4
-                line = plt.Line2D(
-                    (i_th+1-0.1, i_th+1+0.1), (lim, lim), lw=3, color=colour[i_frs], label=label
-                )
-                ax.add_line(line)
-                if i_th == 0: legend_handles.append(line)
-                x_offset = i_frs*0.05 - 0.05
-                ax.annotate(
-                    s='', xy=(i_th+1+x_offset, lim), xytext=(i_th+1+x_offset, lim+3),
-                    arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[i_frs]}
-                )
-            else:
-                print 'No points for DIM {0} FRS {1} NULL {2}!'.format(dim, frs, min_llh)
-    try:
-        yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True))
-        # ax.set_ylim(yranges)
-        ax.set_ylim([-30, -20])
-    except: pass
-
-    ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right',
-              title='dimension {0}'.format(dim))
-    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():
-        ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1)
-
-    for of in outformat:
-        fig.savefig('../images/freq/lim_DIM{0}.'.format(dim)+of, bbox_inches='tight', dpi=150)
diff --git a/bout/plot_corr.py b/bout/plot_corr.py
deleted file mode 100644
index a75fe8a..0000000
--- a/bout/plot_corr.py
+++ /dev/null
@@ -1,193 +0,0 @@
-import os
-
-import numpy as np
-import numpy.ma as ma
-
-import scipy.interpolate as interpolate
-
-import matplotlib as mpl
-mpl.use('Agg')
-from matplotlib import pyplot as plt
-from matplotlib.offsetbox import AnchoredText
-from matplotlib import rc
-
-rc('text', usetex=False)
-rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
-
-fix_sfr_mfr = [
-    (1, 1, 1, 1, 2, 0),
-    # (1, 1, 1, 1, 0, 0),
-    (1, 1, 1, 0, 1, 0),
-]
-
-# FR
-# dimension         = [3, 6]
-dimension         = [3, 6]
-sigma_ratio       = ['0.01']
-energy_dependance = 'spectral'
-spectral_index    = -2
-binning           = [1e4, 1e7, 5]
-fix_mixing        = 'False'
-fix_mixing_almost = 'False'
-scale_region      = "1E10"
-
-# Likelihood
-likelihood = 'golemfit'
-confidence = 4.61 # 90% for 2DOF
-outformat = ['png']
-
-
-def gen_identifier(measured_ratio, source_ratio, dimension, sigma_ratio=0.01):
-    mr = np.array(measured_ratio) / float(np.sum(measured_ratio))
-    sr = np.array(source_ratio) / float(np.sum(source_ratio))
-    si = sigma_ratio
-    out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_single_scale'.format(
-        int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000),
-        int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), dimension
-    )
-    return out
-
-
-def get_units(dimension):
-    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}'
-
-
-def myround(x, base=5, 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))
-    else: int(base * np.round(float(x)/base))
-
-
-colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
-
-labels = [r'$sin^2(2\mathcal{O}_{12})$', r'$sin^2(2\mathcal{O}_{13})$', r'$sin^2(2\mathcal{O}_{23})$']
-for i_dim, dim in enumerate(dimension):
-    for i_frs, frs in enumerate(fix_sfr_mfr):
-        print '== DIM{0}'.format(dim)
-        print '== FRS = {0}'.format(frs)
-        outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}/fix_ifr/0_01/'.format(likelihood, dim, spectral_index)
-        infile = outchain_head + '/angles_corr/fr_co_evidence'+ gen_identifier(frs[:3], frs[-3:], dim) + '.npy'
-        # infile = '../mnrun/fr_co_evidence_033_033_033_0010_sfr_033_066_000_DIM6_single_scale.npy'
-        try:
-            array = ma.masked_invalid(np.load(infile))
-        except IOError:
-            print 'failed to open {0}'.format(infile)
-            continue
-        print 'array', array
-        print 'array', array.shape
-        for i_scen in xrange(len(labels)):
-            d = array[i_scen].reshape(len(array[i_scen])**2, 3)
-            fig = plt.figure(figsize=(7, 5))
-            ax = fig.add_subplot(111)
-            xranges = [np.inf, -np.inf]
-            legend_handles = []
-            ax.set_ylim(0, 1)
-            ax.set_ylabel(labels[i_scen])
-            xlabel = r'${\rm log}_{10} \Lambda' + get_units(dim) + r'$'
-            ax.set_xlabel(xlabel)
-
-            x = d[:,0]
-            y = d[:,1]
-            z = d[:,2]
-
-            print 'x', x
-            print 'y', y
-            print 'z', z
-            null_idx = np.argmin(z)
-            null = z[null_idx]
-            print 'null = {0}, for scale = {1}'.format(null, x[null_idx])
-            z = 2*(z - null)
-            print 'scale', x
-            print 'delta_llh', z
-
-            # x_ = np.linspace(np.min(x), np.max(x), 30)
-            # y_ = np.linspace(np.min(y), np.max(y), 30)
-            # z_ = interpolate.gridddata((x, y), z, (x_[None,:], y_[:,None]), method='nearest')
-
-            data = np.array([x, y, z, np.ones(x.shape)]).T
-            print 'data', data
-            data_clean = []
-            for d in data:
-                if not np.any(np.isnan(d)): data_clean.append(d)
-            data = np.vstack(data_clean)
-            sort_column = 3
-            data_sorted = data[data[:,sort_column].argsort()]
-            uni, c = np.unique(data[:,sort_column], return_counts=True)
-            print uni, c
-            print len(uni)
-            print np.unique(c)
-            
-            n = len(uni)
-            assert len(np.unique(c)) == 1
-            c = c[0]
-            col_array = []
-            for col in data_sorted.T:
-                col_array.append(col.reshape(n, c))
-            col_array = np.stack(col_array)
-            sep_arrays = []
-            for x_i in xrange(n):
-                sep_arrays.append(col_array[:,x_i])
-            
-            print len(sep_arrays)
-            sep_arrays = sep_arrays[0][:3]
-            print sep_arrays
-
-            llh_90_percent = (sep_arrays[2] < confidence)
-            data_90_percent = sep_arrays.T[llh_90_percent].T
-            print 'data_90_percent', data_90_percent
-
-            ax.tick_params(axis='x', labelsize=11)
-            ax.tick_params(axis='y', labelsize=11)
-
-            mini, maxi = np.min(x), np.max(x)
-            ax.set_xlim((mini, maxi))
-            ax.set_ylim(0, 1)
-            ax.grid(b=False)
-            
-            x_v = data_90_percent[0].round(decimals=4)
-            y_v = data_90_percent[1].round(decimals=4)
-            uniques = np.unique(x_v)
-            print 'uniques', uniques
-            if len(uniques) == 1: continue
-            bw = np.min(np.diff(uniques))
-            print 'bw', bw
-            print np.diff(uniques)
-            uni_x_split = np.split(uniques, np.where(np.diff(uniques) > bw*(1e20))[0] + 1)
-            print 'len(uni_x_split)', len(uni_x_split)
-            for uni_x in uni_x_split:
-                p_x_l, p_y_l = [], []
-                p_x_u, p_y_u = [], []
-                for uni in uni_x:
-                    idxes = np.where(x_v == uni)[0]
-                    ymin, ymax = 1, 0
-                    for idx in idxes:
-                        if y_v[idx] < ymin: ymin = y_v[idx]
-                        if y_v[idx] > ymax: ymax = y_v[idx]
-                    p_x_l.append(uni)
-                    p_y_l.append(ymin)
-                    p_x_u.append(uni)
-                    p_y_u.append(ymax)
-                p_x_l, p_y_l = np.array(p_x_l, dtype=np.float64), np.array(p_y_l, dtype=np.float64)
-                p_x_u, p_y_u = np.array(list(reversed(p_x_u)), dtype=np.float64), np.array(list(reversed(p_y_u)), dtype=np.float64)
-                p_x = np.hstack([p_x_l, p_x_u])
-                p_y = np.hstack([p_y_l, p_y_u])
-                p_x = np.r_[p_x, p_x[0]]
-                p_y = np.r_[p_y, p_y[0]]
-                print 'p_x', p_x
-                print 'p_y', p_y
-                try:
-                    tck, u = interpolate.splprep([p_x, p_y], s=1e-5, per=True)
-                    xi, yi = interpolate.splev(np.linspace(0, 1, 1000), tck)
-                except:
-                    xi, yi = p_x, p_y
-                ax.fill(xi, yi, 'r', edgecolor='r', linewidth=1)
-            
-            for of in outformat:
-                plt.savefig('../images/freq/lim_corr_DIM{0}_AN{1}_FRS{2}'.format(dim, i_scen, i_frs)+of, bbox_inches='tight', dpi=150)
-
diff --git a/bout/plot_full.py b/bout/plot_full.py
deleted file mode 100644
index f2e1919..0000000
--- a/bout/plot_full.py
+++ /dev/null
@@ -1,127 +0,0 @@
-import os
-
-import numpy as np
-import numpy.ma as ma
-
-import matplotlib as mpl
-mpl.use('Agg')
-from matplotlib import pyplot as plt
-from matplotlib.offsetbox import AnchoredText
-from matplotlib import rc
-
-rc('text', usetex=False)
-rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
-
-fix_sfr_mfr = [
-    (1, 1, 1, 1, 2, 0),
-    # (1, 1, 1, 1, 0, 0),
-    (1, 1, 1, 0, 1, 0),
-]
-
-# FR
-# dimension         = [3, 6]
-dimension         = [3, 4, 5, 6, 7, 8]
-sigma_ratio       = ['0.01']
-energy_dependance = 'spectral'
-spectral_index    = -2
-binning           = [1e4, 1e7, 5]
-fix_mixing        = 'False'
-fix_mixing_almost = 'False'
-scale_region      = "1E10"
-
-# Likelihood
-likelihood = 'golemfit'
-confidence = 2.71 # 90% for 1DOF
-outformat = ['png']
-
-
-def gen_identifier(measured_ratio, source_ratio, dimension, sigma_ratio=0.01):
-    mr = np.array(measured_ratio) / float(np.sum(measured_ratio))
-    sr = np.array(source_ratio) / float(np.sum(source_ratio))
-    si = sigma_ratio
-    out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_single_scale'.format(
-        int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000),
-        int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), dimension
-    )
-    return out
-
-
-def get_units(dimension):
-    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}'
-
-
-def myround(x, base=5, 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))
-    else: int(base * np.round(float(x)/base))
-
-
-colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
-
-fig = plt.figure(figsize=(7, 5))
-ax = fig.add_subplot(111)
-
-colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
-yranges = [np.inf, -np.inf]
-legend_handles = []
-ax.set_xlim(dimension[0]-1, dimension[-1]+1)
-xticks = [''] + range(dimension[0], dimension[-1]+1) + ['']
-ax.set_xticklabels(xticks)
-ax.set_xlabel(r'BSM operator dimension ' + r'$d$')
-ax.set_ylabel(r'${\rm log}_{10} \Lambda / GeV^{-d+4}$')
-for i_dim, dim in enumerate(dimension):
-    for i_frs, frs in enumerate(fix_sfr_mfr):
-        outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}/fix_ifr/0_01/'.format(likelihood, dim, spectral_index)
-        infile = outchain_head + '/bayes_factor/fr_fr_evidence' + gen_identifier(frs[:3], frs[-3:], dim) + '.npy'
-        try:
-            array = np.load(infile)
-        except IOError:
-            print 'failed to open {0}'.format(infile)
-            continue
-        print 'array', array
-        print 'array', array.shape
-        scale, llhs = array.T
-        print 'scale min', scale[np.argmin(llhs)]
-        null = llhs[np.argmin(llhs)]
-        # null = llhs[0]
-        # TODO(shivesh): negative or not?
-        reduced_ev = 2*(llhs - null)
-        print 'reduced_ev', reduced_ev
-        al = scale[reduced_ev < confidence]
-        if len(al) > 0:
-            label = '[{0}, {1}, {2}]'.format(frs[3], frs[4], frs[5])
-            lim = al[0]
-            print 'frs, dim, lim = ', frs, dim, lim
-            if lim < yranges[0]: yranges[0] = lim
-            if lim > yranges[1]: yranges[1] = lim+4
-            line = plt.Line2D(
-                (dim-0.1, dim+0.1), (lim, lim), lw=3, color=colour[i_frs], label=label
-            )
-            ax.add_line(line)
-            if i_dim == 0: legend_handles.append(line)
-            x_offset = i_frs*0.05 - 0.05
-            ax.annotate(
-                s='', xy=(dim+x_offset, lim), xytext=(dim+x_offset, lim+3),
-                arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[i_frs]}
-            )
-
-        else:
-            print 'No points for DIM {0} FRS {1} NULL {2}!'.format(dim, frs, null)
-            # print 'scales, reduced_ev', np.dstack([scale.data, reduced_ev.data])
-yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True))
-ax.set_ylim(yranges)
-
-ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right')
-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():
-    ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1)
-
-for of in outformat:
-    fig.savefig('../images/freq/full_corr.'+of, bbox_inches='tight', dpi=150)
diff --git a/nufit_u.txt b/misc/nufit_u.txt
index a9c9cb2..a9c9cb2 100644
--- a/nufit_u.txt
+++ b/misc/nufit_u.txt
diff --git a/plot_bf.py b/plot_bf.py
deleted file mode 100755
index 13664b9..0000000
--- a/plot_bf.py
+++ /dev/null
@@ -1,414 +0,0 @@
-#! /usr/bin/env python
-# author : S. Mandalia
-#          s.p.mandalia@qmul.ac.uk
-#
-# date   : April 14, 2018
-
-"""
-HESE BSM flavour ratio sensivity plotting script
-"""
-
-from __future__ import absolute_import, division
-
-import os
-
-import numpy as np
-import numpy.ma as ma
-
-import matplotlib as mpl
-mpl.use('Agg')
-from matplotlib import rc
-from matplotlib import pyplot as plt
-from matplotlib.offsetbox import AnchoredText
-
-import fr
-from utils import misc as misc_utils
-from utils.fr import normalise_fr
-from utils.plot import bayes_factor_plot, myround, get_units
-
-
-rc('text', usetex=False)
-rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
-
-fix_sfr_mfr = [
-    (1, 1, 1, 1, 2, 0),
-    # (1, 1, 1, 1, 0, 0),
-    (1, 1, 1, 0, 1, 0),
-]
-
-# FR
-dimension         = [3, 6]
-# dimension         = [3, 4, 5, 6, 7, 8]
-sigma_ratio       = ['0.01']
-energy_dependance = 'spectral'
-spectral_index    = -2
-binning           = [1e4, 1e7, 5]
-fix_mixing        = 'False'
-fix_mixing_almost = 'False'
-scale_region      = "1E10"
-
-# Likelihood
-likelihood = 'golemfit'
-
-# Nuisance
-convNorm        = 1.
-promptNorm      = 0.
-muonNorm        = 1.
-astroNorm       = 6.9
-astroDeltaGamma = 2.5
-
-# GolemFit
-ast  = 'p2_0'
-data = 'real'
-
-# Bayes Factor
-bayes_bins        = 100
-bayes_live_points = 1000
-bayes_tolerance   = 0.01
-bayes_eval_bin    = True # set to 'all' to run normally
-
-# Plot
-plot_bayes        = False
-plot_angles_limit = True
-plot_angles_corr  = False
-outformat         = ['png']
-# significance      = np.log(10**(1/2.))
-significance      = np.log(10**(3/2.))
-
-
-bayes_array = ma.masked_equal(np.zeros((len(dimension), len(fix_sfr_mfr), bayes_bins, 2)), 0)
-angles_lim_array = np.zeros((len(dimension), len(fix_sfr_mfr), 3, bayes_bins, 2))
-angles_corr_array = np.zeros((len(dimension), len(fix_sfr_mfr), 3, bayes_bins, bayes_bins, 3))
-for i_dim, dim in enumerate(dimension):
-    if energy_dependance == 'mono':
-        outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/{2:.0E}'.format(likelihood, dim, en)
-    elif energy_dependance == 'spectral':
-        outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format(likelihood, dim, spectral_index)
-
-    bayes_output = 'None'
-    angles_lim_output = 'None'
-    angles_corr_output = 'None'
-    for sig in sigma_ratio:
-        for i_frs, frs in enumerate(fix_sfr_mfr):
-            outchains = outchain_head + '/fix_ifr/{0}/'.format(str(sig).replace('.', '_'))
-            if plot_bayes:
-                bayes_output = outchains + '/bayes_factor/'
-            if plot_angles_limit:
-                angles_lim_output = outchains + '/angles_limit/'
-            if plot_angles_corr:
-                angles_corr_output = outchains + '/angles_corr/'
-
-            argstring = '--measured-ratio {0} {1} {2} --fix-source-ratio True --source-ratio {3} {4} {5} --dimension {6} --seed 24 --outfile {7} --run-mcmc False --likelihood {8} --plot-angles False --bayes-output {9} --angles-lim-output {10} --bayes-bins {11} --angles-corr-output {12}'.format(frs[0], frs[1], frs[2], frs[3], frs[4], frs[5], dim, outchains, likelihood, bayes_output, angles_lim_output, bayes_bins, angles_corr_output)
-            args = fr.parse_args(argstring)
-            fr.process_args(args)
-            # misc_utils.print_args(args)
-
-            if plot_bayes:
-                infile = args.bayes_output+'/fr_evidence'+misc_utils.gen_identifier(args)
-            if plot_angles_limit:
-                infile = args.angles_lim_output+'/fr_an_evidence'+misc_utils.gen_identifier(args)
-            if plot_angles_corr:
-                infile = args.angles_corr_output+'/fr_co_evidence' + misc_utils.gen_identifier(args)
-            scan_scales = np.linspace(
-                np.log10(args.scale_region[0]), np.log10(args.scale_region[1]), args.bayes_bins
-            )
-            print 'scan_scales', scan_scales
-            raw = []
-            fail = 0
-            if plot_angles_corr:
-                scan_angles = np.linspace(0, 1, args.bayes_bins)
-            for i_sc, sc in enumerate(scan_scales):
-                if plot_angles_corr:
-                    for i_an, an in enumerate(scan_angles):
-                        idx = i_sc*args.bayes_bins + i_an
-                        infile_s = infile + '_idx_{0}'.format(idx)
-                        try:
-                            lf = np.load(infile_s+'.npy')
-                            print 'lf.shape', lf.shape
-                        except IOError:
-                            fail += 1
-                            print 'failed to open {0}'.format(infile_s)
-                            lf = np.full((3, 1, 1, 3), np.nan)
-                            pass
-                        for x in xrange(len(lf)):
-                            angles_corr_array[i_dim][i_frs][x][i_sc][i_an] = np.array(lf[x])
-                    continue
-                try:
-                    infile_s = infile + '_scale_{0:.0E}'.format(np.power(10, sc))
-                    lf = np.load(infile_s+'.npy')
-                    print lf.shape
-                    if plot_angles_limit:
-                        if len(lf.shape) == 3: lf = lf[:,0,:]
-                    raw.append(lf)
-                except IOError:
-                    fail += 1
-                    print 'failed to open {0}'.format(infile_s)
-                    if plot_bayes:
-                        raw.append([0, 0])
-                    if plot_angles_limit:
-                        raw.append(np.zeros((3, 2)))
-                    pass
-            print 'failed to open {0} files'.format(fail)
-
-            if plot_bayes:
-                raw = np.vstack(raw)
-            if plot_angles_limit:
-                raw = np.vstack(raw).reshape(args.bayes_bins, 3, 2)
-                a = ma.masked_equal(np.zeros((3, args.bayes_bins, 2)), 0)
-                for i_x, x in enumerate(raw):
-                    for i_y, y in enumerate(x):
-                        a[i_y][i_x] = ma.masked_equal(y, 0)
-            if plot_angles_corr:
-                a = angles_corr_array[i_dim][i_frs]
-                a = ma.masked_invalid(a, 0)
-                # for i_sc in xrange(len(scan_scales)):
-                #     for i_a in xrange(len(scan_angles)):
-                #         try:
-                #             bayes_factor_plot(
-                #                 a[i_sc,:,i_a,:][:,(0,2)], './mnrun/corr/test_corr_DIM{0}_FR{1}_AN{2}_SC{3}'.format(dim, i_frs, i_a, i_sc), ['png'], args
-                #             )
-                #         except: pass
-
-            if plot_bayes:
-                bayes_array[i_dim][i_frs] = ma.masked_equal(raw, 0)
-                bayes_factor_plot( bayes_array[i_dim][i_frs], './mnrun/test_full_DIM{0}_FR{1}'.format(dim, i_frs), ['png'], args
-                )
-
-            if plot_angles_limit:
-                angles_lim_array[i_dim][i_frs] = ma.masked_equal(a, 0)
-                for i_a, angle in enumerate(a):
-                    bayes_factor_plot(
-                        angle, './mnrun/test_angles_DIM{0}_FR{1}_AN{2}'.format(i_dim, i_frs, i_a), ['png'], args
-                    )
-
-if plot_bayes:
-    fig = plt.figure(figsize=(7, 5))
-    ax = fig.add_subplot(111)
-
-    colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
-    yranges = [np.inf, -np.inf]
-    legend_handles = []
-    ax.set_xlim(dimension[0]-1, dimension[-1]+1)
-    xticks = [''] + range(dimension[0], dimension[-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}$')
-    for i_dim, dim in enumerate(dimension):
-        for i_frs, frs in enumerate(fix_sfr_mfr):
-            scale, evidences = bayes_array[i_dim][i_frs].T
-            null = evidences[np.argmin(scale)]
-            # TODO(shivesh): negative or not?
-            reduced_ev = -(evidences - null)
-            al = scale[reduced_ev > significance]
-            if len(al) > 0:
-                label = '[{0}, {1}, {2}]'.format(frs[3], frs[4], frs[5])
-                lim = al[0]
-                print 'frs, dim, lim = ', frs, dim, lim
-                if lim < yranges[0]: yranges[0] = lim
-                if lim > yranges[1]: yranges[1] = lim+4
-                line = plt.Line2D(
-                    (dim-0.1, dim+0.1), (lim, lim), lw=3, color=colour[i_frs], label=label
-                )
-                ax.add_line(line)
-                if i_dim == 0: legend_handles.append(line)
-                x_offset = i_frs*0.05 - 0.05
-                ax.annotate(
-                    s='', xy=(dim+x_offset, lim), xytext=(dim+x_offset, lim+3),
-                    arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[i_frs]}
-                )
-
-            else:
-                print 'No points for DIM {0} FRS {1} NULL {2}!'.format(dim, frs, null)
-                # print 'scales, reduced_ev', np.dstack([scale.data, reduced_ev.data])
-    try:
-        yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True))
-        ax.set_ylim(yranges)
-    except: pass
-
-    ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right')
-    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():
-        ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1)
-
-    for of in outformat:
-        fig.savefig('./images/bayes/bayes_factor.'+of, bbox_inches='tight', dpi=150)
-
-if plot_angles_limit:
-    colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
-    for i_dim, dim in enumerate(dimension):
-        fig = plt.figure(figsize=(7, 5))
-        ax = fig.add_subplot(111)
-        yranges = [np.inf, -np.inf]
-        legend_handles = []
-        xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$']
-        ax.set_xlim(0, len(xticks)+1)
-        ax.set_xticklabels([''] + xticks + [''])
-        ax.set_xlabel(r'BSM operator angle')
-        ylabel = r'${\rm log}_{10} \Lambda^{-1}' + get_units(dim) + r'$'
-        ax.set_ylabel(ylabel)
-        for i_th in xrange(len(xticks)):
-            for i_frs, frs in enumerate(fix_sfr_mfr):
-                scale, evidences = angles_lim_array[i_dim][i_frs][i_th].T
-                null = evidences[np.argmin(scale)]
-                # TODO(shivesh): negative or not?
-                reduced_ev = -(evidences - null)
-                al = scale[reduced_ev > significance]
-                # print 'scales, reduced_ev', np.dstack([scale, reduced_ev])
-                if len(al) > 0:
-                    label = '[{0}, {1}, {2}]'.format(frs[3], frs[4], frs[5])
-                    lim = al[0]
-                    print 'frs, dim, lim = ', frs, dim, lim
-                    if lim < yranges[0]: yranges[0] = lim
-                    if lim > yranges[1]: yranges[1] = lim+4
-                    line = plt.Line2D(
-                        (i_th+1-0.1, i_th+1+0.1), (lim, lim), lw=3, color=colour[i_frs], label=label
-                    )
-                    ax.add_line(line)
-                    if i_th == 0: legend_handles.append(line)
-                    x_offset = i_frs*0.05 - 0.05
-                    ax.annotate(
-                        s='', xy=(i_th+1+x_offset, lim), xytext=(i_th+1+x_offset, lim+3),
-                        arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[i_frs]}
-                    )
-        try:
-            yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True))
-            ax.set_ylim(yranges)
-        except: pass
-
-        ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right',
-                  title='dimension {0}'.format(dim))
-
-        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():
-            ax.axvline(x=xmaj, ls='-', color='gray', alpha=0.4, linewidth=1)
-
-        for of in outformat:
-            fig.savefig('./images/bayes/angles_limit_DIM{0}'.format(dim)+'.'+of, bbox_inches='tight', dpi=150)
-
-if plot_angles_corr:
-    colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
-    labels = [r'$sin^2(2\mathcal{O}_{12})$', r'$sin^2(2\mathcal{O}_{13})$', r'$sin^2(2\mathcal{O}_{23})$']
-    for i_dim, dim in enumerate(dimension):
-        for i_frs, frs in enumerate(fix_sfr_mfr):
-            print '== DIM{0}'.format(dim)
-            print '== FRS = {0}'.format(frs)
-            array = angles_corr_array[i_dim][i_frs]
-            print 'array', array
-            print 'array.shape', array.shape
-            for i_scen in xrange(len(labels)):
-                d = array[i_scen].reshape(len(array[i_scen])**2, 3)
-                print 'd.shape', d.shape
-                fig = plt.figure(figsize=(7, 5))
-                ax = fig.add_subplot(111)
-                xranges = [np.inf, -np.inf]
-                legend_handles = []
-                ax.set_ylim(0, 1)
-                ax.set_ylabel(labels[i_scen])
-                xlabel = r'${\rm log}_{10} \Lambda^{-1}' + get_units(dim) + r'$'
-                ax.set_xlabel(xlabel)
-
-                x = d[:,0]
-                y = d[:,1]
-                z = d[:,2]
-
-                data_clean = []
-                for id in d:
-                    if not np.any(np.isnan(id)): data_clean.append(id)
-                d_c = np.vstack(data_clean)
-
-                x = d_c[:,0]
-                y = d_c[:,1]
-                z = d_c[:,2]
-
-                # print 'x', x
-                # print 'y', y
-                # print 'z', z
-                null_idx = np.argmin(x)
-                null = z[null_idx]
-                print 'null = {0}, for scale = {1}'.format(null, x[null_idx])
-                z = -(z - null)
-                print 'scale', x
-                print 'bayes_factor', z
-
-                # x_ = np.linspace(np.min(x), np.max(x), 30)
-                # y_ = np.linspace(np.min(y), np.max(y), 30)
-                # z_ = interpolate.gridddata((x, y), z, (x_[None,:], y_[:,None]), method='nearest')
-
-                data = np.array([x, y, z, np.ones(x.shape)]).T
-                sort_column = 3
-                data_sorted = data[data[:,sort_column].argsort()]
-                uni, c = np.unique(data[:,sort_column], return_counts=True)
-                print uni, c
-                print len(uni)
-                print np.unique(c)
-                
-                n = len(uni)
-                assert len(np.unique(c)) == 1
-                c = c[0]
-                col_array = []
-                for col in data_sorted.T:
-                    col_array.append(col.reshape(n, c))
-                col_array = np.stack(col_array)
-                sep_arrays = []
-                for x_i in xrange(n):
-                    sep_arrays.append(col_array[:,x_i])
-                
-                print len(sep_arrays)
-                sep_arrays = sep_arrays[0][:3]
-                print sep_arrays
-
-                allowed_bf = (sep_arrays[2] < significance) # Shade the excluded region
-                data_allowed_bf = sep_arrays.T[allowed_bf].T
-                print 'data_allowed_bf', data_allowed_bf
-
-                ax.tick_params(axis='x', labelsize=11)
-                ax.tick_params(axis='y', labelsize=11)
-
-                mini, maxi = np.min(scan_scales), np.max(scan_scales)
-                ax.set_xlim((mini, maxi))
-                ax.set_ylim(0, 1)
-                ax.grid(b=False)
-                
-                x_v = data_allowed_bf[0].round(decimals=4)
-                y_v = data_allowed_bf[1].round(decimals=4)
-                uniques = np.unique(x_v)
-                # print 'uniques', uniques
-                if len(uniques) == 1: continue
-                bw = np.min(np.diff(uniques))
-                # print 'bw', bw
-                print np.diff(uniques)
-                uni_x_split = np.split(uniques, np.where(np.diff(uniques) > bw*(1e20))[0] + 1)
-                # print 'len(uni_x_split)', len(uni_x_split)
-                for uni_x in uni_x_split:
-                    p_x_l, p_y_l = [], []
-                    p_x_u, p_y_u = [], []
-                    for uni in uni_x:
-                        idxes = np.where(x_v == uni)[0]
-                        ymin, ymax = 1, 0
-                        for idx in idxes:
-                            if y_v[idx] < ymin: ymin = y_v[idx]
-                            if y_v[idx] > ymax: ymax = y_v[idx]
-                        p_x_l.append(uni)
-                        p_y_l.append(ymin)
-                        p_x_u.append(uni)
-                        p_y_u.append(ymax)
-                    p_x_l, p_y_l = np.array(p_x_l, dtype=np.float64), np.array(p_y_l, dtype=np.float64)
-                    p_x_u, p_y_u = np.array(list(reversed(p_x_u)), dtype=np.float64), np.array(list(reversed(p_y_u)), dtype=np.float64)
-                    p_x = np.hstack([p_x_l, p_x_u])
-                    p_y = np.hstack([p_y_l, p_y_u])
-                    p_x = np.r_[p_x, p_x[0]]
-                    p_y = np.r_[p_y, p_y[0]]
-                    print 'p_x', p_x
-                    print 'p_y', p_y
-                    try:
-                        tck, u = interpolate.splprep([p_x, p_y], s=1e-5, per=True)
-                        xi, yi = interpolate.splev(np.linspace(0, 1, 1000), tck)
-                    except:
-                        xi, yi = p_x, p_y
-                    ax.fill(xi, yi, 'r', edgecolor='r', linewidth=1)
-                
-                for of in outformat:
-                    plt.savefig('./images/bayes/lim_corr_DIM{0}_AN{1}_FRS{2}'.format(dim, i_scen, i_frs)+of, bbox_inches='tight', dpi=150)
-
diff --git a/plot_sens.py b/plot_sens.py
index 6b82ebc..a9704fe 100755
--- a/plot_sens.py
+++ b/plot_sens.py
@@ -349,6 +349,12 @@ def main():
                 args      = args,
             )
         elif args.run_method in fixed_angle_categ:
+            plot_utils.plot_sens_fixed_angle_pretty(
+                data      = data,
+                outfile   = baseoutfile + '/fixed_angle_pretty',
+                outformat = ['png'],
+                args      = args,
+            )
             plot_utils.plot_sens_fixed_angle(
                 data      = data,
                 outfile   = baseoutfile + '/FIXED_ANGLE',
diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py
index 78034d2..6734278 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    = 'corr_angle', # full, fixed_angle, corr_angle
+    run_method    = 'fixed_angle', # full, fixed_angle, corr_angle
     stat_method   = 'bayesian',
-    sens_bins     = 10,
+    sens_bins     = 20,
     seed          = 'None'
 ))
 
@@ -65,7 +65,7 @@ GLOBAL_PARAMS.update(dict(
 # GolemFit
 GLOBAL_PARAMS.update(dict(
     ast  = 'p2_0',
-    data = 'real'
+    data = 'realisation'
 ))
 
 # Plot
diff --git a/utils/enums.py b/utils/enums.py
index c40d681..7fcddae 100644
--- a/utils/enums.py
+++ b/utils/enums.py
@@ -11,8 +11,9 @@ from enum import Enum
 
 
 class DataType(Enum):
-    REAL   = 1
-    ASIMOV = 2
+    REAL        = 1
+    ASIMOV      = 2
+    REALISATION = 3
 
 
 class EnergyDependance(Enum):
diff --git a/utils/gf.py b/utils/gf.py
index 06a6125..0d098f1 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -71,7 +71,7 @@ def steering_params(args):
     return params
 
 
-def set_up_as(fitter, params):
+def setup_asimov(fitter, params):
     print 'Injecting the model', params
     asimov_params = gf.FitParameters(gf.sampleTag.MagicTau)
     for parm in params:
@@ -79,13 +79,23 @@ def set_up_as(fitter, params):
     fitter.SetupAsimov(asimov_params)
 
 
+def setup_realisation(fitter, params):
+    print 'Injecting the model', params
+    asimov_params = gf.FitParameters(gf.sampleTag.MagicTau)
+    for parm in params:
+        asimov_params.__setattr__(parm.name, float(parm.value))
+    fitter.Swallow(fitter.SpitExpectation(asimov_params))
+
+
 def setup_fitter(args, asimov_paramset):
     datapaths = gf.DataPaths()
     sparams = steering_params(args)
     npp = gf.NewPhysicsParams()
     fitter = gf.GolemFit(datapaths, sparams, npp)
     if args.data is DataType.ASIMOV:
-        set_up_as(fitter, asimov_paramset)
+        setup_asimov(fitter, asimov_paramset)
+    elif args.data is DataType.REALISATION:
+        setup_realisation(fitter, asimov_paramset)
     elif args.data is DataType.REAL:
         print 'Using MagicTau DATA'
     return fitter
diff --git a/utils/param.py b/utils/param.py
index 5684e91..f8d64eb 100644
--- a/utils/param.py
+++ b/utils/param.py
@@ -18,7 +18,7 @@ import numpy as np
 
 from utils.plot import get_units
 from utils.fr import fr_to_angles
-from utils.enums import Likelihood, ParamTag, PriorsCateg
+from utils.enums import DataType, Likelihood, ParamTag, PriorsCateg
 
 
 class Param(object):
diff --git a/utils/plot.py b/utils/plot.py
index 4dd8cc4..350fa82 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -10,6 +10,7 @@ Plotting functions for the BSM flavour ratio analysis
 from __future__ import absolute_import, division
 
 import os
+import socket
 from copy import deepcopy
 
 import numpy as np
@@ -20,6 +21,10 @@ mpl.use('Agg')
 from matplotlib import rc
 from matplotlib import pyplot as plt
 from matplotlib.offsetbox import AnchoredText
+import matplotlib.patches as patches
+import matplotlib.gridspec as gridspec
+from matplotlib.lines import Line2D
+from matplotlib.patches import Patch
 
 from getdist import plots, mcsamples
 
@@ -28,8 +33,9 @@ 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)
-rc('font', **{'family':'serif', 'serif':['DejaVu Sans'], 'size':18})
+plt.style.use(os.environ['GOLEMSOURCEPATH']+'/GolemFit/scripts/paper/paper.mplstyle')
+if 'submitter' in socket.gethostname():
+    rc('text', usetex=False)
 
 
 def centers(x):
@@ -37,12 +43,12 @@ def centers(x):
 
 
 def get_units(dimension):
-    if dimension == 3: return r' / \:GeV'
+    if dimension == 3: return r' / \:{\rm 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' / \:{rm GeV}^{-1}'
+    if dimension == 6: return r' / \:{rm GeV}^{-2}'
+    if dimension == 7: return r' / \:{rm GeV}^{-3}'
+    if dimension == 8: return r' / \:{rm GeV}^{-4}'
 
 
 def calc_nbins(x):
@@ -120,13 +126,13 @@ def gen_figtext(args):
     if args.fix_source_ratio:
         sr1, sr2, sr3 = misc_utils.solve_ratio(args.source_ratio)
         t += 'Source flavour ratio = [{0}, {1}, {2}]'.format(sr1, sr2, sr3)
-        if args.data is DataType.ASIMOV:
+        if args.data in [DataType.ASIMOV, DataType.REALISATION]:
             t += '\nIC observed flavour ratio = [{0}, {1}, {2}]'.format(
                 mr1, mr2, mr3
             )
         t += '\nDimension = {0}'.format(args.dimension)
     else:
-        if args.data is DataType.ASIMOV:
+        if args.data in [DataType.ASIMOV, DataType.REALISATION]:
             t += 'IC observed flavour ratio = [{0}, {1}, ' \
                     '{2}]\nDimension = {3}'.format(
                         mr1, mr2, mr3, args.dimension, int(args.energy)
@@ -185,6 +191,13 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset):
         #         '1E{2}'.format(itv[0], itv[2], itv[1])
         #     )
 
+        if args.data is DataType.REAL:
+            fig.text(0.8, 0.7, 'IceCube Preliminary', color='red', fontsize=15,
+                     ha='center', va='center')
+        elif args.data in [DataType.ASIMOV, DataType.REALISATION]:
+            fig.text(0.8, 0.7, 'IceCube Simulation', color='red', fontsize=15,
+                     ha='center', va='center')
+
         for of in outformat:
             g.export(outfile+'_angles.'+of)
 
@@ -229,6 +242,13 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset):
 
         g = plot_Tchain(Tchain, trns_axes_labels, trns_ranges)
 
+        if args.data is DataType.REAL:
+            fig.text(0.8, 0.7, 'IceCube Preliminary', color='red', fontsize=15,
+                     ha='center', va='center')
+        elif args.data in [DataType.ASIMOV, DataType.REALISATION]:
+            fig.text(0.8, 0.7, 'IceCube Simulation', color='red', fontsize=15,
+                     ha='center', va='center')
+
         mpl.pyplot.figtext(0.5, 0.7, fig_text, fontsize=15)
         for of in outformat:
             g.export(outfile+'_elements.'+of)
@@ -285,7 +305,7 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
     if args.data is DataType.REAL:
         fig.text(0.8, 0.14, 'IceCube Preliminary', color='red', fontsize=9,
                  ha='center', va='center')
-    elif args.data is DataType.ASIMOV:
+    elif args.data in [DataType.ASIMOV, DataType.REALISATION]:
         fig.text(0.8, 0.14, 'IceCube Simulation', color='red', fontsize=9,
                  ha='center', va='center')
 
@@ -375,30 +395,210 @@ def plot_sens_full(data, outfile, outformat, args):
         fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150)
 
 
+def plot_sens_fixed_angle_pretty(data, outfile, outformat, args):
+    print 'Making FIXED_ANGLE_PRETTY sensitivity plot'
+    argsc = deepcopy(args)
+    dims = len(data)
+    LV_atmo_90pc_limits = {
+        3: (2E-24, 1E-1),
+        4: (2.7E-28, 3.16E-25),
+        5: (1.5E-32, 1.12E-27),
+        6: (9.1E-37, 2.82E-30),
+        7: (3.6E-41, 1.77E-32),
+        8: (1.4E-45, 1.00E-34)
+    }
+
+    show_data = True
+
+    en = np.log10([1E4, 1E7])
+    bote = {
+        3: (-21-(en[0]+en[0]*0), -21-(en[1]+en[1]*0)),
+        4: (-21-(en[0]+en[0]*1), -21-(en[1]+en[1]*1)),
+        5: (-21-(en[0]+en[0]*2), -21-(en[1]+en[1]*2)),
+        6: (-21-(en[0]+en[0]*3), -21-(en[1]+en[1]*3)),
+        7: (-21-(en[0]+en[0]*4), -21-(en[1]+en[1]*4)),
+        8: (-21-(en[0]+en[0]*5), -21-(en[1]+en[1]*5))
+    }
+
+    colour = {0:'red', 1:'blue', 2:'green'}
+    rgb_co = {0:[1,0,0], 1:[0,0,1], 2:[0.0, 0.5019607843137255, 0.0]}
+    yticks = [r'$\mathcal{O}_{e\mu}$', r'$\mathcal{O}_{e\tau}$', r'$\mathcal{O}_{\mu\tau}$']
+    xlims = (-60, -20)
+
+    fig = plt.figure(figsize=(8, 6))
+    gs = gridspec.GridSpec(dims, 1)
+    gs.update(hspace=0.15)
+
+    first_ax = None
+    legend_log = []
+    legend_elements = []
+
+    for idim in xrange(len(data)):
+        dim = args.dimensions[idim]
+        print '== dim', dim
+        argsc.dimension = dim
+        gs0 = gridspec.GridSpecFromSubplotSpec(
+            len(yticks), 1, subplot_spec=gs[idim], hspace=0
+        )
+
+        for ian in xrange(len(yticks)):
+            print '=== an', ian
+            ax = fig.add_subplot(gs0[ian])
+            if first_ax is None: first_ax = ax
+            ax.set_xlim(xlims)
+            ylim = (0.5, 1.5)
+            ax.set_ylim(ylim)
+            # ax.set_yticks([ylim[0], 1., ylim[1]])
+            # ax.set_yticklabels([''] + [yticks[ian]] + [''], fontsize=13)
+            # ax.yaxis.get_major_ticks()[0].set_visible(False)
+            # ax.yaxis.get_major_ticks()[2].set_visible(False)
+            ax.set_yticks([1.])
+            ax.set_yticklabels([yticks[ian]], fontsize=13)
+            ax.yaxis.tick_right()
+            for xmaj in ax.xaxis.get_majorticklocs():
+                ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.2, linewidth=1)
+            ax.get_xaxis().set_visible(False)
+            linestyle = (5, (10, 10))
+            if ian == 0:
+                ax.spines['bottom'].set_alpha(0.6)
+            elif ian == 1:
+                ax.text(
+                    -0.04, ylim[0], '$d = {0}$'.format(dim), fontsize=16,
+                    rotation='90', verticalalignment='center',
+                    transform=ax.transAxes
+                )
+                dim_label_flag = False
+                ax.spines['top'].set_alpha(0.6)
+                ax.spines['bottom'].set_alpha(0.6)
+            elif ian == 2:
+                ax.spines['top'].set_alpha(0.6)
+
+            for isrc in xrange(len(data[idim])):
+                src = args.source_ratios[isrc]
+                print '== src', src
+                argsc.source_ratio = src
+
+                if show_data:
+                    alpha = 0.03
+                    col = 'black'
+                else:
+                    alpha = 0.07
+                    col = 'blue'
+                ax.add_patch(patches.Rectangle(
+                    (bote[dim][1], ylim[0]), bote[dim][0]-bote[dim][1], np.diff(ylim),
+                    fill=True, facecolor=col, alpha=alpha, linewidth=0
+                ))
+
+                scales, statistic = data[idim][isrc][ian].T
+                tck, u = interpolate.splprep([scales, statistic], s=0)
+                scales, statistic = interpolate.splev(np.linspace(0, 1, 1000), tck)
+                min_idx = np.argmin(scales)
+                null = statistic[min_idx]
+                if args.stat_method is StatCateg.BAYESIAN:
+                    reduced_ev = -(statistic - null)
+                    al = scales[reduced_ev > np.log(10**(3/2.))] # Strong degree of belief
+                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
+                if len(al) == 0:
+                    print 'No points for DIM {0} FRS {1}!'.format(dim, src)
+                    continue
+                if reduced_ev[-1] < np.log(10**(3/2.)) - 0.1:
+                    print 'Peaked contour does not exclude large scales! For ' \
+                        'DIM {0} FRS{1}!'.format(dim, src)
+                    continue
+                lim = al[0]
+                print 'limit = {0}'.format(lim)
+
+                if show_data:
+                    ax.axvline(x=lim, color=colour[isrc], alpha=1., linewidth=1.5)
+                    ax.add_patch(patches.Rectangle(
+                        (lim, ylim[0]), 100, np.diff(ylim), fill=True, facecolor=colour[isrc],
+                        alpha=0.3, linewidth=0
+                    ))
+
+                    if isrc not in legend_log:
+                        legend_log.append(isrc)
+                        label = '{0} : {1} : {2} at source'.format(*misc_utils.solve_ratio(src))
+                        legend_elements.append(
+                            Patch(facecolor=rgb_co[isrc]+[0.3],
+                                  edgecolor=rgb_co[isrc]+[1], label=label)
+                        )
+
+            if ian == 2:
+                LV_lim = np.log10(LV_atmo_90pc_limits[dim])
+                ax.add_patch(patches.Rectangle(
+                    (LV_lim[1], ylim[0]), LV_lim[0]-LV_lim[1], np.diff(ylim),
+                    fill=False, hatch='\\\\'
+                ))
+
+    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)
+
+    legend_elements.append(
+        Patch(facecolor=col, alpha=alpha+0.1, edgecolor='k', label='max sensitivity')
+    )
+    legend_elements.append(
+        Patch(facecolor='none', hatch='\\\\', edgecolor='k', label='IceCube arXiv:1709.03434')
+    )
+
+    legend = first_ax.legend(
+        handles=legend_elements, prop=dict(size=11), loc='upper left',
+        title='Excluded regions', framealpha=1., edgecolor='black',
+        frameon=True
+    )
+    first_ax.set_zorder(10)
+    plt.setp(legend.get_title(), fontsize='11')
+    legend.get_frame().set_linestyle('-')
+
+    if show_data: ybound = 0.65
+    else: ybound = 0.73
+    if args.data is DataType.REAL and show_data:
+        # fig.text(0.295, 0.684, 'IceCube Preliminary', color='red', fontsize=13,
+        fig.text(0.278, ybound, r'\bf IceCube Preliminary', 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)
+
+    misc_utils.make_dir(outfile)
+    for of in outformat:
+        print 'Saving plot as {0}'.format(outfile+'.'+of)
+        fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150)
+
+
 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}_{e\mu}$', r'$\mathcal{O}_{e\tau}$', r'$\mathcal{O}_{\mu\tau}$']
     argsc = deepcopy(args)
+
+    LV_atmo_90pc_limits = dict(np.genfromtxt('./misc/LV_atmo_90pc_limits.txt'))
+    ylims = {
+        3 : (-28, -22), 4 : (-34, -25), 5 : (-42, -28), 6 : (-48, -33),
+        7 : (-54, -37), 8 : (-61, -40)
+    }
+
     for idim in xrange(len(data)):
         dim = args.dimensions[idim]
         print '= dim', dim
         argsc.dimension = dim
 
-        yranges = [np.inf, -np.inf]
+        # yranges = [np.inf, -np.inf]
         legend_handles = []
 
         fig = plt.figure(figsize=(7, 5))
         ax = fig.add_subplot(111)
         ax.set_xlim(0, len(xticks)+1)
-        ax.set_xticklabels([''] + xticks + [''], fontsize=14)
+        ax.set_xticklabels([''] + xticks + [''], fontsize=16)
         ax.set_xlabel(r'BSM operator angle', fontsize=16)
         ax.set_ylabel(r'${\mathrm {log}}_{10} \left (\Lambda^{-1}' + \
-                      get_units(dim) +r'\right )$', fontsize=16)
+                      get_units(dim) +r'\right )$', fontsize=17)
 
-        # ax.tick_params(axis='x', labelsize=11)
-        ax.tick_params(axis='y', labelsize=14)
+        ax.tick_params(axis='y', labelsize=15)
 
         for isrc in xrange(len(data[idim])):
             src = args.source_ratios[isrc]
@@ -418,24 +618,19 @@ def plot_sens_fixed_angle(data, outfile, outformat, args):
                 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
-                print 'reduced_ev', reduced_ev
                 if len(al) == 0:
-                    print 'No points for DIM {0} FRS {1}\nreduced_ev {2}!'.format(
-                        dim, src, reduced_ev
-                    )
+                    print 'No points for DIM {0} FRS {1}!'.format(dim, src)
                     continue
                 if reduced_ev[-1] < np.log(10**(3/2.)) - 0.1:
                     print 'Peaked contour does not exclude large scales! For ' \
-                        'DIM {0} FRS{1} reduced_ev {2}!'.format(
-                        dim, src, reduced_ev
-                    )
+                        'DIM {0} FRS{1}!'.format(dim, src)
                     continue
-                arr_len = dim-2
+                arr_len = 1.5
                 lim = al[0]
                 print 'limit = {0}'.format(lim)
                 label = '{0} : {1} : {2}'.format(*misc_utils.solve_ratio(src))
-                if lim < yranges[0]: yranges[0] = lim-arr_len
-                if lim > yranges[1]: yranges[1] = lim+arr_len+2
+                # if lim < yranges[0]: yranges[0] = lim-arr_len
+                # if lim > yranges[1]: yranges[1] = lim+arr_len+2
                 # if lim > yranges[1]: yranges[1] = lim
                 xoff = 0.15
                 line = plt.Line2D(
@@ -446,19 +641,30 @@ def plot_sens_fixed_angle(data, outfile, outformat, args):
                     legend_handles.append(line)
                 x_offset = isrc*xoff/2. - xoff/2.
                 ax.annotate(
-                    s='', xy=(ian+1+x_offset, lim-0.02), xytext=(ian+1+x_offset, lim+arr_len),
+                    s='', xy=(ian+1+x_offset, lim+0.02), xytext=(ian+1+x_offset, lim-arr_len),
                     arrowprops={'arrowstyle': '<-', 'lw': 1.5, 'color':colour[isrc]}
                 )
-
-        try:
-            yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True))
-            ax.set_ylim(yranges)
-        except: pass
-
-        legend = ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right',
+                if ian == 2:
+                    lim = np.log10(LV_atmo_90pc_limits[dim])
+                    ax.add_patch(patches.Rectangle(
+                        (ian+1-xoff, lim), 2*xoff, 100, fill=True,
+                        facecolor='orange', alpha=0.3, linewidth=1, edgecolor='k'
+                    ))
+                    ax.annotate(s='IC atmospheric', xy=(ian+1, lim+0.13),
+                                color='red', rotation=90, fontsize=7,
+                                horizontalalignment='center',
+                                verticalalignment='bottom')
+
+        # try:
+        #     yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True))
+        #     ax.set_ylim(yranges)
+        # except: pass
+        ax.set_ylim(ylims[dim])
+
+        legend = ax.legend(handles=legend_handles, prop=dict(size=10), loc='lower left',
                            title=r'$\nu_e:\nu_\mu:\nu_\tau{\rm\:\:at\:\:source}$',
                            framealpha=1., edgecolor='black')
-        plt.setp(legend.get_title(), fontsize='8')
+        plt.setp(legend.get_title(), fontsize='10')
         legend.get_frame().set_linestyle('-')
 
         an_text = 'Dimension {0}'.format(dim)
@@ -468,20 +674,22 @@ def plot_sens_fixed_angle(data, outfile, outformat, args):
         at.patch.set_boxstyle("round,pad=0.1,rounding_size=0.5")
         ax.add_artist(at)
 
-        fig.text(0.42, 0.8, r'Excluded', color='red', fontsize=16, ha='center',
+        fig.text(0.52, 0.8, r'Excluded', color='red', fontsize=16, ha='center',
                  va='center', fontweight='bold')
+        # fig.text(0.52, 0.76, r'with strong evidence', color='red', fontsize=11,
+        #          ha='center', va='center')
 
         if args.data is DataType.REAL:
-            fig.text(0.805, 0.14, 'IceCube Preliminary', color='red', fontsize=9,
+            fig.text(0.77, 0.14, 'IceCube Preliminary', color='red', fontsize=10,
                      ha='center', va='center')
-        elif args.data is DataType.ASIMOV:
-            fig.text(0.805, 0.14, 'IceCube Simulation', color='red', fontsize=9,
+        elif args.data in [DataType.ASIMOV, DataType.REALISATION]:
+            fig.text(0.77, 0.14, 'IceCube Simulation', color='red', fontsize=10,
                      ha='center', va='center')
 
         for ymaj in ax.yaxis.get_majorticklocs():
-            ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1)
+            ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.2, linewidth=1)
         for xmaj in ax.xaxis.get_majorticklocs():
-            ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1)
+            ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.2, linewidth=1)
 
         out = outfile + '_DIM{0}'.format(dim)
         misc_utils.make_dir(out)
@@ -489,11 +697,6 @@ def plot_sens_fixed_angle(data, outfile, outformat, args):
             print 'Saving plot as {0}'.format(out+'.'+of)
             fig.savefig(out+'.'+of, bbox_inches='tight', dpi=150)
 
-        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():
-            ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1)
-
 
 def plot_sens_corr_angle(data, outfile, outformat, args):
     print 'Making CORR_ANGLE sensitivity plot'