Fri Apr 13 22:00:22 CDT 2018

5 files changed, 392 insertions, 158 deletions

diff --git a/fr.py b/fr.py
index 689c8e5..dc4f075 100755
--- a/fr.py
+++ b/fr.py
@@ -10,6 +10,7 @@ HESE BSM flavour ratio analysis script
 
 from __future__ import absolute_import, division
 
+import os
 import argparse
 from functools import partial
 
@@ -33,11 +34,11 @@ def define_nuisance():
     """
     tag = ParamTag.NUISANCE
     nuisance = ParamSet(
-        Param(name='convNorm',        value=1., seed=[0.5, 2. ], ranges=[0. , 50.], std=0.3,  tag=tag),
-        Param(name='promptNorm',      value=0., seed=[0. , 6. ], ranges=[0. , 50.], std=0.05, tag=tag),
-        Param(name='muonNorm',        value=1., seed=[0. , 2. ], ranges=[0. , 50.], std=0.1,  tag=tag),
-        Param(name='astroNorm',       value=1., seed=[4. , 10.], ranges=[0. , 50.], std=0.1,  tag=tag),
-        Param(name='astroDeltaGamma', value=2., seed=[2. , 3. ], ranges=[-5., 5. ], std=0.1,  tag=tag)
+        Param(name='convNorm',        value=1.,  seed=[0.5, 2. ], ranges=[0. , 50.], std=0.3,  tag=tag),
+        Param(name='promptNorm',      value=0.,  seed=[0. , 6. ], ranges=[0. , 50.], std=0.05, tag=tag),
+        Param(name='muonNorm',        value=1.,  seed=[0.1, 2. ], ranges=[0. , 50.], std=0.1,  tag=tag),
+        Param(name='astroNorm',       value=6.9, seed=[0.1, 10.], ranges=[0. , 50.], std=0.1,  tag=tag),
+        Param(name='astroDeltaGamma', value=2.5, seed=[1. , 3. ], ranges=[-5., 5. ], std=0.1,  tag=tag)
     )
     return nuisance
 
@@ -73,11 +74,12 @@ def get_paramsets(args):
 
     if not args.fix_mixing and not args.fix_mixing_almost:
         tag = ParamTag.MMANGLES
+        e = 1e-10
         mcmc_paramset.extend([
-            Param(name='s_12^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag),
-            Param(name='c_13^4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag),
-            Param(name='s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag),
-            Param(name='dcp', value=np.pi, ranges=[0., 2*np.pi], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag)
+            Param(name='s_12^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag),
+            Param(name='c_13^4', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag),
+            Param(name='s_23^2', value=0.5, ranges=[0.+e, 1.-e], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag),
+            Param(name='dcp', value=np.pi, ranges=[0.+e, 2*np.pi-e], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag)
         ])
     if args.fix_mixing_almost:
         tag = ParamTag.MMANGLES
@@ -128,17 +130,22 @@ def process_args(args):
         if args.energy_dependance is EnergyDependance.MONO:
             args.scale = np.power(
                 10, np.round(np.log10(MASS_EIGENVALUES[1]/args.energy)) - \
-                np.log10(args.energy**(args.dimension-3))
+                np.log10(args.energy**(args.dimension-3)) + 6
             )
         elif args.energy_dependance is EnergyDependance.SPECTRAL:
             args.scale = np.power(
                 10, np.round(
                     np.log10(MASS_EIGENVALUES[1]/np.power(10, np.average(np.log10(args.binning)))) \
                     - np.log10(np.power(10, np.average(np.log10(args.binning)))**(args.dimension-3))
-                )
+                ) + 6
             )
             """Estimate the scale of when to expect the BSM term to contribute."""
 
+    if args.bayes_eval_bin.lower() == 'all':
+        args.bayes_eval_bin = None
+    else:
+        args.bayes_eval_bin = int(args.bayes_eval_bin)
+
 
 def parse_args():
     """Parse command line arguments"""
@@ -180,10 +187,18 @@ def parse_args():
         help='Number of bins for the Bayes factor plot'
     )
     parser.add_argument(
+        '--bayes-eval-bin', type=str, default='all',
+        help='Which bin to evalaute for Bayes factor plot'
+    )
+    parser.add_argument(
         '--bayes-live-points', type=int, default=400,
         help='Number of live points for MultiNest evaluations'
     )
     parser.add_argument(
+        '--bayes-tolerance', type=float, default=0.5,
+        help='Tolerance for MultiNest'
+    )
+    parser.add_argument(
         '--bayes-output', type=str, default='./mnrun/',
         help='Folder to store MultiNest evaluations'
     )
@@ -212,7 +227,7 @@ def parse_args():
         help='Set the new physics scale'
     )
     parser.add_argument(
-        '--scale-region', type=float, default=1e10,
+        '--scale-region', type=float, default=5e5,
         help='Set the size of the box to scan for the scale'
     )
     parser.add_argument(
@@ -296,38 +311,43 @@ def main():
         mcmc_paramset = mcmc_paramset
     )
 
+    out = args.bayes_output+'/fr_evidence'
+
+    sc_range = mcmc_paramset.from_tag(ParamTag.SCALE)[0].ranges
+    scales = np.linspace(
+        sc_range[0], sc_range[1], args.bayes_bins
+    )
     if args.run_bayes_factor:
-        # TODO(shivesh)
         import pymultinest
-        from pymultinest.solve import solve
-        from pymultinest.watch import ProgressPrinter
 
         if not args.run_mcmc and args.likelihood is Likelihood.GOLEMFIT:
             fitter = gf_utils.setup_fitter(args, asimov_paramset)
         else:
             fitter = None
 
-        sc_range = mcmc_paramset.from_tag(ParamTag.SCALE)[0].ranges
-        scales = np.linspace(
-            sc_range[0], sc_range[1], args.bayes_bins+1
-        )
-
         p = mcmc_paramset.from_tag(ParamTag.SCALE, invert=True)
         n_params = len(p)
-        prior_ranges = p.ranges
+        prior_ranges = p.seeds
 
         def CubePrior(cube, ndim, nparams):
             # default are uniform priors
             return ;
 
         arr = []
-        for s in scales:
+        for s_idx, s in enumerate(scales):
+            if args.bayes_eval_bin is not None:
+                if args.bayes_eval_bin == s_idx:
+                    out += '_scale_{0:.0E}'.format(np.power(10, s))
+                else: continue
+
+            print '== SCALE = {0:.0E}'.format(np.power(10, s))
             theta = np.zeros(n_params)
             def lnProb(cube, ndim, nparams):
                 for i in range(ndim):
                     prange = prior_ranges[i][1] - prior_ranges[i][0]
                     theta[i] = prange*cube[i] + prior_ranges[i][0]
                 theta_ = np.array(theta.tolist() + [s])
+                # print 'mcmc_paramset', mcmc_paramset
                 return llh_utils.triangle_llh(
                     theta=theta_,
                     args=args,
@@ -336,29 +356,35 @@ def main():
                     fitter=fitter
                 )
 
-            prefix = 'mnrun_{0:.0E}'.format(np.power(10, s)) + '_' + outfile[2:]
+            prefix = 'mnrun_{0:.0E}'.format(np.power(10, s)) + \
+                '_' + os.path.basename(outfile) + '_'
             print 'begin running evidence calculation for {0}'.format(prefix)
             result = pymultinest.run(
                 LogLikelihood=lnProb,
                 Prior=CubePrior,
                 n_dims=n_params,
+                importance_nested_sampling=True,
                 n_live_points=args.bayes_live_points,
+                evidence_tolerance=args.bayes_tolerance,
                 outputfiles_basename=prefix,
-                # resume=False
+                resume=False,
+                verbose=True
             )
 
-            analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params)
+            analyzer = pymultinest.Analyzer(
+                outputfiles_basename=prefix, n_params=n_params
+            )
             a_lnZ = analyzer.get_stats()['global evidence']
             print 'Evidence = {0}'.format(a_lnZ)
-
             arr.append([s, a_lnZ])
-        out = args.bayes_output+'fr_evidence_'+outfile[2:]+'.npy'
+
         misc_utils.make_dir(out)
-        np.save(out, np.array(arr))
+        np.save(out+'.npy', np.array(arr))
 
-    b_out = args.bayes_output+'fr_evidence_'+outfile[2:]
+    dirname = os.path.dirname(out)
     plot_utils.bayes_factor_plot(
-        infile=b_out+'.npy', outfile=b_out, outformat=['png'], args=args
+        dirname=dirname, outfile=out, outformat=['png'], args=args,
+        xlim=sc_range
     )
 
     print "DONE!"
diff --git a/sens_bayes.py b/sens_bayes.py
new file mode 100755
index 0000000..b0030d4
--- /dev/null
+++ b/sens_bayes.py
@@ -0,0 +1,175 @@
+#! /usr/bin/env python
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : April 10, 2018
+
+"""
+HESE BSM flavour ratio analysis script
+"""
+
+from __future__ import absolute_import, division
+
+import numpy as np
+import matplotlib as mpl
+mpl.use('Agg')
+from matplotlib import pyplot as plt
+from matplotlib import rc
+
+import GolemFitPy as gf
+import pymultinest
+from pymultinest.solve import solve
+from pymultinest.watch import ProgressPrinter
+
+import fr
+from utils import gf as gf_utils
+from utils import likelihood as llh_utils
+from utils import misc as misc_utils
+from utils.enums import Likelihood, ParamTag
+from utils.plot import plot_BSM_angles_limit
+
+rc('text', usetex=False)
+rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
+
+
+RUN = False
+
+
+z = 0.
+scenarios = [
+    [np.sin(np.pi/2.)**2, z, z, z],
+    [z, np.cos(np.pi/2.)**4, z, z],
+    [z, z, np.sin(np.pi/2.)**2, z],
+]
+xticks = [r'$\mathcal{O}_{12}$', r'$\mathcal{O}_{13}$', r'$\mathcal{O}_{23}$']
+
+def fit_flags(flag_dict):
+    flags = gf.FitParametersFlag()
+    for key in flag_dict.iterkeys():
+        flags.__setattr__(key, flag_dict[key])
+    return flags
+
+default_flags = {
+    # False means it's not fixed in minimization
+    'astroFlavorAngle1'         : True,
+    'astroFlavorAngle2'         : True,
+    # 'astroENorm'                : True,
+    # 'astroMuNorm'               : True,
+    # 'astroTauNorm'              : True,
+    'convNorm'                  : False,
+    'promptNorm'                : False,
+    'muonNorm'                  : False,
+    'astroNorm'                 : False,
+    'astroParticleBalance'      : True,
+    'astroDeltaGamma'           : False,
+    'cutoffEnergy'              : True,
+    'CRDeltaGamma'              : False,
+    'piKRatio'                  : False,
+    'NeutrinoAntineutrinoRatio' : True,
+    'darkNorm'                  : True,
+    'domEfficiency'             : True,
+    'holeiceForward'            : True,
+    'anisotropyScale'           : True,
+    'astroNormSec'              : True,
+    'astroDeltaGammaSec'        : True
+}
+
+
+def main():
+    args = fr.parse_args()
+    fr.process_args(args)
+    misc_utils.print_args(args)
+
+    bins = 10
+
+    asimov_paramset, mcmc_paramset = fr.get_paramsets(args)
+
+    sc_range = mcmc_paramset.from_tag(ParamTag.SCALE)[0].ranges
+    scan_scales = np.linspace(sc_range[0], sc_range[1], bins)
+    print 'scan_scales', scan_scales
+
+    p = mcmc_paramset.from_tag([ParamTag.SCALE, ParamTag.MMANGLES], invert=True)
+    n_params = len(p)
+    prior_ranges = p.seeds
+
+    outfile = './sens'
+    if RUN:
+        if args.likelihood is Likelihood.GOLEMFIT:
+            fitter = gf_utils.setup_fitter(args, asimov_paramset)
+            fitter.SetFitParametersFlag(fit_flags(default_flags))
+        elif args.likelihood is Likelihood.GAUSSIAN:
+            fitter = None
+
+        def CubePrior(cube, ndim, nparams):
+            # default are uniform priors
+            return ;
+
+        data = np.zeros((len(scenarios), bins, 2))
+        mm_angles = mcmc_paramset.from_tag(ParamTag.MMANGLES)
+        sc_angles = mcmc_paramset.from_tag(ParamTag.SCALE)[0]
+        for idx, scen in enumerate(scenarios):
+            scales = []
+            evidences = []
+            for yidx, an in enumerate(mm_angles):
+                an.value = scen[yidx]
+            for sc in scan_scales:
+                sc_angles.value = sc
+                def lnProb(cube, ndim, nparams):
+                    for i in range(ndim):
+                        prange = prior_ranges[i][1] - prior_ranges[i][0]
+                        p[i].value = prange*cube[i] + prior_ranges[i][0]
+                    for name in p.names:
+                        mcmc_paramset[name].value = p[name].value
+                    theta = mcmc_paramset.values
+                    # print 'theta', theta
+                    # print 'mcmc_paramset', mcmc_paramset
+                    return llh_utils.triangle_llh(
+                        theta=theta,
+                        args=args,
+                        asimov_paramset=asimov_paramset,
+                        mcmc_paramset=mcmc_paramset,
+                        fitter=fitter
+                    )
+                # TODO(shivesh)
+                prefix = 'mnrun_{0:.0E}'.format(np.power(10, sc)) + '_' + misc_utils.gen_outfile_name(args)[2:]
+                print 'begin running evidence calculation for {0}'.format(prefix)
+                result = pymultinest.run(
+                    LogLikelihood=lnProb,
+                    Prior=CubePrior,
+                    n_dims=n_params,
+                    importance_nested_sampling=True,
+                    n_live_points=args.bayes_live_points,
+                    evidence_tolerance=args.bayes_tolerance,
+                    outputfiles_basename=prefix,
+                    resume=False,
+                    verbose=True
+                )
+
+                analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params)
+                a_lnZ = analyzer.get_stats()['global evidence']
+                print 'Evidence = {0}'.format(a_lnZ)
+                scales.append(sc)
+                evidences.append(a_lnZ)
+
+            for i, d in enumerate(evidences):
+                data[idx][i][0] = scales[i]
+                data[idx][i][1] = d
+
+        np.save(outfile + '.npy', data)
+
+    plot_BSM_angles_limit(
+        infile=outfile+'.npy',
+        outfile=outfile,
+        xticks=xticks,
+        outformat=['png'],
+        args=args,
+        bayesian=True
+    )
+
+
+main.__doc__ = __doc__
+
+
+if __name__ == '__main__':
+    main()
+
diff --git a/submitter/make_dag.py b/submitter/make_dag.py
index 641e00e..be13ac8 100644
--- a/submitter/make_dag.py
+++ b/submitter/make_dag.py
@@ -12,7 +12,7 @@ f_fr = (2, 1, 0)
 g_fr = (1, 1, 0)
 
 full_scan_mfr = [
-    (1, 1, 1), (1, 1, 0)
+    # (1, 1, 1), (1, 1, 0)
 ]
 
 fix_sfr_mfr = [
@@ -35,11 +35,11 @@ burnin   = 500
 nsteps   = 2000
 nwalkers = 60
 seed     = 24
-threads  = 12
+threads  = 4
 mcmc_seed_type = 'uniform'
 
 # FR
-dimension         = [3, 6]
+dimension         = [6]
 energy            = [1e6]
 likelihood        = 'golemfit'
 no_bsm            = 'False'
@@ -56,30 +56,35 @@ fix_mixing_almost = 'False'
 likelihood = 'golemfit'
 
 # Nuisance
-astroDeltaGamma = 2.
-astroNorm       = 1.
 convNorm        = 1.
-muonNorm        = 1.
 promptNorm      = 0.
+muonNorm        = 1.
+astroNorm       = 6.9
+astroDeltaGamma = 2.5
 
 # GolemFit
 ast  = 'p2_0'
 data = 'real'
 
 # Bayes Factor
-run_bayes_factor  = 'True'
+run_bayes_factor  = 'False'
 bayes_bins        = 10
 bayes_live_points = 200
+bayes_tolerance   = 0.01
+bayes_eval_bin    = True # set to 'all' to run normally
 
 # Plot
 plot_angles   = 'False'
 plot_elements = 'False'
-plot_bayes    = 'True'
+plot_bayes    = 'False'
 
 outfile = 'dagman_FR.submit'
 golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
 condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/submit.sub'
 
+if bayes_eval_bin != 'all': b_runs = bayes_bins
+else: b_runs = 1
+
 with open(outfile, 'w') as f:
     job_number = 1
     for dim in dimension:
@@ -101,104 +106,112 @@ with open(outfile, 'w') as f:
                     if run_bayes_factor == 'True':
                         bayes_output = outchains + '/bayes_factor/'
                     outchains += 'mcmc_chain'
-                    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}\tlikelihood="{1}"\n'.format(job_number, likelihood))
-                    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}\tplot_bayes="{1}"\n'.format(job_number, plot_bayes))
-                    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/'
-                #     outchains += 'mcmc_chain'
-                #     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}\tlikelihood="{1}"\n'.format(job_number, likelihood))
-                #     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}\tplot_bayes="{1}"\n'.format(job_number, plot_bayes))
-                #     job_number += 1
+                    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}\tlikelihood="{1}"\n'.format(job_number, likelihood))
+                        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))
+                        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/'
+                    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}\tlikelihood="{1}"\n'.format(job_number, likelihood))
+                        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))
+                        job_number += 1
diff --git a/submitter/submit.sub b/submitter/submit.sub
index e563924..e9e66bd 100644
--- a/submitter/submit.sub
+++ b/submitter/submit.sub
@@ -1,5 +1,5 @@
 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) --likelihood $(likelihood) --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)"
+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) --likelihood $(likelihood) --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)"
 
 # All logs will go to a single file 
 log    = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).log
@@ -17,7 +17,7 @@ getenv = True
 transfer_output_files = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/metaouts/
 
 request_memory = 30GB
-request_cpus = 12
+request_cpus = 4
 
 Universe        = vanilla 
 Notification    = never 
diff --git a/utils/plot.py b/utils/plot.py
index 4180eb3..0ff89c2 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -9,6 +9,7 @@ Plotting functions for the BSM flavour ratio analysis
 
 from __future__ import absolute_import, division
 
+import os
 import argparse
 
 import numpy as np
@@ -259,13 +260,20 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
             g.export(outfile+'_elements.'+of)
 
 
-def bayes_factor_plot(infile, outfile, outformat, args):
+def bayes_factor_plot(dirname, outfile, outformat, args, xlim):
     """Make Bayes factor plot."""
     if not args.plot_bayes: return
     print "Making Bayes Factor plot"
+    print 'dirname', dirname
     fig_text = gen_figtext(args)
 
-    raw = np.load(infile)
+    raw = []
+    for root, dirs, filenames in os.walk(dirname):
+        for fn in filenames:
+            if fn[-4:] == '.npy':
+                raw.append(np.load(os.path.join(root, fn)))
+    raw = np.sort(np.vstack(raw), axis=0)
+    print 'raw', raw
     scales, evidences = raw.T
     null = evidences[0]
 
@@ -274,15 +282,16 @@ def bayes_factor_plot(infile, outfile, outformat, args):
     fig = plt.figure(figsize=(7, 5))
     ax = fig.add_subplot(111)
 
-    ax.set_xlabel(r'{\rm log}_{10} \Lambda ' + get_units(args.dimension))
+    ax.set_xlim(xlim)
+    ax.set_xlabel(r'${\rm log}_{10} \Lambda ' + get_units(args.dimension) +r'$')
     ax.set_ylabel(r'Bayes Factor')
 
     ax.plot(scales, reduced_ev)
 
     for ymaj in ax.yaxis.get_majorticklocs():
-        ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.7, 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.7, linewidth=1)
+        ax.axvline(x=xmaj, ls='-', color='gray', alpha=0.4, linewidth=1)
 
     at = AnchoredText(
         fig_text, prop=dict(size=7), frameon=True, loc=2
@@ -300,26 +309,37 @@ def myround(x, base=5, up=False, down=False):
     else: int(base * np.round(float(x)/base))
 
 
-def plot_BSM_angles_limit(infile, outfile, xticks, outformat, args):
+def plot_BSM_angles_limit(infile, outfile, xticks, outformat, args, bayesian):
     """Make BSM angles vs scale limit plot."""
     print "Making BSM angles limit plot."""
     fig_text = gen_figtext(args)
 
     raw = np.load(infile)
+    print 'raw', raw
+    print 'raw.shape', raw.shape
     sc_ranges = (
-        myround(np.min(raw[0][:,0]), down=True),
-        myround(np.max(raw[0][:,0]), up=True)
+        myround(np.min(raw[0][:,0]), up=True),
+        myround(np.max(raw[0][:,0]), down=True)
     )
 
     proc = []
-    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]))
+    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
 
     fig = plt.figure(figsize=(7, 5))
     ax = fig.add_subplot(111)
@@ -341,14 +361,14 @@ def plot_BSM_angles_limit(infile, outfile, xticks, outformat, args):
         #     ec='r', lw=2
         # )
         ax.annotate(
-            s='', xy=l, xytext=(l[0], l[1]-1.5),
+            s='', xy=l, xytext=(l[0], l[1]+1.5),
             arrowprops={'arrowstyle': '<-', 'lw': 1.5, 'color':'r'}
         )
 
     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)
 
     for of in outformat:
         fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150)