Sun Apr 22 23:18:44 CDT 2018

23 files changed, 2267 insertions, 1146 deletions

diff --git a/bout/plot.py b/bout/plot.py
new file mode 100644
index 0000000..717cf81
--- /dev/null
+++ b/bout/plot.py
@@ -0,0 +1,129 @@
+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
new file mode 100644
index 0000000..a75fe8a
--- /dev/null
+++ b/bout/plot_corr.py
@@ -0,0 +1,193 @@
+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
new file mode 100644
index 0000000..f2e1919
--- /dev/null
+++ b/bout/plot_full.py
@@ -0,0 +1,127 @@
+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/fr.py b/fr.py
index 906c0bf..363cb8e 100755
--- a/fr.py
+++ b/fr.py
@@ -5,27 +5,25 @@
 # date   : March 17, 2018
 
 """
-HESE BSM flavour ratio analysis script
+HESE BSM flavour ratio MCMC analysis script
 """
 
 from __future__ import absolute_import, division
 
-import os
 import argparse
 from functools import partial
 
 import numpy as np
 
-import GolemFitPy as gf
-
+from utils import fr as fr_utils
 from utils import gf as gf_utils
 from utils import likelihood as llh_utils
 from utils import mcmc as mcmc_utils
 from utils import misc as misc_utils
 from utils import plot as plot_utils
 from utils.enums import EnergyDependance, Likelihood, MCMCSeedType, ParamTag
-from utils.fr import MASS_EIGENVALUES, normalise_fr, fr_to_angles
-from utils.misc import enum_parse, Param, ParamSet
+from utils.fr import estimate_scale, normalise_fr
+from utils.param import Param, ParamSet, get_paramsets
 
 
 def define_nuisance():
@@ -52,56 +50,6 @@ def nuisance_argparse(parser):
         )
 
 
-def get_paramsets(args):
-    """Make the paramsets for generating the Asmimov MC sample and also running
-    the MCMC.
-    """
-    asimov_paramset = []
-    mcmc_paramset = []
-    if args.likelihood == Likelihood.GOLEMFIT:
-        nuisance_paramset = define_nuisance()
-        asimov_paramset.extend(nuisance_paramset.params)
-        mcmc_paramset.extend(nuisance_paramset.params)
-        for parm in nuisance_paramset:
-            parm.value = args.__getattribute__(parm.name)
-    tag = ParamTag.BESTFIT
-    flavour_angles = fr_to_angles(args.measured_ratio)
-    asimov_paramset.extend([
-        Param(name='astroFlavorAngle1', value=flavour_angles[0], ranges=[0., 1.], std=0.2, tag=tag),
-        Param(name='astroFlavorAngle2', value=flavour_angles[1], ranges=[-1., 1.], std=0.2, tag=tag),
-    ])
-    asimov_paramset = ParamSet(asimov_paramset)
-
-    if not args.fix_mixing and not args.fix_mixing_almost:
-        tag = ParamTag.MMANGLES
-        e = 1e-10
-        mcmc_paramset.extend([
-            Param(name='s_12^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag),
-            Param(name='c_13^4', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag),
-            Param(name='s_23^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag),
-            Param(name='dcp', value=np.pi, ranges=[0.+e, 2*np.pi-e], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag)
-        ])
-    if args.fix_mixing_almost:
-        tag = ParamTag.MMANGLES
-        mcmc_paramset.extend([
-            Param(name='s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag)
-        ])
-    if not args.fix_scale:
-        tag = ParamTag.SCALE
-        mcmc_paramset.append(
-            Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3,
-                  tex=r'{\rm log}_{10}\Lambda'+plot_utils.get_units(args.dimension), tag=tag)
-        )
-    if not args.fix_source_ratio:
-        tag = ParamTag.SRCANGLES
-        mcmc_paramset.extend([
-            Param(name='s_phi4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'sin^4(\phi)', tag=tag),
-            Param(name='c_2psi', value=0.5, ranges=[0., 1.], std=0.2, tex=r'cos(2\psi)', tag=tag)
-        ])
-    mcmc_paramset = ParamSet(mcmc_paramset)
-    return asimov_paramset, mcmc_paramset
-
-
 def process_args(args):
     """Process the input args."""
     if args.fix_mixing and args.fix_scale:
@@ -110,10 +58,6 @@ def process_args(args):
         raise NotImplementedError(
             '--fix-mixing and --fix-mixing-almost cannot be used together'
         )
-    if args.run_bayes_factor and args.fix_scale:
-        raise NotImplementedError(
-            '--run-bayes-factor and --fix-scale cannot be used together'
-        )
 
     args.measured_ratio = normalise_fr(args.measured_ratio)
     if args.fix_source_ratio:
@@ -125,26 +69,9 @@ def process_args(args):
         )
 
     if not args.fix_scale:
-        if args.energy_dependance is EnergyDependance.MONO:
-            args.scale = np.power(
-                10, np.round(np.log10(MASS_EIGENVALUES[1]/args.energy)) - \
-                np.log10(args.energy**(args.dimension-3)) + 6
-            )
-        elif args.energy_dependance is EnergyDependance.SPECTRAL:
-            args.scale = np.power(
-                10, np.round(
-                    np.log10(MASS_EIGENVALUES[1]/np.power(10, np.average(np.log10(args.binning)))) \
-                    - np.log10(np.power(10, np.average(np.log10(args.binning)))**(args.dimension-3))
-                ) + 6
-            )
-            """Estimate the scale of when to expect the BSM term to contribute."""
+        args.scale = estimate_scale(args)
         args.scale_region = (args.scale/args.scale_region, args.scale*args.scale_region)
 
-    if args.bayes_eval_bin.lower() == 'all':
-        args.bayes_eval_bin = None
-    else:
-        args.bayes_eval_bin = int(args.bayes_eval_bin)
-
 
 def parse_args(args=None):
     """Parse command line arguments"""
@@ -153,108 +80,7 @@ def parse_args(args=None):
         formatter_class=misc_utils.SortingHelpFormatter,
     )
     parser.add_argument(
-        '--measured-ratio', type=int, nargs=3, default=[1, 1, 1],
-        help='Set the central value for the measured flavour ratio at IceCube'
-    )
-    parser.add_argument(
-        '--sigma-ratio', type=float, default=0.01,
-        help='Set the 1 sigma for the measured flavour ratio for a gaussian LLH'
-    )
-    parser.add_argument(
-        '--fix-source-ratio', type=misc_utils.parse_bool, default='False',
-        help='Fix the source flavour ratio'
-    )
-    parser.add_argument(
-        '--energy-dependance', default='spectral', type=partial(enum_parse, c=EnergyDependance),
-        choices=EnergyDependance,
-        help='Type of energy dependance to use in the BSM fit'
-    )
-    parser.add_argument(
-        '--spectral-index', default=-2, type=int,
-        help='Spectral index for spectral energy dependance'
-    )
-    parser.add_argument(
-        '--binning', default=[1e4, 1e7, 5], type=float, nargs=3,
-        help='Binning for spectral energy dependance'
-    )
-    parser.add_argument(
-        '--run-bayes-factor', type=misc_utils.parse_bool, default='False',
-        help='Make the bayes factor plot for the scale'
-    )
-    parser.add_argument(
-        '--run-angles-limit', type=misc_utils.parse_bool, default='False',
-        help='Make the limit vs BSM angles plot'
-    )
-    parser.add_argument(
-        '--run-angles-correlation', type=misc_utils.parse_bool, default='False',
-        help='Make the limit vs BSM angles plot'
-    )
-    parser.add_argument(
-        '--bayes-bins', type=int, default=10,
-        help='Number of bins for the Bayes factor plot'
-    )
-    parser.add_argument(
-        '--bayes-eval-bin', type=str, default='all',
-        help='Which bin to evalaute for Bayes factor plot'
-    )
-    parser.add_argument(
-        '--bayes-live-points', type=int, default=400,
-        help='Number of live points for MultiNest evaluations'
-    )
-    parser.add_argument(
-        '--bayes-tolerance', type=float, default=0.5,
-        help='Tolerance for MultiNest'
-    )
-    parser.add_argument(
-        '--bayes-output', type=str, default='./mnrun/',
-        help='Folder to store MultiNest evaluations'
-    )
-    parser.add_argument(
-        '--angles-lim-output', type=str, default='./mnrun/',
-        help='Folder to store MultiNest evaluations'
-    )
-    parser.add_argument(
-        '--angles-corr-output', type=str, default='./mnrun/',
-        help='Folder to store MultiNest evaluations'
-    )
-    parser.add_argument(
-        '--source-ratio', type=int, nargs=3, default=[2, 1, 0],
-        help='Set the source flavour ratio for the case when you want to fix it'
-    )
-    parser.add_argument(
-        '--no-bsm', type=misc_utils.parse_bool, default='False',
-        help='Turn off BSM terms'
-    )
-    parser.add_argument(
-        '--fix-mixing', type=misc_utils.parse_bool, default='False',
-        help='Fix all mixing parameters to bi-maximal mixing'
-    )
-    parser.add_argument(
-        '--fix-mixing-almost', type=misc_utils.parse_bool, default='False',
-        help='Fix all mixing parameters except s23'
-    )
-    parser.add_argument(
-        '--fix-scale', type=misc_utils.parse_bool, default='False',
-        help='Fix the new physics scale'
-    )
-    parser.add_argument(
-        '--scale', type=float, default=1e-30,
-        help='Set the new physics scale'
-    )
-    parser.add_argument(
-        '--scale-region', type=float, default=1e7,
-        help='Set the size of the box to scan for the scale'
-    )
-    parser.add_argument(
-        '--dimension', type=int, default=3,
-        help='Set the new physics dimension to consider'
-    )
-    parser.add_argument(
-        '--energy', type=float, default=1000,
-        help='Set the energy scale'
-    )
-    parser.add_argument(
-        '--seed', type=int, default=99,
+        '--seed', type=misc_utils.seed_parse, default='25',
         help='Set the random seed value'
     )
     parser.add_argument(
@@ -263,12 +89,12 @@ def parse_args(args=None):
     )
     parser.add_argument(
         '--outfile', type=str, default='./untitled',
-        help='Path to output chains'
+        help='Path to output results'
     )
+    fr_utils.fr_argparse(parser)
+    gf_utils.gf_argparse(parser)
     llh_utils.likelihood_argparse(parser)
     mcmc_utils.mcmc_argparse(parser)
-    gf_utils.gf_argparse(parser)
-    plot_utils.plot_argparse(parser)
     nuisance_argparse(parser)
     if args is None: return parser.parse_args()
     else: return parser.parse_args(args.split())
@@ -279,31 +105,31 @@ def main():
     process_args(args)
     misc_utils.print_args(args)
 
-    np.random.seed(args.seed)
+    if args.seed is not None:
+        np.random.seed(args.seed)
 
-    asimov_paramset, mcmc_paramset = get_paramsets(args)
+    asimov_paramset, llh_paramset = get_paramsets(args)
     outfile = misc_utils.gen_outfile_name(args)
     print '== {0:<25} = {1}'.format('outfile', outfile)
 
     if args.run_mcmc:
         if args.likelihood is Likelihood.GOLEMFIT:
             fitter = gf_utils.setup_fitter(args, asimov_paramset)
-        else:
-            fitter = None
+        else: fitter = None
 
         ln_prob = partial(
             llh_utils.ln_prob, args=args, fitter=fitter,
-            asimov_paramset=asimov_paramset, mcmc_paramset=mcmc_paramset
+            asimov_paramset=asimov_paramset, llh_paramset=llh_paramset
         )
 
-        ndim = len(mcmc_paramset)
+        ndim = len(llh_paramset)
         if args.mcmc_seed_type == MCMCSeedType.UNIFORM:
             p0 = mcmc_utils.flat_seed(
-                mcmc_paramset, nwalkers=args.nwalkers
+                llh_paramset, nwalkers=args.nwalkers
             )
         elif args.mcmc_seed_type == MCMCSeedType.GAUSSIAN:
             p0 = mcmc_utils.gaussian_seed(
-                mcmc_paramset, nwalkers=args.nwalkers
+                llh_paramset, nwalkers=args.nwalkers
             )
 
         samples = mcmc_utils.mcmc(
@@ -320,275 +146,12 @@ def main():
         mcmc_utils.save_chains(samples, outfile)
 
     plot_utils.chainer_plot(
-        infile        = outfile+'.npy',
-        outfile       = outfile[:5]+outfile[5:].replace('data', 'plots'),
-        outformat     = ['pdf'],
-        args          = args,
-        mcmc_paramset = mcmc_paramset
-    )
-
-    out = args.bayes_output+'/fr_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
-    )
-    if args.run_bayes_factor:
-        import pymultinest
-
-        if not args.run_mcmc and args.likelihood is Likelihood.GOLEMFIT:
-            fitter = gf_utils.setup_fitter(args, asimov_paramset)
-        else: fitter = None
-
-        p = mcmc_paramset.from_tag(ParamTag.SCALE, invert=True)
-        n_params = len(p)
-        prior_ranges = p.seeds
-
-        def CubePrior(cube, ndim, nparams):
-            # default are uniform priors
-            return ;
-
-        arr = []
-        for s_idx, sc in enumerate(scan_scales):
-            if args.bayes_eval_bin is not None:
-                if args.bayes_eval_bin == s_idx:
-                    out += '_scale_{0:.0E}'.format(np.power(10, sc))
-                else: continue
-
-            print '== SCALE = {0:.0E}'.format(np.power(10, sc))
-            theta = np.zeros(n_params)
-            def lnProb(cube, ndim, nparams):
-                for i in range(ndim):
-                    prange = prior_ranges[i][1] - prior_ranges[i][0]
-                    theta[i] = prange*cube[i] + prior_ranges[i][0]
-                theta_ = np.array(theta.tolist() + [sc])
-                # print 'mcmc_paramset', mcmc_paramset
-                return llh_utils.triangle_llh(
-                    theta=theta_,
-                    args=args,
-                    asimov_paramset=asimov_paramset,
-                    mcmc_paramset=mcmc_paramset,
-                    fitter=fitter
-                )
-
-            prefix = 'mnrun/' + os.path.basename(out) + '_'
-            if args.bayes_eval_bin is not None:
-                prefix += '{0}_'.format(s_idx)
-            print 'begin running evidence calculation for {0}'.format(prefix)
-            result = pymultinest.run(
-                LogLikelihood=lnProb,
-                Prior=CubePrior,
-                n_dims=n_params,
-                importance_nested_sampling=True,
-                n_live_points=args.bayes_live_points,
-                evidence_tolerance=args.bayes_tolerance,
-                outputfiles_basename=prefix,
-                resume=False,
-                verbose=True
-            )
-
-            analyzer = pymultinest.Analyzer(
-                outputfiles_basename=prefix, n_params=n_params
-            )
-            a_lnZ = analyzer.get_stats()['global evidence']
-            print 'Evidence = {0}'.format(a_lnZ)
-            arr.append([sc, a_lnZ])
-
-        misc_utils.make_dir(out)
-        np.save(out+'.npy', np.array(arr))
-
-    dirname = os.path.dirname(out)
-    plot_utils.bayes_factor_plot(
-        dirname=dirname, outfile=out, outformat=['png'], args=args
-    )
-
-    out = args.angles_lim_output+'/fr_an_evidence' + misc_utils.gen_identifier(args)
-    if args.run_angles_limit:
-        import pymultinest
-
-        scenarios = [
-            [np.sin(np.pi/2.)**2, 0, 0, 0],
-            [0, np.cos(np.pi/2.)**4, 0, 0],
-            [0, 0, np.sin(np.pi/2.)**2, 0],
-        ]
-        p = mcmc_paramset.from_tag([ParamTag.SCALE, ParamTag.MMANGLES], invert=True)
-        n_params = len(p)
-        prior_ranges = p.seeds
-
-        if not args.run_mcmc and args.likelihood is Likelihood.GOLEMFIT:
-            fitter = gf_utils.setup_fitter(args, asimov_paramset)
-        else: fitter = None
-
-        def CubePrior(cube, ndim, nparams):
-            # default are uniform priors
-            return ;
-
-        if args.bayes_eval_bin is not None:
-            data = np.zeros((len(scenarios), 1, 2))
-        else: data = np.zeros((len(scenarios), args.bayes_bins, 2))
-        mm_angles = mcmc_paramset.from_tag(ParamTag.MMANGLES)
-        sc_angles = mcmc_paramset.from_tag(ParamTag.SCALE)[0]
-        for idx, scen in enumerate(scenarios):
-            scales, evidences = [], []
-            for yidx, an in enumerate(mm_angles):
-                an.value = scen[yidx]
-            for s_idx, sc in enumerate(scan_scales):
-                if args.bayes_eval_bin is not None:
-                    if args.bayes_eval_bin == s_idx:
-                        if idx == 0:
-                            out += '_scale_{0:.0E}'.format(np.power(10, sc))
-                    else: continue
-
-                print '== SCALE = {0:.0E}'.format(np.power(10, sc))
-                sc_angles.value = sc
-                def lnProb(cube, ndim, nparams):
-                    for i in range(ndim):
-                        prange = prior_ranges[i][1] - prior_ranges[i][0]
-                        p[i].value = prange*cube[i] + prior_ranges[i][0]
-                    for name in p.names:
-                        mcmc_paramset[name].value = p[name].value
-                    theta = mcmc_paramset.values
-                    # print 'theta', theta
-                    # print 'mcmc_paramset', mcmc_paramset
-                    return llh_utils.triangle_llh(
-                        theta=theta,
-                        args=args,
-                        asimov_paramset=asimov_paramset,
-                        mcmc_paramset=mcmc_paramset,
-                        fitter=fitter
-                    )
-                prefix = 'mnrun/' + os.path.basename(out) + '_'
-                if args.bayes_eval_bin is not None:
-                    prefix += '{0}_{1}_'.format(idx, s_idx)
-                print 'begin running evidence calculation for {0}'.format(prefix)
-                result = pymultinest.run(
-                    LogLikelihood=lnProb,
-                    Prior=CubePrior,
-                    n_dims=n_params,
-                    importance_nested_sampling=True,
-                    n_live_points=args.bayes_live_points,
-                    evidence_tolerance=args.bayes_tolerance,
-                    outputfiles_basename=prefix,
-                    resume=False,
-                    verbose=True
-                )
-
-                analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params)
-                a_lnZ = analyzer.get_stats()['global evidence']
-                print 'Evidence = {0}'.format(a_lnZ)
-                scales.append(sc)
-                evidences.append(a_lnZ)
-
-            for i, d in enumerate(evidences):
-                data[idx][i][0] = scales[i]
-                data[idx][i][1] = d
-
-        misc_utils.make_dir(out)
-        print 'saving to {0}.npy'.format(out)
-        np.save(out+'.npy', np.array(data))
-
-    dirname = os.path.dirname(out)
-    plot_utils.plot_BSM_angles_limit(
-        dirname=dirname, outfile=outfile, outformat=['png'],
-        args=args, bayesian=True
+        infile       = outfile+'.npy',
+        outfile      = outfile[:5]+outfile[5:].replace('data', 'plots'),
+        outformat    = ['pdf'],
+        args         = args,
+        llh_paramset = llh_paramset
     )
