Fri Apr 27 18:37:45 CDT 2018

15 files changed, 515 insertions, 433 deletions

diff --git a/sens.py b/sens.py
index 980ef3f..50fe0c8 100755
--- a/sens.py
+++ b/sens.py
@@ -25,6 +25,7 @@ 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.misc import DTYPE
 from utils.param import Param, ParamSet, get_paramsets
 
 from utils import multinest as mn_utils
@@ -34,11 +35,13 @@ 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.29, 0.31], 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.998, 1.0], 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.46, 0.61], ranges=[0., 1.],      std=0.069,   tex=r's_{23}^2', prior=g_prior, tag=tag),
-        Param(name='dcp',    value=4.08404,        seed=[0, 2*np.pi], ranges=[0., 2*np.pi], std=2.0,     tex=r'\delta_{CP}', tag=tag),
+        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
@@ -76,9 +79,9 @@ def process_args(args):
         raise NotImplementedError(
             '--fix-mixing and --fix-mixing-almost cannot be used together'
         )
-    if args.mn_run and args.fix_scale:
+    if args.sens_run and args.fix_scale:
         raise NotImplementedError(
-            '--mn-run and --fix-scale cannot be used together'
+            '--sens-run and --fix-scale cannot be used together'
         )
 
     args.measured_ratio = fr_utils.normalise_fr(args.measured_ratio)
@@ -94,13 +97,13 @@ def process_args(args):
         args.scale = fr_utils.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
+    if args.sens_eval_bin.lower() == 'all':
+        args.sens_eval_bin = None
     else:
-        args.mn_eval_bin = int(args.mn_eval_bin)
+        args.sens_eval_bin = int(args.sens_eval_bin)
 
     if args.stat_method is StatCateg.FREQUENTIST and \
-       args.likelihood is Likelihood.GOLEMFIT::
+       args.likelihood is Likelihood.GOLEMFIT:
         args.likelihood = Likelihood.GF_FREQ
 
 
@@ -123,6 +126,10 @@ def parse_args(args=None):
         help='Path to output chains'
     )
     parser.add_argument(
+        '--sens-run', type=misc_utils.parse_bool, default='True',
+        help='Generate sensitivities'
+    )
+    parser.add_argument(
         '--run-method', default='full',
         type=partial(misc_utils.enum_parse, c=SensitivityCateg),
         choices=SensitivityCateg,
@@ -134,7 +141,15 @@ def parse_args(args=None):
         help='Statistical method to employ'
     )
     parser.add_argument(
-        '--plot-statistic', type=parse_bool, default='False',
+        '--sens-bins', type=int, default=10,
+        help='Number of bins for the Bayes factor plot'
+    )
+    parser.add_argument(
+        '--sens-eval-bin', type=str, default='all',
+        help='Which bin to evalaute for Bayes factor plot'
+    )
+    parser.add_argument(
+        '--plot-statistic', type=misc_utils.parse_bool, default='False',
         help='Plot MultiNest evidence or LLH value'
     )
     fr_utils.fr_argparse(parser)
@@ -161,8 +176,10 @@ def main():
     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 is SensitivityCateg.FIXED_ANGLE or \
-            args.run_method is SensitivityCateg.CORR_ANGLE:
+    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):
@@ -171,55 +188,66 @@ def main():
                           if mmangles[x].name != prms.name])
             )
 
+    scan_scales = np.linspace(
+        np.log10(args.scale_region[0]), np.log10(args.scale_region[1]), args.sens_bins
+    )
+
+    if args.sens_eval_bin is None:
+        eval_dim = args.sens_bins
+    else: eval_dim = 1
+
+    if args.run_method is SensitivityCateg.CORR_ANGLE:
+        scan_angles = np.linspace(0+1e-9, 1-1e-9, eval_dim)
+    else: scan_angles = np.array([0])
+    print 'scan_scales', scan_scales
+    print 'scan_angles', scan_angles
+
     out = args.outfile+'/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \
         + misc_utils.gen_identifier(args)
-    if args.mn_run:
-        if args.likelihood is Likelihood.GOLEMFIT:
+    if args.sens_run:
+        if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]:
             fitter = gf_utils.setup_fitter(args, asimov_paramset)
