Mon 8 Apr 15:15:12 CDT 2019

8 files changed, 787 insertions, 7 deletions

diff --git a/fig2.py b/fig2.py
index 54a724d..e60c987 100755
--- a/fig2.py
+++ b/fig2.py
@@ -108,7 +108,7 @@ def main():
         map(fr_utils.angles_to_fr, flavour_angles)
     )
 
-    nbins = 15
+    nbins = 25
 
     fig = plt.figure(figsize=(8, 8))
     ax = fig.add_subplot(111)
diff --git a/plot_command b/plot_command
index 54ec50e..cf733f5 100644
--- a/plot_command
+++ b/plot_command
@@ -1 +1,3 @@
 python plot_sens.py --data real --dimensions 3 4 5 6 7 8 --fix-source-ratio True --source-ratios 1 2 0 1 0 0 0 1 0 --split-jobs True --stat-method bayesian --run-method fixed_angle --infile /data/user/smandalia/flavour_ratio/data --likelihood golemfit --sens-bins 20
+
+python contour.py --data realisation --debug True --injected-ratio 1 1 1 --likelihood golemfit --mn-live-points 500 --mn-tolerance 0.3 --outfile /data/user/smandalia/flavour_ratio/data/contour/golemfit/realisation/ --plot-chains True --plot-triangle True --run-scan False
diff --git a/plot_sens.py b/plot_sens.py
index ba3e923..f91b5f2 100755
--- a/plot_sens.py
+++ b/plot_sens.py
@@ -163,7 +163,7 @@ def parse_args(args=None):
         help='Set the new physics dimensions to consider'
     )
     parser.add_argument(
-        '--source-ratios', type=int, nargs='*', default=[2, 1, 0],
+        '--source-ratios', type=int, nargs='*', default=[1, 2, 0],
         help='Set the source flavour ratios for the case when you want to fix it'
     )
     if args is None: return parser.parse_args()