-
-    out = args.angles_corr_output+'/fr_co_evidence' + misc_utils.gen_identifier(args)
-    if args.run_angles_correlation:
-        if args.bayes_eval_bin is None: assert 0
-        import pymultinest
-
-        scenarios = [
-            [1, 0, 0, 0],
-            [0, 1, 0, 0],
-            [0, 0, 1, 0],
-        ]
-        nuisance = mcmc_paramset.from_tag(ParamTag.NUISANCE)
-        mm_angles = mcmc_paramset.from_tag(ParamTag.MMANGLES)
-        sc_angles = mcmc_paramset.from_tag(ParamTag.SCALE)[0]
-
-        if not args.run_mcmc and args.likelihood is Likelihood.GOLEMFIT:
-            fitter = gf_utils.setup_fitter(args, asimov_paramset)
-        else: fitter = None
-
-        def CubePrior(cube, ndim, nparams):
-            # default are uniform priors
-            return ;
-
-        scan_angles = np.linspace(0, 1, args.bayes_bins)
-
-        if args.bayes_eval_bin is not None:
-            data = np.zeros((len(scenarios), 1, 1, 3))
-        else: data = np.zeros((len(scenarios), scale_bins, angle_bins, 3))
-        for idx, scen in enumerate(scenarios):
-            for an in mm_angles:
-                an.value = 0
-            keep = mcmc_paramset.from_tag(ParamTag.MMANGLES)[idx]
-            q = ParamSet(nuisance.params + [x for x in mm_angles if x.name != keep.name])
-            n_params = len(q)
-            prior_ranges = q.seeds
-
-            scales, angles, evidences = [], [], []
-            for s_idx, sc in enumerate(scan_scales):
-                for a_idx, an in enumerate(scan_angles):
-                    index = s_idx*args.bayes_bins + a_idx
-                    if args.bayes_eval_bin is not None:
-                        if args.bayes_eval_bin == index:
-                            if idx == 0:
-                                out += '_scale_{0:.0E}'.format(np.power(10, sc))
-                        else: continue
-
-                    print '== SCALE = {0:.0E}'.format(np.power(10, sc))
-                    print '== ANGLE = {0:.2f}'.format(an)
-                    sc_angles.value = sc
-                    keep.value = an
-                    def lnProb(cube, ndim, nparams):
-                        for i in range(ndim-1):
-                            prange = prior_ranges[i][1] - prior_ranges[i][0]
-                            q[i].value = prange*cube[i] + prior_ranges[i][0]
-                        for name in q.names:
-                            mcmc_paramset[name].value = q[name].value
-                        theta = mcmc_paramset.values
-                        # print 'theta', theta
-                        # print 'mcmc_paramset', mcmc_paramset
-                        return llh_utils.triangle_llh(
-                            theta=theta,
-                            args=args,
-                            asimov_paramset=asimov_paramset,
-                            mcmc_paramset=mcmc_paramset,
-                            fitter=fitter
-                        )
-                    prefix = 'mnrun/' + os.path.basename(out) + '_'
-                    if args.bayes_eval_bin is not None:
-                        prefix += '{0}_{1}_{2}'.format(idx, s_idx, a_idx)
-
-                    print 'begin running evidence calculation for {0}'.format(prefix)
-                    result = pymultinest.run(
-                        LogLikelihood=lnProb,
-                        Prior=CubePrior,
-                        n_dims=n_params,
-                        importance_nested_sampling=True,
-                        n_live_points=args.bayes_live_points,
-                        evidence_tolerance=args.bayes_tolerance,
-                        outputfiles_basename=prefix,
-                        resume=False,
-                        verbose=True
-                    )
-
-                    analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params)
-                    a_lnZ = analyzer.get_stats()['global evidence']
-                    print 'Evidence = {0}'.format(a_lnZ)
-                    scales.append(sc)
-                    angles.append(an)
-                    evidences.append(a_lnZ)
-
-                for i, d in enumerate(evidences):
-                    data[idx][i][i][0] = scales[i]
-                    data[idx][i][i][1] = angles[i]
-                    data[idx][i][i][2] = d
-
-        misc_utils.make_dir(out)
-        print 'saving to {0}.npy'.format(out)
-        np.save(out+'.npy', np.array(data))
-
     print "DONE!"
 
 