-            if args.StatCateg is StatCateg.FREQUENTIST:
-                fitter.SetFitParametersFlag(gf_utils.fit_flags(llh_paramset)
+            if args.stat_method is StatCateg.FREQUENTIST:
+                flags, gf_nuisance = gf_utils.fit_flags(llh_paramset)
+                llh_paramset = llh_paramset.remove_params(gf_nuisance)
+                asimov_paramset = asimov_paramset.remove_params(gf_nuisance)
+                st_paramset_arr = [x.remove_params(gf_nuisance)
+                                   for x in st_paramset_arr]
+                fitter.SetFitParametersFlag(flags)
         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:
             statistic_arr = np.full((eval_dim, 2), np.nan)
         elif args.run_method is SensitivityCateg.FIXED_ANGLE:
             statistic_arr = np.full((len(st_paramset_arr), eval_dim, 2), np.nan)
         elif args.run_method is SensitivityCateg.CORR_ANGLE:
-            statistic_arr = np.full((len(st_paramset_arr), eval_dim, eval_dim, 3), np.nan)
-
+            statistic_arr = np.full(
+                (len(st_paramset_arr), eval_dim, eval_dim, 3), np.nan
+            )
 
-        for idx_scen, mn_paramset in enumerate(st_paramset_arr):
+        for idx_scen, sens_paramset in enumerate(st_paramset_arr):
             print '|||| SCENARIO = {0}'.format(idx_scen)
-            for x in mmangles: x.value = 0.
-            if args.run_method is SensitivityCateg.FIXED_ANGLE:
+            if args.run_method in [SensitivityCateg.FIXED_ONE_ANGLE, SensitivityCateg.FIXED_ANGLE]:
+                if SensitivityCateg.FIXED_ANGLE:
+                    for x in mmangles: x.value = 0.+1e-9
                 if idx_scen == 0 or idx_scen == 2:
-                    mmangles[idx_scen].value = np.sin(np.pi/2.)**2
+                    mmangles[idx_scen].value = np.sin(np.pi/4., dtype=DTYPE)**2
                     """s_12^2 or s_23^2"""
                 elif idx_scen == 1:
-                    mmangles[idx_scen].value = np.cos(np.pi/2.)**4
+                    mmangles[idx_scen].value = np.cos(np.pi/4., dtype=DTYPE)**4
                     """c_13^4"""
 
             for idx_an, an in enumerate(scan_angles):
-                if args.run_method is SensitivityCateg.CORR_ANGLE:
+                if args.run_method in [SensitivityCateg.CORR_ANGLE,
+                                       SensitivityCateg.CORR_ONE_ANGLE]:
                     print '|||| ANGLE = {0:<04.2}'.format(an)
-                    mmangles[idx_an].value = an
+                    if SensitivityCateg.CORR_ANGLE:
+                        for x in mmangles: x.value = 0.+1e-9
+                    mmangles[idx_scen].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:
+                    if args.sens_eval_bin is not None:
+                        if idx_sc == args.sens_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))
@@ -231,7 +259,7 @@ def main():
                     if args.stat_method is StatCateg.BAYESIAN:
                         try:
                             stat = mn_utils.mn_evidence(
-                                mn_paramset     = mn_paramset,
+                                mn_paramset     = sens_paramset,
                                 llh_paramset    = llh_paramset,
                                 asimov_paramset = asimov_paramset,
                                 args            = args,
@@ -242,26 +270,30 @@ def main():
                             continue
                         print '## Evidence = {0}'.format(stat)
                     elif args.stat_method is StatCateg.FREQUENTIST:
-                        llh_paramset_freq = [x for parm in llh_paramset if
-                                             x.name not in asimov_paramset.names]
                         def fn(x):
-                            for idx, parm in enumerate(llh_paramset_freq):
-                                parm.value = x[idx]
-                            theta = llh_paramset_freq.values
-                            try:
-                                llh = llh_utils.ln_prob(
-                                    theta=theta, args=args, asimov_paramset=asimov_paramset,
-                                    mcmc_paramset=mcmc_paramset_freq, fitter=fitter
-                                )
-                            except:
-                                print 'Failed run, continuing'
-                                return np.inf
+                            pranges = sens_paramset.seeds
+                            for i, name in enumerate(sens_paramset.names):
+                                llh_paramset[name].value = \
+                                    (pranges[i][1]-pranges[i][0])*x[i] + pranges[i][0]
+                            theta = llh_paramset.values
+                            print 'llh_paramset', llh_paramset
+                            llh = llh_utils.ln_prob(
+                                theta=theta, args=args, asimov_paramset=asimov_paramset,
+                                llh_paramset=llh_paramset, fitter=fitter
+                            )
+                            print 'llh', llh
                             return -llh
-                        
-                        x0 = np.average(llh_paramset_freq.values, axis=1)
-                        res = minimize(fun=fn, x0=x0, method='L-BFGS-B',
-                                       bounds=llh_paramset.seed, fitter=fitter)
+
+                        n_params = len(sens_paramset)
+                        x0 = np.full(n_params, 0.7)
+                        bounds = np.dstack([np.zeros(n_params), np.ones(n_params)])[0]
+                        try:
+                            res = minimize(fun=fn, x0=x0, method='L-BFGS-B', bounds=bounds)
+                        except AssertionError:
+                            print 'Failed run, continuing'
+                            continue
                         stat = -fn(res.x)
+                        print '=== final llh', stat
                     if args.run_method is SensitivityCateg.FULL:
                         statistic_arr[idx_sc] = np.array([sc, stat])
                     elif args.run_method is SensitivityCateg.FIXED_ANGLE:
@@ -271,10 +303,10 @@ def main():
 
         misc_utils.make_dir(out)
         print 'Saving to {0}'.format(out+'.npy')
-        np.save(out+'.npy', np.array(statistic_arr))
+        np.save(out+'.npy', statistic_arr)
 
     if args.plot_statistic:
-        if args.mn_run: raw = statistic_arr
+        if args.sens_run: raw = statistic_arr
         else: raw = np.load(out+'.npy')
         data = ma.masked_invalid(raw, 0)
 
@@ -292,9 +324,9 @@ def main():
             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{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}$'
+                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[idx_scen],
                     outfile     = outfile,
@@ -313,14 +345,14 @@ def main():
                 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)
+                    _label = label + r'{0:<04.2}$'.format(an)
                     plot_utils.plot_statistic(
                         data        = data[idx_scen][idx_an][:,1:],
                         outfile     = outfile,
                         outformat   = ['png'],
                         args        = args,
                         scale_param = scale,
-                        label       = label
+                        label       = _label
                     )
 
 
