Sun May 6 21:22:46 CDT 2018

19 files changed, 1111 insertions, 171 deletions

diff --git a/fr.py b/fr.py
index 0e14098..30f9793 100755
--- a/fr.py
+++ b/fr.py
@@ -33,28 +33,28 @@ def define_nuisance():
     g_prior = PriorsCateg.GAUSSIAN
     hg_prior = PriorsCateg.HALFGAUSS
     e = 1e-9
-    # nuisance.extend([
-    #     Param(name='s_12_2', value=0.307,          seed=[0.26, 0.35],     ranges=[0., 1.],      std=0.013,   tex=r's_{12}^2', prior=g_prior,  tag=tag),
-    #     Param(name='c_13_4', value=1-(0.02206)**2, seed=[0.995, 1-e],     ranges=[0., 1.],      std=0.00147, tex=r'c_{13}^4', prior=hg_prior, tag=tag),
-    #     Param(name='s_23_2', value=0.538,          seed=[0.31, 0.75],     ranges=[0., 1.],      std=0.069,   tex=r's_{23}^2', prior=g_prior,  tag=tag),
-    #     Param(name='dcp',    value=4.08404,        seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0,     tex=r'\delta_{CP}', tag=tag),
-    #     Param(
-    #         name='m21_2', value=7.40E-23, seed=[7.2E-23, 7.6E-23], ranges=[6.80E-23, 8.02E-23],
-    #         std=2.1E-24, tex=r'\Delta m_{21}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag
-    #     ),
-    #     Param(
-    #         name='m3x_2', value=2.494E-21, seed=[2.46E-21, 2.53E-21], ranges=[2.399E-21, 2.593E-21],
-    #         std=3.3E-23, tex=r'\Delta m_{3x}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag
-    #     )
-    # ])
+    nuisance.extend([
+        Param(name='s_12_2', value=0.307,            seed=[0.26, 0.35],     ranges=[0., 1.],      std=0.013,   tex=r's_{12}^2', prior=g_prior,  tag=tag),
+        Param(name='c_13_4', value=(1-(0.02206))**2, seed=[0.950, 0.961],   ranges=[0., 1.],      std=0.00147, tex=r'c_{13}^4', prior=g_prior, tag=tag),
+        Param(name='s_23_2', value=0.538,            seed=[0.31, 0.75],     ranges=[0., 1.],      std=0.069,   tex=r's_{23}^2', prior=g_prior,  tag=tag),
+        Param(name='dcp',    value=4.08404,          seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0,     tex=r'\delta_{CP}', tag=tag),
+        Param(
+            name='m21_2', value=7.40E-23, seed=[7.2E-23, 7.6E-23], ranges=[6.80E-23, 8.02E-23],
+            std=2.1E-24, tex=r'\Delta m_{21}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag
+        ),
+        Param(
+            name='m3x_2', value=2.494E-21, seed=[2.46E-21, 2.53E-21], ranges=[2.399E-21, 2.593E-21],
+            std=3.3E-23, tex=r'\Delta m_{3x}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag
+        )
+    ])
     tag = ParamTag.NUISANCE
