refactor sensitivity/limit scripts

8 files changed, 203 insertions, 456 deletions

diff --git a/utils/enums.py b/utils/enums.py
index 441a16f..e85158d 100644
--- a/utils/enums.py
+++ b/utils/enums.py
@@ -20,23 +20,9 @@ class DataType(Enum):
     REALISATION = 3
 
 
-class EnergyDependance(Enum):
-    MONO     = 1
-    SPECTRAL = 2
-
-
 class Likelihood(Enum):
-    FLAT     = 1
-    GAUSSIAN = 2
-    GOLEMFIT = 3
-    GF_FREQ  = 4
-
-
-class MixingScenario(Enum):
-    T12  = 1
-    T13  = 2
-    T23  = 3
-    NONE = 4
+    GOLEMFIT = 1
+    GF_FREQ  = 2
 
 
 class ParamTag(Enum):
@@ -60,14 +46,6 @@ class MCMCSeedType(Enum):
     GAUSSIAN = 2
 
 
-class SensitivityCateg(Enum):
-    FULL            = 1
-    FIXED_ANGLE     = 2
-    CORR_ANGLE      = 3
-    FIXED_ONE_ANGLE = 4
-    CORR_ONE_ANGLE  = 5
-
-
 class StatCateg(Enum):
     BAYESIAN    = 1
     FREQUENTIST = 2
@@ -79,6 +57,7 @@ class SteeringCateg(Enum):
 
 
 class Texture(Enum):
-    OEU = 0
-    OET = 1
-    OUT = 2
+    OEU  = 1
+    OET  = 2
+    OUT  = 3
+    NONE = 4
diff --git a/utils/fr.py b/utils/fr.py
index 1d12fe5..b8eba44 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -13,7 +13,7 @@ from functools import partial
 
 import numpy as np
 
-from utils.enums import EnergyDependance, MixingScenario
+from utils.enums import Texture
 from utils.misc import enum_parse, parse_bool
 
 import mpmath as mp
@@ -40,7 +40,17 @@ ASIN   = mp.asin
 EXP    = mp.exp
 
 MASS_EIGENVALUES = [7.40E-23, 2.515E-21]
-"""SM mass eigenvalues"""
+"""SM mass eigenvalues."""
+
+SCALE_BOUNDARIES = {
+    3 : (-32, -20),
+    4 : (-40, -24),
+    5 : (-48, -27),
+    6 : (-56, -30),
+    7 : (-64, -33),
+    8 : (-72, -36)
+}
+"""Boundaries to scan the NP scale for each dimension."""
 
 
 def determinant(x):
@@ -244,38 +254,13 @@ def normalise_fr(fr):
     return np.array(fr) / float(np.sum(fr))
 
 
-def estimate_scale(args):
-    """Estimate the scale at which new physics will enter."""
-    try: m_eign = args.m3x_2
-    except: m_eign = MASS_EIGENVALUES[1]
-    if hasattr(args, 'scale'):
-        if args.scale != 0:
-            scale = args.scale
-            scale_region = (scale/args.scale_region, scale*args.scale_region)
-            return scale, scale_region
-    if args.energy_dependance is EnergyDependance.MONO:
-        scale = np.power(
-            10, np.round(np.log10(m_eign/args.energy)) - \
-            np.log10(args.energy**(args.dimension-3))
-        )
-        scale_region = (scale/args.scale_region, scale*args.scale_region)
-    elif args.energy_dependance is EnergyDependance.SPECTRAL:
-        lower_s = (m_eign/args.binning[-1]) / (args.binning[-1]**(args.dimension-3))
-        upper_s = (m_eign/args.binning[0]) / (args.binning[0]**(args.dimension-3))
-        scale = np.power(10, np.average(np.log10([lower_s, upper_s])))
-        diff = upper_s / lower_s
-        scale_region = (lower_s/np.power(10, args.dimension), upper_s*diff*np.power(10, args.dimension))
-        scale_region = [np.power(10, np.round(np.log10(x))) for x in scale_region]
-    return scale, scale_region
-
-
 def fr_argparse(parser):
     parser.add_argument(
-        '--measured-ratio', type=float, nargs=3, default=[1, 1, 1],
-        help='Set the central value for the measured flavour ratio at IceCube'
+        '--injected-ratio', type=float, nargs=3, required=False,
+        help='Injected ratio if not using data'
     )
     parser.add_argument(
-        '--source-ratio', type=float, nargs=3, default=[2, 1, 0],
+        '--source-ratio', type=float, nargs=3, default=[1, 2, 0],
         help='Set the source flavour ratio for the case when you want to fix it'
     )
     parser.add_argument(
@@ -287,51 +272,13 @@ def fr_argparse(parser):
         help='Set the new physics dimension to consider'
     )
     parser.add_argument(
-        '--energy-dependance', default='spectral', type=partial(enum_parse, c=EnergyDependance),
-        choices=EnergyDependance,
-        help='Type of energy dependance to use in the BSM fit'
-    )
-    parser.add_argument(
-        '--spectral-index', default=-2, type=float,
-        help='Spectral index for spectral energy dependance'
-    )
-    parser.add_argument(
-        '--fold-index', default='True', type=parse_bool,
-        help='Fold in the spectral index when using GolemFit'
+        '--texture', type=partial(enum_parse, c=Texture),
+        default='none', choices=Texture, help='Set the BSM mixing texture'
     )
     parser.add_argument(
         '--binning', default=[6e4, 1e7, 20], type=float, nargs=3,
         help='Binning for spectral energy dependance'
     )