diff --git a/submitter/make_dag.py b/submitter/make_dag.py
deleted file mode 100644
index 0e41c9a..0000000
--- a/submitter/make_dag.py
+++ /dev/null
@@ -1,241 +0,0 @@
-#! /usr/bin/env python
-
-import os
-import numpy as np
-
-a_fr = (1, 2, 0)
-b_fr = (1, 0, 0)
-c_fr = (0, 1, 0)
-d_fr = (0, 0, 1)
-e_fr = (1, 1, 1)
-f_fr = (2, 1, 0)
-g_fr = (1, 1, 0)
-
-full_scan_mfr = [
-    # (1, 1, 1), (1, 1, 0)
-]
-
-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, 0, 1),
-    # (1, 1, 0, 1, 2, 0),
-    # (1, 1, 0, 1, 0, 0),
-    # (1, 1, 0, 0, 1, 0),
-    # (1, 0, 0, 1, 0, 0),
-    # (0, 1, 0, 0, 1, 0),
-    # (1, 2, 0, 1, 2, 0),
-    # (1, 2, 0, 0, 1, 0),
-]
-
-# MCMC
-run_mcmc = 'True'
-burnin   = 2500
-nsteps   = 10000
-nwalkers = 60
-seed     = 'None'
-threads  = 1
-mcmc_seed_type = 'uniform'
-
-# FR
-dimension         = [3, 6]
-energy            = [1e6]
-no_bsm            = 'False'
-sigma_ratio       = ['0.01']
-scale             = "1E-20 1E-30"
-scale_region      = "1E10"
-energy_dependance = 'spectral'
-spectral_index    = -2
-binning           = [1e4, 1e7, 5]
-fix_mixing        = 'False'
-fix_mixing_almost = 'False'
-
-# Likelihood
-likelihood = 'gaussian'
-
-# Nuisance
-convNorm        = 1.
-promptNorm      = 0.
-muonNorm        = 1.
-astroNorm       = 6.9
-astroDeltaGamma = 2.5
-
-# GolemFit
-ast  = 'p2_0'
-data = 'real'
-
-# Bayes Factor
-run_bayes_factor       = 'False'
-run_angles_limit       = 'False'
-run_angles_correlation = 'False'
-bayes_bins             = 100
-bayes_live_points      = 3000
-bayes_tolerance        = 0.01
-bayes_eval_bin         = 'all' # set to 'all' to run normally
-
-# Plot
-plot_angles       = 'True'
-plot_elements     = 'False'
-plot_bayes        = 'False'
-plot_angles_limit = 'False'
-
-outfile = 'dagman_FR.submit'
-golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
-condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/submit.sub'
-
-if bayes_eval_bin != 'all':
-    if run_angles_correlation == 'True':
-        b_runs = bayes_bins**2
-    else:
-        b_runs = bayes_bins
-else: b_runs = 1
-
-with open(outfile, 'w') as f:
-    job_number = 1
-    for dim in dimension:
-        print 'dimension', dim
-        for en in energy:
-            print 'energy {0:.0E}'.format(en)
-
-            if energy_dependance == 'mono':
-                outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/{2:.0E}'.format(likelihood, dim, en)
-            elif energy_dependance == 'spectral':
-                outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format(likelihood, dim, spectral_index)
-
-            bayes_output = 'None'
-            angles_lim_output = 'None'
-            angles_corr_output = 'None'
-            for sig in sigma_ratio:
-                print 'sigma', sig
-                for frs in fix_sfr_mfr:
-                    print frs
-                    outchains = outchain_head + '/fix_ifr/{0}/'.format(str(sig).replace('.', '_'))
-                    if run_bayes_factor == 'True':
-                        bayes_output = outchains + '/bayes_factor/'
-                    if run_angles_limit == 'True':
-                        angles_lim_output = outchains + '/angles_limit/'
-                    if run_angles_correlation == 'True':
-                        angles_corr_output = outchains + '/angles_corr/'
-                    outchains += 'mcmc_chain'
-                    for r in range(b_runs):
-                        print 'run', r
-                        f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
-                        f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
-                        f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
-                        f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
-                        f.write('VARS\tjob{0}\tsigma_ratio="{1}"\n'.format(job_number, sig))
-                        f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, 'True'))
-                        f.write('VARS\tjob{0}\tsr0="{1}"\n'.format(job_number, frs[3]))
-                        f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, frs[4]))
-                        f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, frs[5]))
-                        f.write('VARS\tjob{0}\tfix_scale="{1}"\n'.format(job_number, 'False'))
-                        f.write('VARS\tjob{0}\tscale="{1}"\n'.format(job_number, 0))
-                        f.write('VARS\tjob{0}\tscale_region="{1}"\n'.format(job_number, scale_region))
-                        f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
-                        f.write('VARS\tjob{0}\tenergy="{1}"\n'.format(job_number, en))
-                        f.write('VARS\tjob{0}\tlikelihood="{1}"\n'.format(job_number, likelihood))
-                        f.write('VARS\tjob{0}\tburnin="{1}"\n'.format(job_number, burnin))
-                        f.write('VARS\tjob{0}\tnwalkers="{1}"\n'.format(job_number, nwalkers))
-                        f.write('VARS\tjob{0}\tnsteps="{1}"\n'.format(job_number, nsteps))
-                        f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outchains))
-                        f.write('VARS\tjob{0}\tfix_mixing="{1}"\n'.format(job_number, fix_mixing))
-                        f.write('VARS\tjob{0}\tno_bsm="{1}"\n'.format(job_number, no_bsm))
-                        f.write('VARS\tjob{0}\trun_mcmc="{1}"\n'.format(job_number, run_mcmc))
-                        f.write('VARS\tjob{0}\tastroDeltaGamma="{1}"\n'.format(job_number, astroDeltaGamma))
-                        f.write('VARS\tjob{0}\tastroNorm="{1}"\n'.format(job_number, astroNorm))
-                        f.write('VARS\tjob{0}\tconvNorm="{1}"\n'.format(job_number, convNorm))
-                        f.write('VARS\tjob{0}\tmuonNorm="{1}"\n'.format(job_number, muonNorm))
-                        f.write('VARS\tjob{0}\tpromptNorm="{1}"\n'.format(job_number, promptNorm))
-                        f.write('VARS\tjob{0}\tdata="{1}"\n'.format(job_number, data))
-                        f.write('VARS\tjob{0}\tast="{1}"\n'.format(job_number, ast))
-                        f.write('VARS\tjob{0}\tplot_angles="{1}"\n'.format(job_number, plot_angles))
-                        f.write('VARS\tjob{0}\tplot_elements="{1}"\n'.format(job_number, plot_elements))
-                        f.write('VARS\tjob{0}\tseed="{1}"\n'.format(job_number, seed))
-                        f.write('VARS\tjob{0}\tthreads="{1}"\n'.format(job_number, threads))
-                        f.write('VARS\tjob{0}\tmcmc_seed_type="{1}"\n'.format(job_number, mcmc_seed_type))
-                        f.write('VARS\tjob{0}\tenergy_dependance="{1}"\n'.format(job_number, energy_dependance))
-                        f.write('VARS\tjob{0}\tspectral_index="{1}"\n'.format(job_number, spectral_index))
-                        f.write('VARS\tjob{0}\tbinning_0="{1}"\n'.format(job_number, binning[0]))
-                        f.write('VARS\tjob{0}\tbinning_1="{1}"\n'.format(job_number, binning[1]))
-                        f.write('VARS\tjob{0}\tbinning_2="{1}"\n'.format(job_number, binning[2]))
-                        f.write('VARS\tjob{0}\tfix_mixing_almost="{1}"\n'.format(job_number, fix_mixing_almost))
-                        # f.write('VARS\tjob{0}\trun_bayes_factor="{1}"\n'.format(job_number, run_bayes_factor))
-                        # f.write('VARS\tjob{0}\tbayes_bins="{1}"\n'.format(job_number, bayes_bins))
-                        # f.write('VARS\tjob{0}\tbayes_output="{1}"\n'.format(job_number, bayes_output))
-                        # f.write('VARS\tjob{0}\tbayes_live_points="{1}"\n'.format(job_number, bayes_live_points))
-                        # f.write('VARS\tjob{0}\tbayes_tolerance="{1}"\n'.format(job_number, bayes_tolerance))
-                        # f.write('VARS\tjob{0}\tplot_bayes="{1}"\n'.format(job_number, plot_bayes))
-                        # f.write('VARS\tjob{0}\tbayes_eval_bin="{1}"\n'.format(job_number, r))
-                        # f.write('VARS\tjob{0}\trun_angles_limit="{1}"\n'.format(job_number, run_angles_limit))
-                        # f.write('VARS\tjob{0}\tangles_lim_output="{1}"\n'.format(job_number, angles_lim_output))
-                        # f.write('VARS\tjob{0}\tplot_angles_limit="{1}"\n'.format(job_number, plot_angles_limit))
-                        # f.write('VARS\tjob{0}\trun_angles_correlation="{1}"\n'.format(job_number, run_angles_correlation))
-                        # f.write('VARS\tjob{0}\tangles_corr_output="{1}"\n'.format(job_number, angles_corr_output))
-                        job_number += 1
-
-                for frs in full_scan_mfr:
-                    print frs
-                    outchains = outchain_head + '/full_scan/{0}'.format(str(sig).replace('.', '_'))
-                    if run_bayes_factor == 'True':
-                        bayes_output = outchains + '/bayes_factor/'
-                    if run_angles_limit == 'True':
-                        angles_lim_output = outchains + '/angles_limit/'
-                    if run_angles_correlation == 'True':
-                        angles_corr_output = outchains + '/angles_corr/'
-                    outchains += 'mcmc_chain'
-                    for r in range(b_runs):
-                        print 'run', r
-                        f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
-                        f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
-                        f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
-                        f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
-                        f.write('VARS\tjob{0}\tsigma_ratio="{1}"\n'.format(job_number, sig))
-                        f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, 'False'))
-                        f.write('VARS\tjob{0}\tsr0="{1}"\n'.format(job_number, 0))
-                        f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, 0))
-                        f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, 0))
-                        f.write('VARS\tjob{0}\tfix_scale="{1}"\n'.format(job_number, 'False'))
-                        f.write('VARS\tjob{0}\tscale="{1}"\n'.format(job_number, 0))
-                        f.write('VARS\tjob{0}\tscale_region="{1}"\n'.format(job_number, scale_region))
-                        f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
-                        f.write('VARS\tjob{0}\tenergy="{1}"\n'.format(job_number, en))
-                        f.write('VARS\tjob{0}\tlikelihood="{1}"\n'.format(job_number, likelihood))
-                        f.write('VARS\tjob{0}\tburnin="{1}"\n'.format(job_number, burnin))
-                        f.write('VARS\tjob{0}\tnwalkers="{1}"\n'.format(job_number, nwalkers))
-                        f.write('VARS\tjob{0}\tnsteps="{1}"\n'.format(job_number, nsteps))
-                        f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outchains))
-                        f.write('VARS\tjob{0}\tfix_mixing="{1}"\n'.format(job_number, fix_mixing))
-                        f.write('VARS\tjob{0}\tno_bsm="{1}"\n'.format(job_number, no_bsm))
-                        f.write('VARS\tjob{0}\trun_mcmc="{1}"\n'.format(job_number, run_mcmc))
-                        f.write('VARS\tjob{0}\tastroDeltaGamma="{1}"\n'.format(job_number, astroDeltaGamma))
-                        f.write('VARS\tjob{0}\tastroNorm="{1}"\n'.format(job_number, astroNorm))
-                        f.write('VARS\tjob{0}\tconvNorm="{1}"\n'.format(job_number, convNorm))
-                        f.write('VARS\tjob{0}\tmuonNorm="{1}"\n'.format(job_number, muonNorm))
-                        f.write('VARS\tjob{0}\tpromptNorm="{1}"\n'.format(job_number, promptNorm))
-                        f.write('VARS\tjob{0}\tdata="{1}"\n'.format(job_number, data))
-                        f.write('VARS\tjob{0}\tast="{1}"\n'.format(job_number, ast))
-                        f.write('VARS\tjob{0}\tplot_angles="{1}"\n'.format(job_number, plot_angles))
-                        f.write('VARS\tjob{0}\tplot_elements="{1}"\n'.format(job_number, plot_elements))
-                        f.write('VARS\tjob{0}\tseed="{1}"\n'.format(job_number, seed))
-                        f.write('VARS\tjob{0}\tthreads="{1}"\n'.format(job_number, threads))
-                        f.write('VARS\tjob{0}\tmcmc_seed_type="{1}"\n'.format(job_number, mcmc_seed_type))
-                        f.write('VARS\tjob{0}\tenergy_dependance="{1}"\n'.format(job_number, energy_dependance))
-                        f.write('VARS\tjob{0}\tspectral_index="{1}"\n'.format(job_number, spectral_index))
-                        f.write('VARS\tjob{0}\tbinning_0="{1}"\n'.format(job_number, binning[0]))
-                        f.write('VARS\tjob{0}\tbinning_1="{1}"\n'.format(job_number, binning[1]))
-                        f.write('VARS\tjob{0}\tbinning_2="{1}"\n'.format(job_number, binning[2]))
-                        f.write('VARS\tjob{0}\tfix_mixing_almost="{1}"\n'.format(job_number, fix_mixing_almost))
-                        f.write('VARS\tjob{0}\trun_bayes_factor="{1}"\n'.format(job_number, run_bayes_factor))
-                        f.write('VARS\tjob{0}\tbayes_bins="{1}"\n'.format(job_number, bayes_bins))
-                        f.write('VARS\tjob{0}\tbayes_output="{1}"\n'.format(job_number, bayes_output))
-                        f.write('VARS\tjob{0}\tbayes_live_points="{1}"\n'.format(job_number, bayes_live_points))
-                        f.write('VARS\tjob{0}\tbayes_tolerance="{1}"\n'.format(job_number, bayes_tolerance))
-                        f.write('VARS\tjob{0}\tplot_bayes="{1}"\n'.format(job_number, plot_bayes))
-                        f.write('VARS\tjob{0}\tbayes_eval_bin="{1}"\n'.format(job_number, r))
-                        f.write('VARS\tjob{0}\trun_angles_limit="{1}"\n'.format(job_number, run_angles_limit))
-                        f.write('VARS\tjob{0}\tangles_lim_output="{1}"\n'.format(job_number, angles_lim_output))
-                        f.write('VARS\tjob{0}\tplot_angles_limit="{1}"\n'.format(job_number, plot_angles_limit))
-                        f.write('VARS\tjob{0}\trun_angles_correlation="{1}"\n'.format(job_number, run_angles_correlation))
-                        f.write('VARS\tjob{0}\tangles_corr_output="{1}"\n'.format(job_number, angles_corr_output))
-                        job_number += 1
diff --git a/submitter/mcmc_dag.py b/submitter/mcmc_dag.py
new file mode 100644
index 0000000..79f9b6d
--- /dev/null
+++ b/submitter/mcmc_dag.py
@@ -0,0 +1,116 @@
+#! /usr/bin/env python
+
+import os
+import numpy as np
+
+full_scan_mfr = [
+    # (1, 1, 1), (1, 0, 0)
+]
+
+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, 0, 1),
+    # (1, 1, 0, 1, 2, 0),
+    # (1, 1, 0, 1, 0, 0),
+    # (1, 1, 0, 0, 1, 0),
+    # (1, 0, 0, 1, 0, 0),
+    # (0, 1, 0, 0, 1, 0),
+    # (1, 2, 0, 1, 2, 0),
+    # (1, 2, 0, 0, 1, 0),
+]
+
+GLOBAL_PARAMS = {}
+
+# MCMC
+GLOBAL_PARAMS.update(dict(
+    run_mcmc = 'True',
+    burnin   = 2500,
+    nsteps   = 10000,
+    nwalkers = 60,
+    seed     = 'None',
+    mcmc_seed_type = 'uniform'
+))
+
+# FR
+dimension         = [3, 6]
+GLOBAL_PARAMS.update(dict(
+    threads           = 1,
+    binning           = '1e4 1e7 5',
+    no_bsm            = 'False',
+    scale_region      = "1E10",
+    energy_dependance = 'spectral',
+    spectral_index    = -2,
+    fix_mixing        = 'False',
+    fix_mixing_almost = 'False'
+))
+
+# Likelihood
+GLOBAL_PARAMS.update(dict(
+    likelihood  = 'gaussian',
+    sigma_ratio = '0.01'
+))
+
+# GolemFit
+GLOBAL_PARAMS.update(dict(
+    ast  = 'p2_0',
+    data = 'real'
+))
+
+# Plot
+GLOBAL_PARAMS.update(dict(
+    plot_angles       = 'True',
+    plot_elements     = 'False',
+))
+
+outfile = 'dagman_FR_MCMC.submit'
+golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
+condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/mcmc_submit.sub'
+
+with open(outfile, 'w') as f:
+    job_number = 1
+    for dim in dimension:
+        print 'dimension', dim
+        outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format(
+            GLOBAL_PARAMS['likelihood'], dim, GLOBAL_PARAMS['spectral_index']
+        )
+        for frs in fix_sfr_mfr:
+            print 'frs', frs
+            outchains = outchain_head + '/fix_ifr/'
+            if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
+                outchains += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
+            outchains += 'mcmc_chain'
+            f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
+            f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
+            f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
+            f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
+            f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
+            f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, True))
+            f.write('VARS\tjob{0}\tsr0="{1}"\n'.format(job_number, frs[3]))
+            f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, frs[4]))
+            f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, frs[5]))
+            for key in GLOBAL_PARAMS.keys():
+                f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, GLOBAL_PARAMS[key]))
+            f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outfile))
+            job_number += 1
+
+        for frs in full_scan_mfr:
+            print 'frs', frs
+            outchains = outchain_head + '/full/'
+            if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
+                outchains += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
+            outchains += 'mcmc_chain'
+            f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
+            f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
+            f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
+            f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
+            f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
+            f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, False))
+            f.write('VARS\tjob{0}\tsr0="{1}"\n'.format(job_number, 0))
+            f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, 0))
+            f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, 0))
+            for key in GLOBAL_PARAMS.keys():
+                f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, GLOBAL_PARAMS[key]))
+            f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outfile))
+            job_number += 1
diff --git a/submitter/mcmc_submit.sub b/submitter/mcmc_submit.sub
new file mode 100644
index 0000000..2032cb6
--- /dev/null
+++ b/submitter/mcmc_submit.sub
@@ -0,0 +1,39 @@
+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)"
+
+# All logs will go to a single file 
+log    = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).log
+output = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).out
+error  = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).err
+
+getenv = True
+# environment = "X509_USER_PROXY=x509up_u14830" 
+
+# Stage user cert to the node (Gridftp-Users is already on CVMFS)
+# transfer_input_files = /tmp/x509up_u14830
+
+# but do not try to copy outputs back (see: https://htcondor-wiki.cs.wisc.edu/index.cgi/tktview?tn=3081) 
+# +TransferOutput="" 
+Transfer_output_files = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/metaouts/
+
+request_memory = 1GB
+request_cpus = 1
+
+Universe        = vanilla 
+Notification    = never 
+
+# run on both SL5 and 6
+# +WantRHEL6 = True
+# +WantSLC6 = False
+
+# # run on OSG
+# +WantGlidein = True
+
+# +TransferOutput=""
+
+# Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6")
+# Requirements = IS_GLIDEIN
+# Requirements = (OpSysMajorVer =?= 6)
+
+# GO! 
+queue
diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py
new file mode 100644
index 0000000..2cecfe9
--- /dev/null
+++ b/submitter/sens_dag.py
@@ -0,0 +1,134 @@
+#! /usr/bin/env python
+
+import os
+import numpy as np
+
+full_scan_mfr = [
+    # (1, 1, 1), (1, 2, 0)
+]
+
+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, 0, 1),
+    # (1, 1, 0, 1, 2, 0),
+    # (1, 1, 0, 1, 0, 0),
+    # (1, 1, 0, 0, 1, 0),
+    # (1, 0, 0, 1, 0, 0),
+    # (0, 1, 0, 0, 1, 0),
+    # (1, 2, 0, 1, 2, 0),
+    # (1, 2, 0, 0, 1, 0),
+]
+
+GLOBAL_PARAMS = {}
+
+# Bayes Factor
+GLOBAL_PARAMS.update(dict(
+    sens_run      = 'True',
+    run_method    = 'full',
+    stat_method   = 'bayesian',
+    sens_bins     = 10,
+    sens_eval_bin = 'all' # set to 'all' to run normally
+))
+
+# MultiNest
+GLOBAL_PARAMS.update(dict(
+    mn_live_points = 400,
+    mn_tolerance   = 0.01,
+    mn_output      = './mnrun'
+))
+
+# FR
+dimension         = [3, 6]
+GLOBAL_PARAMS.update(dict(
+    threads           = 1,
+    binning           = '1e4 1e7 5',
+    no_bsm            = 'False',
+    scale_region      = "1E10",
+    energy_dependance = 'spectral',
+    spectral_index    = -2,
+    fix_mixing        = 'False',
+    fix_mixing_almost = 'False'
+))
+
+# Likelihood
+GLOBAL_PARAMS.update(dict(
+    likelihood  = 'gaussian',
+    sigma_ratio = '0.01'
+))
+
+# GolemFit
+GLOBAL_PARAMS.update(dict(
+    ast  = 'p2_0',
+    data = 'real'
+))
+
+# Plot
+GLOBAL_PARAMS.update(dict(
+    plot_statistic = 'True'
+))
+
+outfile = 'dagman_FR_SENS.submit'
+golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
+condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub'
+
+if GLOBAL_PARAMS['sens_eval_bin'].lower() != 'all':
+    if GLOBAL_PARAMS['run_method'].lower() == 'corr_angle':
+        sens_runs = GLOBAL_PARAMS['sens_bins']**2
+    else:
+        sens_runs = GLOBAL_PARAMS['sens_bins']
+else: sens_runs = 1
+
+with open(outfile, 'w') as f:
+    job_number = 1
+    for dim in dimension:
+        print 'dimension', dim
+        outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format(
+            GLOBAL_PARAMS['likelihood'], dim, GLOBAL_PARAMS['spectral_index']
+        )
+        for frs in fix_sfr_mfr:
+            print 'frs', frs
+            output = outchain_head + '/fix_ifr/'
+            if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
+                output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
+            output += 'fr_stat'
+            for r in xrange(sens_runs):
+                print 'run', r
+                f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
+                f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
+                f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
+                f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
+                f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
+                f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, True))
+                f.write('VARS\tjob{0}\tsr0="{1}"\n'.format(job_number, frs[3]))
+                f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, frs[4]))
+                f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, frs[5]))
+                f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, r))
+                for key in GLOBAL_PARAMS.keys():
+                    f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, GLOBAL_PARAMS[key]))
+                f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outfile))
+                job_number += 1
+
+        for frs in full_scan_mfr:
+            print 'frs', frs
+            output = outchain_head + '/full/'
+            if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
+                output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
+            output += 'fr_stat'
+            for r in xrange(sens_runs):
+                print 'run', r
+                f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
+                f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
+                f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
+                f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
+                f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
+                f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, False))
+                f.write('VARS\tjob{0}\tsr0="{1}"\n'.format(job_number, 0))
+                f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, 0))
+                f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, 0))
+                f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, r))
+                for key in GLOBAL_PARAMS.keys():
+                    f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, GLOBAL_PARAMS[key]))
+                f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outfile))
+                job_number += 1
diff --git a/submitter/sens_submit.sub b/submitter/sens_submit.sub
new file mode 100644
index 0000000..1a02608
--- /dev/null
+++ b/submitter/sens_submit.sub
@@ -0,0 +1,39 @@
+Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/fr.py
+Arguments = "--ast $(ast) --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) --outfile $(outfile) --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) --sens-run $(sens_run) --run-method $(run_method) --stat-method $(stat_method) --sens-bins $(sens_bins) --sens-eval-bin $(sens_eval_bin) --mn-live-points $(mn_live_points) --mn-tolerance $(mn_tolerance) --mn-output $(mn_output) --plot-statistic $(plot_statistic)"
+
+# All logs will go to a single file 
+log    = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).log
+output = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).out
+error  = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).err
+
+getenv = True
+# environment = "X509_USER_PROXY=x509up_u14830" 
+
+# Stage user cert to the node (Gridftp-Users is already on CVMFS)
+# transfer_input_files = /tmp/x509up_u14830
+
+# but do not try to copy outputs back (see: https://htcondor-wiki.cs.wisc.edu/index.cgi/tktview?tn=3081) 
+# +TransferOutput="" 
+Transfer_output_files = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/metaouts/
+
+request_memory = 1GB
+request_cpus = 1
+
+Universe        = vanilla 
+Notification    = never 
+
+# run on both SL5 and 6
+# +WantRHEL6 = True
+# +WantSLC6 = False
+
+# # run on OSG
+# +WantGlidein = True
+
+# +TransferOutput=""
+
+# Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6")
+# Requirements = IS_GLIDEIN
+# Requirements = (OpSysMajorVer =?= 6)
+
+# GO! 
+queue
diff --git a/submitter/submit.sub b/submitter/submit.sub
deleted file mode 100644
index 8e3c5af..0000000
--- a/submitter/submit.sub
+++ /dev/null
@@ -1,41 +0,0 @@
-Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/fr.py
-Arguments = "--ast $(ast) --astroDeltaGamma $(astroDeltaGamma) --astroNorm $(astroNorm) --burnin $(burnin) --convNorm $(convNorm) --data $(data) --dimension $(dimension) --energy $(energy) --fix-mixing $(fix_mixing) --fix-scale $(fix_scale) --fix-source-ratio $(fix_source_ratio) --likelihood $(likelihood) --measured-ratio $(mr0) $(mr1) $(mr2) --muonNorm $(muonNorm) --no-bsm $(no_bsm) --nsteps $(nsteps) --nwalkers $(nwalkers) --outfile $(outfile) --plot-angles $(plot_angles) --plot-elements $(plot_elements) --promptNorm $(promptNorm) --run-mcmc $(run_mcmc) --scale $(scale) --scale-region $(scale_region) --seed $(seed) --sigma-ratio $(sigma_ratio) --source-ratio $(sr0) $(sr1) $(sr2) --threads $(threads) --mcmc-seed-type $(mcmc_seed_type) --energy-dependance $(energy_dependance) --spectral-index $(spectral_index) --binning $(binning_0) $(binning_1) $(binning_2) --fix-mixing-almost $(fix_mixing_almost)"
-
-# --run-bayes-factor $(run_bayes_factor) --bayes-bins $(bayes_bins) --bayes-output $(bayes_output) --bayes-live-points $(bayes_live_points) --plot-bayes $(plot_bayes) --bayes-tolerance $(bayes_tolerance) --bayes-eval-bin $(bayes_eval_bin) --run-angles-limit $(run_angles_limit) --angles-lim-out $(angles_lim_output) --plot-angles-limit $(plot_angles_limit) --run-angles-correlation $(run_angles_correlation) --angles-corr-output $(angles_corr_output)"
-
-# All logs will go to a single file 
-log    = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).log
-output = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).out
-error  = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).err
-
-getenv = True
-# environment = "X509_USER_PROXY=x509up_u14830" 
-
-# Stage user cert to the node (Gridftp-Users is already on CVMFS)
-# transfer_input_files = /tmp/x509up_u14830
-
-# but do not try to copy outputs back (see: https://htcondor-wiki.cs.wisc.edu/index.cgi/tktview?tn=3081) 
-# +TransferOutput="" 
-Transfer_output_files = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/metaouts/
-
-request_memory = 1GB
-request_cpus = 4
-
-Universe        = vanilla 
-Notification    = never 
-
-# run on both SL5 and 6
-# +WantRHEL6 = True
-# +WantSLC6 = False
-
-# # run on OSG
-# +WantGlidein = True
-
-# +TransferOutput=""
-
-# Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6")
-# Requirements = IS_GLIDEIN
-# Requirements = (OpSysMajorVer =?= 6)
-
-# GO! 
-queue
diff --git a/utils/enums.py b/utils/enums.py
index 8a9e868..b80b712 100644
--- a/utils/enums.py
+++ b/utils/enums.py
@@ -39,8 +39,9 @@ class ParamTag(Enum):
 
 
 class PriorsCateg(Enum):