-    # nuisance.extend([
-    #     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)
-    # ])
+    nuisance.extend([
+        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 ParamSet(nuisance)
 
 
diff --git a/plot_sens.py b/plot_sens.py
index e51d55d..0c862eb 100755
--- a/plot_sens.py
+++ b/plot_sens.py
@@ -10,21 +10,355 @@ HESE BSM flavour ratio analysis plotting script
 
 from __future__ import absolute_import, division
 
+import os
+import argparse
+from functools import partial
+from copy import deepcopy
+
 import numpy as np
-import matplotlib as mpl
-mpl.use('Agg')
-from matplotlib import rc
-from matplotlib import pyplot as plt
-from matplotlib.offsetbox import AnchoredText
-
-rc('text', usetex=False)
-rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
-
-FRS = [
-    (1, 1, 1, 1, 2, 0),
-    (1, 1, 1, 0, 1, 0),
-    # (1, 1, 1, 1, 0, 0),
-    # (1, 1, 1, 0, 0, 1),
-]
-
-DIMENSION = [3, 6]
+import numpy.ma as ma
+
+from utils import fr as fr_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 PriorsCateg, SensitivityCateg, StatCateg
+from utils.param import Param, ParamSet, get_paramsets
+
+from utils import multinest as mn_utils
+
+
+def define_nuisance():
+    """Define the nuisance parameters."""
+    tag = ParamTag.SM_ANGLES
+    g_prior = PriorsCateg.GAUSSIAN
+    hg_prior = PriorsCateg.HALFGAUSS
+    e = 1e-9
+    nuisance = [
+        Param(name='s_12_2', value=0.307,          seed=[0.26, 0.35],     ranges=[0., 1.],      std=0.013,   tex=r's_{12}^2', prior=g_prior,  tag=tag),
+        Param(name='c_13_4', value=1-(0.02206)**2, seed=[0.995, 1-e],     ranges=[0., 1.],      std=0.00147, tex=r'c_{13}^4', prior=hg_prior, tag=tag),
+        Param(name='s_23_2', value=0.538,          seed=[0.31, 0.75],     ranges=[0., 1.],      std=0.069,   tex=r's_{23}^2', prior=g_prior,  tag=tag),
+        Param(name='dcp',    value=4.08404,        seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0,     tex=r'\delta_{CP}', tag=tag),
+        Param(
+            name='m21_2', value=7.40E-23, seed=[7.2E-23, 7.6E-23], ranges=[6.80E-23, 8.02E-23],
+            std=2.1E-24, tex=r'\Delta m_{21}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag
+        ),
+        Param(
+            name='m3x_2', value=2.494E-21, seed=[2.46E-21, 2.53E-21], ranges=[2.399E-21, 2.593E-21],
+            std=3.3E-23, tex=r'\Delta m_{3x}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag
+        )
+    ]
+    tag = ParamTag.NUISANCE
+    nuisance.extend([
+        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 ParamSet(nuisance)
+
+
+def nuisance_argparse(parser):
+    nuisance = define_nuisance()
+    for parm in nuisance:
+        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.fix_scale:
+        raise NotImplementedError(
+            '--fix-scale not implemented'
+        )
+
+    args.measured_ratio = fr_utils.normalise_fr(args.measured_ratio)
+    if args.fix_source_ratio:
+        assert len(args.source_ratios) % 3 == 0
+        srs = [args.source_ratios[3*x:3*x+3]
+               for x in range(int(len(args.source_ratios)/3))]
+        args.source_ratios = map(fr_utils.normalise_fr, srs)
+
+    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 args.split_jobs and args.run_method is SensitivityCateg.FULL:
+        raise NotImplementedError(
+            'split_jobs and run_method not implemented'
+        )
+
+    args.dimensions = np.sort(args.dimensions)
+
+    args_copy = deepcopy(args)
+    scale_regions = []
+    for dim in args.dimensions:
+        args_copy.dimension = dim
+        _, scale_region = fr_utils.estimate_scale(args_copy)
+        scale_regions.append(scale_region)
+    args.scale_region = [np.min(scale_regions), np.max(scale_regions)]
+    args.scale = np.power(10., np.average(np.log10(args.scale_region)))
+
+
+def parse_args(args=None):
+    """Parse command line arguments"""
+    parser = argparse.ArgumentParser(
+        description="HESE BSM flavour ratio analysis plotting script",
+        formatter_class=misc_utils.SortingHelpFormatter,
+    )
+    parser.add_argument(
+        '--infile', type=str, default='./untitled',
+        help='Path to input dir'
+    )
+    parser.add_argument(
+        '--run-method', default='full',
+        type=partial(misc_utils.enum_parse, c=SensitivityCateg),
+        choices=SensitivityCateg,
+        help='Choose which type of sensivity plot to make'
+    )
+    parser.add_argument(
+        '--stat-method', default='bayesian',
+        type=partial(misc_utils.enum_parse, c=StatCateg), choices=StatCateg,
+        help='Statistical method to employ'
+    )
+    parser.add_argument(
+        '--sens-bins', type=int, default=10,
+        help='Number of bins for the Bayes factor plot'
+    )
+    parser.add_argument(
+        '--split-jobs', type=misc_utils.parse_bool, default='False',
+        help='Did the jobs get split'
+    )
+    parser.add_argument(
+        '--plot', type=misc_utils.parse_bool, default='True',
+        help='Make sensitivity plots'
+    )
+    parser.add_argument(
+        '--plot-statistic', type=misc_utils.parse_bool, default='False',
+        help='Plot MultiNest evidence or LLH value'
+    )
+    fr_utils.fr_argparse(parser)
+    llh_utils.likelihood_argparse(parser)
+    mn_utils.mn_argparse(parser)
+    nuisance_argparse(parser)
+    misc_utils.remove_option(parser, 'dimension')
+    misc_utils.remove_option(parser, 'source_ratio')
+    misc_utils.remove_option(parser, 'scale')
+    misc_utils.remove_option(parser, 'scale_region')
+    parser.add_argument(
+        '--dimensions', type=int, nargs='*', default=[3, 6],
+        help='Set the new physics dimensions to consider'
+    )
+    parser.add_argument(
+        '--source-ratios', type=int, nargs='*', default=[2, 1, 0],
+        help='Set the source flavour ratios for the case when you want to fix it'
+    )
+    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)
+
+    asimov_paramset, llh_paramset = get_paramsets(args, define_nuisance())
+
+    scale = llh_paramset.from_tag(ParamTag.SCALE)[0]
+    mmangles = llh_paramset.from_tag(ParamTag.MMANGLES)
+    if args.run_method is SensitivityCateg.FULL:
+        st_paramset_arr = [llh_paramset.from_tag(ParamTag.SCALE, invert=True)]
+    elif args.run_method in [SensitivityCateg.FIXED_ANGLE, SensitivityCateg.CORR_ANGLE]:
+        nscale_pmset = llh_paramset.from_tag([ParamTag.SCALE, ParamTag.MMANGLES], invert=True)
+        st_paramset_arr = [nscale_pmset] * 3
+    elif args.run_method in [SensitivityCateg.FIXED_ONE_ANGLE, SensitivityCateg.CORR_ONE_ANGLE]:
+        nscale_pmset = llh_paramset.from_tag(ParamTag.SCALE, invert=True)
+        st_paramset_arr = []
+        for x in xrange(3):
+            st_paramset_arr.append(
+                ParamSet([prms for prms in nscale_pmset
+                          if mmangles[x].name != prms.name])
+            )
+
+    corr_angles_categ = [SensitivityCateg.CORR_ANGLE, SensitivityCateg.CORR_ONE_ANGLE]
+    fixed_angle_categ = [SensitivityCateg.FIXED_ANGLE, SensitivityCateg.FIXED_ONE_ANGLE]
+
+    if args.run_method in corr_angles_categ:
+        scan_angles = np.linspace(0+1e-9, 1-1e-9, args.sens_bins)
+    else: scan_angles = np.array([0])
+    print 'scan_angles', scan_angles
+
+    dims = len(args.dimensions)
+    srcs = len(args.source_ratios)
+    if args.run_method is SensitivityCateg.FULL:
+        statistic_arr = np.full((dims, srcs, args.sens_bins, 2), np.nan)
+    elif args.run_method in fixed_angle_categ:
+        statistic_arr = np.full((dims, srcs, len(st_paramset_arr), args.sens_bins, 2), np.nan)
+    elif args.run_method in corr_angles_categ:
+        statistic_arr = np.full(
+            (dims, srcs, len(st_paramset_arr), args.sens_bins, args.sens_bins, 3), np.nan
+        )
+
+    print 'Loading data'
+    for idim, dim in enumerate(args.dimensions):
+        argsc = deepcopy(args)
+        argsc.dimension = dim
+        _, scale_region = fr_utils.estimate_scale(argsc)
+        argsc.scale_region = scale_region
+        scan_scales = np.linspace(
+            np.log10(scale_region[0]), np.log10(scale_region[1]), args.sens_bins
+        )
+
+        for isrc, src in enumerate(args.source_ratios):
+            argsc.source_ratio = src
+            infile = args.infile
+            if args.likelihood is Likelihood.GOLEMFIT:
+                infile += '/golemfit/'
+            elif args.likelihood is Likelihood.GAUSSIAN:
+                infile += '/gaussian/'
+            infile += '/DIM{0}/fix_ifr/'.format(dim)
+            if args.likelihood is Likelihood.GAUSSIAN:
+                infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_'))
+            infile += 'fr_stat/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \
+                + misc_utils.gen_identifier(argsc)
+            print '== {0:<25} = {1}'.format('infile', infile)
+
+            if args.split_jobs:
+                for idx_an, an in enumerate(scan_angles):
+                    for idx_sc, sc in enumerate(scan_scales):
+                        filename = infile + '_scale_{0:.0E}'.format(np.power(10, sc))
+                        try:
+                            if args.run_method in fixed_angle_categ:
+                                print 'Loading from {0}'.format(filename+'.npy')
+                                statistic_arr[idim][isrc][:,idx_sc] = np.load(filename+'.npy')[:,0]
+                            if args.run_method in corr_angles_categ:
+                                filename += '_angle_{0:<04.2}'.format(an)
+                                print 'Loading from {0}'.format(filename+'.npy')
+                                statistic_arr[idim][isrc][:,idx_an,idx_sc] = np.load(filename+'.npy')[:,0,0]
+                        except:
+                            print 'Unable to load file {0}'.format(filename+'.npy')
+                            continue
+            else:
+                print 'Loading from {0}'.format(infile+'.npy')
+                try:
+                    statistic_arr[idim][isrc] = np.load(infile+'.npy')
+                except:
+                    print 'Unable to load file {0}'.format(infile+'.npy')
+                    continue
+
+    print 'statistic_arr', statistic_arr
+
+    data = ma.masked_invalid(statistic_arr)
+    if args.plot_statistic:
+        print 'Plotting statistic'
+
+        argsc = deepcopy(args)
+        base_infile = args.infile
+        if args.likelihood is Likelihood.GOLEMFIT:
+            base_infile += '/golemfit/'
+        elif args.likelihood is Likelihood.GAUSSIAN:
+            base_infile += '/gaussian/'
+        for idim, dim in enumerate(args.dimensions):
+            argsc.dimension = dim
+            _, scale_region = fr_utils.estimate_scale(argsc)
+            argsc.scale_region = scale_region
+            infile = base_infile + '/DIM{0}/fix_ifr/'.format(dim)
+            if args.likelihood is Likelihood.GAUSSIAN:
+                infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_'))
+            infile += 'fr_stat/'
+
+            for isrc, src in enumerate(args.source_ratios):
+                argsc.source_ratio = src
+                finfile = infile +'/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \
+                    + misc_utils.gen_identifier(argsc)
+                basename = os.path.dirname(finfile) + '/statrun/' + os.path.basename(finfile)
+                baseoutfile = basename[:5]+basename[5:].replace('data', 'plots')
+                if args.run_method is SensitivityCateg.FULL:
+                    outfile = baseoutfile
+                    plot_utils.plot_statistic(
+                        data        = data[idim][isrc],
+                        outfile     = outfile,
+                        outformat   = ['png'],
+                        args        = argsc,
+                        scale_param = scale,
+                    )
+                if args.run_method in fixed_angle_categ:
+                    for idx_scen in xrange(len(st_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{\pi}{4}$'
+                        elif idx_scen == 1: label = r'$\mathcal{O}_{13}=\frac{\pi}{4}$'
+                        elif idx_scen == 2: label = r'$\mathcal{O}_{23}=\frac{\pi}{4}$'
+                        plot_utils.plot_statistic(
+                            data        = data[idim][isrc][idx_scen],
+                            outfile     = outfile,
+                            outformat   = ['png'],
+                            args        = argsc,
+                            scale_param = scale,
+                            label       = label
+                        )
+                elif args.run_method in corr_angles_categ:
+                    for idx_scen in xrange(len(st_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(float(an))
+                            outfile = basescenoutfile + '_ANGLE{0}'.format(idx_an)
+                            _label = label + r'{0:<04.2}$'.format(an)
+                            plot_utils.plot_statistic(
+                                data        = data[idim][isrc][idx_scen][idx_an][:,1:],
+                                outfile     = outfile,
+                                outformat   = ['png'],
+                                args        = argsc,
+                                scale_param = scale,
+                                label       = _label
+                            )
+
+    if args.plot:
+        print 'Plotting sensitivities'
+
+        basename = args.infile + '{0}/{1}'.format(args.likelihood, args.stat_method)
+        baseoutfile = basename[:5]+basename[5:].replace('data', 'plots')
+
+        if args.run_method is SensitivityCateg.FULL:
+            plot_utils.plot_sens_full(
+                data      = data,
+                outfile   = basename + '/FULL',
+                outformat = ['png'],
+                args      = args,
+            )
+        elif args.run_method in fixed_angle_categ:
+            plot_utils.plot_sens_fixed_angle(
+                data      = data,
+                outfile   = basename + '/FIXED_ANGLE',
+                outformat = ['png'],
+                args      = args,
+            )
+        elif args.run_method in corr_angles_categ:
+            plot_utils.plot_sens_corr_angle(
+                data      = data,
+                outfile   = basename + '/CORR_ANGLE',
+                outformat = ['png'],
+                args      = args,
+            )
+
+
+main.__doc__ = __doc__
+
+
+if __name__ == '__main__':
+    main()
diff --git a/sens.py b/sens.py
index 489f126..1b04275 100755
--- a/sens.py
+++ b/sens.py
@@ -34,14 +34,14 @@ from utils import multinest as mn_utils
 def define_nuisance():
     """Define the nuisance parameters."""
     tag = ParamTag.SM_ANGLES
+    nuisance = []
     g_prior = PriorsCateg.GAUSSIAN
-    hg_prior = PriorsCateg.HALFGAUSS
     e = 1e-9
-    nuisance = [
-        Param(name='s_12_2', value=0.307,          seed=[0.26, 0.35],     ranges=[0., 1.],      std=0.013,   tex=r's_{12}^2', prior=g_prior,  tag=tag),
-        Param(name='c_13_4', value=1-(0.02206)**2, seed=[0.995, 1-e],     ranges=[0., 1.],      std=0.00147, tex=r'c_{13}^4', prior=hg_prior, tag=tag),
-        Param(name='s_23_2', value=0.538,          seed=[0.31, 0.75],     ranges=[0., 1.],      std=0.069,   tex=r's_{23}^2', prior=g_prior,  tag=tag),
-        Param(name='dcp',    value=4.08404,        seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0,     tex=r'\delta_{CP}', tag=tag),
+    nuisance.extend([
+        Param(name='s_12_2', value=0.307,            seed=[0.26, 0.35],     ranges=[0., 1.],      std=0.013,   tex=r's_{12}^2', prior=g_prior,  tag=tag),
+        Param(name='c_13_4', value=(1-(0.02206))**2, seed=[0.950, 0.961],   ranges=[0., 1.],      std=0.00147, tex=r'c_{13}^4', prior=g_prior, tag=tag),
+        Param(name='s_23_2', value=0.538,            seed=[0.31, 0.75],     ranges=[0., 1.],      std=0.069,   tex=r's_{23}^2', prior=g_prior,  tag=tag),
+        Param(name='dcp',    value=4.08404,          seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0,     tex=r'\delta_{CP}', tag=tag),
         Param(
             name='m21_2', value=7.40E-23, seed=[7.2E-23, 7.6E-23], ranges=[6.80E-23, 8.02E-23],
             std=2.1E-24, tex=r'\Delta m_{21}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag
@@ -50,7 +50,7 @@ def define_nuisance():
             name='m3x_2', value=2.494E-21, seed=[2.46E-21, 2.53E-21], ranges=[2.399E-21, 2.593E-21],
             std=3.3E-23, tex=r'\Delta m_{3x}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag
         )
-    ]
+    ])
     tag = ParamTag.NUISANCE
     nuisance.extend([
         Param(name='convNorm',        value=1.,  seed=[0.5, 2. ], ranges=[0. , 50.], std=0.3,  tag=tag),
@@ -173,10 +173,8 @@ def main():
         np.random.seed(args.seed)
 
     asimov_paramset, llh_paramset = get_paramsets(args, define_nuisance())
-    scale = llh_paramset.from_tag(ParamTag.SCALE)[0]
-    outfile = misc_utils.gen_outfile_name(args)
-    print '== {0:<25} = {1}'.format('outfile', outfile)
 
+    scale = llh_paramset.from_tag(ParamTag.SCALE)[0]
     mmangles = llh_paramset.from_tag(ParamTag.MMANGLES)
     if args.run_method is SensitivityCateg.FULL:
         st_paramset_arr = [llh_paramset.from_tag(ParamTag.SCALE, invert=True)]
@@ -246,7 +244,13 @@ def main():
             for idx_an, an in enumerate(scan_angles):
                 if args.run_method in corr_angles_categ:
                     for x in mmangles: x.value = 0.+1e-9
-                    mmangles[idx_scen].value = an
+                    angle = np.arcsin(np.sqrt(an))/2.
+                    if idx_scen == 0 or idx_scen == 2:
+                        mmangles[idx_scen].value = np.sin(angle)**2
+                        """s_12^2 or s_23^2"""
+                    elif idx_scen == 1:
+                        mmangles[idx_scen].value = np.cos(angle)**4
+                        """c_13^4"""
 
                 for idx_sc, sc in enumerate(scan_scales):
                     if args.sens_eval_bin is not None:
@@ -275,6 +279,7 @@ def main():
                             )
                         except:
                             print 'Failed run, continuing'
+                            # raise
                             continue
                         print '## Evidence = {0}'.format(stat)
                     elif args.stat_method is StatCateg.FREQUENTIST:
@@ -289,18 +294,27 @@ def main():
                                 llh_paramset=llh_paramset, fitter=fitter
                             )
                             # print 'llh_paramset', llh_paramset
-                            # print 'llh', llh
+                            # print '-llh', -llh
                             return -llh
 
                         n_params = len(sens_paramset)
-                        x0 = np.full(n_params, 0.7)
                         bounds = np.dstack([np.zeros(n_params), np.ones(n_params)])[0]
+                        x0_arr = np.linspace(0+1e-3, 1-1e-3, 3)
+                        stat = -np.inf
                         try:
-                            res = minimize(fun=fn, x0=x0, method='L-BFGS-B', bounds=bounds)
+                            for x0_v in x0_arr:
+                                x0 = np.full(n_params, x0_v)
+                                # res = minimize(fun=fn, x0=x0, method='Powell', tol=1e-3)
+                                res = minimize(fun=fn, x0=x0, method='Nelder-Mead', tol=1e-3)
+                                # res = minimize(fun=fn, x0=x0, method='L-BFGS-B', tol=1e-3, bounds=bounds)
+                                # res = minimize(fun=fn, x0=x0, method='BFGS', tol=1e-3)
+                                s = -fn(res.x)
+                                if s > stat:
+                                    stat = s
                         except AssertionError:
                             print 'Failed run, continuing'
+                            # raise
                             continue
-                        stat = -fn(res.x)
                         print '=== final llh', stat
                     if args.run_method is SensitivityCateg.FULL:
                         statistic_arr[idx_sc] = np.array([sc, stat])
@@ -323,7 +337,7 @@ def main():
         print 'Plotting statistic'
         if args.sens_run: raw = statistic_arr
         else: raw = np.load(out+'.npy')
-        data = ma.masked_invalid(raw, 0)
+        data = ma.masked_invalid(raw)
 
         basename = os.path.dirname(out) + '/statrun/' + os.path.basename(out)
         baseoutfile = basename[:5]+basename[5:].replace('data', 'plots')
diff --git a/submitter/mcmc_dag.py b/submitter/mcmc_dag.py
index 7f6eae6..d06e920 100644
--- a/submitter/mcmc_dag.py
+++ b/submitter/mcmc_dag.py
@@ -35,6 +35,7 @@ GLOBAL_PARAMS.update(dict(
 
 # FR
 dimension         = [3, 6]
+# dimension         = [4, 5, 7, 8]
 GLOBAL_PARAMS.update(dict(
     threads           = 1,
     binning           = '1e4 1e7 5',
@@ -44,7 +45,7 @@ GLOBAL_PARAMS.update(dict(
     spectral_index    = -2,
     fix_mixing        = 'False',
     fix_mixing_almost = 'False',
-    fold_index        = 'False'
+    fold_index        = 'True'
 ))
 
 # Likelihood
@@ -65,7 +66,7 @@ GLOBAL_PARAMS.update(dict(
     plot_elements     = 'False',
 ))
 
-outfile = 'dagman_FR_MCMC.submit'
+outfile = 'dagman_FR_MCMC_{0}.submit'.format(GLOBAL_PARAMS['likelihood'])
 golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
 condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/mcmc_submit.sub'
 
@@ -115,3 +116,5 @@ with open(outfile, 'w') as f:
                 f.write('VARS\tjob{0}\t{1}="{2}"\n'.format(job_number, key, GLOBAL_PARAMS[key]))
             f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outchains))
             job_number += 1
+
+    print 'dag file = {0}'.format(outfile)
diff --git a/submitter/mcmc_submit.sub b/submitter/mcmc_submit.sub
index f006350..c565205 100644
--- a/submitter/mcmc_submit.sub
+++ b/submitter/mcmc_submit.sub
@@ -1,5 +1,5 @@
 Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/fr.py
-Arguments = "--ast $(ast) --burnin $(burnin) --data $(data) --dimension $(dimension) --fix-mixing $(fix_mixing) --fix-source-ratio $(fix_source_ratio) --likelihood $(likelihood) --measured-ratio $(mr0) $(mr1) $(mr2) --no-bsm $(no_bsm) --nsteps $(nsteps) --nwalkers $(nwalkers) --outfile $(outfile) --plot-angles $(plot_angles) --plot-elements $(plot_elements) --run-mcmc $(run_mcmc) --scale-region $(scale_region) --seed $(seed) --sigma-ratio $(sigma_ratio) --source-ratio $(sr0) $(sr1) $(sr2) --threads $(threads) --energy-dependance $(energy_dependance) --spectral-index $(spectral_index) --binning $(binning) --fix-mixing-almost $(fix_mixing_almost) --fold-index $(fold_index)"
+Arguments = "--ast $(ast) --burnin $(burnin) --data $(data) --dimension $(dimension) --fix-mixing $(fix_mixing) --fix-source-ratio $(fix_source_ratio) --likelihood $(likelihood) --measured-ratio $(mr0) $(mr1) $(mr2) --no-bsm $(no_bsm) --nsteps $(nsteps) --nwalkers $(nwalkers) --outfile $(outfile) --plot-angles $(plot_angles) --plot-elements $(plot_elements) --run-mcmc $(run_mcmc) --scale-region $(scale_region) --seed $(seed) --sigma-ratio $(sigma_ratio) --source-ratio $(sr0) $(sr1) $(sr2) --threads $(threads) --energy-dependance $(energy_dependance) --spectral-index $(spectral_index) --binning $(binning) --fix-mixing-almost $(fix_mixing_almost) --fold-index $(fold_index) --mcmc-seed-type $(mcmc_seed_type)"
 
 # All logs will go to a single file 
 log    = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).log
@@ -16,7 +16,7 @@ getenv = True
 # +TransferOutput="" 
 Transfer_output_files = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/metaouts/
 
-request_memory = 2GB
+request_memory = 3GB
 request_cpus = 1
 
 Universe        = vanilla 
@@ -33,7 +33,7 @@ Notification    = never
 
 # Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6")
 # Requirements = IS_GLIDEIN
-# Requirements = (OpSysMajorVer =?= 6)
+Requirements = (OpSysMajorVer =?= 6)
 
 # GO! 
 queue
diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py
index 17b063d..00be770 100644
--- a/submitter/sens_dag.py
+++ b/submitter/sens_dag.py
@@ -27,21 +27,23 @@ GLOBAL_PARAMS = {}
 sens_eval_bin = 'all' # set to 'all' to run normally
 GLOBAL_PARAMS.update(dict(
     sens_run      = 'True',
-    run_method    = 'fixed_angle', # full, fixed_angle, corr_angle
+    run_method    = 'full', # full, fixed_angle, corr_angle
     stat_method   = 'frequentist',
-    sens_bins     = 10,
+    sens_bins     = 40,
     seed          = 'None'
 ))
 
 # MultiNest
 GLOBAL_PARAMS.update(dict(
-    mn_live_points = 400,
+    mn_live_points = 800,
     mn_tolerance   = 0.01,
     mn_output      = './mnrun'
 ))
 
 # FR
-dimension         = [3, 6]
+# dimension         = [3, 6]
+dimension         = [4, 5, 7, 8]
+# dimension         = [3, 4, 5, 6, 7, 8]
 GLOBAL_PARAMS.update(dict(
     threads           = 1,
     binning           = '1e4 1e7 5',
@@ -50,12 +52,13 @@ GLOBAL_PARAMS.update(dict(
     energy_dependance = 'spectral',
     spectral_index    = -2,
     fix_mixing        = 'False',
-    fix_mixing_almost = 'False'
+    fix_mixing_almost = 'False',
+    fold_index        = 'False'
 ))
 
 # Likelihood
 GLOBAL_PARAMS.update(dict(
-    likelihood  = 'golemfit',
+    likelihood  = 'gaussian',
     sigma_ratio = '0.01'
 ))
 
@@ -70,8 +73,8 @@ GLOBAL_PARAMS.update(dict(
     plot_statistic = 'True'
 ))
 
-outfile = 'dagman_FR_SENS_{0}_{1}.submit'.format(
-    GLOBAL_PARAMS['stat_method'], GLOBAL_PARAMS['run_method']
+outfile = 'dagman_FR_SENS_{0}_{1}_{2}.submit'.format(
+    GLOBAL_PARAMS['stat_method'], GLOBAL_PARAMS['run_method'], GLOBAL_PARAMS['likelihood']
 )
 golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
 condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub'
@@ -141,3 +144,5 @@ with open(outfile, 'w') as f:
                     f.write('VARS\tjob{0}\t{1}="{2}"\n'.format(job_number, key, GLOBAL_PARAMS[key]))
                 f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, output))
                 job_number += 1
+
+    print 'dag file = {0}'.format(outfile)
diff --git a/submitter/sens_submit.sub b/submitter/sens_submit.sub
index 44ec0a3..f3a3e4c 100644
--- a/submitter/sens_submit.sub
+++ b/submitter/sens_submit.sub
@@ -16,7 +16,7 @@ getenv = True
 # +TransferOutput="" 
 Transfer_output_files = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/metaouts/
 
-request_memory = 1GB
+request_memory = 3GB
 request_cpus = 1
 
 Universe        = vanilla 
@@ -33,7 +33,7 @@ Notification    = never
 
 # Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6")
 # Requirements = IS_GLIDEIN
-# Requirements = (OpSysMajorVer =?= 6)
+Requirements = (OpSysMajorVer =?= 6)
 
 # GO! 
 queue
diff --git a/test/test_LV.py b/test/test_LV.py
new file mode 100644
index 0000000..72d0a9c
--- /dev/null
+++ b/test/test_LV.py
@@ -0,0 +1,96 @@
+#!/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.py b/test/test_NSI.py
new file mode 100644
index 0000000..d144420
--- /dev/null
+++ b/test/test_NSI.py
@@ -0,0 +1,106 @@
+#!/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()
+npp = gf.NewPhysicsParams()
+
+steer.quiet = False
+steer.fastmode = True
+
+golem = gf.GolemFit(dp, steer, npp)
+
+fit_params = gf.FitParameters(gf.sampleTag.HESE)
+golem.SetupAsimov(fit_params)
+
+fig = plt.figure(figsize=[6, 5])
+ax = fig.add_subplot(111)
+ax.set_xscale('log')
+ax.set_yscale('log')
+
+binning = golem.GetEnergyBinsMC()
+ax.set_xlim(binning[0], binning[-1])
+# ax.set_ylim(binning[0], binning[-1])
+
+print 'NULL min_llh', golem.MinLLH().likelihood
+
+# 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 expectation', exp
+print
+
+npp.type = gf.NewPhysicsType.NonStandardInteraction
+npp.epsilon_mutau = 0.1
+golem.SetNewPhysicsParams(npp)
+
+print '0.1 mutau min_llh', golem.MinLLH().likelihood
+
+# 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 expectation', exp
+print
+
+np.epsilon_mutau = 0.2
+golem.SetNewPhysicsParams(npp)
+
+print '0.2 mutau min_llh', golem.MinLLH().likelihood
+
+# 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 expectation', exp
+print
+
+np.epsilon_mutau = 0.3
+golem.SetNewPhysicsParams(npp)
+
+print '0.3 mutau min_llh', golem.MinLLH().likelihood
+
+# 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 expectation', exp
+print
+
+np.epsilon_mutau = 0.4
+golem.SetNewPhysicsParams(npp)
+
+print '0.4 mutau min_llh', golem.MinLLH().likelihood
+
+# 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 expectation', exp
+print
+
+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/test/test_all_gf.py b/test/test_all_gf.py
new file mode 100644
index 0000000..04b0980
--- /dev/null
+++ b/test/test_all_gf.py
@@ -0,0 +1,79 @@
+import numpy as np
+import matplotlib as mpl
+mpl.use('Agg')
+import matplotlib.pyplot as plt
+
+import GolemFitPy as gf
+
+FASTMODE = True
+PARAMETERS = [
+    # 'astroFlavorAngle1', 'astroFlavorAngle2',
+    'convNorm',
+    # 'promptNorm', 'muonNorm', 'astroNorm'
+]
+DEFAULTS = [
+    # 4/9., 0., 1., 0., 1., 6.9
+    1.
+]
+RANGES = [
+    # (0, 1), (-1, 1), (0.01, 10), (0., 30), (0.01, 10), (0.01, 30)
+    (0.01, 10)
+]
+BINS = 50
+
+def steering_params():
+    steering_categ = 'p2_0'
+    params = gf.SteeringParams(gf.sampleTag.HESE)
+    if FASTMODE:
+        params.fastmode = True
+    params.quiet = False
+    params.simToLoad= steering_categ.lower()
+    return params
+
+def set_up_as(fitter, params):
+    print 'Injecting the model', params
+    asimov_params = gf.FitParameters(gf.sampleTag.HESE)
+    for x in params.iterkeys():
+        asimov_params.__setattr__(x, float(params[x]))
+    fitter.SetupAsimov(asimov_params)
+    priors = gf.Priors()
+    priors.convNormWidth = 9e9
+    fitter.SetFitPriors(priors)
+
+def setup_fitter(asimov_paramset):
+    datapaths = gf.DataPaths()
+    sparams = steering_params()
+    npp = gf.NewPhysicsParams()
+    fitter = gf.GolemFit(datapaths, sparams, npp)
+    set_up_as(fitter, asimov_paramset)
+    return fitter
+
+def get_llh(fitter, params):
+    fitparams = gf.FitParameters(gf.sampleTag.HESE)
+    for x in params.iterkeys():
+        fitparams.__setattr__(x, float(params[x]))
+    llh = -fitter.EvalLLH(fitparams)
+    return llh
+
+for ip, param in enumerate(PARAMETERS):
+    asimov_paramset = {param: DEFAULTS[ip]}
+    print 'injecting', asimov_paramset
+    fitter = setup_fitter(asimov_paramset)
+    binning = np.linspace(RANGES[ip][0], RANGES[ip][1], BINS)
+    llhs = []
+    for b in binning:
+        test_paramset = {param: b}
+        print 'testing', test_paramset
+        llh = get_llh(fitter, test_paramset)
+        print 'llh', llh
+        llhs.append(llh)
+    plt.plot(binning, llhs)
+    plt.axvline(x=DEFAULTS[ip])
+    plt.xlabel(param)
+    plt.ylabel('LLH')
+    outfile = 'llh_profile_noprior_'
+    if FASTMODE:
+        plt.savefig(outfile + 'fastmode_{0}.png'.format(param))
+    else:
+        plt.savefig(outfile + '{0}.png'.format(param))
+    plt.clf()
diff --git a/test/test_gf.py b/test/test_gf.py
index 6f8f463..4b77924 100644
--- a/test/test_gf.py
+++ b/test/test_gf.py
@@ -1,9 +1,14 @@
+import numpy.ma as ma
+
 import GolemFitPy as gf
 
+FASTMODE = True
+
 def steering_params():
     steering_categ = 'p2_0'
     params = gf.SteeringParams(gf.sampleTag.HESE)
-    params.fastmode = True
+    if FASTMODE:
+        params.fastmode = True
     params.quiet = False
     params.simToLoad= steering_categ.lower()
     return params
@@ -50,7 +55,7 @@ test_paramset = {'astroFlavorAngle1': 4/9., 'astroFlavorAngle2': 0}
 print 'testing', test_paramset
 print 'llh', get_llh(fitter, test_paramset)
 
-test_paramset = {'astroFlavorAngle1': 4/9, 'astroFlavorAngle2': 0}
+test_paramset = {'astroFlavorAngle1': 4/9., 'astroFlavorAngle2': 0}
 print 'testing', test_paramset
 print 'llh', get_llh(fitter, test_paramset)
 
@@ -68,9 +73,17 @@ for an1 in angle1_binning:
     llhs = []
     for an2 in angle2_binning:
         test_paramset = {'astroFlavorAngle1': an1, 'astroFlavorAngle2': an2}
-        llhs.append(get_llh(fitter, test_paramset))
+        llh = get_llh(fitter, test_paramset)
+        if abs(llh) > 9e9:
+            llhs.append(np.nan)
+        else:
+            llhs.append(llh)
+    llhs = ma.masked_invalid(llhs)
     plt.plot(angle2_binning, llhs, label='astroFlavorAngle1 = {0}'.format(an1))
 plt.xlabel('astroFlavorAngle2')
 plt.ylabel('LLH')
 plt.legend()
-plt.savefig('llh_profile_fastmode.png'.format(an1))
+if FASTMODE:
+    plt.savefig('llh_profile_fastmode.png'.format(an1))
+else:
+    plt.savefig('llh_profile.png'.format(an1))
diff --git a/test/test_gf_simple.py b/test/test_gf_simple.py
index 51460c1..3a6ebb5 100644
--- a/test/test_gf_simple.py
+++ b/test/test_gf_simple.py
@@ -1,11 +1,20 @@
+import numpy as np
+import matplotlib as mpl
+mpl.use('Agg')
+import matplotlib.pyplot as plt
+
 import GolemFitPy as gf
 
+FASTMODE = True
+
 dp = gf.DataPaths()
 npp = gf.NewPhysicsParams()
 sp = gf.SteeringParams(gf.sampleTag.HESE)
 
 sp.quiet = False
-# sp.fastmode = True
+if FASTMODE:
+    sp.fastmode = True
+sp.frequentist = True
 
 golem = gf.GolemFit(dp, sp, npp)
 
@@ -20,4 +29,33 @@ fp_sh.astroFlavorAngle1 = 0.36
 fp_sh.astroFlavorAngle2 = -0.57
 
 print 'Eval fp = {0}'.format(golem.EvalLLH(fp))
+
+# energy_centers = golem.GetEnergyBinsMC()[:-1]+ np.diff(golem.GetEnergyBinsMC())/2.
+
+# plt.hist(energy_centers,bins=golem.GetEnergyBinsMC(),
+#          weights=np.sum(golem.GetExpectation(fp),axis=(0,1,2,3)),
+#          histtype="step", lw = 2, label='injected')
+
+# data_energy_dist = np.sum(golem.GetDataDistribution(),axis=(0,1,2,3))
+# energy_centers=golem.GetEnergyBinsData()[:-1]+ np.diff(golem.GetEnergyBinsData())/2.
+# plt.errorbar(energy_centers,data_energy_dist,yerr = np.sqrt(data_energy_dist),fmt='o')
+
 print 'Eval fp_sh = {0}'.format(golem.EvalLLH(fp_sh))
+
+# plt.hist(energy_centers,bins=golem.GetEnergyBinsMC(),
+#          weights=np.sum(golem.GetExpectation(fp_sh),axis=(0,1,2,3)),
+#          histtype="step", lw = 2, label='test')
+
+# data_energy_dist = np.sum(golem.GetDataDistribution(),axis=(0,1,2,3))
+# energy_centers=golem.GetEnergyBinsData()[:-1]+ np.diff(golem.GetEnergyBinsData())/2.
+# plt.errorbar(energy_centers,data_energy_dist,yerr = np.sqrt(data_energy_dist),fmt='o')
+
+# plt.loglog(nonposy="clip")
+# plt.xlabel(r"Deposited energy/GeV")
+# plt.ylabel(r"Events")
+
+# outname = 'Expectation'
+# if FASTMODE:
+#     plt.savefig(outname + 'fastmode.png')
+# else:
+#     plt.savefig(outname + '.png')
diff --git a/utils/fr.py b/utils/fr.py
index 61b5a61..5acd6f8 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -212,7 +212,7 @@ def estimate_scale(args):
         upper_s = (m_eign/args.binning[0]) / (args.binning[0]**(args.dimension-3))
         scale = np.power(10, np.average(np.log10([lower_s, upper_s])))
         diff = upper_s / lower_s
-        scale_region = (lower_s/diff, upper_s*diff)
+        scale_region = (lower_s/np.power(10, args.dimension), upper_s*diff*np.power(10, args.dimension))
         scale_region = [np.power(10, np.round(np.log10(x))) for x in scale_region]
     return scale, scale_region
 
diff --git a/utils/gf.py b/utils/gf.py
index ea4f61f..098faa9 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -57,9 +57,14 @@ def steering_params(args):
     steering_categ = args.ast
     params = gf.SteeringParams()
     params.quiet = False
+    params.fastmode = True
     params.simToLoad= steering_categ.name.lower()
     params.evalThreads = args.threads
     # params.evalThreads = thread_factors(args.threads)[1]
+
+    # For Tianlu
+    # params.years = [999]
+
     return params
 
 
diff --git a/utils/likelihood.py b/utils/likelihood.py
index 2e4a22d..c2822a1 100644
--- a/utils/likelihood.py
+++ b/utils/likelihood.py
@@ -167,9 +167,6 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
     for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)):
         param.value = flavour_angles[idx]
 
-    print 'llh_paramset', llh_paramset
-    print 'hypo_paramset', hypo_paramset
-    print 'fr', fr
     if args.likelihood is Likelihood.FLAT:
         llh = 1.
     elif args.likelihood is Likelihood.GAUSSIAN:
@@ -179,7 +176,6 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
         llh = gf_utils.get_llh(fitter, hypo_paramset)
     elif args.likelihood is Likelihood.GF_FREQ:
         lhh = gf_utils.get_llh_freq(fitter, hypo_paramset)
-    print 'llh', llh
     return llh
 
 
diff --git a/utils/misc.py b/utils/misc.py
index cad03bc..5483675 100644
--- a/utils/misc.py
+++ b/utils/misc.py
@@ -12,6 +12,7 @@ from __future__ import absolute_import, division
 import os
 import errno
 import multiprocessing
+from fractions import gcd
 
 import argparse
 from operator import attrgetter
@@ -33,14 +34,18 @@ class SortingHelpFormatter(argparse.HelpFormatter):
         super(SortingHelpFormatter, self).add_arguments(actions)
 
 
+def solve_ratio(fr):
+    denominator = reduce(gcd, fr)
+    return [int(x/denominator) for x in fr]
+
+
 def gen_identifier(args):
     f = '_DIM{0}'.format(args.dimension)
-    mr1, mr2, mr3 = args.measured_ratio
+    mr1, mr2, mr3 = solve_ratio(args.measured_ratio)
     if args.fix_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)
+        sr1, sr2, sr3 = solve_ratio(args.source_ratio)
+        f += '_sfr_{0:G}_{1:G}_{2:G}_mfr_{3:G}_{4:G}_{5:G}'.format(
+            sr1, sr2, sr3, mr1, mr2, mr3
         )
         if args.fix_mixing:
             f += '_fix_mixing'
@@ -128,6 +133,21 @@ def make_dir(outfile):
         else:
             raise
 
+def remove_option(parser, arg):
+    for action in parser._actions:
+        if (vars(action)['option_strings']
+            and vars(action)['option_strings'][0] == arg) \
+                or vars(action)['dest'] == arg:
+            parser._remove_action(action)
+
+    for action in parser._action_groups:
+        vars_action = vars(action)
+        var_group_actions = vars_action['_group_actions']
+        for x in var_group_actions:
+            if x.dest == arg:
+                var_group_actions.remove(x)
+                return
+
 
 def seed_parse(s):
     if s.lower() == 'none':
diff --git a/utils/multinest.py b/utils/multinest.py
index 9dd0742..f3595f9 100644
--- a/utils/multinest.py
+++ b/utils/multinest.py
@@ -16,7 +16,7 @@ import numpy as np
 from pymultinest import analyse, run
 
 from utils import likelihood
-from utils.misc import gen_outfile_name, make_dir
+from utils.misc import gen_identifier, make_dir
 
 
 def CubePrior(cube, ndim, n_params):
@@ -37,13 +37,15 @@ def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset,
         llh_paramset[pm].value = mn_paramset[pm].value
     theta = llh_paramset.values
     # print 'llh_paramset', llh_paramset
-    return likelihood.ln_prob(
+    llh = likelihood.ln_prob(
         theta=theta,
         args=args,
         asimov_paramset=asimov_paramset,
         llh_paramset=llh_paramset,
         fitter=fitter
     )
+    # print 'llh', llh
+    return llh
 
 
 def mn_argparse(parser):
@@ -77,8 +79,8 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter):
         fitter          = fitter
     )
 
-    prefix = './mnrun/DIM{0}/{1}_{2:>019}_'.format(
-        args.dimension, gen_outfile_name(args), np.random.randint(0, 2**63)
+    prefix = './mnrun/DIM{0}/{1}_{2:>010}_'.format(
+        args.dimension, gen_identifier(args), np.random.randint(0, 2**30)
     )
     make_dir(prefix)
     print 'Running evidence calculation for {0}'.format(prefix)
diff --git a/utils/param.py b/utils/param.py
index f861264..25558c0 100644
--- a/utils/param.py
+++ b/utils/param.py
@@ -252,10 +252,16 @@ def get_paramsets(args, nuisance_paramset):
         ])
     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 hasattr(args, 'dimension'):
+            llh_paramset.append(
+                Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3,
+                      tex=r'{\rm log}_{10}\Lambda^{-1}'+get_units(args.dimension), tag=tag)
+            )
+        elif hasattr(args, 'dimensions'):
+            llh_paramset.append(
+                Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3,
+                      tex=r'{\rm log}_{10}\Lambda^{-1} / GeV^{-d+4}', tag=tag)
+            )
     if not args.fix_source_ratio:
         tag = ParamTag.SRCANGLES
         llh_paramset.extend([
diff --git a/utils/plot.py b/utils/plot.py
index 3d94cc1..8f4eba8 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -10,6 +10,7 @@ Plotting functions for the BSM flavour ratio analysis
 from __future__ import absolute_import, division
 
 import os
+from copy import deepcopy
 
 import numpy as np
 import matplotlib as mpl
@@ -112,18 +113,18 @@ def plot_Tchain(Tchain, axes_labels, ranges):
 def gen_figtext(args):
     """Generate the figure text."""
     t = ''
-    mr1, mr2, mr3 = args.measured_ratio
+    mr1, mr2, mr3 = misc_utils.solve_ratio(args.measured_ratio)
     if args.fix_source_ratio:
-        sr1, sr2, sr3 = args.source_ratio
-        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 = misc_utils.solve_ratio(args.source_ratio)
+        t += 'Source flavour ratio = [{0}, {1}, {2}]\nIC ' \
+                'observed flavour ratio = [{3}, {4}, ' \
+                '{5}]\nDimension = {6}'.format(
                     sr1, sr2, sr3, mr1, mr2, mr3, args.dimension,
                     int(args.energy)
                 )
     else:
-        t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \
-                '{2:.2f}]\nDimension = {3}'.format(
+        t += 'IC observed flavour ratio = [{0}, {1}, ' \
+                '{2}]\nDimension = {3}'.format(
                     mr1, mr2, mr3, args.dimension, int(args.energy)
                 )
     if args.fix_scale:
@@ -229,28 +230,33 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset):
             g.export(outfile+'_elements.'+of)
 
 
+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))
+
+
 def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
     """Make MultiNest factor or LLH value plot."""
-    print "Making Statistic plot"
+    print 'Making Statistic plot'
     fig_text = gen_figtext(args)
     if label is not None: fig_text += '\n' + label
 
     print 'data', data
     print 'data.shape', data.shape
     scales, statistic = data.T
+    min_idx = np.argmin(scales)
+    null = statistic[min_idx]
     if args.stat_method is StatCateg.BAYESIAN:
-        min_idx = np.argmin(scales)
-        null = statistic[min_idx]
         reduced_ev = -(statistic - null)
     elif args.stat_method is StatCateg.FREQUENTIST:
-        min_idx = np.argmin(scales)
-        null = statistic[min_idx]
         reduced_ev = -2*(statistic - null)
 
     fig = plt.figure(figsize=(7, 5))
     ax = fig.add_subplot(111)
 
-    ax.set_xlim(scale_param.ranges)
+    ax.set_xlim(np.log10(args.scale_region))
     ax.set_xlabel('$'+scale_param.tex+'$')
     if args.stat_method is StatCateg.BAYESIAN:
         ax.set_ylabel(r'Bayes Factor')
@@ -271,85 +277,302 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
     ax.add_artist(at)
     misc_utils.make_dir(outfile)
     for of in outformat:
+        print 'Saving as {0}'.format(outfile+'.'+of)
         fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150)
 
 
-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))
+def plot_sens_full(data, outfile, outformat, args):
+    print 'Making FULL sensitivity plot'
 