diff --git a/fr_edit.py b/fr_edit.py
new file mode 100755
index 0000000..df3b43b
--- /dev/null
+++ b/fr_edit.py
@@ -0,0 +1,581 @@
+#! /usr/bin/env python
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : March 17, 2018
+
+"""
+HESE BSM flavour ratio analysis script
+"""
+
+from __future__ import absolute_import, division
+
+import os
+import argparse
+from functools import partial
+
+import numpy as np
+
+import GolemFitPy as gf
+
+from utils import gf as gf_utils
+from utils import likelihood as llh_utils
+from utils import mcmc as mcmc_utils
+from utils import misc as misc_utils
+from utils import plot as plot_utils
+from utils.enums import EnergyDependance, Likelihood, MCMCSeedType, ParamTag
+from utils.fr import MASS_EIGENVALUES, normalise_fr, fr_to_angles
+from utils.misc import enum_parse, Param, ParamSet
+
+
+def define_nuisance():
+    """Define the nuisance parameters to scan over with default values,
+    ranges and sigma.
+    """
+    tag = ParamTag.NUISANCE
+    nuisance = ParamSet(
+        Param(name='convNorm',        value=1.,  seed=[0.5, 2. ], ranges=[0. , 50.], std=0.3,  tag=tag),
+        Param(name='promptNorm',      value=0.,  seed=[0. , 6. ], ranges=[0. , 50.], std=0.05, tag=tag),
+        Param(name='muonNorm',        value=1.,  seed=[0.1, 2. ], ranges=[0. , 50.], std=0.1,  tag=tag),
+        Param(name='astroNorm',       value=6.9, seed=[0.1, 10.], ranges=[0. , 50.], std=0.1,  tag=tag),
+        Param(name='astroDeltaGamma', value=2.5, seed=[1. , 3. ], ranges=[-5., 5. ], std=0.1,  tag=tag)
+    )
+    return nuisance
+
+
+def nuisance_argparse(parser):
+    nuisance_paramset = define_nuisance()
+    for parm in nuisance_paramset:
+        parser.add_argument(
+            '--'+parm.name, type=float, default=parm.value,
+            help=parm.name+' to inject'
+        )
+
+
+def get_paramsets(args):
+    """Make the paramsets for generating the Asmimov MC sample and also running
+    the MCMC.
+    """
+    asimov_paramset = []
+    mcmc_paramset = []
+    if args.likelihood == Likelihood.GOLEMFIT:
+        nuisance_paramset = define_nuisance()
+        asimov_paramset.extend(nuisance_paramset.params)
+        mcmc_paramset.extend(nuisance_paramset.params)
+        for parm in nuisance_paramset:
+            parm.value = args.__getattribute__(parm.name)
+    tag = ParamTag.BESTFIT
+    flavour_angles = fr_to_angles(args.measured_ratio)
+    asimov_paramset.extend([
+        Param(name='astroFlavorAngle1', value=flavour_angles[0], ranges=[0., 1.], std=0.2, tag=tag),
+        Param(name='astroFlavorAngle2', value=flavour_angles[1], ranges=[-1., 1.], std=0.2, tag=tag),
+    ])
+    asimov_paramset = ParamSet(asimov_paramset)
+
+    if not args.fix_mixing and not args.fix_mixing_almost:
+        tag = ParamTag.MMANGLES
+        e = 1e-10
+        mcmc_paramset.extend([
+            Param(name='s_12^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag),
+            Param(name='c_13^4', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag),
+            Param(name='s_23^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag),
+            Param(name='dcp', value=np.pi, ranges=[0., 2*np.pi], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag)
+        ])
+    if args.fix_mixing_almost:
+        tag = ParamTag.MMANGLES
+        mcmc_paramset.extend([
+            Param(name='s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag)
+        ])
+    if not args.fix_scale:
+        tag = ParamTag.SCALE
+        mcmc_paramset.append(
+            Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3,
+                  tex=r'{\rm log}_{10}\Lambda'+plot_utils.get_units(args.dimension), tag=tag)
+        )
+    if not args.fix_source_ratio:
+        tag = ParamTag.SRCANGLES
+        mcmc_paramset.extend([
+            Param(name='s_phi4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'sin^4(\phi)', tag=tag),
+            Param(name='c_2psi', value=0.5, ranges=[0., 1.], std=0.2, tex=r'cos(2\psi)', tag=tag)
+        ])
+    mcmc_paramset = ParamSet(mcmc_paramset)
+    return asimov_paramset, mcmc_paramset
+
+
+def process_args(args):
+    """Process the input args."""
+    if args.fix_mixing and args.fix_scale:
+        raise NotImplementedError('Fixed mixing and scale not implemented')
+    if args.fix_mixing and args.fix_mixing_almost:
+        raise NotImplementedError(
+            '--fix-mixing and --fix-mixing-almost cannot be used together'
+        )
+    if args.run_bayes_factor and args.fix_scale:
+        raise NotImplementedError(
+            '--run-bayes-factor and --fix-scale cannot be used together'
+        )
+
+    args.measured_ratio = normalise_fr(args.measured_ratio)
+    if args.fix_source_ratio:
+        args.source_ratio = normalise_fr(args.source_ratio)
+
+    if args.energy_dependance is EnergyDependance.SPECTRAL:
+        args.binning = np.logspace(
+            np.log10(args.binning[0]), np.log10(args.binning[1]), args.binning[2]+1
+        )
+
+    if not args.fix_scale:
+        if args.energy_dependance is EnergyDependance.MONO:
+            args.scale = np.power(
+                10, np.round(np.log10(MASS_EIGENVALUES[1]/args.energy)) - \
+                np.log10(args.energy**(args.dimension-3)) + 6
+            )
+        elif args.energy_dependance is EnergyDependance.SPECTRAL:
+            args.scale = np.power(
+                10, np.round(
+                    np.log10(MASS_EIGENVALUES[1]/np.power(10, np.average(np.log10(args.binning)))) \
+                    - np.log10(np.power(10, np.average(np.log10(args.binning)))**(args.dimension-3))
+                )
+            )
+            """Estimate the scale of when to expect the BSM term to contribute."""
+        args.scale_region = (args.scale/args.scale_region, args.scale*args.scale_region)
+
+    if args.bayes_eval_bin.lower() == 'all':
+        args.bayes_eval_bin = None
+    else:
+        args.bayes_eval_bin = map(int, args.bayes_eval_bin.split())
+
+
+def parse_args(args=None):
+    """Parse command line arguments"""
+    parser = argparse.ArgumentParser(
+        description="BSM flavour ratio analysis",
+        formatter_class=misc_utils.SortingHelpFormatter,
+    )
+    parser.add_argument(
+        '--measured-ratio', type=int, nargs=3, default=[1, 1, 1],
+        help='Set the central value for the measured flavour ratio at IceCube'
+    )
+    parser.add_argument(
+        '--sigma-ratio', type=float, default=0.01,
+        help='Set the 1 sigma for the measured flavour ratio for a gaussian LLH'
+    )
+    parser.add_argument(
+        '--fix-source-ratio', type=misc_utils.parse_bool, default='False',
+        help='Fix the source flavour ratio'
+    )
+    parser.add_argument(
+        '--energy-dependance', default='spectral', type=partial(enum_parse, c=EnergyDependance),
+        choices=EnergyDependance,
+        help='Type of energy dependance to use in the BSM fit'
+    )
+    parser.add_argument(
+        '--spectral-index', default=-2, type=int,
+        help='Spectral index for spectral energy dependance'
+    )
+    parser.add_argument(
+        '--binning', default=[1e4, 1e7, 5], type=float, nargs=3,
+        help='Binning for spectral energy dependance'
+    )
+    parser.add_argument(
+        '--run-bayes-factor', type=misc_utils.parse_bool, default='False',
+        help='Make the bayes factor plot for the scale'
+    )
+    parser.add_argument(
+        '--run-angles-limit', type=misc_utils.parse_bool, default='False',
+        help='Make the limit vs BSM angles plot'
+    )
+    parser.add_argument(
+        '--run-angles-correlation', type=misc_utils.parse_bool, default='False',
+        help='Make the limit vs BSM angles plot'
+    )
+    parser.add_argument(
+        '--bayes-bins', type=int, default=10,
+        help='Number of bins for the Bayes factor plot'
+    )
+    parser.add_argument(
+        '--bayes-eval-bin', type=str, default='all',
+        help='Which bin to evalaute for Bayes factor plot'
+    )
+    parser.add_argument(
+        '--bayes-live-points', type=int, default=400,
+        help='Number of live points for MultiNest evaluations'
+    )
+    parser.add_argument(
+        '--bayes-tolerance', type=float, default=0.5,
+        help='Tolerance for MultiNest'
+    )
+    parser.add_argument(
+        '--bayes-output', type=str, default='./mnrun/',
+        help='Folder to store MultiNest evaluations'
+    )
+    parser.add_argument(
+        '--angles-lim-output', type=str, default='./mnrun/',
+        help='Folder to store MultiNest evaluations'
+    )
+    parser.add_argument(
+        '--angles-corr-output', type=str, default='./mnrun/',
+        help='Folder to store MultiNest evaluations'
+    )
+    parser.add_argument(
+        '--source-ratio', type=int, nargs=3, default=[2, 1, 0],
+        help='Set the source flavour ratio for the case when you want to fix it'
+    )
+    parser.add_argument(
+        '--no-bsm', type=misc_utils.parse_bool, default='False',
+        help='Turn off BSM terms'
+    )
+    parser.add_argument(
+        '--fix-mixing', type=misc_utils.parse_bool, default='False',
+        help='Fix all mixing parameters to bi-maximal mixing'
+    )
+    parser.add_argument(
+        '--fix-mixing-almost', type=misc_utils.parse_bool, default='False',
+        help='Fix all mixing parameters except s23'
+    )
+    parser.add_argument(
+        '--fix-scale', type=misc_utils.parse_bool, default='False',
+        help='Fix the new physics scale'
+    )
+    parser.add_argument(
+        '--scale', type=float, default=1e-30,
+        help='Set the new physics scale'
+    )
+    parser.add_argument(
+        '--scale-region', type=float, default=1e10,
+        help='Set the size of the box to scan for the scale'
+    )
+    parser.add_argument(
+        '--dimension', type=int, default=3,
+        help='Set the new physics dimension to consider'
+    )
+    parser.add_argument(
+        '--energy', type=float, default=1000,
+        help='Set the energy scale'
+    )
+    parser.add_argument(
+        '--seed', type=int, default=99,
+        help='Set the random seed value'
+    )
+    parser.add_argument(
+        '--threads', type=misc_utils.thread_type, default='1',
+        help='Set the number of threads to use (int or "max")'
+    )
+    parser.add_argument(
+        '--outfile', type=str, default='./untitled',
+        help='Path to output chains'
+    )
+    llh_utils.likelihood_argparse(parser)
+    mcmc_utils.mcmc_argparse(parser)
+    gf_utils.gf_argparse(parser)
+    plot_utils.plot_argparse(parser)
+    nuisance_argparse(parser)
+    if args is None: return parser.parse_args()
+    else: return parser.parse_args(args.split())
+
+
+def main():
+    args = parse_args()
+    process_args(args)
+    misc_utils.print_args(args)
+
+    np.random.seed(args.seed)
+
+    asimov_paramset, mcmc_paramset = get_paramsets(args)
+    outfile = misc_utils.gen_outfile_name(args)
+    print '== {0:<25} = {1}'.format('outfile', outfile)
+
+    if args.run_mcmc:
+        if args.likelihood is Likelihood.GOLEMFIT:
+            fitter = gf_utils.setup_fitter(args, asimov_paramset)
+        else:
+            fitter = None
+
+        ln_prob = partial(
+            llh_utils.ln_prob, args=args, fitter=fitter,
+            asimov_paramset=asimov_paramset, mcmc_paramset=mcmc_paramset
+        )
+
+        ndim = len(mcmc_paramset)
+        if args.mcmc_seed_type == MCMCSeedType.UNIFORM:
+            p0 = mcmc_utils.flat_seed(
+                mcmc_paramset, nwalkers=args.nwalkers
+            )
+        elif args.mcmc_seed_type == MCMCSeedType.GAUSSIAN:
+            p0 = mcmc_utils.gaussian_seed(
+                mcmc_paramset, nwalkers=args.nwalkers
+            )
+
+        samples = mcmc_utils.mcmc(
+            p0       = p0,
+            ln_prob  = ln_prob,
+            ndim     = ndim,
+            nwalkers = args.nwalkers,
+            burnin   = args.burnin,
+            nsteps   = args.nsteps,
+            threads  = 1
+            # TODO(shivesh): broken because you cannot pickle a GolemFitPy object
+            # threads      = misc_utils.thread_factors(args.threads)[0]
+        )
+        mcmc_utils.save_chains(samples, outfile)
+
+    plot_utils.chainer_plot(
+        infile        = outfile+'.npy',
+        outfile       = outfile[:5]+outfile[5:].replace('data', 'plots'),
+        outformat     = ['pdf'],
+        args          = args,
+        mcmc_paramset = mcmc_paramset
+    )
+
+    out = args.bayes_output+'/fr_fr_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
+    )
+    if args.run_bayes_factor:
+        from scipy.optimize import minimize
+        def fit_flags(flag_dict):
+            flags = gf.FitParametersFlag()
+            for key in flag_dict.iterkeys():
+                flags.__setattr__(key, flag_dict[key])
+            return flags
+
+        args.likelihood = Likelihood.GF_FREQ
+        default_flags = {
+            # False means it's not fixed in minimization
+            'astroFlavorAngle1'         : True,
+            'astroFlavorAngle2'         : True,
+            'astroENorm'                : False,
+            'astroMuNorm'               : False,
+            'astroTauNorm'              : False,
+            'convNorm'                  : False,
+            'promptNorm'                : False,
+            'muonNorm'                  : False,
+            'astroNorm'                 : False,
+            'astroParticleBalance'      : True,
+            'astroDeltaGamma'           : False,
+            'cutoffEnergy'              : True,
+            'CRDeltaGamma'              : True,
+            'piKRatio'                  : True,
+            'NeutrinoAntineutrinoRatio' : True,
+            'darkNorm'                  : True,
+            'domEfficiency'             : True,
+            'holeiceForward'            : True,
+            'anisotropyScale'           : True,
+            'astroNormSec'              : True,
+            'astroDeltaGammaSec'        : True
+        }
+
+        mcmc_paramset_freq = mcmc_paramset.from_tag(ParamTag.NUISANCE, invert=True)
+        scale = mcmc_paramset.from_tag(ParamTag.SCALE)
+
+        if not args.run_mcmc and args.likelihood is Likelihood.GF_FREQ:
+            fitter = gf_utils.setup_fitter(args, asimov_paramset)
+            fitter.SetFitParametersFlag(fit_flags(default_flags))
+        else: fitter = None
+
+        data = np.zeros((args.bayes_bins, 2))
+        for s_idx, sc in enumerate(scan_scales):
+            print '== SCALE = {0:.0E}'.format(np.power(10, sc))
+            def fn(scen):
+                theta = map(float, scen) + [0, sc]
+                try:
+                    llh = llh_utils.triangle_llh(
+                        theta=theta, args=args, asimov_paramset=asimov_paramset,
+                        mcmc_paramset=mcmc_paramset_freq, fitter=fitter
+                    )
+                    print 'llh', llh
+                except:
+                    print 'forbidden'
+                    return 9e99
+                return -llh
+            res = minimize(fun=fn, x0=[0.1, 0.1, 0.1], method='L-BFGS-B',
+                           bounds=[(0, 1), (0, 1), (0, 1)])
+            llh = fn(res.x)
+            data[s_idx][0] = sc
+            data[s_idx][1] = llh
+
+        misc_utils.make_dir(out)
+        np.save(out+'.npy', np.array(data))
+
+    # dirname = os.path.dirname(out)
+    # plot_utils.bayes_factor_plot(
+    #     dirname=dirname, outfile=out, outformat=['png'], args=args
+    # )
+
+    out = args.angles_lim_output+'/fr_anfr_evidence' + misc_utils.gen_identifier(args)
+    if args.run_angles_limit:
+        def fit_flags(flag_dict):
+            flags = gf.FitParametersFlag()
+            for key in flag_dict.iterkeys():
+                flags.__setattr__(key, flag_dict[key])
+            return flags
+
+        args.likelihood = Likelihood.GF_FREQ
+        default_flags = {
+            # False means it's not fixed in minimization
+            'astroFlavorAngle1'         : True,
+            'astroFlavorAngle2'         : True,
+            # 'astroENorm'                : True,
+            # 'astroMuNorm'               : True,
+            # 'astroTauNorm'              : True,
+            'convNorm'                  : False,
+            'promptNorm'                : False,
+            'muonNorm'                  : False,
+            'astroNorm'                 : False,
+            'astroParticleBalance'      : True,
+            'astroDeltaGamma'           : False,
+            'cutoffEnergy'              : True,
+            'CRDeltaGamma'              : True,
+            'piKRatio'                  : True,
+            'NeutrinoAntineutrinoRatio' : True,
+            'darkNorm'                  : True,
+            'domEfficiency'             : True,
+            'holeiceForward'            : True,
+            'anisotropyScale'           : True,
+            'astroNormSec'              : True,
+            'astroDeltaGammaSec'        : True
+        }
+
+        mcmc_paramset_freq = mcmc_paramset.from_tag(ParamTag.NUISANCE, invert=True)
+        scenarios = [
+            [np.sin(np.pi/2.)**2, 0, 0, 0],
+            [0, np.cos(np.pi/2.)**4, 0, 0],
+            [0, 0, np.sin(np.pi/2.)**2, 0],
+        ]
+
+        if not args.run_mcmc and args.likelihood is Likelihood.GF_FREQ:
+            fitter = gf_utils.setup_fitter(args, asimov_paramset)
+            fitter.SetFitParametersFlag(fit_flags(default_flags))
+        else: fitter = None
+
+        data = np.zeros((len(scenarios), args.bayes_bins, 2))
+        mm_angles = mcmc_paramset_freq.from_tag(ParamTag.MMANGLES)
+        sc_angles = mcmc_paramset_freq.from_tag(ParamTag.SCALE)[0]
+        for idx, scen in enumerate(scenarios):
+            scales, llhs = [], []
+            for yidx, an in enumerate(mm_angles):
+                an.value = scen[yidx]
+            for s_idx, sc in enumerate(scan_scales):
+                print '== SCALE = {0:.0E}'.format(np.power(10, sc))
+                theta = scen + [sc]
+                print 'theta', theta
+                llh = llh_utils.triangle_llh(
+                    theta=theta, args=args, asimov_paramset=asimov_paramset,
+                    mcmc_paramset=mcmc_paramset_freq, fitter=fitter
+                )
+                print 'llh', llh
+                scales.append(sc)
+                llhs.append(llh)
+
+            for i, d in enumerate(llhs):
+                data[idx][i][0] = scales[i]
+                data[idx][i][1] = d
+
+        misc_utils.make_dir(out)
+        print 'saving to {0}.npy'.format(out)
+        np.save(out+'.npy', np.array(data))
+
+    dirname = os.path.dirname(out)
+    plot_utils.plot_BSM_angles_limit(
+        dirname=dirname, outfile=outfile, outformat=['png'],
+        args=args, bayesian=True
+    )
+
+    out = args.angles_corr_output+'/fr_co_evidence' + misc_utils.gen_identifier(args)
+    if args.run_angles_correlation:
+        def fit_flags(flag_dict):
+            flags = gf.FitParametersFlag()
+            for key in flag_dict.iterkeys():
+                flags.__setattr__(key, flag_dict[key])
+            return flags
+
+        args.likelihood = Likelihood.GF_FREQ
+        default_flags = {
+            # False means it's not fixed in minimization
+            'astroFlavorAngle1'         : True,
+            'astroFlavorAngle2'         : True,
+            # 'astroENorm'                : True,
+            # 'astroMuNorm'               : True,
+            # 'astroTauNorm'              : True,
+            'convNorm'                  : False,
+            'promptNorm'                : False,
+            'muonNorm'                  : False,
+            'astroNorm'                 : False,
+            'astroParticleBalance'      : True,
+            'astroDeltaGamma'           : False,
+            'cutoffEnergy'              : True,
+            'CRDeltaGamma'              : True,
+            'piKRatio'                  : True,
+            'NeutrinoAntineutrinoRatio' : True,
+            'darkNorm'                  : True,
+            'domEfficiency'             : True,
+            'holeiceForward'            : True,
+            'anisotropyScale'           : True,
+            'astroNormSec'              : True,
+            'astroDeltaGammaSec'        : True
+        }
+
+        mcmc_paramset_freq = mcmc_paramset.from_tag(ParamTag.NUISANCE, invert=True)
+
+        scenarios = [
+            [1, 0, 0, 0],
+            [0, 1, 0, 0],
+            [0, 0, 1, 0],
+        ]
+        mm_angles = mcmc_paramset_freq.from_tag(ParamTag.MMANGLES)
+        sc_angles = mcmc_paramset_freq.from_tag(ParamTag.SCALE)[0]
+
+        if not args.run_mcmc and args.likelihood is Likelihood.GF_FREQ:
+            fitter = gf_utils.setup_fitter(args, asimov_paramset)
+            fitter.SetFitParametersFlag(fit_flags(default_flags))
+        else: fitter = None
+
+        scan_angles = np.linspace(0, 1, args.bayes_bins)
+
+        data = np.zeros((len(scenarios), args.bayes_bins, args.bayes_bins, 3))
+        for idx, scen in enumerate(scenarios):
+            for an in mm_angles:
+                an.value = 0
+            keep = mcmc_paramset_freq.from_tag(ParamTag.MMANGLES)[idx]
+
+            for s_idx, sc in enumerate(scan_scales):
+                scales, angles, llhs = [], [], []
+                for a_idx, an in enumerate(scan_angles):
+                    print '== SCALE = {0:.0E}'.format(np.power(10, sc))
+                    print '== ANGLE = {0:.2f}'.format(an)
+                    sc_angles.value = sc
+                    s2_2 = an
+                    if s_idx == 0: keep.value = np.sin(np.arcsin(np.sqrt(s2_2))/2.)**2
+                    elif s_idx == 1: keep.value = np.cos(np.arcsin(np.sqrt(s2_2))/2.)**4
+                    elif s_idx == 2: keep.value = np.sin(np.arcsin(np.sqrt(s2_2))/2.)**2
+                    theta = mcmc_paramset_freq.values
+                    print 'mcmc_paramset_freq', mcmc_paramset_freq
+                    try:
+                        llh = llh_utils.triangle_llh(
+                            theta=theta, args=args, asimov_paramset=asimov_paramset,
+                            mcmc_paramset=mcmc_paramset_freq, fitter=fitter
+                        )
+                    except AssertionError:
+                        print 'forbidden by unitarity'
+                        data[idx][s_idx][a_idx][0] = np.nan
+                        data[idx][s_idx][a_idx][1] = np.nan
+                        data[idx][s_idx][a_idx][2] = np.nan
+                        continue
+                    print 'llh', llh
+                    data[idx][s_idx][a_idx][0] = sc
+                    data[idx][s_idx][a_idx][1] = an
+                    data[idx][s_idx][a_idx][2] = llh
+
+        misc_utils.make_dir(out)
+        print 'saving to {0}.npy'.format(out)
+        np.save(out+'.npy', np.array(data))
+
+    print "DONE!"
+
+
+main.__doc__ = __doc__
+
+
+if __name__ == '__main__':
+    main()
+
diff --git a/plot_bf.py b/plot_bf.py
index faefcde..13664b9 100644..100755
--- a/plot_bf.py
+++ b/plot_bf.py
@@ -24,7 +24,7 @@ from matplotlib.offsetbox import AnchoredText
 import fr
 from utils import misc as misc_utils
 from utils.fr import normalise_fr
-from utils.plot import myround, get_units
+from utils.plot import bayes_factor_plot, myround, get_units
 
 
 rc('text', usetex=False)
@@ -33,11 +33,12 @@ 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),
+    (1, 1, 1, 0, 1, 0),
 ]
 
 # FR