-    UNIFORM  = 1
-    GAUSSIAN = 2
+    UNIFORM   = 1
+    GAUSSIAN  = 2
+    HALFGAUSS = 3
 
 
 class MCMCSeedType(Enum):
@@ -49,9 +50,11 @@ class MCMCSeedType(Enum):
 
 
 class SensitivityCateg(Enum):
-    FULL        = 1
-    FIXED_ANGLE = 2
-    CORR_ANGLE  = 3
+    FULL            = 1
+    FIXED_ANGLE     = 2
+    CORR_ANGLE      = 3
+    FIXED_ONE_ANGLE = 4
+    CORR_ONE_ANGLE  = 5
 
 
 class StatCateg(Enum):
diff --git a/utils/fr.py b/utils/fr.py
index 342a848..a82e081 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -13,9 +13,9 @@ import argparse
 from functools import partial
 
 import numpy as np
-from scipy import linalg
 
 from utils.enums import EnergyDependance
+from utils.misc import DTYPE, CDTYPE, PI
 from utils.misc import enum_parse, parse_bool
 
 
@@ -44,11 +44,11 @@ def angles_to_fr(src_angles):
     """
     sphi4, c2psi = src_angles
 
-    psi = (0.5)*np.arccos(c2psi)
+    psi = (0.5)*np.arccos(c2psi, dtype=DTYPE)
 
-    sphi2 = np.sqrt(sphi4)
+    sphi2 = np.sqrt(sphi4, dtype=DTYPE)
     cphi2 = 1. - sphi2
-    spsi2 = np.sin(psi)**2
+    spsi2 = np.sin(psi, dtype=DTYPE)**2
     cspi2 = 1. - spsi2
 
     x = sphi2*cspi2
@@ -80,16 +80,16 @@ def angles_to_u(bsm_angles):
 
     """
     s12_2, c13_4, s23_2, dcp = bsm_angles
