Wed 20 Mar 12:11:23 CDT 2019

11 files changed, 493 insertions, 45 deletions

diff --git a/contour_emcee.py b/contour_emcee.py
new file mode 100755
index 0000000..5dc3c98
--- /dev/null
+++ b/contour_emcee.py
@@ -0,0 +1,229 @@
+#! /usr/bin/env python
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : November 26, 2018
+
+"""
+HESE flavour ratio contour
+"""
+
+from __future__ import absolute_import, division
+
+import os
+import argparse
+from functools import partial
+
+import numpy as np
+
+from utils import fr as fr_utils
+from utils import gf as gf_utils
+from utils import llh as llh_utils
+from utils import misc as misc_utils
+from utils import mcmc as mcmc_utils
+from utils import plot as plot_utils
+from utils.enums import str_enum
+from utils.enums import DataType, Likelihood, MCMCSeedType, ParamTag, PriorsCateg
+from utils.param import Param, ParamSet, get_paramsets
+
+from pymultinest import Analyzer, run
+
+
+def define_nuisance():
+    """Define the nuisance parameters."""
+    nuisance = []
+    tag = ParamTag.NUISANCE
+    lg_prior = PriorsCateg.LIMITEDGAUSS
+    nuisance.extend([
+        # Param(name='convNorm',                  value=1.,   seed=[0.5,  2. ],  ranges=[0.1,   10.],   std=0.4, prior=lg_prior, tag=tag),
+        # Param(name='promptNorm',                value=0.,   seed=[0.,   6. ],  ranges=[0.,    20.],   std=2.4, prior=lg_prior, tag=tag),
+        Param(name='convNorm',                  value=1.,   seed=[0.5,  2. ],  ranges=[0.1,   10.],   std=0.4, tag=tag),
+        Param(name='promptNorm',                value=0.,   seed=[0.,   6. ],  ranges=[0.,    20.],   std=2.4, tag=tag),
+        Param(name='muonNorm',                  value=1.,   seed=[0.1,  2. ],  ranges=[0.,    10.],   std=0.1, tag=tag),
+        Param(name='astroNorm',                 value=6.9,  seed=[0.,   5. ],  ranges=[0.,    20.],   std=1.5, tag=tag),
+        Param(name='astroDeltaGamma',           value=2.5,  seed=[2.4,  3. ],  ranges=[-5.,   5. ],   std=0.1, tag=tag),
+        Param(name='CRDeltaGamma',              value=0.,   seed=[-0.1, 0.1 ], ranges=[-1.,   1. ],   std=0.1, tag=tag),
+        Param(name='NeutrinoAntineutrinoRatio', value=1.,   seed=[0.8,  1.2 ], ranges=[0.,    2. ],   std=0.1, tag=tag),
+        Param(name='anisotropyScale',           value=1.,   seed=[0.8,  1.2 ], ranges=[0.,    2. ],   std=0.1, tag=tag),
+        Param(name='domEfficiency',             value=0.99, seed=[0.8,  1.2 ], ranges=[0.8,   1.2 ],  std=0.1, tag=tag),
+        Param(name='holeiceForward',            value=0.,   seed=[-0.8, 0.8 ], ranges=[-4.42, 1.58 ], std=0.1, tag=tag),
+        Param(name='piKRatio',                  value=1.0,  seed=[0.8,  1.2 ], ranges=[0.,    2. ],   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.likelihood is not Likelihood.GOLEMFIT \
+       and args.likelihood is not Likelihood.GF_FREQ:
+        raise AssertionError(
+            'Likelihood method {0} not supported for this '
+            'script!\nChoose either GOLEMFIT or GF_FREQ'.format(
+                str_enum(args.likelihood)
+            )
+        )
+
+
+def parse_args(args=None):
+    """Parse command line arguments"""
+    parser = argparse.ArgumentParser(
+        description="BSM flavour ratio analysis",
+        formatter_class=misc_utils.SortingHelpFormatter,
+    )
+    parser.add_argument(
+        '--injected-ratio', type=float, nargs=3, default=[1, 1, 1],
+        help='Set the central value for the injected flavour ratio at IceCube'
+    )
+    parser.add_argument(
+        '--seed', type=misc_utils.seed_parse, default='25',
+        help='Set the random seed value'
+    )
+    parser.add_argument(
+        '--threads', type=misc_utils.thread_type, default='1',
+        help='Set the number of threads to use (int or "max")'
+    )
+    parser.add_argument(
+        '--outfile', type=str, default='./untitled',
+        help='Path to output results'
+    )
+    try:
+        gf_utils.gf_argparse(parser)
+    except: pass
+    llh_utils.likelihood_argparse(parser)
+    mcmc_utils.mcmc_argparse(parser)
+    nuisance_argparse(parser)
+    misc_utils.remove_option(parser, 'sigma_ratio')
+    if args is None: return parser.parse_args()
+    else: return parser.parse_args(args.split())
+
+
+def gen_identifier(args):
+    f = '_{0}_{1}'.format(*map(str_enum, (args.likelihood, args.data)))
+    if args.data is not DataType.REAL:
+        ir1, ir2, ir3 = misc_utils.solve_ratio(args.injected_ratio)
+        f += '_INJ_{0:03d}_{1:03d}_{2:03d}'.format(ir1, ir2, ir3)
+    return f
+
+
+def gen_figtext(args, asimov_paramset):
+    f = ''
+    if args.data is DataType.REAL:
+        f += 'IceCube Preliminary'
+    else:
+        ir1, ir2, ir3 = misc_utils.solve_ratio(args.injected_ratio)
+        f += 'Injected ratio = [{0}, {1}, {2}]'.format(ir1, ir2, ir3)
+        for param in asimov_paramset:
+            f += '\nInjected {0:20s} = {1:.3f}'.format(
+                param.name, param.nominal_value
+            )
+    return f
+
+
+def triangle_llh(theta, args, hypo_paramset, fitter):
+    """Log likelihood function for a given theta."""
+    if len(theta) != len(hypo_paramset):
+        raise AssertionError(
+            'Dimensions of scan is not the same as the input '
+            'params\ntheta={0}\nparamset]{1}'.format(theta, hypo_paramset)
+        )
+    for idx, param in enumerate(hypo_paramset):
+        param.value = theta[idx]
+
+    if args.likelihood is Likelihood.GOLEMFIT:
+        llh = gf_utils.get_llh(fitter, hypo_paramset)
+    elif args.likelihood is Likelihood.GF_FREQ:
+        llh = gf_utils.get_llh_freq(fitter, hypo_paramset)
+
+    return llh
+
+
+def ln_prob(theta, args, hypo_paramset, fitter):
+    lp = llh_utils.lnprior(theta, paramset=hypo_paramset)
+    if not np.isfinite(lp):
+        return -np.inf
+    return lp + triangle_llh(
+        theta,
+        args           = args,
+        hypo_paramset  = hypo_paramset,
+        fitter         = fitter
+    )
+
+
+def main():
+    args = parse_args()
+    process_args(args)
+    misc_utils.print_args(args)
+
+    if args.seed is not None:
+        np.random.seed(args.seed)
+
+    asimov_paramset, hypo_paramset = get_paramsets(args, define_nuisance())
+    hypo_paramset.extend(asimov_paramset.from_tag(ParamTag.BESTFIT))
+    outfile = args.outfile + gen_identifier(args)
+    print '== {0:<25} = {1}'.format('outfile', outfile)
+
+    n_params = len(hypo_paramset)
+    outfile = outfile + '_emcee_'
+
+    print 'asimov_paramset', asimov_paramset
+    print 'hypo_paramset', hypo_paramset
+
+    if args.run_mcmc:
+        fitter = gf_utils.setup_fitter(args, asimov_paramset)
+
+        ln_prob_eval = partial(
+            ln_prob,
+            hypo_paramset  = hypo_paramset,
+            args           = args,
+            fitter         = fitter
+        )
+
+        if args.mcmc_seed_type == MCMCSeedType.UNIFORM:
+            p0 = mcmc_utils.flat_seed(
+                hypo_paramset, nwalkers=args.nwalkers
+            )
+        elif args.mcmc_seed_type == MCMCSeedType.GAUSSIAN:
+            p0 = mcmc_utils.gaussian_seed(
+                hypo_paramset, nwalkers=args.nwalkers
+            )
+
+        samples = mcmc_utils.mcmc(
+            p0       = p0,
+            ln_prob  = ln_prob_eval,
+            ndim     = n_params,
+            nwalkers = args.nwalkers,
+            burnin   = args.burnin,
+            nsteps   = args.nsteps,
+            args     = args,
+            threads  = 1
+            # TODO(shivesh): broken because you cannot pickle a GolemFitPy object
+            # threads      = misc_utils.thread_factors(args.threads)[0]
+        )
+        mcmc_utils.save_chains(samples, outfile)
+
+    of = outfile[:5]+outfile[5:].replace('data', 'plots')+'_posterior'
+    plot_utils.chainer_plot(
+        infile       = outfile+'.npy',
+        outfile      = of,
+        outformat    = ['png'],
+        args         = args,
+        llh_paramset = hypo_paramset,
+        fig_text     = gen_figtext(args, hypo_paramset)
+    )
+
+    print "DONE!"
+
+
+main.__doc__ = __doc__
+
+
+if __name__ == '__main__':
+    main()
diff --git a/fig2.py b/fig2.py
index 90aca90..54a724d 100755
--- a/fig2.py
+++ b/fig2.py
@@ -99,41 +99,16 @@ def main():
     n_params = len(paramset)
     print n_params
 
-    # Data
-    data_path = '{0}/{1}/real'.format(
-        args.contour_dir, str_enum(args.likelihood).lower()
-    )
-    prefix = '{0}/_{1}_REAL_mn_'.format(data_path, str_enum(args.likelihood))
-    analyser = Analyzer(
-        outputfiles_basename=prefix, n_params=n_params
-    )
-    print analyser
-
-    pranges = paramset.ranges
-
-    bf = analyser.get_best_fit()['parameters']
-    for i in xrange(len(bf)):
-        bf[i] = (pranges[i][1]-pranges[i][0])*bf[i] + pranges[i][0]
-    print 'bestfit = ', bf
-    print 'bestfit log_likelihood', analyser.get_best_fit()['log_likelihood']
-
-    print
-    print '{0:50} = {1}'.format('global evidence', analyser.get_stats()['global evidence'])
-    print
-
-    chains = analyser.get_data()[:,2:]
-    for x in chains:
-        for i in xrange(len(x)):
-            x[i] = (pranges[i][1]-pranges[i][0])*x[i] + pranges[i][0]
-
-    llh = -0.5 * analyser.get_data()[:,1]
+    chains = np.load('/data/user/smandalia/flavour_ratio/data/contour_emcee/golemfit/real/_GOLEMFIT_REAL_emcee_.npy')
+    # chains = np.load('/data/user/smandalia/flavour_ratio/data/contour_emcee/golemfit/real/more_sys_flat/_GOLEMFIT_REAL_emcee_.npy')
 
+    print chains
     flavour_angles = chains[:,-2:]
     flavour_ratios = np.array(
         map(fr_utils.angles_to_fr, flavour_angles)
     )
 
-    nbins = 25
+    nbins = 15
 
     fig = plt.figure(figsize=(8, 8))
     ax = fig.add_subplot(111)
@@ -168,7 +143,7 @@ def main():
 
     ax.legend()
 
-    fig.savefig('test.png', bbox_inches='tight', dpi=150)
+    fig.savefig('test_emcee_moresys_flat.png', bbox_inches='tight', dpi=150)
 
     print "DONE!"
 
diff --git a/fig2_austin.py b/fig2_austin.py
index 8a4cc01..cddbed1 100755
--- a/fig2_austin.py
+++ b/fig2_austin.py
@@ -101,7 +101,8 @@ def main():
     print n_params
 
     # Data
-    data_path = '/home/aschneider/programs/GOLEMSPACE/sources/GolemFit/scripts/diffuse/mcmcs/results/dpl_numu_prior_flavor_20190302-162221-a747f528-8aa6-4488-8c80-059572c099fe.json'
+    # data_path = '/home/aschneider/programs/GOLEMSPACE/sources/GolemFit/scripts/diffuse/mcmcs/results/dpl_numu_prior_flavor_20190302-162221-a747f528-8aa6-4488-8c80-059572c099fe.json'
+    data_path = '/home/aschneider/programs/GOLEMSPACE/sources/GolemFit/scripts/diffuse/mcmcs/results/spl_flavor_20190311-170924-5297d736-3c6e-447f-8de7-4a0653a51bb6.json'
     with open(data_path) as f:
         d_json = json.load(f)
 
@@ -110,7 +111,7 @@ def main():
     print 'names', names
     print 'chains.shape', chains.shape
 
-    flavour_angles = chains[:,5:7]
+    flavour_angles = chains[:,4:6]
     flavour_ratios = np.array(
         map(fr_utils.angles_to_fr, flavour_angles)
     )
diff --git a/plot_sens.py b/plot_sens.py
index b12389f..ba3e923 100755
--- a/plot_sens.py
+++ b/plot_sens.py
@@ -363,7 +363,7 @@ def main():
         elif args.run_method in fixed_angle_categ:
             plot_utils.plot_sens_fixed_angle_pretty(
                 data      = data,
-                outfile   = baseoutfile + '/fixed_angle_pretty',
+                outfile   = baseoutfile + '/fixed_angle_pretty_substantial',
                 outformat = ['png', 'pdf'],
                 args      = args,
             )
diff --git a/submitter/contour_dag.py b/submitter/contour_dag.py
new file mode 100644
index 0000000..634801b
--- /dev/null
+++ b/submitter/contour_dag.py
@@ -0,0 +1,73 @@
+#! /usr/bin/env python
+
+import os
+import numpy as np
+
+gfsource      = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
+condor_script = gfsource + '/scripts/flavour_ratio/submitter/contour_submit.sub'
+
+injected_ratios = [
+    (1, 1, 1),
+    (1, 0, 0),
+    (0, 1, 0),
+    (0, 0, 1)
+]
+
+GLOBAL_PARAMS = {}
+
+GLOBAL_PARAMS.update(dict(
+    threads            = 1,
+    save_measured_fr   = 'False',
+    output_measured_fr = './frs/',
+    seed               = None
+))
+
+# MultiNest
+GLOBAL_PARAMS.update(dict(
+    mn_live_points = 5000,
+    mn_tolerance   = 0.3,
+))
+
+# Likelihood
+GLOBAL_PARAMS.update(dict(
+    likelihood  = 'golemfit',
+))
+
+# GolemFit
+GLOBAL_PARAMS.update(dict(
+    ast  = 'p2_0',
+    # data = 'realisation'
+    # data = 'asimov'
+    data = 'real'
+))
+
+# Plot
+GLOBAL_PARAMS.update(dict(
+    plot_chains   = 'False',
+    plot_triangle = 'False'
+))
+
+outfile = 'dagman_FR_CONTOUR_{0}'.format(GLOBAL_PARAMS['data'])
+outfile += '.submit'
+output  = '/data/user/smandalia/flavour_ratio/data/contour/{0}/{1}/'.format(
+    GLOBAL_PARAMS['likelihood'], GLOBAL_PARAMS['data']
+)
+# output += 'nosyst/'
+# output += 'noprompt/'
+# output += 'strictpriors/'
+
+with open(outfile, 'w') as f:
+    job_number = 1
+    for inj in injected_ratios:
+        print 'inj', inj
+        f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
+        f.write('VARS\tjob{0}\tir0="{1}"\n'.format(job_number, inj[0]))
+        f.write('VARS\tjob{0}\tir1="{1}"\n'.format(job_number, inj[1]))
+        f.write('VARS\tjob{0}\tir2="{1}"\n'.format(job_number, inj[2]))
+        for key in GLOBAL_PARAMS.iterkeys():
+            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
+        if GLOBAL_PARAMS['data'] == 'real': break
+
+print 'dag file = {0}'.format(outfile)
diff --git a/submitter/contour_emcee_dag.py b/submitter/contour_emcee_dag.py
new file mode 100644
index 0000000..b16312a
--- /dev/null
+++ b/submitter/contour_emcee_dag.py
@@ -0,0 +1,77 @@
+#! /usr/bin/env python
+
+import os
+import numpy as np
+
+gfsource      = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
+condor_script = gfsource + '/scripts/flavour_ratio/submitter/contour_emcee_submit.sub'
+
+injected_ratios = [
+    (1, 1, 1),
+    (1, 0, 0),
+    (0, 1, 0),
+    (0, 0, 1)
+]
+
+GLOBAL_PARAMS = {}
+
+GLOBAL_PARAMS.update(dict(
+    threads = 1,
+))
+
+# Emcee
+GLOBAL_PARAMS.update(dict(
+    run_mcmc = 'True',
+    burnin   = 250,
+    nsteps   = 500,
+    nwalkers = 60,
+    seed     = 25,
+    mcmc_seed_type = 'uniform'
+))
+
+# Likelihood
+GLOBAL_PARAMS.update(dict(
+    likelihood  = 'golemfit',
+))
+
+# GolemFit
+GLOBAL_PARAMS.update(dict(
+    ast  = 'p2_0',
+    # data = 'realisation'
+    # data = 'asimov'
+    data = 'real'
+))
+
+# Plot
+GLOBAL_PARAMS.update(dict(
+    plot_angles   = 'False',
+    plot_elements = 'False',
+))
+
+outfile = 'dagman_FR_CONTOUR_EMCEE_{0}'.format(GLOBAL_PARAMS['data'])
+outfile += 'more_sys_flat'
+outfile += '.submit'
+
+output  = '/data/user/smandalia/flavour_ratio/data/contour_emcee/{0}/{1}/'.format(
+    GLOBAL_PARAMS['likelihood'], GLOBAL_PARAMS['data']
+)
+# output += 'more_sys/'
+output += 'more_sys_flat/'
+# output += 'noprompt/'
+# output += 'strictpriors/'
+
+with open(outfile, 'w') as f:
+    job_number = 1
+    for inj in injected_ratios:
+        print 'inj', inj
+        f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
+        f.write('VARS\tjob{0}\tir0="{1}"\n'.format(job_number, inj[0]))
+        f.write('VARS\tjob{0}\tir1="{1}"\n'.format(job_number, inj[1]))
+        f.write('VARS\tjob{0}\tir2="{1}"\n'.format(job_number, inj[2]))
+        for key in GLOBAL_PARAMS.iterkeys():
+            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
+        if GLOBAL_PARAMS['data'] == 'real': break
+
+print 'dag file = {0}'.format(outfile)
diff --git a/submitter/contour_emcee_submit.sub b/submitter/contour_emcee_submit.sub
new file mode 100644
index 0000000..df47cb7
--- /dev/null
+++ b/submitter/contour_emcee_submit.sub
@@ -0,0 +1,42 @@
+Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/contour_emcee.py
+Arguments = "--ast $(ast) --data $(data) --likelihood $(likelihood) --injected-ratio $(ir0) $(ir1) $(ir2) --outfile $(outfile) --seed $(seed) --threads $(threads) --run-mcmc $(run_mcmc) --burnin $(burnin) --nsteps $(nsteps) --nwalkers $(nwalkers) --seed $(seed) --mcmc-seed-type $(mcmc_seed_type) --plot-angles $(plot_angles) --plot-elements $(plot_elements)"
+
+# 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 = 3GB
+request_cpus = 1
+
+Universe        = vanilla 
+Notification    = never 
+
+# +AccountingGroup="sanctioned.$ENV(USER)"
+# run on both SL5 and 6
+# +WantRHEL6 = True
+# +WantSLC6 = False
+
+# # run on OSG
+# +WantGlidein = True
+
+# +TransferOutput=""
+
++NATIVE_OS = True
+# 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/contour_submit.sub b/submitter/contour_submit.sub
new file mode 100644
index 0000000..f4e13e9
--- /dev/null
+++ b/submitter/contour_submit.sub
@@ -0,0 +1,42 @@
+Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/contour.py
+Arguments = "--ast $(ast) --data $(data) --likelihood $(likelihood) --injected-ratio $(ir0) $(ir1) $(ir2) --outfile $(outfile) --seed $(seed) --threads $(threads) --mn-live-points $(mn_live_points) --mn-tolerance $(mn_tolerance) --plot-chains $(plot_chains) --plot-triangle $(plot_triangle) --save-measured-fr $(save_measured_fr) --output-measured-fr=$(output_measured_fr)"
+
+# 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 = 3GB
+request_cpus = 1
+
+Universe        = vanilla 
+Notification    = never 
+
+# +AccountingGroup="sanctioned.$ENV(USER)"
+# run on both SL5 and 6
+# +WantRHEL6 = True
+# +WantSLC6 = False
+
+# # run on OSG
+# +WantGlidein = True
+
+# +TransferOutput=""
+
++NATIVE_OS = True
+# 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
index 7423d34..2046efb 100644
--- a/submitter/sens_dag.py
+++ b/submitter/sens_dag.py
@@ -35,8 +35,8 @@ GLOBAL_PARAMS.update(dict(
 
 # MultiNest
 GLOBAL_PARAMS.update(dict(
-    mn_live_points = 1000,
-    # mn_live_points = 500,
+    # mn_live_points = 1000,
+    mn_live_points = 600,
     # mn_live_points = 300,
     # mn_tolerance   = 0.1,
     mn_tolerance   = 0.3,
diff --git a/utils/gf.py b/utils/gf.py
index 2c794d3..13d5728 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -87,8 +87,10 @@ def steering_params(args):
 
     params.sampleToLoad = gf.sampleTag.MagicTau
     params.use_legacy_selfveto_calculation = False
-    params.spline_hole_ice = False
-    params.spline_dom_efficiency = False
+    # params.spline_hole_ice = False
+    # params.spline_dom_efficiency = False
+    params.spline_hole_ice = True
+    params.spline_dom_efficiency = True
 
     return params
 
diff --git a/utils/plot.py b/utils/plot.py
index 91b8b4e..f81f9da 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -41,6 +41,12 @@ from utils.enums import Likelihood, MixingScenario, ParamTag, StatCateg
 from utils.fr import angles_to_u, angles_to_fr
 
 
+bayes_K = 1.   # Substantial degree of belief.
+# bayes_K = 3/2. # Strong degree of belief.
+# bayes_K = 2.   # Very strong degree of belief
+# bayes_K = 5/2. # Decisive degree of belief
+
+
 if os.path.isfile('./plot_llh/paper.mplstyle'):
     plt.style.use('./plot_llh/paper.mplstyle')
 elif os.environ.get('GOLEMSOURCEPATH') is not None:
@@ -218,6 +224,7 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset, fig_text=None):
                      ha='center', va='center')
 
         for of in outformat:
+            print 'Saving', outfile+'_angles.'+of
             g.export(outfile+'_angles.'+of)
 
     if not hasattr(args, 'plot_elements'):
@@ -327,7 +334,7 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
 
     ax.plot(scales_rm, reduced_ev)
 
-    ax.axhline(y=np.log(10**(3/2.)), color='red', alpha=1., linewidth=1.3)
+    ax.axhline(y=np.log(10**(bayes_K)), color='red', alpha=1., linewidth=1.3)
 
     for ymaj in ax.yaxis.get_majorticklocs():
         ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.3, linewidth=1)
@@ -384,7 +391,7 @@ def plot_sens_full(data, outfile, outformat, args):
             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 > np.log(10**(bayes_K))] # Strong degree of belief
                 # al = scales[reduced_ev > 0.4] # Testing
             elif args.stat_method is StatCateg.FREQUENTIST:
                 reduced_ev = -2*(statistic - null)
@@ -468,7 +475,7 @@ def plot_sens_fixed_angle_pretty(data, outfile, outformat, args):
         8: (-21-(en[0]*6), -21-(en[1]+en[1]*6))
     }
 
-    angles = 2
+    angles = 3
 
     colour = {0:'red', 1:'blue', 2:'green'}
     rgb_co = {0:[1,0,0], 1:[0,0,1], 2:[0.0, 0.5019607843137255, 0.0]}
@@ -588,14 +595,14 @@ def plot_sens_fixed_angle_pretty(data, outfile, outformat, args):
                 null = statistic[min_idx]
                 if args.stat_method is StatCateg.BAYESIAN:
                     reduced_ev = -(statistic_rm - null)
-                    al = scales_rm[reduced_ev > np.log(10**(3/2.))] # Strong degree of belief
+                    al = scales_rm[reduced_ev > np.log(10**(bayes_K))] # Strong degree of belief
                 elif args.stat_method is StatCateg.FREQUENTIST:
                     reduced_ev = -2*(statistic_rm - null)
                     al = scales_rm[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks
                 if len(al) == 0:
                     print 'No points for DIM {0} FRS {1}!'.format(dim, src)
                     continue
-                if reduced_ev[-1] < np.log(10**(3/2.)) - 0.1:
+                if reduced_ev[-1] < np.log(10**(bayes_K)) - 0.1:
                     print 'Peaked contour does not exclude large scales! For ' \
                         'DIM {0} FRS{1}!'.format(dim, src)
                     continue
@@ -729,14 +736,14 @@ def plot_sens_fixed_angle(data, outfile, outformat, args):
                 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 > np.log(10**(bayes_K))] # Strong degree of belief
                 elif args.stat_method is StatCateg.FREQUENTIST:
                     reduced_ev = -2*(statistic - null)
                     al = scales[reduced_ev > 2.71] # 90% CL for 1 DOF via Wilks
                 if len(al) == 0:
                     print 'No points for DIM {0} FRS {1}!'.format(dim, src)
                     continue
-                if reduced_ev[-1] < np.log(10**(3/2.)) - 0.1:
+                if reduced_ev[-1] < np.log(10**(bayes_K)) - 0.1:
                     print 'Peaked contour does not exclude large scales! For ' \
                         'DIM {0} FRS{1}!'.format(dim, src)
                     continue
@@ -893,7 +900,7 @@ def plot_sens_corr_angle(data, outfile, outformat, args):
                 print sep_arrays
 
                 if args.stat_method is StatCateg.BAYESIAN:
-                    reduced_pdat_mask = (sep_arrays[2] > np.log(10**(3/2.))) # Strong degree of belief
+                    reduced_pdat_mask = (sep_arrays[2] > np.log(10**(bayes_K))) # 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