refactor

6 files changed, 923 insertions, 0 deletions

diff --git a/utils/__init__.py b/utils/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/utils/__init__.py
diff --git a/utils/enums.py b/utils/enums.py
new file mode 100644
index 0000000..31885de
--- /dev/null
+++ b/utils/enums.py
@@ -0,0 +1,68 @@
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : March 17, 2018
+
+"""
+Define Enums for the BSM flavour ratio analysis
+"""
+
+from enum import Enum
+
+
+class DataType(Enum):
+    REAL    = 1
+    FAKE    = 2
+    ASMIMOV = 3
+
+
+class FitCateg(Enum):
+    HESESPL  = 1
+    HESEDPL  = 2
+    ZPSPL    = 3
+    ZPDPL    = 4
+    NUNUBAR2 = 5
+
+
+class FitFlagCateg(Enum):
+    DEFAULT = 1
+    XS      = 2
+
+
+class FitPriorsCateg(Enum):
+    DEFAULT = 1
+    XS      = 2
+
+
+class Likelihood(Enum):
+    FLAT     = 1
+    GAUSSIAN = 2
+    GOLEMFIT = 3
+
+
+class Priors(Enum):
+    UNIFORM    = 1
+    LOGUNIFORM = 2
+    JEFFREYS   = 3
+
+
+class SteeringCateg(Enum):
+    BASELINE              = 1
+    HOLEICE               = 2
+    ABSORPTION            = 3
+    SCATTERING            = 4
+    SCATTERING_ABSORPTION = 5
+    STD                   = 6
+    STD_HALF1             = 7
+    STD_HALF2             = 8
+    LONGLIFE              = 9
+    DPL                   = 10
+
+
+class XSCateg(Enum):
+    HALF    = 1
+    NOM     = 2
+    TWICE   = 3
+    TWICE01 = 4
+    TWICE02 = 5
+
diff --git a/utils/fr.py b/utils/fr.py
new file mode 100644
index 0000000..ddcb5d2
--- /dev/null
+++ b/utils/fr.py
@@ -0,0 +1,353 @@
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : March 17, 2018
+
+"""
+Useful functions for the BSM flavour ratio analysis
+"""
+
+from __future__ import absolute_import, division
+
+import sys
+
+import numpy as np
+from scipy import linalg
+
+
+MASS_EIGENVALUES = [7.40E-23, 2.515E-21]
+"""SM mass eigenvalues"""
+
+
+def angles_to_fr(src_angles):
+    """Convert angular projection of the source flavour ratio back into the
+    flavour ratio.
+
+    Parameters
+    ----------
+    src_angles : list, length = 2
+        sin(phi)^4 and cos(psi)^2
+
+    Returns
+    ----------
+    flavour ratios (nue, numu, nutau)
+
+    Examples
+    ----------
+    >>> print angles_to_fr((0.3, 0.4))
+    (0.38340579025361626, 0.16431676725154978, 0.45227744249483393)
+
+    """
+    sphi4, c2psi = src_angles
+
+    psi = (0.5)*np.arccos(c2psi)
+
+    sphi2 = np.sqrt(sphi4)
+    cphi2 = 1. - sphi2
+    spsi2 = np.sin(psi)**2
+    cspi2 = 1. - spsi2
+
+    x = sphi2*cspi2
+    y = sphi2*spsi2
+    z = cphi2
+    return x, y, z
+
+
+def angles_to_u(bsm_angles):
+    """Convert angular projection of the mixing matrix elements back into the
+    mixing matrix elements.
+
+    Parameters
+    ----------
+    bsm_angles : list, length = 4
+        sin(12)^2, cos(13)^4, sin(23)^2 and deltacp
+
+    Returns
+    ----------
+    unitary numpy ndarray of shape (3, 3)
+
+    Examples
+    ----------
+    >>> from fr import angles_to_u
+    >>> print angles_to_u((0.2, 0.3, 0.5, 1.5))
+    array([[ 0.66195018+0.j        ,  0.33097509+0.j        ,  0.04757188-0.6708311j ],
+           [-0.34631487-0.42427084j,  0.61741198-0.21213542j,  0.52331757+0.j        ],
+           [ 0.28614067-0.42427084j, -0.64749908-0.21213542j,  0.52331757+0.j        ]])
+
+    """
+    s12_2, c13_4, s23_2, dcp = bsm_angles
+    dcp = np.complex128(dcp)
+
+    c12_2 = 1. - s12_2
+    c13_2 = np.sqrt(c13_4)
+    s13_2 = 1. - c13_2
+    c23_2 = 1. - s23_2
+
+    t12 = np.arcsin(np.sqrt(s12_2))
+    t13 = np.arccos(np.sqrt(c13_2))
+    t23 = np.arcsin(np.sqrt(s23_2))
+
+    c12 = np.cos(t12)
+    s12 = np.sin(t12)
+    c13 = np.cos(t13)
+    s13 = np.sin(t13)
+    c23 = np.cos(t23)
+    s23 = np.sin(t23)
+
+    p1 = np.array([[1   , 0   , 0]                   , [0    , c23 , s23] , [0                   , -s23 , c23]] , dtype=np.complex128)
+    p2 = np.array([[c13 , 0   , s13*np.exp(-1j*dcp)] , [0    , 1   , 0]   , [-s13*np.exp(1j*dcp) , 0    , c13]] , dtype=np.complex128)
+    p3 = np.array([[c12 , s12 , 0]                   , [-s12 , c12 , 0]   , [0                   , 0    , 1]]   , dtype=np.complex128)
+
+    u = np.dot(np.dot(p1, p2), p3)
+    return u
+
+
+def cardano_eqn(ham):
+    """Diagonalise the effective Hamiltonian 3x3 matrix into the form
+    h_{eff} = UE_{eff}U^{dagger} using the procedure in PRD91, 052003 (2015).
+
+    Parameters
+    ----------
+    ham : numpy ndarray of shape (3, 3)
+        sin(12)^2, cos(13)^4, sin(23)^2 and deltacp
+
+    Returns
+    ----------
+    unitary numpy ndarray of shape (3, 3)
+
+    Examples
+    ----------
+    >>> import numpy as np
+    >>> from fr import cardano_eqn
+    >>> ham = np.array(
+    >>>     [[ 0.66195018+0.j        ,  0.33097509+0.j        ,  0.04757188-0.6708311j ],
+    >>>      [-0.34631487-0.42427084j,  0.61741198-0.21213542j,  0.52331757+0.j        ],
+    >>>      [ 0.28614067-0.42427084j, -0.64749908-0.21213542j,  0.52331757+0.j        ]]
+    >>> )
+    >>> print cardano_eqn(ham)
+    array([[-0.11143379-0.58863683j, -0.09067747-0.48219068j, 0.34276625-0.08686465j],
+           [ 0.14835519+0.47511473j, -0.18299305+0.40777481j, 0.31906300+0.82514223j],
+           [-0.62298966+0.07231745j, -0.61407815-0.42709603j, 0.03660313+0.30160428j]])
+
+    """
+    if np.shape(ham) != (3, 3):
+        raise ValueError(
+            'Input matrix should be a square and dimension 3, '
+            'got\n{0}'.format(ham)
+        )
+
+    a = -np.trace(ham)
+    b = (0.5) * ((np.trace(ham))**2 - np.trace(np.dot(ham, ham)))
+    c = -linalg.det(ham)
+
+    Q = (1/9.) * (a**2 - 3*b)
+    R = (1/54.) * (2*a**3 - 9*a*b + 27*c)
+    theta = np.arccos(R / np.sqrt(Q**3))
+
+    E1 = -2 * np.sqrt(Q) * np.cos(theta/3.) - (1/3.)*a
+    E2 = -2 * np.sqrt(Q) * np.cos((theta - 2.*np.pi)/3.) - (1/3.)*a
+    E3 = -2 * np.sqrt(Q) * np.cos((theta + 2.*np.pi)/3.) - (1/3.)*a
+
+    A1 = ham[1][2] * (ham[0][0] - E1) - ham[1][0]*ham[0][2]
+    A2 = ham[1][2] * (ham[0][0] - E2) - ham[1][0]*ham[0][2]
+    A3 = ham[1][2] * (ham[0][0] - E3) - ham[1][0]*ham[0][2]
+
+    B1 = ham[2][0] * (ham[1][1] - E1) - ham[2][1]*ham[1][0]
+    B2 = ham[2][0] * (ham[1][1] - E2) - ham[2][1]*ham[1][0]
+    B3 = ham[2][0] * (ham[1][1] - E3) - ham[2][1]*ham[1][0]
+
+    C1 = ham[1][0] * (ham[2][2] - E1) - ham[1][2]*ham[2][0]
+    C2 = ham[1][0] * (ham[2][2] - E2) - ham[1][2]*ham[2][0]
+    C3 = ham[1][0] * (ham[2][2] - E3) - ham[1][2]*ham[2][0]
+
+    N1 = np.sqrt(abs(A1*B1)**2 + abs(A1*C1)**2 + abs(B1*C1)**2)
+    N2 = np.sqrt(abs(A2*B2)**2 + abs(A2*C2)**2 + abs(B2*C2)**2)
+    N3 = np.sqrt(abs(A3*B3)**2 + abs(A3*C3)**2 + abs(B3*C3)**2)
+
+    mm = np.array([
+        [np.conjugate(B1)*C1 / N1, np.conjugate(B2)*C2 / N2, np.conjugate(B3)*C3 / N3],
+        [A1*C1 / N1, A2*C2 / N2, A3*C3 / N3],
+        [A1*B1 / N1, A2*B2 / N2, A3*B3 / N3]
+    ])
+    return mm
+
+
+def normalise_fr(fr):
+    """Normalise an input flavour combination to a flavour ratio.
+
+    Parameters
+    ----------
+    fr : list, length = 3
+        flavour combination
+
+    Returns
+    ----------
+    numpy ndarray flavour ratio
+
+    Examples
+    ----------
+    >>> from fr import normalise_fr
+    >>> print normalise_fr((1, 2, 3))
+    array([ 0.16666667,  0.33333333,  0.5       ])
+
+    """
+    return np.array(fr) / float(np.sum(fr))
+
+
+NUFIT_U = angles_to_u((0.307, (1-0.02195)**2, 0.565, 3.97935))
+"""NuFIT mixing matrix (s_12^2, c_13^4, s_23^2, dcp)"""
+
+
+def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
+                   nufit_u=NUFIT_U, no_bsm=False, fix_mixing=False,
+                   fix_scale=False, scale=None, check_uni=True):
+    """Construct the BSM mixing matrix from the BSM parameters.
+
+    Parameters
+    ----------
+    theta : list, length > 3
+        BSM parameters
+
+    dim : int
+        Dimension of BSM physics
+
+    energy : float
+        Energy in GeV
+
+    mass_eigenvalues : list, length = 2
+        SM mass eigenvalues
+
+    nufit_u : numpy ndarray, dimension 3
+        SM NuFIT mixing matrix
+
+    no_bsm : bool
+        Turn off BSM behaviour
+
+    fix_mixing : bool
+        Fix the BSM mixing angles
+
+    fix_scale : bool
+        Fix the BSM scale
+
+    scale : float
+        Used with fix_scale - scale at which to fix
+
+    check_uni : bool
+        Check the resulting BSM mixing matrix is unitary
+
+    Returns
+    ----------
+    unitary numpy ndarray of shape (3, 3)
+
+    Examples
+    ----------
+    >>> from fr import params_to_BSMu
+    >>> print params_to_BSMu((0.2, 0.3, 0.5, 1.5, -20), dim=3, energy=1000)
+    array([[ 0.18658169 -6.34190523e-01j, -0.26460391 +2.01884200e-01j, 0.67247096 -9.86808417e-07j],
+           [-0.50419832 +2.14420570e-01j, -0.36013768 +5.44254868e-01j, 0.03700961 +5.22039894e-01j],
+           [-0.32561308 -3.95946524e-01j,  0.64294909 -2.23453580e-01j, 0.03700830 +5.22032403e-01j]])
+
+    """
+    if np.shape(nufit_u) != (3, 3):
+        raise ValueError(
+            'Input matrix should be a square and dimension 3, '
+            'got\n{0}'.format(ham)
+        )
+
+    if fix_mixing:
+        s12_2, c13_4, s23_2, dcp, sc2 = 0.5, 1.0-1E-6, 0.5, 0., theta
+    elif fix_scale:
+        s12_2, c13_4, s23_2, dcp = theta
+        sc2 = np.log10(scale)
+    else:
+        s12_2, c13_4, s23_2, dcp, sc2 = theta
+    sc2 = np.power(10., sc2)
+    sc1 = sc2 / 100.
+
+    mass_matrix = np.array(
+        [[0, 0, 0], [0, mass_eigenvalues[0], 0], [0, 0, mass_eigenvalues[1]]]
+    )
+    sm_ham = (1./(2*energy))*np.dot(nufit_u, np.dot(mass_matrix, nufit_u.conj().T))
+    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(
+            [[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_ham = sm_ham + bsm_term
+        eg_vector = cardano_eqn(bsm_ham)
+
+    if check_uni:
+        tu = test_unitarity(eg_vector)
+        if not abs(np.trace(tu) - 3.) < 1e-5 or \
+           not abs(np.sum(tu) - 3.) < 1e-5:
+            raise AssertionError(
+                'Matrix is not unitary!\neg_vector\n{0}\ntest '
+                'u\n{1}'.format(eg_vector, tu)
+            )
+    return eg_vector
+
+
+def test_unitarity(x, prnt=False):
+    """Test the unitarity of a matrix.
+
+    Parameters
+    ----------
+    x : numpy ndarray
+        Matrix to evaluate
+
+    prnt : bool
+        Print the result
+
+    Returns
+    ----------
+    numpy ndarray
+
+    Examples
+    ----------
+    >>> from fr import test_unitarity
+    >>> x = np.identity(3)
+    >>> print test_unitarity(x)
+    array([[ 1.,  0.,  0.],
+           [ 0.,  1.,  0.],
+           [ 0.,  0.,  1.]])
+
+    """
+    f = abs(np.dot(x, x.conj().T))
+    if prnt:
+        print 'Unitarity test:\n{0}'.format(f)
+    return f
+
+
+def u_to_fr(source_fr, matrix):
+    """Compute the observed flavour ratio assuming decoherence.
+
+    Parameters
+    ----------
+    source_fr : list, length = 3
+        Source flavour ratio components
+
+    matrix : numpy ndarray, dimension 3
+        Mixing matrix
+
+    Returns
+    ----------
+    Measured flavour ratio
+
+    Examples
+    ----------
+    >>> from fr import params_to_BSMu, u_to_fr
+    >>> print u_to_fr((1, 2, 0), params_to_BSMu((0.2, 0.3, 0.5, 1.5, -20), 3, 1000))
+        array([ 0.33740075,  0.33176584,  0.33083341])
+
+    """
+    # TODO(shivesh): energy dependence
+    composition = np.einsum(
+        'ai, bi, a -> b', abs(matrix)**2, abs(matrix)**2, source_fr
+    )
+    ratio = composition / np.sum(source_fr)
+    return ratio
+
diff --git a/utils/gf.py b/utils/gf.py
new file mode 100644
index 0000000..766c161
--- /dev/null
+++ b/utils/gf.py
@@ -0,0 +1,105 @@
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : March 17, 2018
+
+"""
+Useful GolemFit wrappers for the BSM flavour ratio analysis
+"""
+
+from __future__ import absolute_import, division
+
+import argparse
+from functools import partial
+
+import GolemFitPy as gf
+
+from utils.enums import *
+from utils.misc import enum_keys, enum_parse
+
+
+def data_distributions(fitter):
+    hdat = fitter.GetDataDistribution()
+    binedges = np.asarray([fitter.GetZenithBinsData(), fitter.GetEnergyBinsData()])
+    return hdat, binedges
+
+
+def fit_flags(fitflag_categ):
+    flags = gf.FitParametersFlag()
+    if fitflag_categ is FitFlagCateg.xs:
+        # False means it's not fixed in minimization
+        flags.NeutrinoAntineutrinoRatio = False
+    return flags
+
+
+def fit_params(fit_categ):
+    params = gf.FitParameters()
+    params.astroNorm = 7.5
+    params.astroDeltaGamma = 0.9
+    if fit_categ is FitCateg.hesespl:
+        params.astroNormSec = 0
+    elif fit_categ is FitCateg.hesedpl:
+        params.astroNormSec=7.0
+    elif fit_categ is FitCateg.zpspl:
+        # zero prompt, single powerlaw
+        params.promptNorm = 0
+        params.astroNormSec = 0
+    elif fit_categ is FitCateg.zpdpl:
+        # zero prompt, double powerlaw
+        params.promptNorm = 0
+        params.astroNormSec=7.0
+    elif fit_categ is FitCateg.nunubar2:
+        params.NeutrinoAntineutrinoRatio = 2
+
+
+def fit_priors(fitpriors_categ):
+    priors = gf.Priors()
+    if fitpriors_categ == FitPriorsCateg.xs:
+        priors.promptNormCenter = 1
+        priors.promptNormWidth = 3
+        priors.astroDeltaGammaCenter = 0
+        priors.astroDeltaGammaWidth = 1
+    return priors
+
+
+def gen_steering_params(steering_categ, quiet=False):
+    params = gf.SteeringParams()
+    if quiet: params.quiet = True
+    params.fastmode = False
+    params.do_HESE_reshuffle = False
+    params.numc_tag = steering_categ.name
+    params.baseline_astro_spectral_index = -2.
+    if steering_categ is SteeringCateg.LONGLIFE:
+        params.years = [999]
+        params.numc_tag = 'std_half1'
+    if steering_categ is SteeringCateg.DPL:
+        params.diffuse_fit_type = gf.DiffuseFitType.DoublePowerLaw
+        params.numc_tag = 'std_half1'
+    return params
+
+
+def gf_argparse(parser):
+    parser.add_argument(
+        '--data', default='real', type=partial(enum_parse, c=DataType),
+        choices=enum_keys(DataType), help='select datatype'
+    )
+    parser.add_argument(
+        '--ast', default='baseline', type=partial(enum_parse, c=SteeringCateg),
+        choices=enum_keys(SteeringCateg),
+        help='use asimov/fake dataset with specific steering'
+    )
+    parser.add_argument(
+        '--aft', default='hesespl', type=partial(enum_parse, c=FitCateg),
+        choices=enum_keys(FitCateg),
+        help='use asimov/fake dataset with specific Fit'
+    )
+    parser.add_argument(
+        '--axs', default='nom', type=partial(enum_parse, c=XSCateg),
+        choices=enum_keys(XSCateg),
+        help='use asimov/fake dataset with xs scaling'
+    )
+    parser.add_argument(
+        '--priors', default='uniform', type=partial(enum_parse, c=Priors),
+        choices=enum_keys(Priors), help='Bayesian priors'
+    )
+
diff --git a/utils/mcmc.py b/utils/mcmc.py
new file mode 100644
index 0000000..d91764f
--- /dev/null
+++ b/utils/mcmc.py
@@ -0,0 +1,167 @@
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : March 17, 2018
+
+"""
+Useful functions to use an MCMC for the BSM flavour ratio analysis
+"""
+
+from __future__ import absolute_import, division
+
+import os
+import errno
+
+import argparse
+from functools import partial
+
+import emcee
+import tqdm
+
+import numpy as np
+from scipy.stats import multivariate_normal
+
+from utils import fr as fr_utils
+from utils.enums import Likelihood
+from utils.misc import enum_keys, enum_parse, parse_bool
+
+
+def lnprior(theta, paramset):
+    """Priors on theta."""
+    ranges = paramset.ranges
+    for value, range in zip(theta, ranges):
+        if range[0] <= value <= range[1]:
+            pass
+        else: return -np.inf
+    return 0.
+
+
+def mcmc(p0, triangle_llh, lnprior, ndim, nwalkers, burnin, nsteps, ntemps=1, threads=1):
+    """Run the MCMC."""
+    sampler = emcee.PTSampler(
+        ntemps, nwalkers, ndim, triangle_llh, lnprior, threads=threads
+    )
+
+    print "Running burn-in"
+    for result in tqdm.tqdm(sampler.sample(p0, iterations=burnin), total=burnin):
+        pos, prob, state = result
+    sampler.reset()
+    print "Finished burn-in"
+
+    print "Running"
+    for _ in tqdm.tqdm(sampler.sample(pos, iterations=nsteps), total=nsteps):
+        pass
+    print "Finished"
+
+    samples = sampler.chain[0, :, :, :].reshape((-1, ndim))
+    print 'acceptance fraction', sampler.acceptance_fraction
+    print 'sum of acceptance fraction', np.sum(sampler.acceptance_fraction)
+    print 'np.unique(samples[:,0]).shape', np.unique(samples[:,0]).shape
+    try:
+        print 'autocorrelation', sampler.acor
+    except:
+        print 'WARNING : NEED TO RUN MORE SAMPLES'
+
+    return samples
+
+
+def mcmc_argparse(parser):
+    parser.add_argument(
+        '--likelihood', default='gaussian', type=partial(enum_parse, c=Likelihood),
+        choices=enum_keys(Likelihood), help='likelihood contour'
+    )
+    parser.add_argument(
+        '--run-mcmc', type=parse_bool, default='True',
+        help='Run the MCMC'
+    )
+    parser.add_argument(
+        '--burnin', type=int, default=100,
+        help='Amount to burnin'
+    )
+    parser.add_argument(
+        '--nwalkers', type=int, default=100,
+        help='Number of walkers'
+    )
+    parser.add_argument(
+        '--nsteps', type=int, default=2000,
+        help='Number of steps to run'
+    )
+
+
+def gaussian_llh(fr, fr_bf, sigma):
+    """Multivariate gaussian likelihood."""
+    cov_fr = np.identity(3) * sigma
+    return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr))
+
+
+def gaussian_seed(paramset, ntemps, nwalkers):
+    """Get gaussian seed values for the MCMC."""
+    ndim = len(paramset)
+    p0 = np.random.normal(
+        paramset.values, paramset.stds, size=[ntemps, nwalkers, ndim]
+    )
+    return p0
+
+
+def save_chains(chains, outfile):
+    """Save the chains.
+
+    Parameters
+    ----------
+    chains : numpy ndarray
+        MCMC chains to save
+
+    outfile : str
+        Output file location of chains
+
+    """
+    try:
+        os.makedirs(outfile[:-len(os.path.basename(outfile))])
+    except OSError as exc:  # Python >2.5
+        if exc.errno == errno.EEXIST and os.path.isdir(outfile[:-len(os.path.basename(outfile))]):
+            pass
+        else:
+            raise
+    print 'Saving chains to location {0}'.format(outfile+'.npy')
+    np.save(outfile+'.npy', chains)
+
+
+def triangle_llh(theta, args):
+    """-Log likelihood function for a given theta."""
+    if args.fix_source_ratio:
+        fr1, fr2, fr3 = args.source_ratio
+        bsm_angles = theta[-5:]
+    else:
+        fr1, fr2, fr3 = fr_utils.angles_to_fr(theta[-2:])
+        bsm_angles = theta[-7:-2]
+
+    u = fr_utils.params_to_BSMu(
+        theta      = bsm_angles,
+        dim        = args.dimension,
+        energy     = args.energy,
+        no_bsm     = args.no_bsm,
+        fix_mixing = args.fix_mixing,
+        fix_scale  = args.fix_scale,
+        scale      = args.scale
+    )
+
+    fr = fr_utils.u_to_fr((fr1, fr2, fr3), u)
+    fr_bf = args.measured_ratio
+    if args.likelihood is Likelihood.FLAT:
+        return 1.
+    elif args.likelihood is Likelihood.GAUSSIAN:
+        return gaussian_llh(fr, fr_bf, args.sigma_ratio)
+    elif args.likelihood is Likelihood.GOLEMFIT:
+        raise NotImplementedError('TODO')
+        import GolemFitPy as gf
+        from collections import namedtuple
+        datapaths = gf.DataPaths()
+        IceModels = namedtuple('IceModels', 'std_half2')
+        fitters = IceModels(*[
+            gf.GolemFit(datapaths,
+                        gf.gen_steering_params(SteeringCateg.__members__[ice]),
+                        xs_params(XSCateg.nom)) for ice in IceModels._fields])
+        for fitter in fitters:
+            fitter.SetFitParametersFlag(fit_flags(FitFlagCateg.xs))
+            fitter.SetFitPriors(fit_priors(FitPriorsCateg.xs))
+
diff --git a/utils/misc.py b/utils/misc.py
new file mode 100644
index 0000000..5c3eb2e
--- /dev/null
+++ b/utils/misc.py
@@ -0,0 +1,230 @@
+# author : S. Mandalia
+#          s.p.mandalia@qmul.ac.uk
+#
+# date   : March 17, 2018
+
+"""
+Misc functions for the BSM flavour ratio analysis
+"""
+
+from __future__ import absolute_import, division
+
+import os
+import errno
+from collections import Sequence
+import multiprocessing
+
+import numpy as np
+
+from utils.enums import Likelihood
+
+
+class Param(object):
+    """Parameter class to store parameters.
+    """
+    def __init__(self, name, value, ranges, std=None, tex=None):
+        self._ranges = None
+        self._tex = None
+
+        self.name = name
+        self.value = value
+        self.ranges = ranges
+        self.std = std
+        self.tex = tex
+
+    @property
+    def ranges(self):
+        return tuple(self._ranges)
+
+    @ranges.setter
+    def ranges(self, values):
+        self._ranges = [val for val in values]
+
+    @property
+    def tex(self):
+        return r'{0}'.format(self.tex)
+
+    @tex.setter
+    def tex(self, t):
+        self._tex = t if t is not None else r'{\rm %s}' % self.name
+
+
+class ParamSet(Sequence):
+    """Container class for a set of parameters.
+    """
+    def __init__(self, *args):
+        param_sequence = []
+        for arg in args:
+            try:
+                param_sequence.extend(arg)
+            except TypeError:
+                param_sequence.append(arg)
+
+        # Disallow duplicated params
+        all_names = [p.name for p in param_sequence]
+        unique_names = set(all_names)
+        if len(unique_names) != len(all_names):
+            duplicates = set([x for x in all_names if all_names.count(x) > 1])
+            raise ValueError('Duplicate definitions found for param(s): ' +
+                             ', '.join(str(e) for e in duplicates))
+
+        # Elements of list must be Param type
+        assert all([isinstance(x, Param) for x in param_sequence]), \
+                'All params must be of type "Param"'
+
+        self._params = param_sequence
+
+    def __len__(self):
+        return len(self._params)
+
+    def __getitem__(self, i):
+        if isinstance(i, int):
+            return self._params[i]
+        elif isinstance(i, basestring):
+            return self._by_name[i]
+
+    def __getattr__(self, attr):
+        try:
+            return super(ParamSet, self).__getattribute__(attr)
+        except AttributeError:
+            t, v, tb = sys.exc_info()
+            try:
+                return self[attr]
+            except KeyError:
+                raise t, v, tb
+
+    def __iter__(self):
+        return iter(self._params)
+
+    @property
+    def _by_name(self):
+        return {obj.name: obj for obj in self._params}
+
+    @property
+    def names(self):
+        return tuple([obj.name for obj in self._params])
+
+    @property
+    def values(self):
+        return tuple([obj.value for obj in self._params])
+
+    @property
+    def ranges(self):
+        return tuple([obj.ranges for obj in self._params])
+
+    @property
+    def stds(self):
+        return tuple([obj.std for obj in self._params])
+
+    @property
+    def params(self):
+        return self._params
+
+    def to_dict(self):
+        return {obj.name: obj.value for obj in self._params}
+
+
+def gen_outfile_name(args):
+    """Generate a name for the output file based on the input args.
+
+    Parameters
+    ----------
+    args : argparse
+        argparse object to print
+
+    """
+    mr = args.measured_ratio
+    si = args.sigma_ratio
+    if args.fix_source_ratio:
+        sr = args.source_ratio
+        if args.fix_mixing:
+            outfile = args.outfile+'_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_fix_mixing'.format(
+                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000),
+                int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), args.dimension
+            )
+        elif args.fix_scale:
+            outfile = args.outfile+'_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_fixed_scale_{8}'.format(
+                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000),
+                int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), args.dimension,
+                args.scale
+            )
+        else:
+            outfile = args.outfile+'_{0:03d}_{1:03d}_{2:03d}_{3:04d}_sfr_{4:03d}_{5:03d}_{6:03d}_DIM{7}_single_scale'.format(
+                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100), int(si*1000),
+                int(sr[0]*100), int(sr[1]*100), int(sr[2]*100), args.dimension
+            )
+    else:
+        if args.fix_mixing:
+            outfile = args.outfile+'_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fix_mixing'.format(
+                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100),
+                int(si*1000), args.dimension
+            )
+        elif args.fix_scale:
+            outfile = args.outfile+'_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fixed_scale_{5}'.format(
+                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100),
+                int(si*1000), args.dimension, args.scale
+            )
+        else:
+            outfile = args.outfile+'_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}'.format(
+                int(mr[0]*100), int(mr[1]*100), int(mr[2]*100),
+                int(si*1000), args.dimension
+            )
+    if args.likelihood is Likelihood.FLAT: outfile += '_flat'
+    return outfile
+
+
+def parse_bool(s):
+    """Parse a string to a boolean.
+
+    Parameters
+    ----------
+    s : str
+        String to parse
+
+    Returns
+    ----------
+    bool
+
+    Examples
+    ----------
+    >>> from misc import parse_bool
+    >>> print parse_bool('true')
+    True
+
+    """
+    if s.lower() == 'true':
+        return True
+    elif s.lower() == 'false':
+        return False
+    else:
+        raise ValueError
+
+
+def print_args(args):
+    """Print the input arguments.
+
+    Parameters
+    ----------
+    args : argparse
+        argparse object to print
+
+    """
+    arg_vars = vars(args)
+    for key in arg_vars.iterkeys():
+        print '== {0:<25} = {1}'.format(key, arg_vars[key])
+
+
+def enum_keys(e):
+    return e.__members__.keys()
+
+
+def enum_parse(s, c):
+    return c[s.upper()]
+
+
+def thread_type(t):
+    if t.lower() == 'max':
+        return multiprocessing.cpu_count()
+    else:
+        return int(t)
+