-    dcp = np.complex128(dcp)
+    dcp = CDTYPE(dcp)
 
     c12_2 = 1. - s12_2
-    c13_2 = np.sqrt(c13_4)
+    c13_2 = np.sqrt(c13_4, dtype=DTYPE)
     s13_2 = 1. - c13_2
     c23_2 = 1. - s23_2
 
-    t12 = np.arcsin(np.sqrt(s12_2))
-    t13 = np.arccos(np.sqrt(c13_2))
-    t23 = np.arcsin(np.sqrt(s23_2))
+    t12 = np.arcsin(np.sqrt(s12_2, dtype=DTYPE))
+    t13 = np.arccos(np.sqrt(c13_2, dtype=DTYPE))
+    t23 = np.arcsin(np.sqrt(s23_2, dtype=DTYPE))
 
     c12 = np.cos(t12)
     s12 = np.sin(t12)
@@ -98,9 +98,9 @@ def angles_to_u(bsm_angles):
     c23 = np.cos(t23)
     s23 = np.sin(t23)
 
-    p1 = np.array([[1   , 0   , 0]                   , [0    , c23 , s23] , [0                   , -s23 , c23]] , dtype=np.complex128)
-    p2 = np.array([[c13 , 0   , s13*np.exp(-1j*dcp)] , [0    , 1   , 0]   , [-s13*np.exp(1j*dcp) , 0    , c13]] , dtype=np.complex128)
-    p3 = np.array([[c12 , s12 , 0]                   , [-s12 , c12 , 0]   , [0                   , 0    , 1]]   , dtype=np.complex128)
+    p1 = np.array([[1   , 0   , 0]                   , [0    , c23 , s23] , [0                   , -s23 , c23]] , dtype=CDTYPE)
+    p2 = np.array([[c13 , 0   , s13*np.exp(-1j*dcp)] , [0    , 1   , 0]   , [-s13*np.exp(1j*dcp) , 0    , c13]] , dtype=CDTYPE)
+    p3 = np.array([[c12 , s12 , 0]                   , [-s12 , c12 , 0]   , [0                   , 0    , 1]]   , dtype=CDTYPE)
 
     u = np.dot(np.dot(p1, p2), p3)
     return u