-
-def plot_BSM_angles_limit(dirname, outfile, outformat, args, bayesian):
-    """Make BSM angles vs scale limit plot."""
-    if not args.plot_angles_limit: return
-    print "Making BSM angles limit plot."""
-    fig_text = gen_figtext(args)
-    xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$']
-
-    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.vstack(raw)
-    raw = raw[np.argsort(raw[:,0])]
-    print 'raw', raw
-    print 'raw.shape', raw.shape
-    sc_ranges = (
-        myround(np.min(raw[0][:,0]), up=True),
-        myround(np.max(raw[0][:,0]), down=True)
-    )
-
-    proc = []
-    if bayesian:
-        for idx, theta in enumerate(raw):
-            scale, evidences = theta.T
-            null = evidences[0]
-            reduced_ev = -(evidences - null)
-            al = scale[reduced_ev > np.log(10**(1/2.))]
-            proc.append((idx+1, al[0]))
-    else:
-        for idx, theta in enumerate(raw):
-            scale, llh = theta.T
-            delta_llh = -2 * (llh - np.max(llh))
-            # 90% CL for 1 dof
-            al = scale[delta_llh > 2.71]
-            proc.append((idx+1, al[0]))
-
-    limits = np.array(proc)
-    print 'limits',  limits
+    colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
+    xticks = ['{0}'.format(x) for x in range(np.min(args.dimensions),
+                                             np.max(args.dimensions)+1)]
+    yranges = [np.inf, -np.inf]
+    legend_handles = []
 
     fig = plt.figure(figsize=(7, 5))
     ax = fig.add_subplot(111)
 