diff --git a/plot_sens_sourcescan.py b/plot_sens_sourcescan.py
new file mode 100755
index 0000000..130817d
--- /dev/null
+++ b/plot_sens_sourcescan.py
@@ -0,0 +1,411 @@
+#! /usr/bin/env python
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : April 28, 2018
+
+"""
+HESE BSM flavour ratio analysis plotting script
+"""
+
+from __future__ import absolute_import, division
+
+import os
+import argparse
+from functools import partial
+from copy import deepcopy
+
+import numpy as np
+import numpy.ma as ma
+
+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 plot as plot_utils
+from utils.enums import EnergyDependance, Likelihood, MixingScenario, ParamTag
+from utils.enums import PriorsCateg, SensitivityCateg, StatCateg
+from utils.param import Param, ParamSet, get_paramsets
+
+from utils import mn as mn_utils
+
+
+def define_nuisance():
+    """Define the nuisance parameters."""
+    tag = ParamTag.SM_ANGLES
+    g_prior = PriorsCateg.GAUSSIAN
+    lg_prior = PriorsCateg.LIMITEDGAUSS
+    e = 1e-9
+    nuisance = [
+        Param(name='s_12_2', value=0.307,            seed=[0.26, 0.35],     ranges=[0., 1.],      std=0.013,   tex=r's_{12}^2', prior=lg_prior,  tag=tag),
+        Param(name='c_13_4', value=(1-(0.02206))**2, seed=[0.950, 0.961],   ranges=[0., 1.],      std=0.00147, tex=r'c_{13}^4', prior=lg_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=lg_prior,  tag=tag),
+        Param(name='dcp',    value=4.08404,          seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0,     tex=r'\delta_{CP}', tag=tag),
+        Param(
+            name='m21_2', value=7.40E-23, seed=[7.2E-23, 7.6E-23], ranges=[6.80E-23, 8.02E-23],
+            std=2.1E-24, tex=r'\Delta m_{21}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag
+        ),
+        Param(
+            name='m3x_2', value=2.494E-21, seed=[2.46E-21, 2.53E-21], ranges=[2.399E-21, 2.593E-21],
+            std=3.3E-23, tex=r'\Delta m_{3x}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag
+        )
+    ]
+    tag = ParamTag.NUISANCE
+    nuisance.extend([
+        Param(name='convNorm',        value=1.,  seed=[0.5, 2. ], ranges=[0.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='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)
+    ])
+    return ParamSet(nuisance)
+
+
+def nuisance_argparse(parser):
+    nuisance = define_nuisance()
+    for parm in nuisance:
+        parser.add_argument(
+            '--'+parm.name, type=float, default=parm.value,
+            help=parm.name+' to inject'
+        )
+
+
+def process_args(args):
+    """Process the input args."""
+    if args.fix_mixing is not MixingScenario.NONE and args.fix_scale:
+        raise NotImplementedError('Fixed mixing and scale not implemented')
+    if args.fix_mixing is not MixingScenario.NONE and args.fix_mixing_almost:
+        raise NotImplementedError(
+            '--fix-mixing and --fix-mixing-almost cannot be used together'
+        )
+    if args.fix_scale:
+        raise NotImplementedError(
+            '--fix-scale not implemented'
+        )
+
+    args.measured_ratio = fr_utils.normalise_fr(args.measured_ratio)
+    # if args.fix_source_ratio:
+    #     assert len(args.source_ratios) % 3 == 0
+    #     srs = [args.source_ratios[3*x:3*x+3]
+    #            for x in range(int(len(args.source_ratios)/3))]
+    #     args.source_ratios = map(fr_utils.normalise_fr, srs)
+
+    if args.energy_dependance is EnergyDependance.SPECTRAL:
+        args.binning = np.logspace(
+            np.log10(args.binning[0]), np.log10(args.binning[1]), args.binning[2]+1
+        )
+
+    if args.split_jobs and args.run_method is SensitivityCateg.FULL:
+        raise NotImplementedError(
+            'split_jobs and run_method not implemented'
+        )
+
+    args.dimensions = np.sort(args.dimensions)
+
+    args_copy = deepcopy(args)
+    scale_regions = []
+    for dim in args.dimensions:
+        args_copy.dimension = dim
+        _, scale_region = fr_utils.estimate_scale(args_copy)
+        scale_regions.append(scale_region)
+    args.scale_region = [np.min(scale_regions), np.max(scale_regions)]
+    args.scale = np.power(10., np.average(np.log10(args.scale_region)))
+
+
+def parse_args(args=None):
+    """Parse command line arguments"""
+    parser = argparse.ArgumentParser(
+        description="HESE BSM flavour ratio analysis plotting script",
+        formatter_class=misc_utils.SortingHelpFormatter,
+    )
+    parser.add_argument(
+        '--infile', type=str, default='./untitled',
+        help='Path to input dir'
+    )
+    parser.add_argument(
+        '--run-method', default='full',
+        type=partial(misc_utils.enum_parse, c=SensitivityCateg),
+        choices=SensitivityCateg,
+        help='Choose which type of sensivity plot to make'
+    )
+    parser.add_argument(
+        '--stat-method', default='bayesian',
+        type=partial(misc_utils.enum_parse, c=StatCateg), choices=StatCateg,
+        help='Statistical method to employ'
+    )
+    parser.add_argument(
+        '--sens-bins', type=int, default=10,
+        help='Number of bins for the Bayes factor plot'
+    )
+    parser.add_argument(
+        '--split-jobs', type=misc_utils.parse_bool, default='False',
+        help='Did the jobs get split'
+    )
+    parser.add_argument(
+        '--plot', type=misc_utils.parse_bool, default='True',
+        help='Make sensitivity plots'
+    )
+    parser.add_argument(
+        '--plot-statistic', type=misc_utils.parse_bool, default='False',
+        help='Plot MultiNest evidence or LLH value'
+    )
+    fr_utils.fr_argparse(parser)
+    gf_utils.gf_argparse(parser)
+    llh_utils.likelihood_argparse(parser)
+    mn_utils.mn_argparse(parser)
+    nuisance_argparse(parser)
+    misc_utils.remove_option(parser, 'dimension')
+    misc_utils.remove_option(parser, 'source_ratio')
+    misc_utils.remove_option(parser, 'scale')
+    misc_utils.remove_option(parser, 'scale_region')
+    parser.add_argument(
+        '--dimensions', type=int, nargs='*', default=[3, 6],
+        help='Set the new physics dimensions to consider'
+    )
+    parser.add_argument(
+        '--source-bins', type=int, default=5,
+        help='Binning in source flavour space'
+    )
+    if args is None: return parser.parse_args()
+    else: return parser.parse_args(args.split())
+
+
+def main():
+    args = parse_args()
+    process_args(args)
+    args.scale = 0
+    misc_utils.print_args(args)
+
+    asimov_paramset, llh_paramset = get_paramsets(args, define_nuisance())
+
+    scale = llh_paramset.from_tag(ParamTag.SCALE)[0]
+    mmangles = llh_paramset.from_tag(ParamTag.MMANGLES)
+    if args.run_method is SensitivityCateg.FULL:
+        st_paramset_arr = [llh_paramset.from_tag(ParamTag.SCALE, invert=True)]
+    elif args.run_method in [SensitivityCateg.FIXED_ANGLE, SensitivityCateg.CORR_ANGLE]:
+        nscale_pmset = llh_paramset.from_tag([ParamTag.SCALE, ParamTag.MMANGLES], invert=True)
+        st_paramset_arr = [nscale_pmset] * 3
+    elif args.run_method in [SensitivityCateg.FIXED_ONE_ANGLE, SensitivityCateg.CORR_ONE_ANGLE]:
+        nscale_pmset = llh_paramset.from_tag(ParamTag.SCALE, invert=True)
+        st_paramset_arr = []
+        for x in xrange(3):
+            st_paramset_arr.append(
+                ParamSet([prms for prms in nscale_pmset
+                          if mmangles[x].name != prms.name])
+            )
+
+    corr_angles_categ = [SensitivityCateg.CORR_ANGLE, SensitivityCateg.CORR_ONE_ANGLE]
+    fixed_angle_categ = [SensitivityCateg.FIXED_ANGLE, SensitivityCateg.FIXED_ONE_ANGLE]
+
+    if args.run_method in corr_angles_categ:
+        scan_angles = np.linspace(0+1e-9, 1-1e-9, args.sens_bins)
+    else: scan_angles = np.array([0])
+    print 'scan_angles', scan_angles
+
+    dims = len(args.dimensions)
+    binning = np.linspace(0, 1, args.source_bins)
+    grid = np.dstack(np.meshgrid(binning, binning)).reshape(
+        args.source_bins*args.source_bins, 2
+    )
+    source_ratios = []
+    for x in grid:
+        if x[0]+x[1] > 1:
+            continue
+        source_ratios.append([x[0], x[1], 1-x[0]-x[1]])
+    args.source_ratios = source_ratios
+    n_source_ratios = map(fr_utils.normalise_fr, source_ratios)
+
+    srcs = len(n_source_ratios)
+    if args.run_method is SensitivityCateg.FULL:
+        statistic_arr = np.full((dims, srcs, args.sens_bins, 2), np.nan)
+    elif args.run_method in fixed_angle_categ:
+        statistic_arr = np.full((dims, srcs, len(st_paramset_arr), args.sens_bins, 2), np.nan)
+    elif args.run_method in corr_angles_categ:
+        statistic_arr = np.full(
+            (dims, srcs, len(st_paramset_arr), args.sens_bins, args.sens_bins, 3), np.nan
+        )
+
+    print 'Loading data'
+    for idim, dim in enumerate(args.dimensions):
+        argsc = deepcopy(args)
+        argsc.dimension = dim
+        _, scale_region = fr_utils.estimate_scale(argsc)
+        argsc.scale_region = scale_region
+        scan_scales = np.linspace(
+            np.log10(scale_region[0]), np.log10(scale_region[1]), args.sens_bins
+        )
+        scan_scales = np.concatenate([[-100.], scan_scales])
+
+        for isrc, src in enumerate(n_source_ratios):
+            argsc.source_ratio = src
+            infile = args.infile
+            if args.likelihood is Likelihood.GOLEMFIT:
+                infile += '/golemfit/'
+            elif args.likelihood is Likelihood.GAUSSIAN:
+                infile += '/gaussian/'
+            if args.likelihood is Likelihood.GAUSSIAN:
+                infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_'))
+            # infile += '/DIM{0}/fix_ifr/{1}/{2}/{3}/fr_stat'.format(
+            infile += '/DIM{0}/fix_ifr/sourcescan/{1}/{2}/{3}/fr_stat'.format(
+            # infile += '/DIM{0}/fix_ifr/prior/{1}/{2}/{3}/fr_stat'.format(
+            # infile += '/DIM{0}/fix_ifr/{1}/{2}/{3}/old/fr_stat'.format(
+            # infile += '/DIM{0}/fix_ifr/seed2/{1}/{2}/{3}/fr_stat'.format(
+            # infile += '/DIM{0}/fix_ifr/100TeV/{1}/{2}/{3}/fr_stat'.format(
+            # infile += '/DIM{0}/fix_ifr/strictprior/{1}/{2}/{3}/fr_stat'.format(
+            # infile += '/DIM{0}/fix_ifr/noprior/{1}/{2}/{3}/fr_stat'.format(
+                dim, *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data])
+            ) + misc_utils.gen_identifier(argsc)
+            print '== {0:<25} = {1}'.format('infile', infile)
+
+            if args.split_jobs:
+                for idx_an, an in enumerate(scan_angles):
+                    for idx_sc, sc in enumerate(scan_scales):
+                        filename = infile + '_scale_{0:.0E}'.format(np.power(10, sc))
+                        try:
+                            if args.run_method in fixed_angle_categ:
+                                print 'Loading from {0}'.format(filename+'.npy')
+                                statistic_arr[idim][isrc][:,idx_sc] = np.load(filename+'.npy')[:,0]
+                            if args.run_method in corr_angles_categ:
+                                filename += '_angle_{0:<04.2}'.format(an)
+                                print 'Loading from {0}'.format(filename+'.npy')
+                                statistic_arr[idim][isrc][:,idx_an,idx_sc] = np.load(filename+'.npy')[:,0,0]
+                        except:
+                            print 'Unable to load file {0}'.format(filename+'.npy')
+                            continue
+            else:
+                print 'Loading from {0}'.format(infile+'.npy')
+                try:
+                    statistic_arr[idim][isrc] = np.load(infile+'.npy')
+                except:
+                    print 'Unable to load file {0}'.format(infile+'.npy')
+                    continue
+
+    data = ma.masked_invalid(statistic_arr)
+
+    print 'data', data
+    print 'data.shape', data.shape
+    if args.plot_statistic:
+        assert 0
+        print 'Plotting statistic'
+
+        argsc = deepcopy(args)
+        for idim, dim in enumerate(args.dimensions):
+            argsc.dimension = dim
+            _, scale_region = fr_utils.estimate_scale(argsc)
+            argsc.scale_region = scale_region
+            base_infile = args.infile
+            if args.likelihood is Likelihood.GOLEMFIT:
+                base_infile += '/golemfit/'
+            elif args.likelihood is Likelihood.GAUSSIAN:
+                base_infile += '/gaussian/'
+            if args.likelihood is Likelihood.GAUSSIAN:
+                base_infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_'))
+            # base_infile += '/DIM{0}/fix_ifr'.format(dim)
+            base_infile += '/DIM{0}/fix_ifr/sourcescan'.format(dim)
+            # base_infile += '/DIM{0}/fix_ifr/prior'.format(dim)
+            # base_infile += '/DIM{0}/fix_ifr/seed2'.format(dim)
+            # base_infile += '/DIM{0}/fix_ifr/100TeV'.format(dim)
+            # base_infile += '/DIM{0}/fix_ifr/strictprior'.format(dim)
+            # base_infile += '/DIM{0}/fix_ifr/noprior'.format(dim)
+
+            for isrc, src in enumerate(n_source_ratios):
+                argsc.source_ratio = src
+                infile = base_infile +'/{0}/{1}/{2}/fr_stat'.format(
+                # infile = base_infile +'/{0}/{1}/{2}/old/fr_stat'.format(
+                    *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data])
+                ) + misc_utils.gen_identifier(argsc)
+                basename = os.path.dirname(infile)
+                baseoutfile = basename[:5]+basename[5:].replace('data', 'plots')
+                baseoutfile += '/' + os.path.basename(infile)
+                if args.run_method is SensitivityCateg.FULL:
+                    outfile = baseoutfile
+                    plot_utils.plot_statistic(
+                        data        = data[idim][isrc],
+                        outfile     = outfile,
+                        outformat   = ['png'],
+                        args        = argsc,
+                        scale_param = scale,
+                    )
+                if args.run_method in fixed_angle_categ:
+                    for idx_scen in xrange(len(st_paramset_arr)):
+                        print '|||| SCENARIO = {0}'.format(idx_scen)
+                        outfile = baseoutfile + '_SCEN{0}'.format(idx_scen)
+                        if idx_scen == 0: label = r'$\mathcal{O}_{12}=\pi/4$'
+                        elif idx_scen == 1: label = r'$\mathcal{O}_{13}=\pi/4$'
+                        elif idx_scen == 2: label = r'$\mathcal{O}_{23}=\pi/4$'
+                        plot_utils.plot_statistic(
+                            data        = data[idim][isrc][idx_scen],
+                            outfile     = outfile,
+                            outformat   = ['png'],
+                            args        = argsc,
+                            scale_param = scale,
+                            label       = label
+                        )
+                elif args.run_method in corr_angles_categ:
+                    for idx_scen in xrange(len(st_paramset_arr)):
+                        print '|||| SCENARIO = {0}'.format(idx_scen)
+                        basescenoutfile = baseoutfile + '_SCEN{0}'.format(idx_scen)
+                        if idx_scen == 0: label = r'$\mathcal{O}_{12}='
+                        elif idx_scen == 1: label = r'$\mathcal{O}_{13}='
+                        elif idx_scen == 2: label = r'$\mathcal{O}_{23}='
+                        for idx_an, an in enumerate(scan_angles):
+                            print '|||| ANGLE = {0:<04.2}'.format(float(an))
+                            outfile = basescenoutfile + '_ANGLE{0}'.format(idx_an)
+                            _label = label + r'{0:<04.2}$'.format(an)
+                            plot_utils.plot_statistic(
+                                data        = data[idim][isrc][idx_scen][idx_an][:,1:],
+                                outfile     = outfile,
+                                outformat   = ['png'],
+                                args        = argsc,
+                                scale_param = scale,
+                                label       = _label
+                            )
+
+    if args.plot:
+        print 'Plotting sensitivities'
+
+        basename = args.infile[:5]+args.infile[5:].replace('data', 'plots')
+        baseoutfile = basename + '/{0}/{1}/{2}/'.format(
+            *map(misc_utils.parse_enum, [args.likelihood, args.stat_method, args.data])
+        )
+
+        if args.run_method is SensitivityCateg.FULL:
+            plot_utils.plot_sens_full(
+                data      = data,
+                outfile   = baseoutfile + '/FULL',
+                outformat = ['png', 'pdf'],
+                args      = args,
+            )
+        elif args.run_method in fixed_angle_categ:
+            # plot_utils.plot_sens_fixed_angle_pretty(
+            #     data      = data,
+            #     outfile   = baseoutfile + '/fixed_angle_pretty_substantial',
+            #     outformat = ['png', 'pdf'],
+            #     args      = args,
+            # )
+            # plot_utils.plot_sens_fixed_angle(
+            #     data      = data,
+            #     outfile   = baseoutfile + '/FIXED_ANGLE',
+            #     outformat = ['png'],
+            #     args      = args,
+            # )
+            plot_utils.plot_source_ternary_1D(
+                data      = data,
+                outfile   = baseoutfile + '/source_ternary',
+                outformat = ['png'],
+                args      = args,
+            )
+        elif args.run_method in corr_angles_categ:
+            plot_utils.plot_sens_corr_angle(
+                data      = data,
+                outfile   = baseoutfile + '/CORR_ANGLE',
+                outformat = ['png', 'pdf'],
+                args      = args,
+            )
+
+
+main.__doc__ = __doc__
+
+
+if __name__ == '__main__':
+    main()
diff --git a/submitter/out b/submitter/out
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/submitter/out
diff --git a/submitter/sens_dag_source.py b/submitter/sens_dag_source.py
new file mode 100644
index 0000000..bdb5924
--- /dev/null
+++ b/submitter/sens_dag_source.py
@@ -0,0 +1,183 @@
+#! /usr/bin/env python
+
+import os
+import numpy as np
+
+full_scan_mfr = [
+    # (1, 1, 1), (1, 2, 0)
+]
+
+bins = 5
+binning = np.linspace(0, 1, bins)
+grid = np.dstack(np.meshgrid(binning, binning)).reshape(bins*bins, 2)
+sources = []
+for x in grid:
+    if x[0]+x[1] > 1:
+        continue
+    sources.append([x[0], x[1], 1-x[0]-x[1]])
+
+# 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),
+# ]
+fix_sfr_mfr = []
+for s in sources:
+    fix_sfr_mfr.append((1, 1, 1, s[0], s[1], s[2]))
+print 'fix_sfr_mfr', fix_sfr_mfr
+print 'len(fix_sfr_mfr)', len(fix_sfr_mfr)
+
+GLOBAL_PARAMS = {}
+
+# Bayes Factor
+sens_eval_bin = 'true' # set to 'all' to run normally
+GLOBAL_PARAMS.update(dict(
+    sens_run      = 'True',
+    run_method    = 'fixed_angle', # full, fixed_angle, corr_angle
+    stat_method   = 'bayesian',
+    sens_bins     = 10,
+    seed          = None
+))
+
+# MultiNest
+GLOBAL_PARAMS.update(dict(
+    # mn_live_points = 1000,
+    # mn_live_points = 600,
+    mn_live_points = 100,
+    # mn_tolerance   = 0.1,
+    mn_tolerance   = 0.3,
+    mn_output      = './mnrun'
+))
+
+# FR
+dimension         = [6]
+# dimension         = [3, 6]
+# dimension         = [3, 4, 5, 6, 7, 8]
+GLOBAL_PARAMS.update(dict(
+    threads            = 1,
+    binning            = '6e4 1e7 20',
+    no_bsm             = 'False',
+    scale_region       = "1E10",
+    energy_dependance  = 'spectral',
+    spectral_index     = -2,
+    fix_mixing         = 'None',
+    fix_mixing_almost  = 'False',
+    fold_index         = 'True',
+    save_measured_fr   = 'False',
+    output_measured_fr = './frs/'
+))
+
+# Likelihood
+GLOBAL_PARAMS.update(dict(
+    likelihood  = 'golemfit',
+    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_{0}_{1}_{2}_{3}'.format(
+    GLOBAL_PARAMS['stat_method'], GLOBAL_PARAMS['run_method'],
+    GLOBAL_PARAMS['likelihood'], GLOBAL_PARAMS['data']
+)
+# outfile += '_seed2'
+# outfile += '_tol03'
+# outfile += '_NULL'
+# outfile += '_prior'
+# outfile += '_strictprior'
+# outfile += '_noprior'
+outfile += '_sourcescan'
+outfile += '.submit'
+golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
+condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub'
+
+if sens_eval_bin.lower() != 'all':
+    if GLOBAL_PARAMS['run_method'].lower() == 'corr_angle':
+        raise NotImplementedError
+        sens_runs = GLOBAL_PARAMS['sens_bins']**2
+    else:
+        sens_runs = GLOBAL_PARAMS['sens_bins'] + 1
+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}'.format(
+            GLOBAL_PARAMS['likelihood'], dim
+        )
+        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 += 'seed2/'
+            # output += 'mn_noverlap/'
+            # output += 'tol_03/'
+            # output += 'prior/'
+            # output += 'strictprior/'
+            # output += 'noprior/'
+            output += 'sourcescan/'
+            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]))
+                if sens_eval_bin.lower() != 'all':
+                    f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, r))
+                else:
+                    f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, 'all'))
+                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
+                # break
+
+        # 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('.', '_'))
+        #     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))
+        #         if sens_eval_bin.lower() != 'all':
+        #             f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, r))
+        #         else:
+        #             f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, 'all'))
+        #         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
+
+    print 'dag file = {0}'.format(outfile)
diff --git a/utils/enums.py b/utils/enums.py
index a7982f8..441a16f 100644
--- a/utils/enums.py
+++ b/utils/enums.py
@@ -76,3 +76,9 @@ class StatCateg(Enum):
 class SteeringCateg(Enum):
     P2_0 = 1
     P2_1 = 2