@@ -142,15 +142,15 @@ def cardano_eqn(ham):
 
     a = -np.trace(ham)
     b = (0.5) * ((np.trace(ham))**2 - np.trace(np.dot(ham, ham)))
-    c = -linalg.det(ham)
+    c = -np.linalg.det(ham)
 
     Q = (1/9.) * (a**2 - 3*b)
     R = (1/54.) * (2*a**3 - 9*a*b + 27*c)
     theta = np.arccos(R / np.sqrt(Q**3))
 
     E1 = -2 * np.sqrt(Q) * np.cos(theta/3.) - (1/3.)*a
-    E2 = -2 * np.sqrt(Q) * np.cos((theta - 2.*np.pi)/3.) - (1/3.)*a
-    E3 = -2 * np.sqrt(Q) * np.cos((theta + 2.*np.pi)/3.) - (1/3.)*a
+    E2 = -2 * np.sqrt(Q) * np.cos((theta - 2.*PI)/3.) - (1/3.)*a
+    E3 = -2 * np.sqrt(Q) * np.cos((theta + 2.*PI)/3.) - (1/3.)*a
 
     A1 = ham[1][2] * (ham[0][0] - E1) - ham[1][0]*ham[0][2]
     A2 = ham[1][2] * (ham[0][0] - E2) - ham[1][0]*ham[0][2]