-dimension         = [3, 4, 5, 6, 7, 8]
+dimension         = [3, 6]
+# dimension         = [3, 4, 5, 6, 7, 8]
 sigma_ratio       = ['0.01']
 energy_dependance = 'spectral'
 spectral_index    = -2
@@ -68,15 +69,16 @@ bayes_eval_bin    = True # set to 'all' to run normally
 
 # Plot
 plot_bayes        = False
-plot_angles_limit = False
-plot_angles_corr  = True
+plot_angles_limit = True
+plot_angles_corr  = False
 outformat         = ['png']
-significance      = np.log(10**(1/2.))
-# significance      = np.log(10**(3/2.))
+# 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)
@@ -96,7 +98,7 @@ for i_dim, dim in enumerate(dimension):
             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'.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)
+            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)
@@ -110,18 +112,33 @@ for i_dim, dim in enumerate(dimension):
             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
-            for sc in scan_scales:
+            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,:]
-                    if plot_angles_corr:
-                        # TODO(shivesh)
-                        assert 0
-                        if len(lf.shape) == 3: lf = lf[:,0,:]
                     raw.append(lf)
                 except IOError:
                     fail += 1
@@ -141,13 +158,28 @@ for i_dim, dim in enumerate(dimension):
                 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:
-                # TODO(shivesh)
                 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))
@@ -160,7 +192,7 @@ if plot_bayes:
     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}$')
+    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
@@ -187,8 +219,11 @@ if plot_bayes:
 
             else:
                 print 'No points for DIM {0} FRS {1} NULL {2}!'.format(dim, frs, null)
-    yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True))
-    ax.set_ylim(yranges)
+                # 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():
@@ -197,7 +232,7 @@ if plot_bayes:
         ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1)
 
     for of in outformat:
-        fig.savefig('./images/bayes_factor.'+of, bbox_inches='tight', dpi=150)
+        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'}
@@ -210,15 +245,16 @@ if plot_angles_limit:
         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'$'
+        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[0]
+                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]
@@ -249,4 +285,130 @@ if plot_angles_limit:
             ax.axvline(x=xmaj, ls='-', color='gray', alpha=0.4, linewidth=1)
 
         for of in outformat:
-            fig.savefig('./images/angles_limit_DIM{0}'.format(dim)+'.'+of, bbox_inches='tight', dpi=150)
+            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/sens.py b/sens.py
new file mode 100755
index 0000000..578528c
--- /dev/null
+++ b/sens.py
@@ -0,0 +1,480 @@
+#! /usr/bin/env python
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : March 17, 2018
+
+"""
+HESE BSM flavour ratio analysis script
+"""
+
+from __future__ import absolute_import, division
+
+import os
+import argparse
+from functools import partial
+
+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
+from utils.enums import EnergyDependance, Likelihood, ParamTag
+from utils.enums import SensitivityCateg, StatCateg
+from utils.fr import estimate_scale, normalise_fr
+from utils.misc import enum_parse, parse_bool
+from utils.param import Param, ParamSet, get_paramsets
+
+from utils import multinest as mn_utils
+
+
+def define_nuisance():
+    """Define the nuisance parameters to scan over with default values,
+    ranges and sigma.
+    """
+    tag = ParamTag.NUISANCE
+    nuisance = ParamSet(
+        Param(name='convNorm',        value=1.,  seed=[0.5, 2. ], ranges=[0. , 50.], std=0.3,  tag=tag),
+        Param(name='promptNorm',      value=0.,  seed=[0. , 6. ], ranges=[0. , 50.], std=0.05, tag=tag),
+        Param(name='muonNorm',        value=1.,  seed=[0.1, 2. ], ranges=[0. , 50.], std=0.1,  tag=tag),
+        Param(name='astroNorm',       value=6.9, seed=[0.1, 10.], ranges=[0. , 50.], std=0.1,  tag=tag),
+        Param(name='astroDeltaGamma', value=2.5, seed=[1. , 3. ], ranges=[-5., 5. ], std=0.1,  tag=tag)
+    )
+    return nuisance
+
+
+def nuisance_argparse(parser):
+    nuisance_paramset = define_nuisance()
+    for parm in nuisance_paramset:
+        parser.add_argument(
+            '--'+parm.name, type=float, default=parm.value,
+            help=parm.name+' to inject'
+        )
+
+
+def process_args(args):
+    """Process the input args."""
+    if args.fix_mixing and args.fix_scale:
+        raise NotImplementedError('Fixed mixing and scale not implemented')
+    if args.fix_mixing and args.fix_mixing_almost:
+        raise NotImplementedError(
+            '--fix-mixing and --fix-mixing-almost cannot be used together'
+        )
+    if args.mn_run and args.fix_scale:
+        raise NotImplementedError(
+            '--mn-run and --fix-scale cannot be used together'
+        )
+
+    args.measured_ratio = normalise_fr(args.measured_ratio)
+    if args.fix_source_ratio:
+        args.source_ratio = normalise_fr(args.source_ratio)
+
+    if args.energy_dependance is EnergyDependance.SPECTRAL:
+        args.binning = np.logspace(
+            np.log10(args.binning[0]), np.log10(args.binning[1]), args.binning[2]+1
+        )
+
+    if not args.fix_scale:
+        args.scale = estimate_scale(args)
+        args.scale_region = (args.scale/args.scale_region, args.scale*args.scale_region)
+
+    if args.mn_eval_bin.lower() == 'all':
+        args.mn_eval_bin = None
+    else:
+        args.mn_eval_bin = int(args.mn_eval_bin)
+
+
+def parse_args(args=None):
+    """Parse command line arguments"""
+    parser = argparse.ArgumentParser(
+        description="BSM flavour ratio analysis",
+        formatter_class=misc_utils.SortingHelpFormatter,
+    )
+    parser.add_argument(
+        '--seed', type=misc_utils.seed_parse, default='25',
+        help='Set the random seed value'
+    )
+    parser.add_argument(
+        '--threads', type=misc_utils.thread_type, default='1',
+        help='Set the number of threads to use (int or "max")'
+    )
+    parser.add_argument(
+        '--outfile', type=str, default='./untitled',
+        help='Path to output chains'
+    )
+    parser.add_argument(
+        '--run-method', default='full',
+        type=partial(enum_parse, c=SensitivityCateg), choices=SensitivityCateg,
+        help='Choose which type of sensivity plot to make'
+    )
+    parser.add_argument(
+        '--stat-method', default='bayesian',
+        type=partial(enum_parse, c=StatCateg), choices=StatCateg,
+        help='Statistical method to employ'
+    )
+    fr_utils.fr_argparse(parser)
+    gf_utils.gf_argparse(parser)
+    llh_utils.likelihood_argparse(parser)
+    mn_utils.mn_argparse(parser)
+    nuisance_argparse(parser)
+    if args is None: return parser.parse_args()
+    else: return parser.parse_args(args.split())
+
+def main():
+    args = parse_args()
+    process_args(args)
+    misc_utils.print_args(args)
+
+    if args.seed is not None:
+        np.random.seed(args.seed)
+
+    asimov_paramset, llh_paramset = get_paramsets(args)
+    scale = llh_paramset.from_tag(ParamTag.SCALE)[0]
+    outfile = misc_utils.gen_outfile_name(args)
+    print '== {0:<25} = {1}'.format('outfile', outfile)
+
+    mmangles = llh_paramset.from_tag(ParamTag.MMANGLES)
+    if args.run_method is SensitivityCateg.FULL:
+        mn_paramset_arr = [llh_paramset.from_tag(ParamTag.SCALE, invert=True)]
+    elif args.run_method is SensitivityCateg.FIXED_ANGLE or \
+            args.run_method is SensitivityCateg.CORR_ANGLE:
+        nscale_pmset = llh_paramset.from_tag(ParamTag.SCALE, invert=True)
+        mn_paramset_arr = []
+        for x in xrange(3):
+            mn_paramset_arr.append(
+                ParamSet([prms for prms in nscale_pmset
+                          if mmangles[x].name != prms.name])
+            )
+
+    out = args.outfile+'/{0}/{1}/fr_evidence'.format(args.stat_method, args.run_method) \
+        + misc_utils.gen_identifier(args)
+    if args.mn_run:
+        if args.likelihood is Likelihood.GOLEMFIT:
+            fitter = gf_utils.setup_fitter(args, asimov_paramset)
+        else: fitter = None
+
+        scan_scales = np.linspace(
+            np.log10(args.scale_region[0]), np.log10(args.scale_region[1]), args.mn_bins
+        )
+        if args.run_method is SensitivityCateg.CORR_ANGLE:
+            scan_angles = np.linspace(0, 1, eval_dim)
+        else: scan_angles = np.array([0])
+        print 'scan_scales', scan_scales
+        print 'scan_angles', scan_angles
+
+        if args.mn_eval_bin is None:
+            eval_dim = args.mn_bins
+        else: eval_dim = 1
+
+        if args.run_method is SensitivityCateg.FULL:
+            evidence_arr = np.full((eval_dim, 2), np.nan)
+        elif args.run_method is SensitivityCateg.FIXED_ANGLE:
+            evidence_arr = np.full((len(mn_paramset_arr), eval_dim, 2), np.nan)
+        elif args.run_method is SensitivityCateg.CORR_ANGLE:
+            evidence_arr = np.full((len(mn_paramset_arr), eval_dim, eval_dim, 3), np.nan)
+
+
+        for idx_scen, mn_paramset in enumerate(mn_paramset_arr):
+            print '|||| SCENARIO = {0}'.format(idx_scen)
+            for x in mmangles: x.value = 0.
+            if args.run_method is SensitivityCateg.FIXED_ANGLE:
+                if idx_scen == 0 or idx_scen == 2:
+                    mmangles[idx_scen].value = np.sin(np.pi/2.)**2
+                    """s_12^2 or s_23^2"""
+                elif idx_scen == 1:
+                    mmangles[idx_scen].value = np.cos(np.pi/2.)**4
+                    """c_13^4"""
+
+            for idx_an, an in enumerate(scan_angles):
+                if args.run_method is SensitivityCateg.CORR_ANGLE:
+                    print '|||| ANGLE = {0:<04.2}'.format(an)
+                    mmangles[idx_an].value = an
+
+                for idx_sc, sc in enumerate(scan_scales):
+                    if args.mn_eval_bin is not None:
+                        if idx_sc == args.mn_eval_bin:
+                            out += '_scale_{0:.0E}'.format(np.power(10, sc))
+                            if args.run_method is SensitivityCateg.CORR_ANGLE:
+                                out += '_angle_{0:>03}'.format(int(an*100))
+                            break
+                        else: continue
+
+                    print '|||| SCALE = {0:.0E}'.format(np.power(10, sc))
+                    scale.value = sc
+                    try:
+                        a_lnZ = mn_utils.mn_evidence(
+                            mn_paramset     = mn_paramset,
+                            llh_paramset    = llh_paramset,
+                            asimov_paramset = asimov_paramset,
+                            args            = args,
+                            fitter          = fitter
+                        )
+                    except:
+                        print 'Failed run, continuing'
+                        continue
+                    print '## Evidence = {0}'.format(a_lnZ)
+                    if args.run_method is SensitivityCateg.FULL:
+                        evidence_arr[idx_sc] = np.array([sc, a_lnZ])
+                    elif args.run_method is SensitivityCateg.FIXED_ANGLE:
+                        evidence_arr[idx_scen][idx_sc] = np.array([sc, a_lnZ])
+                    elif args.run_method is SensitivityCateg.CORR_ANGLE:
+                        evidence_arr[idx_scen][idx_an][idx_sc] = np.array([an, sc, a_lnZ])
+
+        misc_utils.make_dir(out)
+        print 'Saving to {0}'.format(out+'.npy')
+        np.save(out+'.npy', np.array(evidence_arr))
+
+    if args.plot_multinest:
+        if args.mn_run: raw = evidence_arr
+        else: raw = np.load(out+'.npy')
+        data = ma.masked_invalid(raw, 0)
+
+        basename = os.path.dirname(out) + '/mnrun/' + os.path.basename(out)
+        baseoutfile = basename[:5]+basename[5:].replace('data', 'plots')
+        if args.run_method is SensitivityCateg.FULL:
+            plot_utils.plot_multinest(
+                data        = data,
+                outfile     = baseoutfile,
+                outformat   = ['png'],
+                args        = args,
+                scale_param = scale
+            )
+        elif args.run_method is SensitivityCateg.FIXED_ANGLE:
+            for idx_scen in xrange(len(mn_paramset_arr)):
+                print '|||| SCENARIO = {0}'.format(idx_scen)
+                outfile = baseoutfile + '_SCEN{0}'.format(idx_scen)
+                if idx_scen == 0: label = r'$\mathcal{O}_{12}=\frac{1}{2}$'
+                elif idx_scen == 1: label = r'$\mathcal{O}_{13}=\frac{1}{2}$'
+                elif idx_scen == 2: label = r'$\mathcal{O}_{23}=\frac{1}{2}$'
+                plot_utils.plot_multinest(
+                    data        = data[idx_scen],
+                    outfile     = outfile,
+                    outformat   = ['png'],
+                    args        = args,
+                    scale_param = scale,
+                    label       = label
+                )
+        elif args.run_method is SensitivityCateg.CORR_ANGLE:
+            for idx_scen in xrange(len(mn_paramset_arr)):
+                print '|||| SCENARIO = {0}'.format(idx_scen)
+                basescenoutfile = baseoutfile + '_SCEN{0}'.format(idx_scen)
+                if idx_scen == 0: label = r'$\mathcal{O}_{12}='
+                elif idx_scen == 1: label = r'$\mathcal{O}_{13}='
+                elif idx_scen == 2: label = r'$\mathcal{O}_{23}='
+                for idx_an, an in enumerate(scan_angles):
+                    print '|||| ANGLE = {0:<04.2}'.format(an)
+                    outfile = basescenoutfile + '_ANGLE{0}'.format(idx_an)
+                    label += r'{0:<04.2}$'.format(an)
+                    plot_utils.plot_multinest(
+                        data        = data[idx_scen][idx_an][:,1:],
+                        outfile     = outfile,
+                        outformat   = ['png'],
+                        args        = args,
+                        scale_param = scale,
+                        label       = label
+                    )
+
+    # dirname = os.path.dirname(out)
+    # plot_utils.bayes_factor_plot(
+    #     dirname=dirname, outfile=out, outformat=['png'], args=args
+    # )
+
+    # out = args.angles_lim_output+'/fr_an_evidence' + misc_utils.gen_identifier(args)
+    # if args.run_angles_limit:
+    #     import pymultinest
+
+    #     scenarios = [
+    #         [np.sin(np.pi/2.)**2, 0, 0, 0],
+    #         [0, np.cos(np.pi/2.)**4, 0, 0],
+    #         [0, 0, np.sin(np.pi/2.)**2, 0],
+    #     ]
+    #     p = llh_paramset.from_tag([ParamTag.SCALE, ParamTag.MMANGLES], invert=True)
+    #     n_params = len(p)
+    #     prior_ranges = p.seeds
+
+    #     if not args.run_llh and args.likelihood is Likelihood.GOLEMFIT:
+    #         fitter = gf_utils.setup_fitter(args, asimov_paramset)
+    #     else: fitter = None
+
+    #     def CubePrior(cube, ndim, nparams):
+    #         # default are uniform priors
+    #         return ;
+
+    #     if args.bayes_eval_bin is not None:
+    #         data = np.zeros((len(scenarios), 1, 2))
+    #     else: data = np.zeros((len(scenarios), args.bayes_bins, 2))
+    #     mm_angles = llh_paramset.from_tag(ParamTag.MMANGLES)
+    #     sc_angles = llh_paramset.from_tag(ParamTag.SCALE)[0]
+    #     for idx, scen in enumerate(scenarios):
+    #         scales, evidences = [], []
+    #         for yidx, an in enumerate(mm_angles):
+    #             an.value = scen[yidx]
+    #         for s_idx, sc in enumerate(scan_scales):
+    #             if args.bayes_eval_bin is not None:
+    #                 if s_idx in args.bayes_eval_bin:
+    #                     if idx == 0:
+    #                         out += '_scale_{0:.0E}'.format(np.power(10, sc))
+    #                 else: continue
+
+    #             print '== SCALE = {0:.0E}'.format(np.power(10, sc))
+    #             sc_angles.value = sc
+    #             def lnProb(cube, ndim, nparams):
+    #                 for i in range(ndim):
+    #                     prange = prior_ranges[i][1] - prior_ranges[i][0]
+    #                     p[i].value = prange*cube[i] + prior_ranges[i][0]
+    #                 for name in p.names:
+    #                     mcmc_paramset[name].value = p[name].value
+    #                 theta = mcmc_paramset.values
+    #                 # print 'theta', theta
+    #                 # print 'mcmc_paramset', mcmc_paramset
+    #                 return llh_utils.triangle_llh(
+    #                     theta=theta,
+    #                     args=args,
+    #                     asimov_paramset=asimov_paramset,
+    #                     mcmc_paramset=mcmc_paramset,
+    #                     fitter=fitter
+    #                 )
+    #             prefix = 'mnrun/' + os.path.basename(out) + '_'
+    #             if args.bayes_eval_bin is not None:
+    #                 prefix += '{0}_{1}_'.format(idx, s_idx)
+    #             print 'begin running evidence calculation for {0}'.format(prefix)
+    #             result = pymultinest.run(
+    #                 LogLikelihood=lnProb,
+    #                 Prior=CubePrior,
+    #                 n_dims=n_params,
+    #                 importance_nested_sampling=True,
+    #                 n_live_points=args.bayes_live_points,
+    #                 evidence_tolerance=args.bayes_tolerance,
+    #                 outputfiles_basename=prefix,
+    #                 resume=False,
+    #                 verbose=True
+    #             )
+
+    #             analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params)
+    #             a_lnZ = analyzer.get_stats()['global evidence']
+    #             print 'Evidence = {0}'.format(a_lnZ)
+    #             scales.append(sc)
+    #             evidences.append(a_lnZ)
+
+    #         for i, d in enumerate(evidences):
+    #             data[idx][i][0] = scales[i]
+    #             data[idx][i][1] = d
+
+    #     misc_utils.make_dir(out)
+    #     print 'saving to {0}.npy'.format(out)
+    #     np.save(out+'.npy', np.array(data))
+
+    # dirname = os.path.dirname(out)
+    # plot_utils.plot_BSM_angles_limit(
+    #     dirname=dirname, outfile=outfile, outformat=['png'],
+    #     args=args, bayesian=True
+    # )
+
+    # out = args.angles_corr_output+'/fr_co_evidence' + misc_utils.gen_identifier(args)
+    # if args.run_angles_correlation:
+    #     if args.bayes_eval_bin is None: assert 0
+    #     import pymultinest
+
+    #     scenarios = [
+    #         [1, 0, 0, 0],
+    #         [0, 1, 0, 0],
+    #         [0, 0, 1, 0],
+    #     ]
+    #     nuisance = mcmc_paramset.from_tag(ParamTag.NUISANCE)
+    #     mm_angles = mcmc_paramset.from_tag(ParamTag.MMANGLES)
+    #     sc_angles = mcmc_paramset.from_tag(ParamTag.SCALE)[0]
+
+    #     if not args.run_mcmc and args.likelihood is Likelihood.GOLEMFIT:
+    #         fitter = gf_utils.setup_fitter(args, asimov_paramset)
+    #     else: fitter = None
+
+    #     def CubePrior(cube, ndim, nparams):
+    #         # default are uniform priors
+    #         return ;
+
+    #     scan_angles = np.linspace(0, 1, args.bayes_bins)
+
+    #     if args.bayes_eval_bin is not None:
+    #         data = np.zeros((len(scenarios), 1, 1, 3))
+    #     else: data = np.zeros((len(scenarios), args.bayes_bins, args.bayes_bins, 3))
+    #     for idx, scen in enumerate(scenarios):
+    #         for an in mm_angles:
+    #             an.value = 0
+    #         keep = mcmc_paramset.from_tag(ParamTag.MMANGLES)[idx]
+    #         q = ParamSet(nuisance.params + [x for x in mm_angles if x.name != keep.name])
+    #         n_params = len(q)
+    #         prior_ranges = q.seeds
+
+    #         scales, angles, evidences = [], [], []
+    #         for s_idx, sc in enumerate(scan_scales):
+    #             for a_idx, an in enumerate(scan_angles):
+    #                 index = s_idx*args.bayes_bins + a_idx
+    #                 if args.bayes_eval_bin is not None:
+    #                     if index in args.bayes_eval_bin:
+    #                         if idx == 0:
+    #                             out += '_idx_{0}'.format(index)
+    #                     else: continue
+
+    #                 print '== SCALE = {0:.0E}'.format(np.power(10, sc))
+    #                 print '== ANGLE = {0:.2f}'.format(an)
+    #                 sc_angles.value = sc
+    #                 keep.value = an
+    #                 def lnProb(cube, ndim, nparams):
+    #                     for i in range(ndim-1):
+    #                         prange = prior_ranges[i][1] - prior_ranges[i][0]
+    #                         q[i].value = prange*cube[i] + prior_ranges[i][0]
+    #                     for name in q.names:
+    #                         mcmc_paramset[name].value = q[name].value
+    #                     theta = mcmc_paramset.values
+    #                     # print 'theta', theta
+    #                     # print 'mcmc_paramset', mcmc_paramset
+    #                     return llh_utils.triangle_llh(
+    #                         theta=theta,
+    #                         args=args,
+    #                         asimov_paramset=asimov_paramset,
+    #                         mcmc_paramset=mcmc_paramset,
+    #                         fitter=fitter
+    #                     )
+    #                 prefix = 'mnrun/' + os.path.basename(out) + '_'
+    #                 if args.bayes_eval_bin is not None:
+    #                     prefix += '{0}_{1}_{2}'.format(idx, s_idx, a_idx)
+
+    #                 print 'begin running evidence calculation for {0}'.format(prefix)
+    #                 result = pymultinest.run(
+    #                     LogLikelihood=lnProb,
+    #                     Prior=CubePrior,
+    #                     n_dims=n_params,
+    #                     importance_nested_sampling=True,
+    #                     n_live_points=args.bayes_live_points,
+    #                     evidence_tolerance=args.bayes_tolerance,
+    #                     outputfiles_basename=prefix,
+    #                     resume=False,
+    #                     verbose=True
+    #                 )
+
+    #                 analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params)
+    #                 a_lnZ = analyzer.get_stats()['global evidence']
+    #                 print 'Evidence = {0}'.format(a_lnZ)
+    #                 scales.append(sc)
+    #                 angles.append(an)
+    #                 evidences.append(a_lnZ)
+
+    #             for i, d in enumerate(evidences):
+    #                 data[idx][i][i][0] = scales[i]
+    #                 data[idx][i][i][1] = angles[i]
+    #                 data[idx][i][i][2] = d
+
+    #     misc_utils.make_dir(out)
+    #     print 'saving to {0}.npy'.format(out)
+    #     np.save(out+'.npy', np.array(data))
+
+
+main.__doc__ = __doc__
+
+
+if __name__ == '__main__':
+    main()
+
diff --git a/submitter/make_dag.py b/submitter/make_dag.py
index 78b7bff..c6e7400 100644
--- a/submitter/make_dag.py
+++ b/submitter/make_dag.py
@@ -19,7 +19,7 @@ fix_sfr_mfr = [
     (1, 1, 1, 1, 2, 0),
     # (1, 1, 0, 1, 2, 0),
     # (1, 2, 0, 1, 2, 0),
-    (1, 1, 1, 1, 0, 0),
+    # (1, 1, 1, 1, 0, 0),
     # (1, 1, 0, 1, 0, 0),
     # (1, 0, 0, 1, 0, 0),
     (1, 1, 1, 0, 1, 0),
@@ -34,12 +34,12 @@ run_mcmc = 'False'
 burnin   = 2500
 nsteps   = 10000
 nwalkers = 60
-seed     = 24
+seed     = 'None'
 threads  = 4
 mcmc_seed_type = 'uniform'
 
 # FR
-dimension         = [6]
+dimension         = [4, 5, 7, 8]
 energy            = [1e6]
 no_bsm            = 'False'
 sigma_ratio       = ['0.01']
@@ -66,13 +66,13 @@ ast  = 'p2_0'
 data = 'real'
 
 # Bayes Factor
-run_bayes_factor       = 'False'
+run_bayes_factor       = 'True'
 run_angles_limit       = 'False'
-run_angles_correlation = 'True'
-bayes_bins             = 20
-bayes_live_points      = 1000
+run_angles_correlation = 'False'
+bayes_bins             = 100
+bayes_live_points      = 3000
 bayes_tolerance        = 0.01
-bayes_eval_bin         = 'None' # set to 'all' to run normally
+bayes_eval_bin         = 'all' # set to 'all' to run normally
 
 # Plot
 plot_angles       = 'False'
@@ -80,8 +80,8 @@ plot_elements     = 'False'
 plot_bayes        = 'False'
 plot_angles_limit = 'False'
 
-# outfile = 'dagman_FR.submit'.format(dimension[0])
-outfile = 'dagman_FR_angles_correlation_DIM{0}.submit'.format(dimension[0])
+outfile = 'dagman_FR_freq_fullscan_otherdims.submit'
+# outfile = 'dagman_FR_bayes_freq.submit'
 golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
 condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/submit.sub'
 
@@ -106,6 +106,7 @@ with open(outfile, 'w') as f:
 
             bayes_output = 'None'
             angles_lim_output = 'None'
+            angles_corr_output = 'None'
             for sig in sigma_ratio:
                 print 'sigma', sig
                 for frs in fix_sfr_mfr:
diff --git a/test/LV.out b/test/LV.out
deleted file mode 100644
index e9f2b7e..0000000
--- a/test/LV.out
+++ /dev/null
@@ -1,74 +0,0 @@
-test_LV.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 2u, std::allocator<double> > already registered; second conversion method ignored.
-  import GolemFitPy as gf
-test_LV.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 5u, std::allocator<double> > already registered; second conversion method ignored.
-  import GolemFitPy as gf
-reset_steering: 1
-reset_data: 1
-reset_npp: 1
-GolemFit constructor: checking paths
-Constructing DiffuseWeightMaker
-Loading data
-Loading HESE data.
-Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/data/HESEData.txt
-Loaded HESE 80 events.
-Loading MC
-Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 into main simulation
-Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5
-Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 now main simulation size is 420660
-Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 into main simulation
-Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5
-Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 now main simulation size is 616803
-Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 into main simulation
-Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5
-Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 now main simulation size is 740041
-Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 into main simulation
-Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5
-Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 now main simulation size is 740161
-Loaded 740161 events in main simulation set
-Loading Flux weighter
-Loading nusquids objects
-Weighting MC
-Initializing simulation weights
-Applying selfveto correction.
-HESE reshuffle has not been performed.
-Doing HESE reshuffle
-Making data hist
-Making sim hist
-Constructing likelihood problem with default settings
-Setting up Asimov set with parameters
-Conventional normalization 1
-Prompt normalization 1
-Atmospheric muon normalization 1
-Astro component overall normalization 8
-Astro E component 0.333333
-Astro Mu component 0.333333
-Astro Tau component 0.333333
-Astroparticle balance 1
-Astro gamma 2.5
-Astro cutoff 10
-CR gamma 0
-Conv pi/k ratio 1
-Conv particle balance 1
-Conv zenith correction 0
-Dark normalization 0
-DOM efficiency 0.99
-Astro component overall normalization second 0
-Astro gamma second 2
-Holeice forward 0
-Anisotropy scale1
-
-Remaking data hist
-Reconstrucing likelihood problem
-NULL min_llh 835.848062048
-NULL expectation [11.96478761 11.90169848 10.78840154  9.20158936  7.23605538  5.59952977
-  4.39850777  3.35653018  2.58491686  1.92834114  1.46813946  1.11581282
-  0.8492433   0.65154743  0.46727371  0.36211812  0.32718588  0.69746002
-  0.29475116  0.05150809]
-
-Applying new flavor physics to MC weights.
-Entering ApplyNewPhysics
-Entering ConstructNuSQuIDSObjectsForLV.
-Rotated Unitary Transform
-In for loop
-Entering inner for loop
-Leaving inner for loop
diff --git a/test/NSI.out b/test/NSI.out
deleted file mode 100644
index 657a9f3..0000000
--- a/test/NSI.out
+++ /dev/null
@@ -1,75 +0,0 @@
-test_NSI.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 2u, std::allocator<double> > already registered; second conversion method ignored.
-  import GolemFitPy as gf
-test_NSI.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 5u, std::allocator<double> > already registered; second conversion method ignored.
-  import GolemFitPy as gf
-reset_steering: 1
-reset_data: 1
-reset_npp: 1
-GolemFit constructor: checking paths
-Constructing DiffuseWeightMaker
-Loading data
-Loading HESE data.
-Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/data/HESEData.txt
-Loaded HESE 80 events.
-Loading MC
-Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 into main simulation
-Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5
-Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 now main simulation size is 420660
-Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 into main simulation
-Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5
-Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 now main simulation size is 616803
-Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 into main simulation
-Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5
-Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 now main simulation size is 740041
-Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 into main simulation
-Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5
-Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 now main simulation size is 740161
-Loaded 740161 events in main simulation set
-Loading Flux weighter
-Loading nusquids objects
-Weighting MC
-Initializing simulation weights
-Applying selfveto correction.
-HESE reshuffle has not been performed.
-Doing HESE reshuffle
-Making data hist
-Making sim hist
-Constructing likelihood problem with default settings
-Setting up Asimov set with parameters
-Conventional normalization 1
-Prompt normalization 1
-Atmospheric muon normalization 1
-Astro component overall normalization 8
-Astro E component 0.333333
-Astro Mu component 0.333333
-Astro Tau component 0.333333
-Astroparticle balance 1
-Astro gamma 2.5
-Astro cutoff 10
-CR gamma 0
-Conv pi/k ratio 1
-Conv particle balance 1
-Conv zenith correction 0
-Dark normalization 0
-DOM efficiency 0.99
-Astro component overall normalization second 0
-Astro gamma second 2
-Holeice forward 0
-Anisotropy scale1
-
-Remaking data hist
-Reconstrucing likelihood problem
-NULL min_llh 835.848062048
-NULL expectation [11.96478761 11.90169848 10.78840154  9.20158936  7.23605538  5.59952977
-  4.39850777  3.35653018  2.58491686  1.92834114  1.46813946  1.11581282
-  0.8492433   0.65154743  0.46727371  0.36211812  0.32718588  0.69746002
-  0.29475116  0.05150809]
-
-Applying new flavor physics to MC weights.
-Applying New Physics in normal mode to NonStandardInteraction
-HESE reshuffle has been performed.
-0.1 mutau min_llh 836.352921667
-0.1 mutau expectation [16.63858821 16.04604012 14.28034835 11.80228198  9.06453127  6.8815312
-  5.32704912  4.02384379  3.07832186  2.28036523  1.73329938  1.30341591
-  0.99239868  0.75237435  0.53312845  0.4093463   0.35723267  0.71433954
-  0.30173276  0.05276177]
diff --git a/test/test.png b/test/test.png
deleted file mode 100644
index 8aac9e3..0000000
--- a/test/test.png
+++ /dev/null
diff --git a/test/test_LV.py b/test/test_LV.py
deleted file mode 100644
index 96a1863..0000000
--- a/test/test_LV.py
+++ /dev/null
@@ -1,95 +0,0 @@
-#!/usr/bin/env python
-
-from __future__ import absolute_import, division
-
-import numpy as np
-import matplotlib
-matplotlib.use('Agg')
-import matplotlib.pyplot as plt
-from matplotlib import rc
-
-import GolemFitPy as gf
-
-rc('text', usetex=True)
-rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
-
-dp = gf.DataPaths()
-steer = gf.SteeringParams()
-
-fig = plt.figure(figsize=[12, 8])
-ax = fig.add_subplot(111)
-ax.set_xscale('log')
-ax.set_yscale('log')
-
-npp = gf.NewPhysicsParams()
-npp.type = gf.NewPhysicsType.None
-# npp.n_lv = 1
-# npp.lambda_1 = 1.e100
-# npp.lambda_2 = 1.e100
-
-golem = gf.GolemFit(dp, steer, npp)
-
-binning = golem.GetEnergyBinsMC()
-ax.set_xlim(binning[0], binning[-1])
-# ax.set_ylim(binning[0], binning[-1])
-
-fit_params = gf.FitParameters(gf.sampleTag.HESE)
-golem.SetupAsimov(fit_params)
-
-exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
-ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
-        drawstyle='steps-pre', label='NULL', linestyle='--')
-
-print 'NULL min_llh', golem.MinLLH().likelihood
-print 'NULL expectation', exp
-print
-
-npp = gf.NewPhysicsParams()
-npp.type = gf.NewPhysicsType.LorentzViolation
-npp.n_lv = 1
-npp.lambda_1 = 1.e20
-npp.lambda_2 = 1.e20
-
-golem.SetNewPhysicsParams(npp)
-
-exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
-ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
-        drawstyle='steps-pre', label='1e-20 LV', linestyle='--')
-
-print '1e20 LV min_llh', golem.MinLLH().likelihood
-print '1e20 LV expectation', exp
-
-npp = gf.NewPhysicsParams()
-npp.type = gf.NewPhysicsType.LorentzViolation
-npp.n_lv = 1
-npp.lambda_1 = 1.e10
-npp.lambda_2 = 1.e10
-
-golem.SetNewPhysicsParams(npp)
-
-exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
-ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
-        drawstyle='steps-pre', label='1e-10 LV', linestyle='--')
-
-print '1e10 LV min_llh', golem.MinLLH().likelihood
-print '1e10 LV expectation', exp
-
-npp = gf.NewPhysicsParams()
-npp.type = gf.NewPhysicsType.LorentzViolation
-npp.n_lv = 1
-npp.lambda_1 = 1.e-20
-npp.lambda_2 = 1.e-20
-
-golem.SetNewPhysicsParams(npp)
-
-ax.tick_params(axis='x', labelsize=12)
-ax.tick_params(axis='y', labelsize=12)
-ax.set_xlabel(r'Deposited energy / GeV')
-ax.set_ylabel(r'Events')
-for xmaj in ax.xaxis.get_majorticklocs():
-    ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.7, linewidth=1)
-for ymaj in ax.yaxis.get_majorticklocs():
-    ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.7, linewidth=1)
-
-legend = ax.legend(prop=dict(size=12))
-fig.savefig('test.png', bbox_inches='tight', dpi=250)
diff --git a/test/test_NSI.png b/test/test_NSI.png
deleted file mode 100644
index 000ea61..0000000
--- a/test/test_NSI.png
+++ /dev/null
diff --git a/test/test_NSI.py b/test/test_NSI.py
deleted file mode 100644
index 9b49c93..0000000
--- a/test/test_NSI.py
+++ /dev/null
@@ -1,103 +0,0 @@
-#!/usr/bin/env python
-
-from __future__ import absolute_import, division
-
-import numpy as np
-import matplotlib
-matplotlib.use('Agg')
-import matplotlib.pyplot as plt
-from matplotlib import rc
-
-import GolemFitPy as gf
-
-rc('text', usetex=True)
-rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
-
-dp = gf.DataPaths()
-steer = gf.SteeringParams()
-
-fig = plt.figure(figsize=[12, 8])
-ax = fig.add_subplot(111)
-ax.set_xscale('log')
-ax.set_yscale('log')
-
-npp = gf.NewPhysicsParams()
-npp.type = gf.NewPhysicsType.None
-# npp.epsilon_mutau = 0
-# npp.epsilon_prime = 0
-
-golem = gf.GolemFit(dp, steer, npp)
-
-binning = golem.GetEnergyBinsMC()
-ax.set_xlim(binning[0], binning[-1])
-# ax.set_ylim(binning[0], binning[-1])
-
-fit_params = gf.FitParameters(gf.sampleTag.HESE)
-golem.SetupAsimov(fit_params)
-
-exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
-ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
-        drawstyle='steps-pre', label='NULL', linestyle='--')
-
-print 'NULL min_llh', golem.MinLLH().likelihood
-print 'NULL expectation', exp
-print
-
-npp = gf.NewPhysicsParams()
-npp.type = gf.NewPhysicsType.NonStandardInteraction
-npp.epsilon_mutau = 0.1
-# npp.epsilon_prime = 0
-
-golem.SetNewPhysicsParams(npp)
-
-exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
-ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
-        drawstyle='steps-pre', label='0.1 mutau', linestyle='--')
-
-print '0.1 mutau min_llh', golem.MinLLH().likelihood
-print '0.1 mutau expectation', exp
-
-np.epsilon_mutau = 0.2
-
-golem.SetNewPhysicsParams(npp)
-
-exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
-ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
-        drawstyle='steps-pre', label='0.2 mutau', linestyle='--')
-
-print '0.2 mutau min_llh', golem.MinLLH().likelihood
-print '0.2 mutau expectation', exp
-
-np.epsilon_mutau = 0.3
-
-golem.SetNewPhysicsParams(npp)
-
-exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
-ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
-        drawstyle='steps-pre', label='0.3 mutau', linestyle='--')
-
-print '0.3 mutau min_llh', golem.MinLLH().likelihood
-print '0.3 mutau expectation', exp
-
-np.epsilon_mutau = 0.4
-
-golem.SetNewPhysicsParams(npp)
-
-exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
-ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
-        drawstyle='steps-pre', label='0.4 mutau', linestyle='--')
-
-print '0.4 mutau min_llh', golem.MinLLH().likelihood
-print '0.4 mutau expectation', exp
-
-ax.tick_params(axis='x', labelsize=12)
-ax.tick_params(axis='y', labelsize=12)
-ax.set_xlabel(r'Deposited energy / GeV')
-ax.set_ylabel(r'Events')
-for xmaj in ax.xaxis.get_majorticklocs():
-    ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.7, linewidth=1)
-for ymaj in ax.yaxis.get_majorticklocs():
-    ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.7, linewidth=1)
-
-legend = ax.legend(prop=dict(size=12))
-fig.savefig('test_NSI.png', bbox_inches='tight', dpi=250)
diff --git a/utils/enums.py b/utils/enums.py
index 90ca17d..359e104 100644
--- a/utils/enums.py
+++ b/utils/enums.py
@@ -47,6 +47,17 @@ class MCMCSeedType(Enum):
     GAUSSIAN = 2
 
 