-    ax.set_xticklabels(['']+xticks+[''])
-    ax.set_xlim(0, len(xticks)+1)
-    ax.set_ylim(sc_ranges[0], sc_ranges[-1])
-    ax.set_xlabel(r'BSM angle')
-    ylabel = r'${\rm log}_{10} \Lambda' + get_units(args.dimension) + r'$'
-    ax.set_ylabel(ylabel)
+    ax.set_xlim(args.dimensions[0]-1, args.dimensions[-1]+1)
+    ax.set_xticklabels([''] + xticks + [''])
+    ax.set_xlabel(r'BSM operator dimension ' + r'$d$')
+    ax.set_ylabel(r'${\rm log}_{10} \Lambda^{-1} / GeV^{-d+4}$')
+
+    argsc = deepcopy(args)
+    for idim in xrange(len(data)):
+        dim = args.dimensions[idim]
+        print 'dim', dim
+        argsc.dimension = dim
+        for isrc in xrange(len(data[idim])):
+            src = args.source_ratios[isrc]
+            argsc.source_ratio = src
+
+            # fig_text = gen_figtext(argsc)
+            scales, statistic = data[idim][isrc].T
+            min_idx = np.argmin(scales)
+            null = statistic[min_idx]
+            if args.stat_method is StatCateg.BAYESIAN:
+                reduced_ev = -(statistic - null)
+                # al = scales[reduced_ev > np.log(10**(3/2.))] # Strong degree of belief
+                al = scales[reduced_ev > 0.4] # Testing
+            elif args.stat_method is StatCateg.FREQUENTIST:
+                reduced_ev = -2*(statistic - null)
+                al = scales[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks
+            if len(al) == 0:
+                print 'No points for DIM {0} FRS {1} NULL {2}!'.format(
+                    dim, misc_utils.solve_ratio(src), null
+                )
+                print 'Reduced EV {0}'.format(reduced_ev)
+                continue
+            lim = al[0]
+            label = '[{0}, {1}, {2}]'.format(*misc_utils.solve_ratio(src))
+            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[isrc], label=label
+            )
+            ax.add_line(line)
+            if idim == 0: legend_handles.append(line)
+            x_offset = isrc*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[isrc]}
+            )
 