+
+
+class Texture(Enum):
+    OEU = 0
+    OET = 1
+    OUT = 2
diff --git a/utils/plot.py b/utils/plot.py
index f81f9da..5a3d1f7 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -32,12 +32,15 @@ from getdist import plots, mcsamples
 
 import ternary
 from ternary.heatmapping import polygon_generator
+# print ternary.__file__
+# assert 0
 
 import shapely.geometry as geometry
 
 from utils import misc as misc_utils
-from utils.enums import DataType, EnergyDependance
+from utils.enums import DataType, EnergyDependance, str_enum
 from utils.enums import Likelihood, MixingScenario, ParamTag, StatCateg
+from utils.enums import Texture
 from utils.fr import angles_to_u, angles_to_fr
 
 
@@ -595,7 +598,7 @@ 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**(bayes_K))] # Strong degree of belief
+                    al = scales_rm[reduced_ev > np.log(10**(bayes_K))]
                 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
@@ -997,9 +1000,10 @@ def get_tax(ax, scale):
     tax.clear_matplotlib_ticks()
 
     # Set ticks
-    ticks = np.linspace(0, 1, 6)
-    tax.ticks(ticks=ticks, locations=ticks*scale, axis='blr', linewidth=1,
-              offset=0.03, fontsize=fontsize, tick_formats='%.1f')
+    # ticks = np.linspace(0, 1, 6)
+    # tax.ticks(ticks=ticks, locations=ticks*scale, axis='blr', linewidth=1,
+    #           offset=0.03, fontsize=fontsize, tick_formats='%.1f')
+    tax.ticks()
 
     tax._redraw_labels()
 