+class SensitivityCateg(Enum):
+    FULL        = 1
+    FIXED_ANGLE = 2
+    CORR_ANGLE  = 3
+
+
+class StatCateg(Enum):
+    BAYESIAN    = 1
+    FREQUENTIST = 2
+
+
 class SteeringCateg(Enum):
     P2_0 = 1
     P2_1 = 2
diff --git a/utils/fr.py b/utils/fr.py
index 8f71f78..6a405a7 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -9,11 +9,15 @@ Useful functions for the BSM flavour ratio analysis
 
 from __future__ import absolute_import, division
 
-import sys
+import argparse
+from functools import partial
 
 import numpy as np
 from scipy import linalg
 
+from utils.enums import EnergyDependance
+from utils.misc import enum_parse, parse_bool
+
 
 MASS_EIGENVALUES = [7.40E-23, 2.515E-21]
 """SM mass eigenvalues"""
@@ -194,6 +198,83 @@ def normalise_fr(fr):
     return np.array(fr) / float(np.sum(fr))
 
 
+def estimate_scale(args, mass_eigenvalues=MASS_EIGENVALUES):
+    """Estimate the scale at which new physics will enter."""
+    if args.energy_dependance is EnergyDependance.MONO:
+        scale = np.power(
+            10, np.round(np.log10(MASS_EIGENVALUES[1]/args.energy)) - \
+            np.log10(args.energy**(args.dimension-3))
+        )
+    elif args.energy_dependance is EnergyDependance.SPECTRAL:
+        scale = np.power(
+            10, np.round(
+                np.log10(MASS_EIGENVALUES[1]/np.power(10, np.average(np.log10(args.binning)))) \
+                - np.log10(np.power(10, np.average(np.log10(args.binning)))**(args.dimension-3))
+            )
+        )
+    return scale
+
+
+def fr_argparse(parser):
+    parser.add_argument(
+        '--measured-ratio', type=int, nargs=3, default=[1, 1, 1],
+        help='Set the central value for the measured flavour ratio at IceCube'
+    )
+    parser.add_argument(
+        '--source-ratio', type=int, nargs=3, default=[2, 1, 0],
+        help='Set the source flavour ratio for the case when you want to fix it'
+    )
+    parser.add_argument(
+        '--no-bsm', type=parse_bool, default='False',
+        help='Turn off BSM terms'
+    )
+    parser.add_argument(
+        '--dimension', type=int, default=3,
+        help='Set the new physics dimension to consider'
+    )
+    parser.add_argument(
+        '--energy-dependance', default='spectral', type=partial(enum_parse, c=EnergyDependance),
+        choices=EnergyDependance,
+        help='Type of energy dependance to use in the BSM fit'
+    )
+    parser.add_argument(
+        '--spectral-index', default=-2, type=int,
+        help='Spectral index for spectral energy dependance'
+    )
+    parser.add_argument(
+        '--binning', default=[1e4, 1e7, 5], type=float, nargs=3,
+        help='Binning for spectral energy dependance'
+    )
+    parser.add_argument(
+        '--fix-source-ratio', type=parse_bool, default='False',
+        help='Fix the source flavour ratio'
+    )
+    parser.add_argument(
+        '--fix-mixing', type=parse_bool, default='False',
+        help='Fix all mixing parameters to bi-maximal mixing'
+    )
+    parser.add_argument(
+        '--fix-mixing-almost', type=parse_bool, default='False',
+        help='Fix all mixing parameters except s23'
+    )
+    parser.add_argument(
+        '--fix-scale', type=parse_bool, default='False',
+        help='Fix the new physics scale'
+    )
+    parser.add_argument(
+        '--scale', type=float, default=1e-30,
+        help='Set the new physics scale'
+    )
+    parser.add_argument(
+        '--scale-region', type=float, default=1e10,
+        help='Set the size of the box to scan for the scale'
+    )
+    parser.add_argument(
+        '--energy', type=float, default=1000,
+        help='Set the energy scale'
+    )
+
+
 def fr_to_angles(ratios):
     """Convert from flavour ratio into the angular projection of the flavour
     ratios.
@@ -218,7 +299,7 @@ NUFIT_U = angles_to_u((0.307, (1-0.02195)**2, 0.565, 3.97935))
 def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
                    nufit_u=NUFIT_U, no_bsm=False, fix_mixing=False,
                    fix_mixing_almost=False, fix_scale=False, scale=None,
-                   check_uni=True):
+                   check_uni=True, epsilon=1e-9):
     """Construct the BSM mixing matrix from the BSM parameters.
 
     Parameters
@@ -311,8 +392,8 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
 
     if check_uni:
         tu = test_unitarity(eg_vector)
-        if not abs(np.trace(tu) - 3.) < 1e-5 or \
-           not abs(np.sum(tu) - 3.) < 1e-5:
+        if not abs(np.trace(tu) - 3.) < epsilon or \
+           not abs(np.sum(tu) - 3.) < epsilon:
             raise AssertionError(
                 'Matrix is not unitary!\neg_vector\n{0}\ntest '
                 'u\n{1}'.format(eg_vector, tu)
@@ -378,4 +459,3 @@ def u_to_fr(source_fr, matrix):
     )
     ratio = composition / np.sum(source_fr)
     return ratio
-
diff --git a/utils/gf.py b/utils/gf.py
index 3f770eb..20afc75 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -9,7 +9,6 @@ Useful GolemFit wrappers for the BSM flavour ratio analysis
 
 from __future__ import absolute_import, division
 
-import argparse
 import socket
 from functools import partial
 