@@ -164,9 +164,9 @@ def cardano_eqn(ham):
     C2 = ham[1][0] * (ham[2][2] - E2) - ham[1][2]*ham[2][0]
     C3 = ham[1][0] * (ham[2][2] - E3) - ham[1][2]*ham[2][0]
 
-    N1 = np.sqrt(abs(A1*B1)**2 + abs(A1*C1)**2 + abs(B1*C1)**2)
-    N2 = np.sqrt(abs(A2*B2)**2 + abs(A2*C2)**2 + abs(B2*C2)**2)
-    N3 = np.sqrt(abs(A3*B3)**2 + abs(A3*C3)**2 + abs(B3*C3)**2)
+    N1 = np.sqrt(np.abs(A1*B1, dtype=DTYPE)**2 + np.abs(A1*C1, dtype=DTYPE)**2 + np.abs(B1*C1, dtype=DTYPE)**2)
+    N2 = np.sqrt(np.abs(A2*B2, dtype=DTYPE)**2 + np.abs(A2*C2, dtype=DTYPE)**2 + np.abs(B2*C2, dtype=DTYPE)**2)
+    N3 = np.sqrt(np.abs(A3*B3, dtype=DTYPE)**2 + np.abs(A3*C3, dtype=DTYPE)**2 + np.abs(B3*C3, dtype=DTYPE)**2)
 
     mm = np.array([
         [np.conjugate(B1)*C1 / N1, np.conjugate(B2)*C2 / N2, np.conjugate(B3)*C3 / N3],
@@ -195,7 +195,7 @@ def normalise_fr(fr):
     array([ 0.16666667,  0.33333333,  0.5       ])
 
     """
-    return np.array(fr) / float(np.sum(fr))
+    return np.array(fr) / np.sum(fr)
 
 
 def estimate_scale(args):
@@ -300,7 +300,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,
                    sm_u=NUFIT_U, no_bsm=False, fix_mixing=False,
                    fix_mixing_almost=False, fix_scale=False, scale=None,
-                   check_uni=True, epsilon=1e-9):
+                   check_uni=True, epsilon=1e-7):
     """Construct the BSM mixing matrix from the BSM parameters.
 
     Parameters
@@ -393,8 +393,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.) < epsilon or \
-           not abs(np.sum(tu) - 3.) < epsilon:
+        if not np.abs(np.trace(tu) - 3.) < epsilon or \
+           not np.abs(np.sum(tu) - 3.) < epsilon:
             raise AssertionError(
                 'Matrix is not unitary!\neg_vector\n{0}\ntest '
                 'u\n{1}'.format(eg_vector, tu)
@@ -427,7 +427,7 @@ def test_unitarity(x, prnt=False):
            [ 0.,  0.,  1.]])
 
     """
-    f = abs(np.dot(x, x.conj().T))
+    f = np.abs(np.dot(x, x.conj().T), dtype=DTYPE)
     if prnt:
         print 'Unitarity test:\n{0}'.format(f)
     return f
@@ -456,7 +456,8 @@ def u_to_fr(source_fr, matrix):
 
     """
     composition = np.einsum(
-        'ai, bi, a -> b', abs(matrix)**2, abs(matrix)**2, source_fr
+        'ai, bi, a -> b',
+        np.abs(matrix)**2, np.abs(matrix)**2, source_fr,
     )
     ratio = composition / np.sum(source_fr)
     return ratio
diff --git a/utils/gf.py b/utils/gf.py
index 59d1033..17ac029 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -16,6 +16,7 @@ import GolemFitPy as gf
 
 from utils.enums import DataType, SteeringCateg
 from utils.misc import enum_parse, thread_factors
+from utils.param import ParamSet
 
 
 def fit_flags(llh_paramset):
@@ -23,9 +24,6 @@ def fit_flags(llh_paramset):
         # False means it's not fixed in minimization
         'astroFlavorAngle1'         : True,
         'astroFlavorAngle2'         : True,
-        'astroENorm'                : True,
-        'astroMuNorm'               : True,
-        'astroTauNorm'              : True,
         'convNorm'                  : True,
         'promptNorm'                : True,
         'muonNorm'                  : True,
@@ -44,9 +42,14 @@ def fit_flags(llh_paramset):
         'astroDeltaGammaSec'        : True
     }
     flags = gf.FitParametersFlag()
+    gf_nuisance = []
     for param in llh_paramset:
-        flags.__setattr__(param.name, False)
-    return flags
+        if param.name in default_flags:
+            print 'Setting param {0:<15} to float in the ' \
+                'minimisation'.format(param.name)
+            flags.__setattr__(param.name, False)
+            gf_nuisance.append(param)
+    return flags, ParamSet(gf_nuisance)
 
 
 def steering_params(args):
@@ -68,7 +71,7 @@ def set_up_as(fitter, params):
     print 'Injecting the model', params
     asimov_params = gf.FitParameters(gf.sampleTag.HESE)
     for parm in params:
-        asimov_params.__setattr__(parm.name, parm.value)
+        asimov_params.__setattr__(parm.name, float(parm.value))
     fitter.SetupAsimov(asimov_params)
 
 
@@ -84,7 +87,7 @@ def setup_fitter(args, asimov_paramset):
 def get_llh(fitter, params):
     fitparams = gf.FitParameters(gf.sampleTag.HESE)
     for parm in params:
-        fitparams.__setattr__(parm.name, parm.value)
+        fitparams.__setattr__(parm.name, float(parm.value))
     llh = -fitter.EvalLLH(fitparams)
     return llh
 
@@ -93,9 +96,10 @@ def get_llh_freq(fitter, params):
     print 'setting to {0}'.format(params)
     fitparams = gf.FitParameters(gf.sampleTag.HESE)
     for parm in params:
-        fitparams.__setattr__(parm.name, parm.value)
+        fitparams.__setattr__(parm.name, float(parm.value))
     fitter.SetFitParametersSeed(fitparams)