@@ -1156,3 +1160,177 @@ def flavour_contour(frs, ax, nbins, coverage, linewidth=2, color='black'):
     )
     ax.scatter(points.T[0], points.T[1], marker='o', s=2, alpha=1, color=color,
                zorder=3)
+
+def plot_source_ternary(data, outfile, outformat, args):
+    """Ternary plot in the source flavour space for each operator texture."""
+    r_data = np.full((data.shape[0], data.shape[2], data.shape[1], data.shape[3], data.shape[4]), np.nan)
+    for idim in xrange(data.shape[0]):
+        for isrc in xrange(data.shape[1]):
+            for isce in xrange(data.shape[2]):
+                r_data[idim][isce][isrc] = data[idim][isrc][isce]
+    r_data = ma.masked_invalid(r_data)
+    print r_data.shape, 'r_data.shape'
+    nsrcs = int(len(args.source_ratios) / 3)
+
+    for idim, dim in enumerate(args.dimensions):
+        print '|||| DIM = {0}, {1}'.format(idim, dim)
+        for isce in xrange(r_data.shape[1]):
+            print '|||| SCEN = {0}'.format(str_enum(Texture(isce)))
+            fig = plt.figure(figsize=(8, 8))
+            ax  = fig.add_subplot(111)
+            tax = get_tax(ax, scale=nsrcs)
+            interp_dict = {}
+            for isrc, src in enumerate(args.source_ratios):
+                src_sc = tuple(np.array(src)*(nsrcs-1))
+                print '|||| SRC = {0}'.format(src)
+                scales, statistic = ma.compress_rows(r_data[idim][isce][isrc]).T
+                print 'scales', scales
+                print 'statistic', statistic
+
+                try:
+                    tck, u = splprep([scales, statistic], s=0)
+                except:
+                    interp_dict[src_sc] = -60
+                    continue
+                # sc, st = splev(np.linspace(0, 1, 10000), tck)
+                sc, st = splev(np.linspace(0, 1, 20), tck)
+                scales_rm = sc[sc >= scales[1]]
+                statistic_rm = st[sc >= scales[1]]
+                # scales_rm = sc
+                # statistic_rm = st
+
+                min_idx = np.argmin(scales)
+                null = statistic[min_idx]
+                if args.stat_method is StatCateg.BAYESIAN:
+                    reduced_ev = -(statistic_rm - null)
+                    print 'reduced_ev', reduced_ev
+                    al = scales_rm[reduced_ev > np.log(10**(bayes_K))]
+                else:
+                    assert 0
+                if len(al) == 0:
+                    print 'No points for DIM {0} FRS {1}!'.format(dim, src)
+                    interp_dict[src_sc] = -60
+                    continue
+                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)
+                    interp_dict[src_sc] = -60
+                    continue
+                lim = al[0]
+                print 'limit = {0}'.format(lim)
+
+                interp_dict[src_sc] = lim
+            print 'interp_dict', interp_dict
+            print
+            print 'vertices', heatmap(interp_dict, nsrcs)
+            print
+            tax.heatmap(interp_dict, scale=nsrcs, vmin=-60, vmax=-30)
+            misc_utils.make_dir(outfile)
+            for of in outformat:
+                print 'Saving plot as {0}'.format(outfile+'_SCEN{0}.'.format(isce)+of)
+                fig.savefig(outfile+'_SCEN{0}.'.format(isce)+of, bbox_inches='tight', dpi=150)
+            print 'nsrcs', nsrcs
+            assert 0
+
+
+def texture_label(x):
+    if x == Texture.OEU:
+        return r'$\mathcal{O}_{e\mu}$'
+    elif x == Texture.OET:
+        return r'$\mathcal{O}_{e\tau}$'
+    elif x == Texture.OUT:
+        return r'$\mathcal{O}_{\mu\tau}$'
+    else:
+        raise AssertionError
+
+
+def plot_source_ternary_1D(data, outfile, outformat, args):
+    """Ternary plot in the source flavour space for each operator texture."""
+    sources = args.source_ratios
+    x_1D = []
+    i_src_1D = []
+    for i, s in enumerate(sources):
+        if s[2] != 0: continue
+        x_1D.append(s[0])
+        i_src_1D.append(i)
+    i_src_1D = list(reversed(i_src_1D))
+    x_1D = list(reversed(x_1D))
+    print 'x_1D', x_1D
+    def get_edges_from_cen(bincen):
+        hwidth = 0.5*(bincen[1] - bincen[0])
+        return np.append([bincen[0]-hwidth], bincen[:]+hwidth)
+    be = get_edges_from_cen(x_1D)
+
+    r_data = np.full((data.shape[0], data.shape[2], len(i_src_1D), data.shape[3], data.shape[4]), np.nan)
+    for idim in xrange(data.shape[0]):
+        for i1D, isrc in enumerate(i_src_1D):
+            for isce in xrange(data.shape[2]):
+                r_data[idim][isce][i1D] = data[idim][isrc][isce]
+    r_data = ma.masked_invalid(r_data)
+    print r_data.shape, 'r_data.shape'
+
+    for idim, dim in enumerate(args.dimensions):
+        print '|||| DIM = {0}, {1}'.format(idim, dim)
+        fig = plt.figure(figsize=(4, 4))
+        ax  = fig.add_subplot(111)
+        ax.tick_params(axis='x', labelsize=12)
+        ax.tick_params(axis='y', labelsize=12)
+        ax.set_xlabel(r'$x$', fontsize=18)
+        ax.set_ylabel(r'${\rm New\:Physics\:Scale}\:[\:{\rm log}_{10} (\Lambda_{d=6}^{-1}\:/\:{\rm GeV}^{-2})\: ]$', fontsize=18)
+        ax.set_xlim(0, 1)
+        for isce in xrange(r_data.shape[1]):
+            print '|||| SCEN = {0}'.format(str_enum(Texture(isce)))
+            H = np.full(len(x_1D), np.nan)
+            for ix, x in enumerate(x_1D):
+                print '|||| X = {0}'.format(x)
+                scales, statistic = ma.compress_rows(r_data[idim][isce][ix]).T
+                print 'scales', scales
+                print 'statistic', statistic
+
+                try:
+                    tck, u = splprep([scales, statistic], s=0)
+                except:
+                    continue
+                # sc, st = splev(np.linspace(0, 1, 10000), tck)
+                sc, st = splev(np.linspace(0, 1, 20), tck)
+                scales_rm = sc[sc >= scales[1]]
+                statistic_rm = st[sc >= scales[1]]
+                # scales_rm = sc
+                # statistic_rm = st
+
+                min_idx = np.argmin(scales)
+                null = statistic[min_idx]
+                if args.stat_method is StatCateg.BAYESIAN:
+                    reduced_ev = -(statistic_rm - null)
+                    print 'reduced_ev', reduced_ev
+                    al = scales_rm[reduced_ev > np.log(10**(bayes_K))]
+                else:
+                    assert 0
+                if len(al) == 0:
+                    print 'No points for DIM {0} X {1}!'.format(dim, x)
+                    continue
+                if reduced_ev[-1] < np.log(10**(bayes_K)) - 0.1:
+                    print 'Peaked contour does not exclude large scales! For ' \
+                        'DIM {0} X {1}!'.format(dim, x)
+                    continue
+                lim = al[0]
+                print 'limit = {0}'.format(lim)
+
+                H[ix] = lim
+            H = ma.masked_invalid(H)
+            H_0 = np.concatenate([[H[0]], H])
+            ax.step(be, H_0, linewidth=2,
+                    label=texture_label(Texture(isce)), linestyle='-',
+                    drawstyle='steps-pre')
+            print 'H', H
+            print
+        for ymaj in ax.yaxis.get_majorticklocs():
+            ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.3, linewidth=1)
+        for xmaj in be:
+            ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.3, linewidth=1)
+        ax.legend()
+        misc_utils.make_dir(outfile)
+        for of in outformat:
+            print 'Saving plot as {0}'.format(outfile+'_DIM{0}.'.format(dim)+of)
+            fig.savefig(outfile+'_DIM{0}.'.format(dim)+of, bbox_inches='tight', dpi=150)
+