diff --git a/utils/likelihood.py b/utils/likelihood.py
index c1208ab..1981c72 100644
--- a/utils/likelihood.py
+++ b/utils/likelihood.py
@@ -9,7 +9,6 @@ Likelihood functions for the BSM flavour ratio analysis
 
 from __future__ import absolute_import, division
 
-import argparse
 from functools import partial
 
 import numpy as np
@@ -34,6 +33,10 @@ def likelihood_argparse(parser):
         '--likelihood', default='gaussian', type=partial(enum_parse, c=Likelihood),
         choices=Likelihood, help='likelihood contour'
     )
+    parser.add_argument(
+        '--sigma-ratio', type=float, default=0.01,
+        help='Set the 1 sigma for the measured flavour ratio for a gaussian LLH'
+    )
 
 
 def lnprior(theta, paramset):
@@ -46,17 +49,17 @@ def lnprior(theta, paramset):
     return 0.
 
 
-def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter):
+def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
     """-Log likelihood function for a given theta."""
-    if len(theta) != len(mcmc_paramset):
+    if len(theta) != len(llh_paramset):
         raise AssertionError(
             'Length of MCMC scan is not the same as the input '
-            'params\ntheta={0}\nmcmc_paramset]{1}'.format(theta, mcmc_paramset)
+            'params\ntheta={0}\nmcmc_paramset]{1}'.format(theta, llh_paramset)
         )
-    for idx, param in enumerate(mcmc_paramset):
+    for idx, param in enumerate(llh_paramset):
         param.value = theta[idx]
     hypo_paramset = asimov_paramset
-    for param in mcmc_paramset.from_tag(ParamTag.NUISANCE):
+    for param in llh_paramset.from_tag(ParamTag.NUISANCE):
         hypo_paramset[param.name].value = param.value
 
     if args.energy_dependance is EnergyDependance.SPECTRAL:
@@ -67,14 +70,14 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter):
         if args.energy_dependance is EnergyDependance.MONO:
             source_flux = args.source_ratio
         elif args.energy_dependance is EnergyDependance.SPECTRAL:
-	    source_flux = np.array(
-		[fr * np.power(bin_centers, args.spectral_index)
-		 for fr in args.source_ratio]
-	    ).T
+            source_flux = np.array(
+                [fr * np.power(bin_centers, args.spectral_index)
+                 for fr in args.source_ratio]
+            ).T
     else:
         if args.energy_dependance is EnergyDependance.MONO:
             source_flux = fr_utils.angles_to_fr(
-                mcmc_paramset.from_tag(ParamTag.SRCANGLES, values=True)
+                llh_paramset.from_tag(ParamTag.SRCANGLES, values=True)
             )
         elif args.energy_dependance is EnergyDependance.SPECTRAL:
             source_flux = np.array(
@@ -82,7 +85,7 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter):
                  for fr in fr_utils.angles_to_fr(theta[-2:])]
             ).T
 
-    bsm_angles = mcmc_paramset.from_tag(
+    bsm_angles = llh_paramset.from_tag(
         [ParamTag.SCALE, ParamTag.MMANGLES], values=True
     )
 
@@ -134,11 +137,11 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter):
         return gf_utils.get_llh_freq(fitter, hypo_paramset)
 
 
-def ln_prob(theta, args, fitter, asimov_paramset, mcmc_paramset):
-    lp = lnprior(theta, paramset=mcmc_paramset)
+def ln_prob(theta, args, fitter, asimov_paramset, llh_paramset):
+    lp = lnprior(theta, paramset=llh_paramset)
     if not np.isfinite(lp):
         return -np.inf
     return lp + triangle_llh(
         theta, args=args, asimov_paramset=asimov_paramset,
-        mcmc_paramset=mcmc_paramset, fitter=fitter
+        llh_paramset=llh_paramset, fitter=fitter
     )
diff --git a/utils/mcmc.py b/utils/mcmc.py
index c712c3a..97b77dd 100644
--- a/utils/mcmc.py
+++ b/utils/mcmc.py
@@ -9,7 +9,6 @@ Useful functions to use an MCMC for the BSM flavour ratio analysis
 
 from __future__ import absolute_import, division
 
-import argparse
 from functools import partial
 
 import emcee
@@ -72,6 +71,14 @@ def mcmc_argparse(parser):
         type=partial(enum_parse, c=MCMCSeedType), choices=MCMCSeedType,
         help='Type of distrbution to make the initial MCMC seed'
     )
+    parser.add_argument(
+        '--plot-angles', type=misc_utils.parse_bool, default='True',
+        help='Plot MCMC triangle in the angles space'
+    )
+    parser.add_argument(
+        '--plot-elements', type=misc_utils.parse_bool, default='False',
+        help='Plot MCMC triangle in the mixing elements space'
+    )
 
 
 def flat_seed(paramset, nwalkers):
diff --git a/utils/misc.py b/utils/misc.py
index f0c1ad4..59c5edb 100644
--- a/utils/misc.py
+++ b/utils/misc.py
@@ -9,7 +9,7 @@ Misc functions for the BSM flavour ratio analysis
 
 from __future__ import absolute_import, division
 
-import os, sys
+import os
 import errno
 import multiprocessing
 
@@ -19,174 +19,7 @@ from operator import attrgetter
 
 import numpy as np
 
-from utils.enums import Likelihood, ParamTag
-
-
-class Param(object):
-    """Parameter class to store parameters.
-    """
-    def __init__(self, name, value, ranges, seed=None, std=None, tex=None, tag=None):
-        self._seed = None
-        self._ranges = None
-        self._tex = None
-        self._tag = None
-
-        self.name = name
-        self.value = value
-        self.seed = seed
-        self.ranges = ranges
-        self.std = std
-        self.tex = tex
-        self.tag = tag
-
-    @property
-    def ranges(self):
-        return tuple(self._ranges)
-
-    @ranges.setter
-    def ranges(self, values):
-        self._ranges = [val for val in values]
-
-    @property
-    def seed(self):
-        if self._seed is None: return self.ranges
-        return tuple(self._seed)
-
-    @seed.setter
-    def seed(self, values):
-        if values is None: return
-        self._seed = [val for val in values]
-
-    @property
-    def tex(self):
-        return r'{0}'.format(self._tex)
-
-    @tex.setter
-    def tex(self, t):
-        self._tex = t if t is not None else r'{\rm %s}' % self.name
-
-    @property
-    def tag(self):
-        return self._tag
-
-    @tag.setter
-    def tag(self, t):
-        if t is None: self._tag = ParamTag.NONE
-        else:
-            assert t in ParamTag
-            self._tag = t
-
-
-class ParamSet(Sequence):
-    """Container class for a set of parameters.
-    """
-    def __init__(self, *args):
-        param_sequence = []
-        for arg in args:
-            try:
-                param_sequence.extend(arg)
-            except TypeError:
-                param_sequence.append(arg)
-
-        if len(param_sequence) != 0:
-            # Disallow duplicated params
-            all_names = [p.name for p in param_sequence]
-            unique_names = set(all_names)
-            if len(unique_names) != len(all_names):
-                duplicates = set([x for x in all_names if all_names.count(x) > 1])
-                raise ValueError('Duplicate definitions found for param(s): ' +
-                                 ', '.join(str(e) for e in duplicates))
-
-        # Elements of list must be Param type
-        assert all([isinstance(x, Param) for x in param_sequence]), \
-                'All params must be of type "Param"'
-
-        self._params = param_sequence
-
-    def __len__(self):
-        return len(self._params)
-
-    def __getitem__(self, i):
-        if isinstance(i, int):
-            return self._params[i]
-        elif isinstance(i, basestring):
-            return self._by_name[i]
-
-    def __getattr__(self, attr):
-        try:
-            return super(ParamSet, self).__getattribute__(attr)
-        except AttributeError:
-            t, v, tb = sys.exc_info()
-            try:
-                return self[attr]
-            except KeyError:
-                raise t, v, tb
-
-    def __iter__(self):
-        return iter(self._params)
-
-    def __str__(self):
-        o = '\n'
-        for obj in self._params:
-            o += '== {0:<15} = {1:<15}, tag={2:<15}\n'.format(
-                obj.name, obj.value, obj.tag
-            )
-        return o
-
-    @property
-    def _by_name(self):
-        return {obj.name: obj for obj in self._params}
-
-    @property
-    def names(self):
-        return tuple([obj.name for obj in self._params])
-
-    @property
-    def labels(self):
-        return tuple([obj.tex for obj in self._params])
-
-    @property
-    def values(self):
-        return tuple([obj.value for obj in self._params])
-
-    @property
-    def seeds(self):
-        return tuple([obj.seed for obj in self._params])
-
-    @property
-    def ranges(self):
-        return tuple([obj.ranges for obj in self._params])
-
-    @property
-    def stds(self):
-        return tuple([obj.std for obj in self._params])
-
-    @property
-    def tags(self):
-        return tuple([obj.tag for obj in self._params])
-
-    @property
-    def params(self):
-        return self._params
-
-    def to_dict(self):
-        return {obj.name: obj.value for obj in self._params}
-
-    def from_tag(self, tag, values=False, index=False, invert=False):
-        if values and index: assert 0
-        tag = np.atleast_1d(tag)
-        if not invert:
-            ps = [(idx, obj) for idx, obj in enumerate(self._params)
-                  if obj.tag in tag]
-        else:
-            ps = [(idx, obj) for idx, obj in enumerate(self._params)
-                  if obj.tag not in tag]
-        if values:
-            return tuple([io[1].value for io in ps])
-        elif index:
-            return tuple([io[0] for io in ps])
-        else:
-            return ParamSet([io[1] for io in ps])
+from utils.enums import Likelihood
 
 
 class SortingHelpFormatter(argparse.HelpFormatter):
@@ -197,54 +30,32 @@ class SortingHelpFormatter(argparse.HelpFormatter):
 
 
 def gen_identifier(args):
-    mr = args.measured_ratio
-    si = args.sigma_ratio
+    f = '_DIM{0}'.format(args.dimension)
+    mr1, mr2, mr3 = args.measured_ratio
     if args.fix_source_ratio:
-        sr = args.source_ratio
+        sr1, sr2, sr3 = args.source_ratio
+        f += '_sfr_{0:03d}_{1:03d}_{2:03d}_mfr_{3:03d}_{4:03d}_{5:03d}'.format(
+            int(sr1*100), int(sr2*100), int(sr3*100),
+            int(mr1*100), int(mr2*100), int(mr3*100)
+        )
         if args.fix_mixing:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_fix_mixing'.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), args.dimension
-            )
+            f += '_fix_mixing'
         elif args.fix_mixing_almost:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_fix_mixing_almost'.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), args.dimension
-            )
+            f += '_fix_mixing_almost'
         elif args.fix_scale:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_fixed_scale_{8}'.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), args.dimension,
-                args.scale
-            )
-        else:
-            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), args.dimension
-            )
+            f += '_fix_scale_{0}'.format(args.scale)
     else:
-        if args.fix_mixing:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fix_mixing'.format(
-                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100),
-                int(si*1000), args.dimension
-            )
-        elif args.fix_mixing_almost:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fix_mixing_almost'.format(
-                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100),
-                int(si*1000), args.dimension
-            )
-        elif args.fix_scale:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fixed_scale_{5}'.format(
-                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100),
-                int(si*1000), args.dimension, args.scale
-            )
-        else:
-            out = '_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}'.format(
-                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100),
-                int(si*1000), args.dimension
-            )
-    if args.likelihood is Likelihood.FLAT: out += '_flat'
-    return out
+        f += '_mfr_{3:03d}_{4:03d}_{5:03d}'.format(mr1, mr2, mr3)
+    if args.fix_mixing:
+        f += '_fix_mixing'
+    elif args.fix_mixing_almost:
+        f += '_fix_mixing_almost'
+    elif args.fix_scale:
+        f += '_fix_scale_{0}'.format(args.scale)
+    if args.likelihood is Likelihood.FLAT: f += '_flat'
+    elif args.likelihood is Likelihood.GAUSSIAN:
+        f += '_sigma_{0:03d}'.format(int(args.sigma_ratio*1000))
+    return f
 
 
 def gen_outfile_name(args):
@@ -314,6 +125,13 @@ def make_dir(outfile):
             raise
 
 
+def seed_parse(s):
+    if s.lower() == 'none':
+        return None
+    else:
+        return int(s)
+
+
 def thread_type(t):
     if t.lower() == 'max':
         return multiprocessing.cpu_count()