-    return -fitter.MinLLH().likelihood
+    llh = -fitter.MinLLH().likelihood
+    return llh
 
 
 def data_distributions(fitter):
diff --git a/utils/likelihood.py b/utils/likelihood.py
index 04828e8..6387a1e 100644
--- a/utils/likelihood.py
+++ b/utils/likelihood.py
@@ -12,7 +12,7 @@ from __future__ import absolute_import, division
 from functools import partial
 
 import numpy as np
-from scipy.stats import multivariate_normal
+from scipy.stats import multivariate_normal, rv_continuous
 
 import GolemFitPy as gf
 
@@ -22,15 +22,16 @@ from utils.enums import EnergyDependance, Likelihood, ParamTag, PriorsCateg
 from utils.misc import enum_parse
 
 
-def multi_gaussian(fr, fr_bf, sigma):
-    """Multivariate gaussian likelihood."""
-    cov_fr = np.identity(3) * sigma
-    return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr))
+class Gaussian(rv_continuous):
+    """Gaussian for one dimension."""
+    def _pdf(self, x, mu, sig):
+        return (1./np.sqrt(2*np.pi*sig**2))*np.exp(-((x-mu)**2)/(2*sig**2))
 
 
-def log_gaussian(x, mu, sig):
-    """Log gaussian for one dimension."""
-    return np.log((1./np.sqrt(2*np.pi*sig**2))*np.exp(-((x-mu)**2)/(2*sig**2)))
+def multi_gaussian(fr, fr_bf, sigma, offset=-320):
+    """Multivariate gaussian likelihood."""
+    cov_fr = np.identity(3) * sigma
+    return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr)) + offset
 
 
 def likelihood_argparse(parser):
@@ -46,6 +47,11 @@ def likelihood_argparse(parser):
 
 def lnprior(theta, paramset):
     """Priors on theta."""
+    if len(theta) != len(paramset):
+        raise AssertionError(
+            'Length of MCMC scan is not the same as the input '
+            'params\ntheta={0}\nparamset={1}'.format(theta, paramset)
+        )
     for idx, param in enumerate(paramset):
         param.value = theta[idx]
     ranges = paramset.ranges
@@ -57,9 +63,13 @@ def lnprior(theta, paramset):
     prior = 0
     for param in paramset:
         if param.prior is PriorsCateg.GAUSSIAN:
-            prior += log_gaussian(
+            prior += Gaussian().logpdf(
                 param.nominal_value, param.value, param.std
             )
+        elif param.prior is PriorsCateg.HALFGAUSS:
+            prior += Gaussian().logpdf(
+                param.nominal_value, param.value, param.std
+            ) + Gaussian().logcdf(1, param.value, param.std)
     return prior
 
 
@@ -68,7 +78,7 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
     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, llh_paramset)
+            'params\ntheta={0}\nparamset]{1}'.format(theta, llh_paramset)
         )
     for idx, param in enumerate(llh_paramset):
         param.value = theta[idx]
diff --git a/utils/misc.py b/utils/misc.py
index 59c5edb..970c693 100644
--- a/utils/misc.py
+++ b/utils/misc.py
@@ -22,6 +22,11 @@ import numpy as np
 from utils.enums import Likelihood
 
 
+DTYPE  = np.float128
+CDTYPE = np.complex128
+PI     = np.arccos(DTYPE(-1))
+
+
 class SortingHelpFormatter(argparse.HelpFormatter):
     """Sort argparse help options alphabetically."""
     def add_arguments(self, actions):
diff --git a/utils/multinest.py b/utils/multinest.py
index 426a951..005a43a 100644
--- a/utils/multinest.py
+++ b/utils/multinest.py
@@ -16,26 +16,11 @@ import numpy as np
 from pymultinest import analyse, run
 
 from utils import likelihood
-from utils.misc import make_dir, parse_bool
+from utils.misc import gen_outfile_name, make_dir, parse_bool
 
 
-def CubePrior(cube, ndim, n_params, mn_paramset):
-    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]
-    prior = 0
-    for parm in mn_paramset:
-        if parm.prior is PriorsCateg.GAUSSIAN:
-            prior_penatly += multivariate_normal(
-                parm.nominal_value, mean=parm.value, cov=parm.std
-            )
-    print 'prior', prior
-    return prior
+def CubePrior(cube, ndim, n_params):
+    pass
 
 
 def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset,
@@ -63,18 +48,6 @@ def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset,
 
 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'
     )
@@ -104,8 +77,8 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter):
         fitter          = fitter
     )
 
-    prefix = './mnrun/DIM{0}/{1:>019}_'.format(
-        args.dimension, np.random.randint(0, 2**63)
+    prefix = './mnrun/DIM{0}/{1}_{2:>019}_'.format(
+        args.dimension, gen_outfile_name(args), np.random.randint(0, 2**63)
     )
     make_dir(prefix)
     print 'Running evidence calculation for {0}'.format(prefix)
diff --git a/utils/param.py b/utils/param.py
index 13f1a63..4d8106e 100644
--- a/utils/param.py
+++ b/utils/param.py
@@ -205,6 +205,13 @@ class ParamSet(Sequence):
         else:
             return ParamSet([io[1] for io in ps])
 
+    def remove_params(self, params):
+        rm_paramset = []
+        for parm in self.params:
+            if parm.name not in params.names:
+                rm_paramset.append(parm)
+        return ParamSet(rm_paramset)
+
 
 def get_paramsets(args, nuisance_paramset):
     """Make the paramsets for generating the Asmimov MC sample and also running
@@ -216,7 +223,7 @@ def get_paramsets(args, nuisance_paramset):
     llh_paramset.extend(
         [x for x in nuisance_paramset.from_tag(ParamTag.SM_ANGLES)]
     )
-    if args.likelihood is Likelihood.GOLEMFIT:
+    if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]:
         gf_nuisance = [x for x in nuisance_paramset.from_tag(ParamTag.NUISANCE)]
         asimov_paramset.extend(gf_nuisance)
         llh_paramset.extend(gf_nuisance)
diff --git a/utils/plot.py b/utils/plot.py
index 6ae8dc2..0c431cf 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -233,10 +233,11 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset):
 
 def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
     """Make MultiNest factor or LLH value plot."""
-    print "Making MultiNest Factor 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
     if args.stat_method is StatCateg.BAYESIAN:
@@ -244,9 +245,9 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
         null = statistic[min_idx]
         reduced_ev = -(statistic - null)
     elif args.stat_method is StatCateg.FREQUENTIST:
-        min_idx = np.argmin(statistic)
+        min_idx = np.argmin(scales)
         null = statistic[min_idx]
-        reduced_ev = 2*(statistic - null)
+        reduced_ev = -2*(statistic - null)
 
     fig = plt.figure(figsize=(7, 5))
     ax = fig.add_subplot(111)
@@ -256,7 +257,7 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
     if args.stat_method is StatCateg.BAYESIAN:
         ax.set_ylabel(r'Bayes Factor')
     elif args.stat_method is StatCateg.FREQUENTIST:
-        ax.set_ylabel(r'$\Delta {\rm LLH}$')
+        ax.set_ylabel(r'$-2\Delta {\rm LLH}$')
 
     ax.plot(scales, reduced_ev)