-    parser.add_argument(
-        '--fix-source-ratio', type=parse_bool, default='False',
-        help='Fix the source flavour ratio'
-    )
-    parser.add_argument(
-        '--fix-mixing', type=partial(enum_parse, c=MixingScenario),
-        default='None', choices=MixingScenario,
-        help='Fix all mixing parameters to choice of maximal mixing'
-    )
-    parser.add_argument(
-        '--fix-mixing-almost', type=parse_bool, default='False',
-        help='Fix all mixing parameters except s23'
-    )
-    parser.add_argument(
-        '--fix-scale', type=parse_bool, default='False',
-        help='Fix the new physics scale'
-    )
-    parser.add_argument(
-        '--scale', type=float, default=0,
-        help='Set the new physics scale'
-    )
-    parser.add_argument(
-        '--scale-region', type=float, default=1e10,
-        help='Set the size of the box to scan for the scale'
-    )
-    parser.add_argument(
-        '--energy', type=float, default=1000,
-        help='Set the energy scale'
-    )
 
 
 def fr_to_angles(ratios):
@@ -363,8 +310,7 @@ NUFIT_U = angles_to_u((0.307, (1-0.02195)**2, 0.565, 3.97935))
 
 
 def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
-                   sm_u=NUFIT_U, no_bsm=False, fix_mixing=MixingScenario.NONE,
-                   fix_mixing_almost=False, fix_scale=False, scale=None,
+                   sm_u=NUFIT_U, no_bsm=False, texture=Texture.NONE,
                    check_uni=True, epsilon=1e-7):
     """Construct the BSM mixing matrix from the BSM parameters.
 
@@ -388,17 +334,8 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
     no_bsm : bool
         Turn off BSM behaviour
 
-    fix_mixing : MixingScenario
-        Fix the BSM mixing angles
-
-    fix_mixing_almost : bool
-        Fix the BSM mixing angles except one
-
-    fix_scale : bool
-        Fix the BSM scale
-
-    scale : float
-        Used with fix_scale - scale at which to fix
+    texture : Texture
+        BSM mixing texture
 
     check_uni : bool
         Check the resulting BSM mixing matrix is unitary
@@ -422,33 +359,20 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
             'got\n{0}'.format(sm_u)
         )
 
-    if fix_mixing is not MixingScenario.NONE and fix_mixing_almost:
-        raise NotImplementedError(
-            '--fix-mixing and --fix-mixing-almost cannot be used together'
-        )
-
     if not isinstance(theta, (list, tuple)):
         theta = [theta]
 
-    if fix_mixing is MixingScenario.T12:
-        s12_2, c13_4, s23_2, dcp, sc2 = 0.5, 1.0, 0.0, 0., theta
-    elif fix_mixing is MixingScenario.T13:
-        s12_2, c13_4, s23_2, dcp, sc2 = 0.0, 0.25, 0.0, 0., theta
-    elif fix_mixing is MixingScenario.T23:
-        s12_2, c13_4, s23_2, dcp, sc2 = 0.0, 1.0, 0.5, 0., theta
-    elif fix_mixing_almost:
-        s12_2, c13_4, dcp = 0.5, 1.0-1E-6, 0.
-        s23_2, sc2 = theta
-    elif fix_scale:
-        s12_2, c13_4, s23_2, dcp = theta
-        sc2 = scale
+    z = 0.+1e-9
+    if texture is Texture.OEU:
+        np_s12_2, np_c13_4, np_s23_2, np_dcp, sc2 = 0.5, 1.0, z, z, theta
+    elif texture is Texture.OET:
+        np_s12_2, np_c13_4, np_s23_2, np_dcp, sc2 = z, 0.25,  z, z, theta
+    elif texture is Texture.OUT:
+        np_s12_2, np_c13_4, np_s23_2, np_dcp, sc2 = z, 1.0, 0.5, z, theta
     else:
-        s12_2, c13_4, s23_2, dcp, sc2 = theta
+        np_s12_2, np_c13_4, np_s23_2, np_dcp, sc2 = theta
 
-    if len(theta) != 0:
-        sc2 = np.power(10., sc2)
-    else:
-        sc2 = scale
+    sc2 = np.power(10., sc2)
     sc1 = sc2 / 100.
 
     mass_matrix = np.array(
@@ -458,11 +382,11 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
     if no_bsm:
         eg_vector = cardano_eqn(sm_ham)
     else:
-        new_physics_u = angles_to_u((s12_2, c13_4, s23_2, dcp))
-        scale_matrix = np.array(
+        NP_U = angles_to_u((np_s12_2, np_c13_4, np_s23_2, np_dcp))
+        SC_U = np.array(
             [[0, 0, 0], [0, sc1, 0], [0, 0, sc2]]
         )
-        bsm_term = (energy**(dim-3)) * np.dot(new_physics_u, np.dot(scale_matrix, new_physics_u.conj().T))
+        bsm_term = (energy**(dim-3)) * np.dot(NP_U, np.dot(SC_U, NP_U.conj().T))
         bsm_ham = sm_ham + bsm_term
         eg_vector = cardano_eqn(bsm_ham)
 
diff --git a/utils/gf.py b/utils/gf.py
index 13d5728..bc004bd 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -24,6 +24,9 @@ from utils.misc import enum_parse, parse_bool, thread_factors
 from utils.param import ParamSet
 
 
+FITTER = None
+
+
 def fit_flags(llh_paramset):
     default_flags = {
         # False means it's not fixed in minimization
@@ -74,9 +77,6 @@ def steering_params(args):
     elif args.likelihood is Likelihood.GF_FREQ:
         params.frequentist = True;
 
-    # For Tianlu
-    # params.years = [999]
-
     if hasattr(args, 'binning'):
         params.minFitEnergy = args.binning[0]  # GeV
         params.maxFitEnergy = args.binning[-1] # GeV
@@ -95,59 +95,61 @@ def steering_params(args):
     return params
 
 
-def setup_asimov(fitter, params):
+def setup_asimov(params):
     print 'Injecting the model', params
     asimov_params = gf.FitParameters(gf.sampleTag.MagicTau)
     for parm in params:
         asimov_params.__setattr__(parm.name, float(parm.value))
-    fitter.SetupAsimov(asimov_params)
+    FITTER.SetupAsimov(asimov_params)
 
 
-def setup_realisation(fitter, params, seed):
+def setup_realisation(params, seed):
     print 'Injecting the model', params
     asimov_params = gf.FitParameters(gf.sampleTag.MagicTau)
     for parm in params:
         asimov_params.__setattr__(parm.name, float(parm.value))
-    fitter.Swallow(fitter.SpitRealization(asimov_params, seed))
+    FITTER.Swallow(FITTER.SpitRealization(asimov_params, seed))
 
 
 def setup_fitter(args, asimov_paramset):
+    global FITTER
     datapaths = gf.DataPaths()
     sparams = steering_params(args)
     npp = gf.NewPhysicsParams()
-    fitter = gf.GolemFit(datapaths, sparams, npp)
+    FITTER = gf.GolemFit(datapaths, sparams, npp)
     if args.data is DataType.ASIMOV:
-        setup_asimov(fitter, asimov_paramset)
+        setup_asimov(FITTER, asimov_paramset)
     elif args.data is DataType.REALISATION:
         seed = args.seed if args.seed is not None else 1
-        setup_realisation(fitter, asimov_paramset, seed)
+        setup_realisation(FITTER, asimov_paramset, seed)
     elif args.data is DataType.REAL:
         print 'Using MagicTau DATA'
-    return fitter
 
 
-def get_llh(fitter, params):
+def get_llh(params):
     # print 'params', params
     fitparams = gf.FitParameters(gf.sampleTag.MagicTau)
     for parm in params:
         fitparams.__setattr__(parm.name, float(parm.value))
-    llh = -fitter.EvalLLH(fitparams)
+    llh = -FITTER.EvalLLH(fitparams)
     return llh
 
 
-def get_llh_freq(fitter, params):
+def get_llh_freq(params):
     print 'setting to {0}'.format(params)
     fitparams = gf.FitParameters(gf.sampleTag.MagicTau)
     for parm in params:
         fitparams.__setattr__(parm.name, float(parm.value))
-    fitter.SetFitParametersSeed(fitparams)
-    llh = -fitter.MinLLH().likelihood
+    FITTER.SetFitParametersSeed(fitparams)
+    llh = -FITTER.MinLLH().likelihood
     return llh
 
 
-def data_distributions(fitter):
-    hdat = fitter.GetDataDistribution()
-    binedges = np.asarray([fitter.GetZenithBinsData(), fitter.GetEnergyBinsData()])
+def data_distributions():
+    hdat = FITTER.GetDataDistribution()
+    binedges = np.asarray(
+        [FITTER.GetZenithBinsData(), FITTER.GetEnergyBinsData()]
+    )
     return hdat, binedges
 
 
@@ -164,4 +166,3 @@ def gf_argparse(parser):
         choices=SteeringCateg,
         help='use asimov/fake dataset with specific steering'
     )
-
diff --git a/utils/llh.py b/utils/llh.py
index 0c1e97d..93587b9 100644
--- a/utils/llh.py
+++ b/utils/llh.py
@@ -17,7 +17,7 @@ from scipy.stats import multivariate_normal, truncnorm
 
 from utils import fr as fr_utils
 from utils import gf as gf_utils
-from utils.enums import EnergyDependance, Likelihood, ParamTag, PriorsCateg
+from utils.enums import Likelihood, ParamTag, PriorsCateg
 from utils.misc import enum_parse, gen_identifier, parse_bool
 
 
@@ -34,22 +34,11 @@ def multi_gaussian(fr, fr_bf, sigma, offset=-320):
     return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr)) + offset
 
 
-def likelihood_argparse(parser):
+def llh_argparse(parser):
     parser.add_argument(
-        '--likelihood', default='gaussian', type=partial(enum_parse, c=Likelihood),
-        choices=Likelihood, help='likelihood contour'
-    )
-    parser.add_argument(
-        '--sigma-ratio', type=float, default=0.01,
-        help='Set the 1 sigma for the measured flavour ratio for a gaussian LLH'
-    )
-    parser.add_argument(
-        '--save-measured-fr', type=parse_bool, default='False',
-        help='Output the measured flavour ratios'
-    )
-    parser.add_argument(
-        '--output-measured-fr', type=str, default='./frs',
-        help='Output of the measured flavour ratios'
+        '--stat-method', default='bayesian',
+        type=partial(misc_utils.enum_parse, c=StatCateg), choices=StatCateg,
+        help='Statistical method to employ'
     )
 
 
@@ -82,7 +71,7 @@ def lnprior(theta, paramset):
     return prior
 
 
-def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
+def triangle_llh(theta, args, asimov_paramset, llh_paramset):
     """Log likelihood function for a given theta."""
     if len(theta) != len(llh_paramset):
         raise AssertionError(
@@ -95,31 +84,14 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
     for param in llh_paramset.from_tag(ParamTag.NUISANCE):
         hypo_paramset[param.name].value = param.value
 
-    if args.energy_dependance is EnergyDependance.SPECTRAL:
-        bin_centers = np.sqrt(args.binning[:-1]*args.binning[1:])
-        bin_width = np.abs(np.diff(args.binning))
-        if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ] \
-           and args.fold_index:
-            args.spectral_index = -hypo_paramset['astroDeltaGamma'].value
-
-    if args.fix_source_ratio:
-        if args.energy_dependance is EnergyDependance.MONO:
-            source_flux = args.source_ratio
-        elif args.energy_dependance is EnergyDependance.SPECTRAL:
-            source_flux = np.array(
-                [fr * np.power(bin_centers, args.spectral_index)
-                 for fr in args.source_ratio]
-            ).T
-    else:
-        if args.energy_dependance is EnergyDependance.MONO:
-            source_flux = fr_utils.angles_to_fr(
-                llh_paramset.from_tag(ParamTag.SRCANGLES, values=True)
-            )
-        elif args.energy_dependance is EnergyDependance.SPECTRAL:
-            source_flux = np.array(
-                [fr * np.power(bin_centers, args.spectral_index)
-                 for fr in fr_utils.angles_to_fr(theta[-2:])]
-            ).T
+    bin_centers = np.sqrt(args.binning[:-1]*args.binning[1:])
+    bin_width = np.abs(np.diff(args.binning))
+    spectral_index = -hypo_paramset['astroDeltaGamma'].value
+
+    source_flux = np.array(
+        [fr * np.power(bin_centers, spectral_index)
+         for fr in args.source_ratio]
+    ).T
 
     bsm_angles = llh_paramset.from_tag(
         [ParamTag.SCALE, ParamTag.MMANGLES], values=True
@@ -139,21 +111,7 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
 
     if args.no_bsm:
         fr = fr_utils.u_to_fr(source_flux, np.array(sm_u, dtype=np.complex256))
-    elif args.energy_dependance is EnergyDependance.MONO:
-        u = fr_utils.params_to_BSMu(
-            theta             = bsm_angles,
-            dim               = args.dimension,
-            energy            = args.energy,
-            mass_eigenvalues  = mass_eigenvalues,
-            sm_u              = sm_u,
-            no_bsm            = args.no_bsm,
-            fix_mixing        = args.fix_mixing,
-            fix_mixing_almost = args.fix_mixing_almost,
-            fix_scale         = args.fix_scale,
-            scale             = args.scale
-        )
-        fr = fr_utils.u_to_fr(source_flux, u)
-    elif args.energy_dependance is EnergyDependance.SPECTRAL:
+    else:
         mf_perbin = []
         for i_sf, sf_perbin in enumerate(source_flux):
             u = fr_utils.params_to_BSMu(
@@ -163,10 +121,7 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
                 mass_eigenvalues  = mass_eigenvalues,
                 sm_u              = sm_u,
                 no_bsm            = args.no_bsm,
-                fix_mixing        = args.fix_mixing,
-                fix_mixing_almost = args.fix_mixing_almost,
-                fix_scale         = args.fix_scale,
-                scale             = args.scale
+                texture           = args.texture,
             )
             fr = fr_utils.u_to_fr(sf_perbin, u)
             mf_perbin.append(fr)
@@ -181,35 +136,18 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
     for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)):
         param.value = flavour_angles[idx]
 
-    if args.likelihood is Likelihood.FLAT:
-        llh = 1.
-    elif args.likelihood is Likelihood.GAUSSIAN:
-        fr_bf = args.measured_ratio
-        llh = multi_gaussian(fr, fr_bf, args.sigma_ratio)
-    elif args.likelihood is Likelihood.GOLEMFIT:
-        llh = gf_utils.get_llh(fitter, hypo_paramset)
+    if args.likelihood is Likelihood.GOLEMFIT:
+        llh = gf_utils.get_llh(hypo_paramset)
     elif args.likelihood is Likelihood.GF_FREQ:
-        llh = gf_utils.get_llh_freq(fitter, hypo_paramset)
-
-    if args.save_measured_fr and args.burnin is False:
-        n = gen_identifier(args) + '.txt'
-        with open(args.output_measured_fr + n, 'a') as f:
-            f.write(r'{0:.3f} {1:.3f} {2:.3f}'.format(
-                float(fr[0]), float(fr[1]), float(fr[2])
-            ))
-            for p in llh_paramset:
-                f.write(r' {0:.3f}'.format(p.value))
-            f.write('     LLH = {0:.3f}'.format(llh))
-            f.write('\n')
-
+        llh = gf_utils.get_llh_freq(hypo_paramset)
     return llh
 
 
-def ln_prob(theta, args, asimov_paramset, llh_paramset, fitter):
+def ln_prob(theta, args, asimov_paramset, llh_paramset):
     lp = lnprior(theta, paramset=llh_paramset)
     if not np.isfinite(lp):
         return -np.inf
     return lp + triangle_llh(
         theta, args=args, asimov_paramset=asimov_paramset,
-        llh_paramset=llh_paramset, fitter=fitter
+        llh_paramset=llh_paramset
     )
diff --git a/utils/misc.py b/utils/misc.py
index 26bd828..36d5330 100644
--- a/utils/misc.py
+++ b/utils/misc.py
@@ -19,7 +19,8 @@ from operator import attrgetter
 
 import numpy as np
 
-from utils.enums import Likelihood, MixingScenario
+from utils.enums import str_enum
+from utils.enums import Likelihood, Texture
 
 
 class SortingHelpFormatter(argparse.HelpFormatter):
@@ -31,30 +32,20 @@ class SortingHelpFormatter(argparse.HelpFormatter):
 
 def solve_ratio(fr):
     denominator = reduce(gcd, fr)
-    return [int(x/denominator) for x in fr]
+    f = [int(x/denominator) for x in fr]
+    if f[0] > 1E3 or f[1] > 1E3 or f[2] > 1E3:
+        return '{0:.2f}_{1:.2f}_{2:.2f}'.format(fr[0], fr[1], fr[2])
+    else:
+        return '{0}_{1}_{2}'.format(f[0], f[1], f[2])
 
 
 def gen_identifier(args):
     f = '_DIM{0}'.format(args.dimension)
-    mr1, mr2, mr3 = solve_ratio(args.measured_ratio)
-    if args.fix_source_ratio:
-        sr1, sr2, sr3 = solve_ratio(args.source_ratio)
-        f += '_sfr_{0:G}_{1:G}_{2:G}_mfr_{3:G}_{4:G}_{5:G}'.format(
-            sr1, sr2, sr3, mr1, mr2, mr3
-        )
-    else:
-        f += '_mfr_{3:03d}_{4:03d}_{5:03d}'.format(mr1, mr2, mr3)
-
-    if args.fix_mixing is not MixingScenario.NONE:
-        f += '_{0}'.format(args.fix_mixing)
-    elif args.fix_mixing_almost:
-        f += '_fix_mixing_almost'
-    elif args.fix_scale:
-        f += '_fix_scale_{0}'.format(args.scale)
-
-    if args.likelihood is Likelihood.FLAT: f += '_flat'
-    elif args.likelihood is Likelihood.GAUSSIAN:
-        f += '_sigma_{0:03d}'.format(int(args.sigma_ratio*1000))
+    f += '_sfr_' + solve_ratio(args.source_ratio)
+    if args.data in [DataType.ASIMOV, DataType.REALISATION]:
+        f += '_mfr_' + solve_ratio(args.injected_ratio)
+    if args.Texture is not Texture.NONE:
+        f += '_{0}'.format(str_enum(args.texture))
     return f
 
 
@@ -163,3 +154,71 @@ def thread_factors(t):
         if t%x == 0:
             return (x, int(t/x))
 
+
+def centers(x):
+    return (x[:-1]+x[1:])*0.5
+
+
+def get_units(dimension):
+    if dimension == 3: return r' / \:{\rm GeV}'
+    if dimension == 4: return r''
+    if dimension == 5: return r' / \:{rm GeV}^{-1}'
+    if dimension == 6: return r' / \:{rm GeV}^{-2}'
+    if dimension == 7: return r' / \:{rm GeV}^{-3}'
+    if dimension == 8: return r' / \:{rm GeV}^{-4}'
+
+
+def calc_nbins(x):
+    n =  (np.max(x) - np.min(x)) / (2 * len(x)**(-1./3) * (np.percentile(x, 75) - np.percentile(x, 25)))
+    return np.floor(n)
+
+
+def calc_bins(x):
+    nbins = calc_nbins(x)
+    return np.linspace(np.min(x), np.max(x)+2, num=nbins+1)
+
+
+def most_likely(arr):
+    """Return the densest region given a 1D array of data."""
+    binning = calc_bins(arr)
+    harr = np.histogram(arr, binning)[0]
+    return centers(binning)[np.argmax(harr)]
+
+
+def interval(arr, percentile=68.):
+    """Returns the *percentile* shortest interval around the mode."""
+    center = most_likely(arr)
+    sarr = sorted(arr)
+    delta = np.abs(sarr - center)
+    curr_low = np.argmin(delta)
+    curr_up = curr_low
+    npoints = len(sarr)
+    while curr_up - curr_low < percentile/100.*npoints:
+        if curr_low == 0:
+            curr_up += 1
+        elif curr_up == npoints-1:
+            curr_low -= 1
+        elif sarr[curr_up]-sarr[curr_low-1] < sarr[curr_up+1]-sarr[curr_low]:
+            curr_low -= 1
+        elif sarr[curr_up]-sarr[curr_low-1] > sarr[curr_up+1]-sarr[curr_low]:
+            curr_up += 1
+        elif (curr_up - curr_low) % 2:
+            # they are equal so step half of the time up and down
+            curr_low -= 1
+        else:
+            curr_up += 1
+    return sarr[curr_low], center, sarr[curr_up]
+
+
+def flat_angles_to_u(x):
+    """Convert from angles to mixing elements."""
+    return abs(angles_to_u(x)).astype(np.float32).flatten().tolist()
+
+
+def myround(x, base=2, up=False, down=False):
+    if up == down and up is True: assert 0
+    if up: return int(base * np.round(float(x)/base-0.5))
+    elif down: return int(base * np.round(float(x)/base+0.5))
+    else: int(base * np.round(float(x)/base))
+
+
diff --git a/utils/mn.py b/utils/mn.py
index 6582c80..c60d316 100644
--- a/utils/mn.py
+++ b/utils/mn.py
@@ -24,7 +24,7 @@ def CubePrior(cube, ndim, n_params):
 
 
 def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset,
-           args, fitter):
+           args):
     if ndim != len(mn_paramset):
         raise AssertionError(
             'Length of MultiNest scan paramset is not the same as the input '
@@ -41,7 +41,6 @@ def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset,
         args=args,
         asimov_paramset=asimov_paramset,
         llh_paramset=llh_paramset,
-        fitter=fitter
     )
     return llh
 
@@ -61,13 +60,14 @@ def mn_argparse(parser):
     )
 
 
-def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter,
+def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args,
                 identifier='mn'):
     """Run the MultiNest algorithm to calculate the evidence."""
     n_params = len(mn_paramset)
 
     for n in mn_paramset.names:
         llh_paramset[n].value = mn_paramset[n].value
+    print 'llh_paramset', llh_paramset
 
     lnProbEval = partial(
         lnProb,
@@ -75,14 +75,13 @@ def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter,
         llh_paramset    = llh_paramset,
         asimov_paramset = asimov_paramset,
         args            = args,
-        fitter          = fitter
     )
 
     llh = '{0}'.format(args.likelihood).split('.')[1]
     data = '{0}'.format(args.data).split('.')[1]
-    sr1, sr2, sr3 = solve_ratio(args.source_ratio)
-    prefix = './mnrun/DIM{0}/{1}/{2}/s{3}{4}{5}/{6}'.format(
-        args.dimension, data, llh, sr1, sr2, sr3, identifier
+    src_string = solve_ratio(args.source_ratio)
+    prefix = './mnrun/DIM{0}/{1}/{2}/s{3}/{4}'.format(
+        args.dimension, data, llh, src_string, identifier
     )
     make_dir(prefix)
     print 'Running evidence calculation for {0}'.format(prefix)
diff --git a/utils/param.py b/utils/param.py
index 572b65a..558018e 100644
--- a/utils/param.py
+++ b/utils/param.py
@@ -16,9 +16,8 @@ from copy import deepcopy
 
 import numpy as np
 
-from utils.plot import get_units
 from utils.fr import fr_to_angles
-from utils.enums import DataType, Likelihood, MixingScenario, ParamTag, PriorsCateg
+from utils.enums import DataType, Likelihood, ParamTag, PriorsCateg
 
 
 class Param(object):
@@ -219,66 +218,3 @@ class ParamSet(Sequence):
         elif isinstance(p, ParamSet):
             param_sequence.extend(p.params)
         return ParamSet(param_sequence)
-
-
-def get_paramsets(args, nuisance_paramset):
-    """Make the paramsets for generating the Asmimov MC sample and also running
-    the MCMC.
-    """
-    asimov_paramset = []
-    llh_paramset = []
-
-    llh_paramset.extend(
-        [x for x in nuisance_paramset.from_tag(ParamTag.SM_ANGLES)]
-    )
-    if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]:
-        gf_nuisance = [x for x in nuisance_paramset.from_tag(ParamTag.NUISANCE)]
-        asimov_paramset.extend(gf_nuisance)
-        llh_paramset.extend(gf_nuisance)
-    for parm in llh_paramset:
-        parm.value = args.__getattribute__(parm.name)
-    tag = ParamTag.BESTFIT
-    try:
-        flavour_angles = fr_to_angles(args.measured_ratio)
-    except:
-        flavour_angles = fr_to_angles(args.injected_ratio)
-    asimov_paramset.extend([
-        Param(name='astroFlavorAngle1', value=flavour_angles[0], ranges=[0., 1.], std=0.2, tag=tag),
-        Param(name='astroFlavorAngle2', value=flavour_angles[1], ranges=[-1., 1.], std=0.2, tag=tag),
-    ])
-    asimov_paramset = ParamSet(asimov_paramset)
-
-    if hasattr(args, 'fix_source_ratio'):
-        if args.fix_mixing is MixingScenario.NONE and not args.fix_mixing_almost:
-            tag = ParamTag.MMANGLES
-            llh_paramset.extend([
-                Param(name='np_s_12^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag),
-                Param(name='np_c_13^4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag),
-                Param(name='np_s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^2', tag=tag),
-                Param(name='np_dcp', value=np.pi, ranges=[0., 2*np.pi], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag)
-            ])
-        if args.fix_mixing_almost:
-            tag = ParamTag.MMANGLES
-            llh_paramset.extend([
-                Param(name='np_s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag)
-            ])
-        if not args.fix_scale:
-            tag = ParamTag.SCALE
-            if hasattr(args, 'dimension'):
-                llh_paramset.append(
-                    Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3,
-                          tex=r'{\rm log}_{10}\left (\Lambda^{-1}'+get_units(args.dimension)+r'\right )', tag=tag)
-                )
-            elif hasattr(args, 'dimensions'):
-                llh_paramset.append(
-                    Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3,
-                          tex=r'{\rm log}_{10}\left (\Lambda^{-1} / GeV^{-d+4}\right )', tag=tag)
-                )
-        if not args.fix_source_ratio:
-            tag = ParamTag.SRCANGLES
-            llh_paramset.extend([
-                Param(name='s_phi4', value=0.5, ranges=[0., 1.],  std=0.2, tex=r'sin^4(\phi)', tag=tag),
-                Param(name='c_2psi', value=0.5, ranges=[-1., 1.], std=0.2, tex=r'cos(2\psi)', tag=tag)
-            ])
-    llh_paramset = ParamSet(llh_paramset)
-    return asimov_paramset, llh_paramset
diff --git a/utils/plot.py b/utils/plot.py
index 6161cfb..29ef5dc 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -19,12 +19,13 @@ from scipy.interpolate import splev, splprep
 from scipy.ndimage.filters import gaussian_filter
 
 import matplotlib as mpl
+import matplotlib.patches as patches
+import matplotlib.gridspec as gridspec
 mpl.use('Agg')
+
 from matplotlib import rc
 from matplotlib import pyplot as plt
 from matplotlib.offsetbox import AnchoredText
-import matplotlib.patches as patches
-import matplotlib.gridspec as gridspec
 from matplotlib.lines import Line2D
 from matplotlib.patches import Patch
 
@@ -32,22 +33,19 @@ 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, str_enum
-from utils.enums import Likelihood, MixingScenario, ParamTag, StatCateg
-from utils.enums import Texture
+from utils.enums import Likelihood, ParamTag, StatCateg, Texture
 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
+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'):
@@ -58,64 +56,16 @@ if 'submitter' in socket.gethostname():
     rc('text', usetex=False)
 
 
-def centers(x):
-    return (x[:-1]+x[1:])*0.5
-
-
-def get_units(dimension):
-    if dimension == 3: return r' / \:{\rm GeV}'
-    if dimension == 4: return r''
-    if dimension == 5: return r' / \:{rm GeV}^{-1}'
-    if dimension == 6: return r' / \:{rm GeV}^{-2}'
-    if dimension == 7: return r' / \:{rm GeV}^{-3}'
-    if dimension == 8: return r' / \:{rm GeV}^{-4}'
-
-
-def calc_nbins(x):
-    n =  (np.max(x) - np.min(x)) / (2 * len(x)**(-1./3) * (np.percentile(x, 75) - np.percentile(x, 25)))
-    return np.floor(n)
-
-
-def calc_bins(x):
-    nbins = calc_nbins(x)
-    return np.linspace(np.min(x), np.max(x)+2, num=nbins+1)
-
-
-def most_likely(arr):
-    """Return the densest region given a 1D array of data."""
-    binning = calc_bins(arr)
-    harr = np.histogram(arr, binning)[0]
-    return centers(binning)[np.argmax(harr)]
-
-
-def interval(arr, percentile=68.):
-    """Returns the *percentile* shortest interval around the mode."""
-    center = most_likely(arr)
-    sarr = sorted(arr)
-    delta = np.abs(sarr - center)
-    curr_low = np.argmin(delta)
-    curr_up = curr_low
-    npoints = len(sarr)
-    while curr_up - curr_low < percentile/100.*npoints:
-        if curr_low == 0:
-            curr_up += 1
-        elif curr_up == npoints-1:
-            curr_low -= 1
-        elif sarr[curr_up]-sarr[curr_low-1] < sarr[curr_up+1]-sarr[curr_low]:
-            curr_low -= 1
-        elif sarr[curr_up]-sarr[curr_low-1] > sarr[curr_up+1]-sarr[curr_low]:
-            curr_up += 1
-        elif (curr_up - curr_low) % 2:
-            # they are equal so step half of the time up and down
-            curr_low -= 1
-        else:
-            curr_up += 1
-    return sarr[curr_low], center, sarr[curr_up]
-
-
-def flat_angles_to_u(x):
-    """Convert from angles to mixing elements."""
-    return abs(angles_to_u(x)).astype(np.float32).flatten().tolist()
+def gen_figtext(args):
+    """Generate the figure text."""
+    t = ''
+    t += 'Source flavour ratio = [{0}]'.format(solve_ratio(args.source_ratio))
+    if args.data in [DataType.ASIMOV, DataType.REALISATION]:
+        t += '\nIC injected flavour ratio = [{0}]'.format(
+            solve_ratio(args.injected_ratio)
+        )
+    t += '\nDimension = {0}'.format(args.dimension)
+    return t
 
 
 def plot_Tchain(Tchain, axes_labels, ranges):
@@ -139,39 +89,6 @@ def plot_Tchain(Tchain, axes_labels, ranges):
     return g
 
 
-def gen_figtext(args):
-    """Generate the figure text."""
-    t = ''
-    mr1, mr2, mr3 = misc_utils.solve_ratio(args.measured_ratio)
-    if args.fix_source_ratio:
-        sr1, sr2, sr3 = misc_utils.solve_ratio(args.source_ratio)
-        t += 'Source flavour ratio = [{0}, {1}, {2}]'.format(sr1, sr2, sr3)
-        if args.data in [DataType.ASIMOV, DataType.REALISATION]:
-            t += '\nIC observed flavour ratio = [{0}, {1}, {2}]'.format(
-                mr1, mr2, mr3
-            )
-        t += '\nDimension = {0}'.format(args.dimension)
-    else:
-        if args.data in [DataType.ASIMOV, DataType.REALISATION]:
-            t += 'IC observed flavour ratio = [{0}, {1}, ' \
-                    '{2}]\nDimension = {3}'.format(
-                        mr1, mr2, mr3, args.dimension, int(args.energy)
-                    )
-    if args.fix_scale:
-        t += 'Scale = {0}'.format(args.scale)
-    if args.likelihood is Likelihood.GAUSSIAN:
-        t += '\nSigma = {0:.3f}'.format(args.sigma_ratio)
-    if args.energy_dependance is EnergyDependance.SPECTRAL:
-        if not args.fold_index:
-            t += '\nSpectral Index = {0}'.format(int(args.spectral_index))
-        t += '\nBinning = [{0}, {1}] TeV - {2} bins'.format(
-            int(args.binning[0]/1e3), int(args.binning[-1]/1e3), len(args.binning)-1
-        )
-    elif args.energy_dependance is EnergyDependance.MONO:
-        t += '\nEnergy = {0} TeV'.format(int(args.energy/1e3))
-    return t
-
-
 def chainer_plot(infile, outfile, outformat, args, llh_paramset, fig_text=None):
     """Make the triangle plot."""
     if hasattr(args, 'plot_elements'):
@@ -287,13 +204,6 @@ def chainer_plot(infile, outfile, outformat, args, llh_paramset, fig_text=None):
             g.export(outfile+'_elements.'+of)
 
 
-def myround(x, base=2, up=False, down=False):
-    if up == down and up is True: assert 0
-    if up: return int(base * np.round(float(x)/base-0.5))
-    elif down: return int(base * np.round(float(x)/base+0.5))
-    else: int(base * np.round(float(x)/base))
-
-
 def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
     """Make MultiNest factor or LLH value plot."""
     print 'Making Statistic plot'
@@ -337,7 +247,7 @@ def plot_statistic(data, outfile, outformat, args, scale_param, label=None):
 
     ax.plot(scales_rm, reduced_ev)
 
-    ax.axhline(y=np.log(10**(bayes_K)), 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)
@@ -394,7 +304,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**(bayes_K))] # 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)
@@ -437,8 +347,8 @@ def plot_sens_full(data, outfile, outformat, args):
         fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150)
 
 
-def plot_sens_fixed_angle_pretty(data, outfile, outformat, args):
-    print 'Making FIXED_ANGLE_PRETTY sensitivity plot'
+def plot_table_sens(data, outfile, outformat, args):
+    print 'Making TABLE sensitivity plot'
     argsc = deepcopy(args)
     dims = len(data)
     LV_atmo_90pc_limits = {
@@ -598,14 +508,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**(bayes_K))]
+                    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
                 if len(al) == 0:
                     print 'No points for DIM {0} FRS {1}!'.format(dim, src)
                     continue
-                if reduced_ev[-1] < np.log(10**(bayes_K)) - 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
@@ -739,14 +649,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**(bayes_K))] # 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**(bayes_K)) - 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
@@ -903,7 +813,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**(bayes_K))) # 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
@@ -1175,7 +1085,7 @@ def plot_source_ternary(data, outfile, outformat, args):
     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)))
+            print '|||| SCEN = {0}'.format(str_enum(MixingScenario(isce+1)))
             fig = plt.figure(figsize=(8, 8))
             ax  = fig.add_subplot(111)
             tax = get_tax(ax, scale=nsrcs)
@@ -1204,14 +1114,14 @@ def plot_source_ternary(data, outfile, outformat, args):
                 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))]
+                    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:
+                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
@@ -1244,8 +1154,9 @@ def texture_label(x):
         raise AssertionError
 
 
-def plot_source_ternary_1D(data, outfile, outformat, args):
+def plot_x(data, outfile, outformat, args):
     """Ternary plot in the source flavour space for each operator texture."""
+    print 'Making X sensitivity plot'
     sources = args.source_ratios
     x_1D = []
     i_src_1D = []
@@ -1279,7 +1190,7 @@ def plot_source_ternary_1D(data, outfile, outformat, args):
         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)))
+            print '|||| SCEN = {0}'.format(str_enum(MixingScenario(isce+1)))
             H = np.full(len(x_1D), np.nan)
             for ix, x in enumerate(x_1D):
                 print '|||| X = {0}'.format(x)
@@ -1303,13 +1214,13 @@ def plot_source_ternary_1D(data, outfile, outformat, args):
                 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))]
+                    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:
+                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
@@ -1320,7 +1231,7 @@ def plot_source_ternary_1D(data, outfile, outformat, args):
             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='-',
+                    label=texture_label(MixingScenario(isce+1)), linestyle='-',
                     drawstyle='steps-pre')
             print 'H', H
             print