-    for l in limits:
-        line = plt.Line2D(
-            (l[0]-0.1, l[0]+0.1), (l[1], l[1]), lw=3, color='r'
-        )
-        ax.add_line(line)
-        # ax.arrow(
-        #     l[0], l[1], 0, -1.5, head_width=0.05, head_length=0.2, fc='r',
-        #     ec='r', lw=2
-        # )
-        ax.annotate(
-            s='', xy=l, xytext=(l[0], l[1]+1.5),
-            arrowprops={'arrowstyle': '<-', 'lw': 1.5, 'color':'r'}
-        )
+    try:
+        yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True))
+        ax.set_ylim(yranges)
+    except: pass
 
+    ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right')
     for ymaj in ax.yaxis.get_majorticklocs():
-        ax.axhline(y=ymaj, ls='-', color='gray', alpha=0.4, linewidth=1)
+        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)
+        ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1)
 
+    misc_utils.make_dir(outfile)
     for of in outformat:
+        print 'Saving plot as {0}'.format(outfile+'.'+of)
         fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150)
 
+
+def plot_sens_fixed_angle(data, outfile, outformat, args):
+    print 'Making FIXED_ANGLE sensitivity plot'
+
+    colour = {0:'red', 1:'blue', 2:'green', 3:'purple', 4:'orange', 5:'black'}
+    xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$']
+    argsc = deepcopy(args)
+    for idim in xrange(len(data)):
+        dim = args.dimensions[idim]
+        print '= dim', dim
+        argsc.dimension = dim
+
+        yranges = [np.inf, -np.inf]
+        legend_handles = []
+
+        fig = plt.figure(figsize=(7, 5))
+        ax = fig.add_subplot(111)
+        ax.set_xlim(0, len(xticks)+1)
+        ax.set_xticklabels([''] + xticks + [''])
+        ax.set_xlabel(r'BSM operator angle')
+        ax.set_ylabel(r'${\rm log}_{10} \Lambda^{-1} / GeV^{-d+4}$')
+
+        for isrc in xrange(len(data[idim])):
+            src = args.source_ratios[isrc]
+            argsc.source_ratio = src
+            print '== src', src
+
+            for ian in xrange(len(data[idim][isrc])):
+                print '=== an', ian
+                scales, statistic = data[idim][isrc][ian].T
+                min_idx = np.argmin(scales)
+                null = statistic[min_idx]
+                if args.stat_method is StatCateg.BAYESIAN:
+                    reduced_ev = -(statistic - null)
+                    al = scales[reduced_ev > np.log(10**(3/2.))] # Strong degree of belief
+                elif args.stat_method is StatCateg.FREQUENTIST:
+                    reduced_ev = -2*(statistic - null)
+                    al = scales[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks
+                print 'reduced_ev', reduced_ev
+                if len(al) == 0:
+                    print 'No points for DIM {0} FRS {1}\nreduced_ev {2}!'.format(
+                        dim, src, reduced_ev
+                    )
+                    continue
+                lim = al[0]
+                print 'limit = {0}'.format(lim)
+                label = '[{0}, {1}, {2}]'.format(*misc_utils.solve_ratio(src))
+                if lim < yranges[0]: yranges[0] = lim
+                if lim > yranges[1]: yranges[1] = lim+4
+                line = plt.Line2D(
+                    (ian+1-0.1, ian+1+0.1), (lim, lim), lw=3, color=colour[isrc], label=label
+                )
+                ax.add_line(line)
+                if len(legend_handles) < isrc+1:
+                    legend_handles.append(line)
+                x_offset = isrc*0.05 - 0.05
+                ax.annotate(
+                    s='', xy=(ian+1+x_offset, lim), xytext=(ian+1+x_offset, lim+3),
+                    arrowprops={'arrowstyle': '<-', 'lw': 1.2, 'color':colour[isrc]}
+                )
+
+        try:
+            yranges = (myround(yranges[0], up=True), myround(yranges[1], down=True))
+            ax.set_ylim(yranges)
+        except: pass
+
+        ax.legend(handles=legend_handles, prop=dict(size=8), loc='upper right')
+        for ymaj in ax.yaxis.get_majorticklocs():
+            ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1)
+        for xmaj in ax.xaxis.get_majorticklocs():
+            ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1)
+
+        out = outfile + '_DIM{0}'.format(dim)
+        misc_utils.make_dir(out)
+        for of in outformat:
+            print 'Saving plot as {0}'.format(out+'.'+of)
+            fig.savefig(out+'.'+of, bbox_inches='tight', dpi=150)
+
+
+        for ymaj in ax.yaxis.get_majorticklocs():
+            ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1)
+        for xmaj in ax.xaxis.get_majorticklocs():
+            ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1)
+
+
+def plot_sens_corr_angle(data, outfile, outformat, args):
+    print 'Making CORR_ANGLE sensitivity plot'
+
+    labels = [r'$sin^2(2\mathcal{O}_{12})$',
+              r'$sin^2(2\mathcal{O}_{13})$',
+              r'$sin^2(2\mathcal{O}_{23})$']
+
+    argsc = deepcopy(args)
+    for idim in xrange(len(data)):
+        dim = args.dimensions[idim]
+        print '= dim', dim
+        argsc.dimension = dim
+        for isrc in xrange(len(data[idim])):
+            src = args.source_ratios[isrc]
+            argsc.source_ratio = src
+            print '== src', src
+            for ian in xrange(len(data[idim][isrc])):
+                print '=== an', ian
+
+                d = data[idim][isrc][ian].reshape(
+                    len(data[idim][isrc][ian])**2, 3
+                )
+
+                fig = plt.figure(figsize=(7, 5))
+                ax = fig.add_subplot(111)
+                ax.set_ylim(0, 1)
+                ax.set_ylabel(labels[ian])
+                ax.set_xlabel(r'${\rm log}_{10} \Lambda^{-1}'+get_units(dim)+r'$')
+
+                xranges = [np.inf, -np.inf]
+                legend_handles = []
+
+                y = d[:,0]
+                x = d[:,1]
+                z = d[:,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])
+
+                if args.stat_method is StatCateg.BAYESIAN:
+                    z_r = -(z - null)
+                elif args.stat_method is StatCateg.FREQUENTIST:
+                    z_r = -2*(z - null)
+                print 'scale', x
+                print 'reduced ev', z_r
+
+                pdat = np.array([x, y, z_r, np.ones(x.shape)]).T
+                print 'pdat', pdat
+                pdat_clean = []
+                for d in pdat:
+                    if not np.any(np.isnan(d)): pdat_clean.append(d)
+                pdat = np.vstack(pdat_clean)
+                sort_column = 3
+                pdat_sorted = pdat[pdat[:,sort_column].argsort()]
+                uni, c = np.unique(pdat[:,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 pdat_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
+
+                if args.stat_method is StatCateg.BAYESIAN:
+                    reduced_pdat_mask = (sep_arrays[2] > log(10**(3/2.))) # Strong degree of belief
+                elif args.stat_method is StatCateg.FREQUENTIST:
+                    reduced_pdat_mask = (sep_arrays[2] > 4.61) # 90% CL for 2 DOFS via Wilks
+                reduced_pdat = sep_arrays.T[reduced_pdat_mask].T
+                print 'reduced_pdat', reduced_pdat
+
+                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 = reduced_pdat[0].round(decimals=4)
+                y_v = reduced_pdat[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)
+
+                out = outfile + '_DIM{0}_SRC{1}_AN{2}'.format(dim, isrc, ian)
+                misc_utils.make_dir(out)
+                for of in outformat:
+                    print 'Saving plot as {0}'.format(out+'.'+of)
+                    fig.savefig(out+'.'+of, bbox_inches='tight', dpi=150)