diff --git a/utils/multinest.py b/utils/multinest.py
new file mode 100644
index 0000000..3a7ab2d
--- /dev/null
+++ b/utils/multinest.py
@@ -0,0 +1,117 @@
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : April 19, 2018
+
+"""
+Useful functions to use MultiNest for the BSM flavour ratio analysis
+"""
+
+from __future__ import absolute_import, division
+
+from functools import partial
+
+import numpy as np
+
+from pymultinest import analyse, run
+
+from utils.likelihood import triangle_llh
+from utils.misc import make_dir, parse_bool
+
+
+def CubePrior(cube, ndim, n_params):
+    # default are uniform priors
+    return ;
+
+
+def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset,
+           args, fitter):
+    if ndim != len(mn_paramset):
+        raise AssertionError(
+            'Length of MultiNest scan paramset is not the same as the input '
+            'params\ncube={0}\nmn_paramset]{1}'.format(cube, mn_paramset)
+        )
+    pranges = mn_paramset.seeds
+    for i in xrange(ndim):
+        mn_paramset[i].value = (pranges[i][1]-pranges[i][0])*cube[i] + pranges[i][0]
+    for pm in mn_paramset.names:
+        llh_paramset[pm].value = mn_paramset[pm].value
+    theta = llh_paramset.values
+    # print 'llh_paramset', llh_paramset
+    return triangle_llh(
+        theta=theta,
+        args=args,
+        asimov_paramset=asimov_paramset,
+        llh_paramset=llh_paramset,
+        fitter=fitter
+    )
+
+
+def mn_argparse(parser):
+    parser.add_argument(
+        '--mn-run', type=parse_bool, default='True',
+        help='Run MultiNest'
+    )
+    parser.add_argument(
+        '--mn-bins', type=int, default=10,
+        help='Number of bins for the Bayes factor plot'
+    )
+    parser.add_argument(
+        '--mn-eval-bin', type=str, default='all',
+        help='Which bin to evalaute for Bayes factor plot'
+    )
+    parser.add_argument(
+        '--mn-live-points', type=int, default=3000,
+        help='Number of live points for MultiNest evaluations'
+    )
+    parser.add_argument(
+        '--mn-tolerance', type=float, default=0.01,
+        help='Tolerance for MultiNest'
+    )
+    parser.add_argument(
+        '--mn-output', type=str, default='./mnrun/',
+        help='Folder to store MultiNest evaluations'
+    )
+    parser.add_argument(
+        '--plot-multinest', type=parse_bool, default='False',
+        help='Plot MultiNest evidence'
+    )
+
+
+def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter):
+    """Run the MultiNest algorithm to calculate the evidence."""
+    n_params = len(mn_paramset)
+
+    for n in mn_paramset.names:
+        llh_paramset[n].value = mn_paramset[n].value
+
+    lnProbEval = partial(
+        lnProb,
+        mn_paramset     = mn_paramset,
+        llh_paramset    = llh_paramset,
+        asimov_paramset = asimov_paramset,
+        args            = args,
+        fitter          = fitter
+    )
+
+    prefix = './mnrun/DIM{0}/{1:>019}_'.format(
+        args.dimension, np.random.randint(0, 2**63)
+    )
+    make_dir(prefix)
+    print 'Running evidence calculation for {0}'.format(prefix)
+    result = 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
+    )
+
+    analyser = analyse.Analyzer(
+        outputfiles_basename=prefix, n_params=n_params
+    )
+    return analyser.get_stats()['global evidence']
diff --git a/utils/param.py b/utils/param.py
new file mode 100644
index 0000000..bf230df
--- /dev/null
+++ b/utils/param.py
@@ -0,0 +1,235 @@
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : April 19, 2018
+
+"""
+Param class and functions for the BSM flavour ratio analysis
+"""
+
+from __future__ import absolute_import, division
+
+import sys
+
+from collections import Sequence
+
+import numpy as np
+
+from utils.plot import get_units
+from utils.fr import fr_to_angles
+from utils.enums import Likelihood, ParamTag
+
+
+class Param(object):
+    """Parameter class to store parameters."""
+    def __init__(self, name, value, ranges, seed=None, std=None, tex=None, tag=None):
+        self._seed = None
+        self._ranges = None
+        self._tex = None
+        self._tag = None
+
+        self.name = name
+        self.value = value
+        self.seed = seed
+        self.ranges = ranges
+        self.std = std
+        self.tex = tex
+        self.tag = tag
+
+    @property
+    def ranges(self):
+        return tuple(self._ranges)
+
+    @ranges.setter
+    def ranges(self, values):
+        self._ranges = [val for val in values]
+
+    @property
+    def seed(self):
+        if self._seed is None: return self.ranges
+        return tuple(self._seed)
+
+    @seed.setter
+    def seed(self, values):
+        if values is None: return
+        self._seed = [val for val in values]
+
+    @property
+    def tex(self):
+        return r'{0}'.format(self._tex)
+
+    @tex.setter
+    def tex(self, t):
+        self._tex = t if t is not None else r'{\rm %s}' % self.name
+
+    @property
+    def tag(self):
+        return self._tag
+
+    @tag.setter
+    def tag(self, t):
+        if t is None: self._tag = ParamTag.NONE
+        else:
+            assert t in ParamTag
+            self._tag = t
+
+
+class ParamSet(Sequence):
+    """Container class for a set of parameters."""
+    def __init__(self, *args):
+        param_sequence = []
+        for arg in args:
+            try:
+                param_sequence.extend(arg)
+            except TypeError:
+                param_sequence.append(arg)
+
+        if len(param_sequence) != 0:
+            # Disallow duplicated params
+            all_names = [p.name for p in param_sequence]
+            unique_names = set(all_names)
+            if len(unique_names) != len(all_names):
+                duplicates = set([x for x in all_names if all_names.count(x) > 1])
+                raise ValueError('Duplicate definitions found for param(s): ' +
+                                 ', '.join(str(e) for e in duplicates))
+
+        # Elements of list must be Param type
+        assert all([isinstance(x, Param) for x in param_sequence]), \
+                'All params must be of type "Param"'
+
+        self._params = param_sequence
+
+    def __len__(self):
+        return len(self._params)
+
+    def __getitem__(self, i):
+        if isinstance(i, int):
+            return self._params[i]
+        elif isinstance(i, basestring):
+            return self._by_name[i]
+
+    def __getattr__(self, attr):
+        try:
+            return super(ParamSet, self).__getattribute__(attr)
+        except AttributeError:
+            t, v, tb = sys.exc_info()
+            try:
+                return self[attr]
+            except KeyError:
+                raise t, v, tb
+
+    def __iter__(self):
+        return iter(self._params)
+
+    def __str__(self):
+        o = '\n'
+        for obj in self._params:
+            o += '== {0:<15} = {1:<15}, tag={2:<15}\n'.format(
+                obj.name, obj.value, obj.tag
+            )
+        return o
+
+    @property
+    def _by_name(self):
+        return {obj.name: obj for obj in self._params}
+
+    @property
+    def names(self):
+        return tuple([obj.name for obj in self._params])
+
+    @property
+    def labels(self):
+        return tuple([obj.tex for obj in self._params])
+
+    @property
+    def values(self):
+        return tuple([obj.value for obj in self._params])
+
+    @property
+    def seeds(self):
+        return tuple([obj.seed for obj in self._params])
+
+    @property
+    def ranges(self):
+        return tuple([obj.ranges for obj in self._params])
+
+    @property
+    def stds(self):
+        return tuple([obj.std for obj in self._params])
+
+    @property
+    def tags(self):
+        return tuple([obj.tag for obj in self._params])
+
+    @property
+    def params(self):
+        return self._params
+
+    def to_dict(self):
+        return {obj.name: obj.value for obj in self._params}
+
+    def from_tag(self, tag, values=False, index=False, invert=False):
+        if values and index: assert 0
+        tag = np.atleast_1d(tag)
+        if not invert:
+            ps = [(idx, obj) for idx, obj in enumerate(self._params)
+                  if obj.tag in tag]
+        else:
+            ps = [(idx, obj) for idx, obj in enumerate(self._params)
+                  if obj.tag not in tag]
+        if values:
+            return tuple([io[1].value for io in ps])
+        elif index:
+            return tuple([io[0] for io in ps])
+        else:
+            return ParamSet([io[1] for io in ps])
+
+
+def get_paramsets(args):
+    """Make the paramsets for generating the Asmimov MC sample and also running
+    the MCMC.
+    """
+    asimov_paramset = []
+    llh_paramset = []
+    if args.likelihood == Likelihood.GOLEMFIT:
+        nuisance_paramset = define_nuisance()
+        asimov_paramset.extend(nuisance_paramset.params)
+        llh_paramset.extend(nuisance_paramset.params)
+        for parm in nuisance_paramset:
+            parm.value = args.__getattribute__(parm.name)
+    tag = ParamTag.BESTFIT
+    flavour_angles = fr_to_angles(args.measured_ratio)
+    asimov_paramset.extend([
+        Param(name='astroFlavorAngle1', value=flavour_angles[0], ranges=[0., 1.], std=0.2, tag=tag),
+        Param(name='astroFlavorAngle2', value=flavour_angles[1], ranges=[-1., 1.], std=0.2, tag=tag),
+    ])
+    asimov_paramset = ParamSet(asimov_paramset)
+
+    if not args.fix_mixing and not args.fix_mixing_almost:
+        tag = ParamTag.MMANGLES
+        e = 1e-10
+        llh_paramset.extend([
+            Param(name='s_12^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag),
+            Param(name='c_13^4', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag),
+            Param(name='s_23^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag),
+            Param(name='dcp', value=np.pi, ranges=[0.+e, 2*np.pi-e], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag)
+        ])
+    if args.fix_mixing_almost:
+        tag = ParamTag.MMANGLES
+        llh_paramset.extend([
+            Param(name='s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag)
+        ])
+    if not args.fix_scale:
+        tag = ParamTag.SCALE
+        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)
+        )
+    if not args.fix_source_ratio:
+        tag = ParamTag.SRCANGLES
+        llh_paramset.extend([
+            Param(name='s_phi4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'sin^4(\phi)', tag=tag),
+            Param(name='c_2psi', value=0.5, ranges=[0., 1.], std=0.2, tex=r'cos(2\psi)', tag=tag)
+        ])
+    llh_paramset = ParamSet(llh_paramset)
+    return asimov_paramset, llh_paramset
diff --git a/utils/plot.py b/utils/plot.py
index 0b82675..66c888c 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -86,26 +86,6 @@ def interval(arr, percentile=68.):
     return sarr[curr_low], center, sarr[curr_up]
 
 
-def plot_argparse(parser):
-    """Arguments for plotting."""
-    parser.add_argument(
-        '--plot-angles', type=misc_utils.parse_bool, default='True',
-        help='Plot MCMC triangle in the angles space'
-    )
-    parser.add_argument(
-        '--plot-elements', type=misc_utils.parse_bool, default='False',
-        help='Plot MCMC triangle in the mixing elements space'
-    )
-    parser.add_argument(
-        '--plot-bayes', type=misc_utils.parse_bool, default='False',
-        help='Plot Bayes factor'
-    )
-    parser.add_argument(
-        '--plot-angles-limit', type=misc_utils.parse_bool, default='False',
-        help='Plot limit vs BSM angles'
-    )
-
-
 def flat_angles_to_u(x):
     """Convert from angles to mixing elements."""
     return abs(angles_to_u(x)).astype(np.float32).flatten().tolist()
@@ -138,32 +118,19 @@ def gen_figtext(args):
     mr1, mr2, mr3 = args.measured_ratio
     if args.fix_source_ratio:
         sr1, sr2, sr3 = args.source_ratio
-        if args.fix_scale:
-            t += 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \
-                    'observed flavour ratio = [{3:.2f}, {4:.2f}, ' \
-                    '{5:.2f}]\nDimension = {6}\nScale = {7}'.format(
-                        sr1, sr2, sr3, mr1, mr2, mr3, args.dimension,
-                        int(args.energy), args.scale
-                    )
-        else:
-            t += 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \
-                    'observed flavour ratio = [{3:.2f}, {4:.2f}, ' \
-                    '{5:.2f}]\nDimension = {6}'.format(
-                        sr1, sr2, sr3, mr1, mr2, mr3, args.dimension,
-                        int(args.energy)
-                    )
+        t += 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \
+                'observed flavour ratio = [{3:.2f}, {4:.2f}, ' \
+                '{5:.2f}]\nDimension = {6}'.format(
+                    sr1, sr2, sr3, mr1, mr2, mr3, args.dimension,
+                    int(args.energy)
+                )
     else:
-        if args.fix_scale:
-            t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \
-                    '{2:.2f}]\nDimension = {3}\nScale = {4}'.format(
-                        mr1, mr2, mr3, args.dimension, int(args.energy),
-                        args.scale
-                    )
-        else:
-            t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \
-                    '{2:.2f}]\nDimension = {3}'.format(
-                        mr1, mr2, mr3, args.dimension, int(args.energy)
-                    )
+        t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \
+                '{2:.2f}]\nDimension = {3}'.format(
+                    mr1, mr2, mr3, args.dimension, int(args.energy)
+                )
+    if args.fix_scale:
+        t += 'Scale = {0}'.format(args.scale)
     if args.likelihood is Likelihood.GAUSSIAN:
         t += '\nSigma = {0:.3f}'.format(args.sigma_ratio)
     if args.energy_dependance is EnergyDependance.SPECTRAL:
@@ -176,7 +143,7 @@ def gen_figtext(args):
     return t
 
 
-def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
+def chainer_plot(infile, outfile, outformat, args, llh_paramset):
     """Make the triangle plot."""
     if not args.plot_angles and not args.plot_elements:
         return
@@ -187,8 +154,8 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
     misc_utils.make_dir(outfile)
     fig_text = gen_figtext(args)
 
-    axes_labels = mcmc_paramset.labels
-    ranges = mcmc_paramset.ranges
+    axes_labels = llh_paramset.labels
+    ranges = llh_paramset.ranges
 
     if args.plot_angles:
         print "Making triangle plots"
@@ -208,7 +175,7 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
                         ), fontsize=10)
 
         # if not args.fix_mixing:
-        #     sc_index = mcmc_paramset.from_tag(ParamTag.SCALE, index=True)
+        #     sc_index = llh_paramset.from_tag(ParamTag.SCALE, index=True)
         #     itv = interval(Tchain[:,sc_index], percentile=90.)
         #     mpl.pyplot.figtext(
         #         0.5, 0.3, 'Scale 90% Interval = [1E{0}, 1E{1}], Center = '
@@ -222,11 +189,11 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
         print "Making triangle plots"
         if args.fix_mixing_almost:
             raise NotImplementedError
-        nu_index = mcmc_paramset.from_tag(ParamTag.NUISANCE, index=True)
-        fr_index = mcmc_paramset.from_tag(ParamTag.MMANGLES, index=True)
-        sc_index = mcmc_paramset.from_tag(ParamTag.SCALE, index=True)
+        nu_index = llh_paramset.from_tag(ParamTag.NUISANCE, index=True)
+        fr_index = llh_paramset.from_tag(ParamTag.MMANGLES, index=True)
+        sc_index = llh_paramset.from_tag(ParamTag.SCALE, index=True)
         if not args.fix_source_ratio:
-            sr_index = mcmc_paramset.from_tag(ParamTag.SRCANGLES, index=True)
+            sr_index = llh_paramset.from_tag(ParamTag.SRCANGLES, index=True)
 
         nu_elements = raw[:,nu_index]
         fr_elements = np.array(map(flat_angles_to_u, raw[:,fr_index]))
@@ -264,45 +231,39 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
             g.export(outfile+'_elements.'+of)
 
 
-def bayes_factor_plot(dirname, outfile, outformat, args):
-    """Make Bayes factor plot."""
-    if not args.plot_bayes: return
-    print "Making Bayes Factor plot"
-    print 'dirname', dirname
+def plot_multinest(data, outfile, outformat, args, scale_param, label=None):
+    """Make MultiNest factor plot."""
+    print "Making MultiNest Factor plot"
     fig_text = gen_figtext(args)
+    if label is not None: fig_text += '\n' + label
 
-    raw = []
-    for root, dirs, filenames in os.walk(dirname):
-        for fn in filenames:
-            if fn[-4:] == '.npy':
-                raw.append(np.load(os.path.join(root, fn)))
-    raw = np.sort(np.vstack(raw), axis=0)
-    print 'raw', raw
-    print 'raw.shape', raw.shape
-    scales, evidences = raw.T
-    null = evidences[0]
+    print 'data.shape', data.shape
+    scales, evidences = data.T
+    min_idx = np.argmin(scales)
+    null = evidences[min_idx]
 
     reduced_ev = -(evidences - null)
 
     fig = plt.figure(figsize=(7, 5))
     ax = fig.add_subplot(111)
 
-    ax.set_xlim(np.log10(args.scale_region))
-    ax.set_xlabel(r'${\rm log}_{10} \Lambda ' + get_units(args.dimension) +r'$')
+    ax.set_xlim(scale_param.ranges)
+    ax.set_xlabel('$'+scale_param.tex+'$')
     ax.set_ylabel(r'Bayes Factor')
 
     ax.plot(scales, reduced_ev)
 
     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.3, 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.3, linewidth=1)
 
     at = AnchoredText(
-        fig_text, prop=dict(size=7), frameon=True, loc=2
+        '\n'+fig_text, prop=dict(size=7), frameon=True, loc=2
     )
-    at.patch.set_boxstyle("round,pad=0.,rounding_size=0.5")
+    at.patch.set_boxstyle("round,pad=0.1,rounding_size=0.5")
     ax.add_artist(at)
+    misc_utils.make_dir(outfile)
     for of in outformat:
         fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150)
 
@@ -326,7 +287,8 @@ def plot_BSM_angles_limit(dirname, outfile, outformat, args, bayesian):
         for fn in filenames:
             if fn[-4:] == '.npy':
                 raw.append(np.load(os.path.join(root, fn)))
-    raw = np.sort(np.vstack(raw), axis=0)
+    raw = np.vstack(raw)
+    raw = raw[np.argsort(raw[:,0])]
     print 'raw', raw
     print 'raw.shape', raw.shape
     sc_ranges = (