Initial Commit

5 files changed, 1278 insertions, 0 deletions

diff --git a/plot_llh/LVCPT.py b/plot_llh/LVCPT.py
new file mode 100644
index 0000000..f5b02e6
--- /dev/null
+++ b/plot_llh/LVCPT.py
@@ -0,0 +1,431 @@
+
+# coding: utf-8
+
+## The Theory
+
+import numpy
+import MinimalTools as MT
+import PhysConst as PC
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib.collections as mco
+import matplotlib as mpl
+import scipy.interpolate as interpolate
+import scipy.integrate as integrate
+import scipy as sp
+from numpy import linalg as LA
+
+use_cython = False
+
+if use_cython:
+    import cython.cLVCPT as clv
+
+mpl.rc('font', family='serif', size=20)
+
+pc = PC.PhysicsConstants()
+
+degree = np.pi/180.0
+pc.th12 = 33.36*degree#33.48*degree
+pc.th23 = 45.*degree#42.3*degree
+pc.th13 = 8.66*degree#8.5*degree
+pc.delta1 = 0.0#300.0*degree#306.*degree # perhaps better just set to 0.
+pc.dm21sq = 7.5e-5
+pc.dm31sq = 2.47e-3#2.457e-3
+pc.Refresh()
+
+MT.calcU(pc)
+DELTAM2 = MT.flavorM2(pc)
+
+def Hamiltonian(Enu, LVATERM = np.zeros((3,3), dtype=numpy.complex),
+                LVCTERM = np.zeros((3,3), dtype=numpy.complex)):
+    return DELTAM2/(2.0*Enu) + LVATERM + Enu*LVCTERM 
+
+def OscProbFromMixingMatrix(alpha, beta, MixMatrix):
+    return sum([(np.absolute(MixMatrix[i][alpha])*np.absolute(MixMatrix[i][beta]))**2 for i in range(pc.numneu)] )
+    #return sum([(np.absolute(MixMatrix[i][alpha]))**2*(np.absolute(MixMatrix[i][beta]))**2 for i in range(pc.numneu)] )
+    #prob = 0.0;
+    #for i in range(pc.numneu) :
+    #    prob += (np.absolute(MixMatrix[i][alpha]))**2*(np.absolute(MixMatrix[i][beta]))**2
+    #return prob
+
+def OscProb(alpha, Enu, LVATERM = np.zeros((3,3), dtype=numpy.complex),
+                LVCTERM = np.zeros((3,3), dtype=numpy.complex)):
+    eigvals, eigvec = MT.eigenvectors(Hamiltonian(Enu, LVATERM=LVATERM, LVCTERM=LVCTERM))
+    #print eigvec.dtype
+    if use_cython:
+        return [ clv.OscProbFromMixingMatrix(alpha,beta,eigvec) for beta in range(pc.numneu)]
+    else:
+        return [ OscProbFromMixingMatrix(alpha,beta,eigvec) for beta in range(pc.numneu)]
+
+def FlavorRatio(initial_flavor_ratio, Enu, LVATERM = np.zeros((3,3), dtype=numpy.complex),
+                LVCTERM = np.zeros((3,3), dtype=numpy.complex)):
+    final_flavor_ratio = [0.0]*pc.numneu
+    osc_prob_array = [OscProb(beta,Enu,LVATERM=LVATERM,LVCTERM=LVCTERM) for beta in range(pc.numneu)]
+
+    for alpha in range(pc.numneu):
+        for beta,phi in enumerate(initial_flavor_ratio):
+            final_flavor_ratio[alpha] += osc_prob_array[beta][alpha]*phi
+    return final_flavor_ratio
+
+def RRR(initial_flavor_ratio, Enu, LVATERM = np.zeros((3,3), dtype=numpy.complex),
+                LVCTERM = np.zeros((3,3), dtype=numpy.complex)):
+    ffr = FlavorRatio(initial_flavor_ratio,Enu,LVATERM=LVATERM,LVCTERM=LVCTERM)
+    return ffr[1]/ffr[0]
+def SSS(initial_flavor_ratio, Enu, LVATERM = np.zeros((3,3), dtype=numpy.complex),
+                LVCTERM = np.zeros((3,3), dtype=numpy.complex)):
+    ffr = FlavorRatio(initial_flavor_ratio,Enu,LVATERM=LVATERM,LVCTERM=LVCTERM)
+    return ffr[2]/ffr[1]
+
+def PointToList(p1,p2):
+    return [[p1[0],p2[0]],[p1[1],p2[1]]]
+
+def PointFromFlavor(origin,scale,flavor_ratio_list):
+    nu_e_vec = np.array([1.,0.])*scale
+    nu_mu_vec = np.array([1./2.,np.sqrt(3.)/2.])*scale
+    nu_tau_vec = np.array([-1./2.,np.sqrt(3.)/2.])*scale
+    fpos = origin + flavor_ratio_list[0]*nu_e_vec + flavor_ratio_list[1]*nu_mu_vec
+    return [fpos[0],fpos[1]]
+
+def MakeFlavorTriangle(list_of_flavor_ratios, scale = 8,
+                       p = np.array([0.,0.]), save_file = False, PlotPoints = False, PlotTrayectories = False, figure = None, alpha = 1.0,
+                       filename = "triangle",icolor = "green", icolormap = "Greens", divisions = 5, initial_flavor_ratio = [1,0,0],
+                       term = "a", subdivisions = False, triangle_collection = None, color_scale = "lin", return_fig = True, addtext = "",
+                       add_default_text = True, ilw = 1., slw = 0.75, output_format = "eps", inner_line_color = "k", plot_color_bar = False):
+    # i will be nice ...
+    list_of_flavor_ratios = np.array(list_of_flavor_ratios)
+
+    if figure == None:
+        fig = plt.figure(figsize=(scale,scale), frameon = False)
+    else:
+        fig = figure
+
+    ax = fig.add_axes([0, 0, 1, 1])
+    ax.axis('off')
+
+    # delete extra lines
+    frame = plt.gca()
+    frame.axes.get_xaxis().set_visible(False)
+    frame.axes.get_yaxis().set_visible(False)
+
+    s0 = np.array([1.,0.])*scale
+    s1 = np.array([1./2.,np.sqrt(3.)/2.])*scale
+    s2 = np.array([1./2.,-np.sqrt(3.)/2.])*scale
+
+    # make triangle outer frame
+
+    plt.plot(*PointToList(p, p+s0), color = "k", lw = 4)
+    plt.plot(*PointToList(p, p+s1), color = "k", lw = 2)
+    plt.plot(*PointToList(p+s0, p+s1), color = "k", lw = 2)
+
+    # put outer triangle labels
+
+    # ax.text((p+s0*0.5+s0*0.025)[0], (p+s0*0.5-np.array([0,0.15*scale]))[1],r"$\alpha^{\oplus}_{e}$",
+    #             horizontalalignment="center",fontsize = 50, zorder = 10)
+    # ax.text((p+s1*0.5-0.2*s0)[0], (p+s1*0.5+0.1*s0)[1]+scale*0.1,r"$\alpha^{\oplus}_{\tau}$",
+    #             horizontalalignment="center",fontsize = 50, zorder = 10, rotation = 60.)
+    # ax.text((p+s1*0.5 + 0.7*s0)[0], (p+s1*0.5 + 0.6*s0)[1]+0.05*scale,r"$\alpha^{\oplus}_{\mu}$",
+    #             horizontalalignment="center",fontsize = 50, zorder = 10, rotation = -60)
+
+    ax.text((p+s0*0.5+s0*0.025)[0], (p+s0*0.5-np.array([0,0.15*scale]))[1],r"$f_{e,\oplus}$",
+                horizontalalignment="center",fontsize = 50, zorder = 10)
+    ax.text((p+s1*0.5-0.2*s0)[0], (p+s1*0.5+0.1*s0)[1]+scale*0.1,r"$f_{\tau, \oplus}$",
+                horizontalalignment="center",fontsize = 50, zorder = 10, rotation = 60.)
+    ax.text((p+s1*0.5 + 0.7*s0)[0], (p+s1*0.5 + 0.6*s0)[1]+0.05*scale,r"$f_{\mu, \oplus}$",
+                horizontalalignment="center",fontsize = 50, zorder = 10, rotation = -60)
+
+    # construct triangle grid
+    fsl = 35
+    for i in range(divisions+1):
+        subsize = 1./float(divisions)
+
+        ax.text((p+s0*subsize*float(i))[0], (p+s0*subsize*float(i)-np.array([0,0.05*scale]))[1],str(i*subsize),
+                horizontalalignment="center",fontsize = fsl)
+        plt.plot(*PointToList(p+s0*subsize*float(i), p+s1+s2*subsize*float(i)), color = inner_line_color, lw = ilw, ls = "dashed", zorder = -1)
+        ax.text((p+s1-s1*subsize*float(i)-np.array([0.06*scale,0.0]))[0], (p+s1-s1*subsize*float(i))[1],str(i*subsize),
+                horizontalalignment="center",fontsize = fsl)
+        plt.plot(*PointToList(p+s0*subsize*float(divisions-i), p+s1-s1*subsize*float(i)), color = inner_line_color, lw = ilw, ls = "dashed", zorder = -1)
+
+        ax.text((p+s1+s2*subsize*float(i)+np.array([0.05*scale,0.0]))[0], (p+s1+s2*subsize*float(i))[1],str((divisions-i)*subsize),
+                horizontalalignment="center",fontsize = fsl)
+        plt.plot(*PointToList(p+s1*subsize*float(divisions-i), p+s1+s2*subsize*float(i)), color = inner_line_color, lw = ilw, ls = "dashed", zorder = -1)
+
+        if subdivisions and i < divisions:
+            plt.plot(*PointToList(p+s0*subsize*float(i+0.5), p+s1+s2*subsize*float(i+0.5)), color = inner_line_color, lw = slw, ls = "dotted", zorder = -1)
+        if subdivisions and i > 0:
+            plt.plot(*PointToList(p+s0*subsize*float(divisions-(i-0.5)), p+s1-s1*subsize*float(i-0.5)), color = inner_line_color, lw = slw, ls = "dotted", zorder = -1)
+            plt.plot(*PointToList(p+s1*subsize*float(divisions-(i-0.5)), p+s1+s2*subsize*float(i-0.5)), color = inner_line_color, lw = slw, ls = "dotted", zorder = -1)
+
+
+    # plot triangle collection
+    if (triangle_collection != None):
+        # get total number of points
+        total_points = float(sum([ triangle.number_of_points for triangle in triangle_collection]))
+        max_points = float(max([ triangle.number_of_points for triangle in triangle_collection]))
+        color_map = plt.get_cmap(icolormap)
+        for triangle in triangle_collection:
+            if triangle.number_of_points > 0:
+                xx,yy = zip(*triangle.coordinates)
+                if color_scale == "lin":
+                    plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(0.75), alpha = np.sqrt(float(triangle.number_of_points)/max_points))
+                    #plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(float(triangle.number_of_points)/max_points), alpha = alpha)
+                elif color_scale == "log":
+                    plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(0.75), alpha = (np.log10(float(triangle.number_of_points))/np.log10(max_points)))
+                    #plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(0.7), alpha = (np.log10(float(triangle.number_of_points))/np.log10(max_points))**(2./3.))
+                    #plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(np.log10(float(triangle.number_of_points))/np.log10(max_points)), alpha = alpha)
+                    #plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(np.log10(float(triangle.number_of_points))))
+                else : 
+                    raise NameError('Error. Love CA.')
+
+        if(plot_color_bar):
+            # the color bar magic
+            # location set on 0 to 1 scales.
+            left = 0.1
+            bottom = -0.25
+            width = 0.8
+            height = 0.025
+            cbaxes = fig.add_axes([left,bottom,width,height])
+            if color_scale == "lin":
+                norm = mpl.colors.Normalize(vmin = 0., vmax = max_points)
+            elif color_scale == "log":
+                norm = mpl.colors.Normalize(vmin = 0., vmax = 1.0)
+            else :
+                raise NameError('Error. Love CA.')
+            mpl.rcParams.update({'font.size': 10})
+            triangle_colorbar = mpl.colorbar.ColorbarBase(cbaxes, cmap = color_map, norm = norm,
+                                                          orientation = "horizontal", spacing = "proportional",
+							 # )
+                                                          format ='%.0e')
+            cbaxes.set_xlabel("Likelihood", fontsize = 12)
+
+    # plot flavor ratio points
+    if PlotTrayectories :
+        if len(list_of_flavor_ratios.shape) == 3 :
+            for flavor_ratio_l in list_of_flavor_ratios:
+                flv_ratio_coords = map(lambda f:PointFromFlavor(p,scale,np.array(f)),flavor_ratio_l)
+                xc, yc = zip(*flv_ratio_coords)
+                plt.plot(xc,yc, color = icolor,
+                        ms = 10, linewidth = 4, zorder = 0)
+        elif len(list_of_flavor_ratios.shape) == 2 :
+            flv_ratio_coords = map(lambda f:PointFromFlavor(p,scale,np.array(f)),list_of_flavor_ratios)
+            xc, yc = zip(*flv_ratio_coords)
+
+            plt.plot(xc,yc, color = icolor,
+                    ms = 10, linewidth = 4, zorder = 0)
+        else:
+            raise NameError('Check your input flavor list array and the joined flag. Love CA.')
+    elif PlotPoints:
+        if len(list_of_flavor_ratios.shape) !=2 :
+            print len(list_of_flavor_ratios.shape)
+            raise NameError('Check your input flavor list array and the joined flag. Love CA.')
+        for i,flavor_ratio in enumerate(list_of_flavor_ratios):
+            if len(icolor) != len(list_of_flavor_ratios):
+                icolor_ = icolor
+            else:
+                icolor_ = icolor[i]
+            plt.plot(*PointFromFlavor(p,scale,np.array(flavor_ratio)), color = icolor_,
+                     marker = 'o', ms = 10, linewidth = 0,
+                     markeredgecolor=None,markeredgewidth=0.1, zorder = 1000)
+
+    # put back color scale axis
+    if add_default_text:
+        ax.text((s0/5.+0.9*s1)[0],(s0/5.+0.9*s1)[1],
+                "LV "+term+"-term scan with\n $\ \phi_e:\ \phi_\\mu:\ \phi_\\tau = "+str(initial_flavor_ratio[0])+":\ "+str(initial_flavor_ratio[1])+":\ "+str(initial_flavor_ratio[2])+"$"+" \n "+addtext,
+                fontsize = 20)
+
+    if(save_file):
+        # save figure
+        plt.savefig("./plots/"+filename+"."+output_format, dpi = 600, bbox_inches='tight')
+    else:
+        # show figure
+        plt.show()
+    if return_fig:
+        return fig
+
+
+def s_bario(p,p0,p1,p2):
+    return (p0[1]*p2[0] - p0[0]*p2[1] + (p2[1] - p0[1])*p[0] + (p0[0] - p2[0])*p[1])
+
+def t_bario(p,p0,p1,p2):
+    return (p0[0]*p1[1] - p0[1]*p1[0] + (p0[1] - p1[1])*p[0] + (p1[0] - p0[0])*p[1])
+
+def IsInTriangle(p,p0,p1,p2,area):
+    s = s_bario(p,p0,p1,p2)
+    t = t_bario(p,p0,p1,p2)
+    #print s,t,2.0*area - s - t
+    return  s >= -1.e-15 and  t >= -1.0e-15 and s+t <= 2.0*area
+
+
+class Triangle:
+    coordinates = []
+    area = 0.0
+    number_of_points = 0.0
+    n_t = 0
+    i = 0 
+    j = 0
+    orientation = ""
+    
+    def IsPointIn(self,point):
+        p0 = self.coordinates[0]
+        p1 = self.coordinates[1]
+        p2 = self.coordinates[2]
+        return IsInTriangle(point,p0,p1,p2,self.area)
+
+
+def GenerateTriangles(scale, divisions, p = np.array([0.,0.])):
+    s0 = np.array([1.,0.])*scale/float(divisions)
+    s1 = np.array([1./2.,np.sqrt(3.)/2.])*scale/float(divisions)
+    s2 = np.array([1./2.,-np.sqrt(3.)/2.])*scale/float(divisions)
+    
+    area = np.sqrt(3)*(LA.norm(s0)/2.0)**2
+    
+    n_t = 0
+    
+    triangle_collection = []
+    for i in range(divisions):
+        for j in range(divisions-i):
+            lower_triangle = Triangle()
+
+            p0_l = p + i*s0 + j*s1 
+            p1_l = p0_l + s0
+            p2_l = p0_l + s1
+            
+            lower_triangle.coordinates = [p0_l,p1_l,p2_l]
+            lower_triangle.n_t = n_t
+            lower_triangle.i = i
+            lower_triangle.j = j
+            lower_triangle.orientation = "L"
+            lower_triangle.area = area
+            
+            n_t += 1
+            # append to triangle collection
+            triangle_collection.append(lower_triangle)
+            
+            upper_triangle = Triangle()
+
+            p0_u = p2_l
+            p1_u = p1_l
+            p2_u = p1_l + s1
+            
+            upper_triangle.coordinates = [p0_u,p1_u,p2_u]
+            upper_triangle.n_t = n_t
+            upper_triangle.i = i
+            upper_triangle.j = j
+            upper_triangle.orientation = "U"
+            upper_triangle.area = area
+            
+            n_t += 1
+            # append to triangle collection
+            triangle_collection.append(upper_triangle)
+    return triangle_collection
+
+def AddPointToTriangleCollectionLegacy(flavor_ratio, triangle_collection,
+                  p = np.array([0.,0.]), scale = 8, divisions = 10):
+    point = PointFromFlavor(p,scale,np.array(flavor_ratio))
+    electron = 0; tau = 2;
+    # the silly way
+    for triangle in triangle_collection:
+        if(triangle.IsPointIn(point)):
+            triangle.number_of_points += 1.
+
+def AddPointToTriangleCollection(flavor_ratio, triangle_collection,
+                  p = np.array([0.,0.]), scale = 8, divisions = 10):
+    point = PointFromFlavor(p,scale,np.array(flavor_ratio))
+    electron = 0; muon = 1; tau = 2;
+    # the smart way
+    u_i = int(flavor_ratio[electron]*float(divisions))
+    u_j = int(flavor_ratio[muon]*float(divisions))
+    index = u_i*(2*divisions-u_i+1) + 2*u_j
+    if triangle_collection[index].IsPointIn(point):
+        triangle_collection[index].number_of_points += 1.
+    else:
+        triangle_collection[index+1].number_of_points += 1.
+# legacy
+    #elif triangle_collection[index+1].IsPointIn(point):
+    #    triangle_collection[index+1].number_of_points += 1.
+    #else:
+    #    print "Math error."
+    #    print point, "\n",u_i, u_j, "\n", triangle_collection[index].coordinates, "\n", triangle_collection[index+1].coordinates
+    #    raise NameError("Error triangle location math")
+
+class AnarchySampling:
+    def __init__(self, n_sample, LV_scale_1, LV_scale_2, term):
+        self.n_sample = n_sample
+        self.th12_sample = np.arcsin(np.sqrt(np.random.uniform(0.,1., size=n_sample)))
+        self.th13_sample = np.arccos(np.sqrt(np.sqrt(np.random.uniform(0.,1., size=n_sample))))
+        self.th23_sample = np.arcsin(np.sqrt(np.random.uniform(0.,1., size=n_sample)))
+        self.delta_sample = np.random.uniform(0.,2.*np.pi, size=n_sample)
+
+        self.LV_scale_1 = LV_scale_1
+        self.LV_scale_2 = LV_scale_2
+
+        self.term = term
+
+def GenerateFlavorRatioPoints(Initial_Flavor_Ratio, SamplingObject, gamma = 2.0,
+                              Log10Emax = 7., Log10Emin = 4.0, Epoints = 30,
+                              save_list = False, save_avg = True):
+    flavor_tray_list = []
+    flavor_avg_list = []
+
+    # energy things
+
+    Erange = np.logspace(Log10Emin,Log10Emax,Epoints) # in GeV
+    Emin = Erange[0]
+    Emax = Erange[-1]
+
+    if gamma == 1 or gamma == 1.0:
+        spectral_normalization = np.log(Emax)-np.log(Emin)
+    else:
+        spectral_normalization = (Emax**(1.-gamma) - Emin**(1.-gamma))/(1.-gamma)
+
+    spectral_function = lambda Enu: Enu**(-gamma)/spectral_normalization
+
+    # loop over random parameters
+    for i in range(SamplingObject.n_sample):
+        lv_term = MT.LVP()
+
+        lv_term.th12 = SamplingObject.th12_sample[i]
+        lv_term.th13 = SamplingObject.th13_sample[i]
+        lv_term.th23 = SamplingObject.th23_sample[i]
+        lv_term.delta1 = SamplingObject.delta_sample[i]
+
+        lv_term.LVS21 = SamplingObject.LV_scale_1
+        lv_term.LVS31 = SamplingObject.LV_scale_2
+
+        lv_term.Refresh()
+
+        LVTERM = MT.LVTerm(lv_term);
+
+        if SamplingObject.term == "a":
+            flavor_ratio_list = np.array(map(lambda Enu : FlavorRatio(Initial_Flavor_Ratio, Enu*pc.GeV, LVATERM = LVTERM), Erange))
+        elif SamplingObject.term == "c":
+            flavor_ratio_list = np.array(map(lambda Enu : FlavorRatio(Initial_Flavor_Ratio, Enu*pc.GeV, LVCTERM = LVTERM), Erange))
+        else :
+            raise NameError('Only a or c term.'+ str(term))
+
+        if save_avg: 
+            if Epoints != 1:
+                flavor_avg = [0.]*lv_term.numneu
+                for alpha in range(lv_term.numneu):
+                    #inter = interpolate.interp1d(Erange,flavor_ratio_list[:,alpha])
+                    inter = interpolate.UnivariateSpline(Erange,flavor_ratio_list[:,alpha])
+                    flavor_avg[alpha] = integrate.quad(lambda Enu : inter(Enu)*spectral_function(Enu),
+                                                       Emin,Emax, limit = 75, epsrel = 1e-2, epsabs = 1.0e-2)[0]
+                    #flavor_avg[alpha] = integrate.quadrature(lambda Enu : inter(Enu)*spectral_function(Enu),
+                    #                                         Emin,Emax, maxiter = 75, rtol = 1e-3, tol = 1.e-3)[0]
+                flavor_avg_list.append(flavor_avg)
+            else:
+                flavor_avg = flavor_ratio_list[0]
+                flavor_avg_list.append(flavor_avg)
+
+        if save_list:
+            flavor_tray_list.append(flavor_ratio_list)
+
+    if save_list and save_avg:
+        return flavor_tray_list, flavor_avg_list
+    elif save_list:
+        return flavor_tray_list
+    elif save_avg:
+        return flavor_avg_list
+    else :
+        print "Math is broken."
+        return None
diff --git a/plot_llh/MinimalTools.py b/plot_llh/MinimalTools.py
new file mode 100644
index 0000000..4ae8360
--- /dev/null
+++ b/plot_llh/MinimalTools.py
@@ -0,0 +1,166 @@
+import numpy as np
+from PhysConst import PhysicsConstants
+
+def eigenvectors(M):
+        """ Calculates the eigenvectors and eigenvalues ordered by eigenvalue size
+
+        @type  M     : matrix
+        @param M     : matrix M
+
+        @rtype       : list
+        @return      : [eigenvalues list, eigenvector list]
+        """
+        D,V = np.linalg.eig(M)
+        DV = []
+        VT = V.T
+        for i,eigenvalue in enumerate(D):
+            DV.append([eigenvalue,VT[i]]) 
+
+        DV = sorted(DV,key = lambda x : x[0].real)#np.abs(x[0].real))
+
+        V2 = []
+        D2 = []
+        for e in DV:
+            V2.append(e[1])
+            D2.append(e[0])
+        return D2,V2
+
+# General Rotation Matrix
+def R(i,j,cp,param):
+        """ Rotation Matrix
+        Calculates the R_ij rotations. Also incorporates CP-phases when necesary.
+        @type   i       :   int
+        @param  i       :   i-column.
+        @type   j       :   int
+        @param  j       :   j-row.
+        @type   cp      :   int
+        @param  cp      :   if cp = 0 : no CP-phase. else CP-phase = CP_array[cp]
+
+        @rtype          :   numpy.array
+        @return         :   returns the R_ij rotation matrix.
+        """
+        # if cp = 0 -> no complex phase
+        # R_ij, i<j
+        if(j<i):
+            # no funny business
+            l = i
+            i = j
+            j = l
+
+        # rotation matrix generator
+        R = np.zeros([param.numneu,param.numneu],complex)
+        # diagonal terms
+        for k in np.arange(0,param.numneu,1):
+            if(k != i-1 and k != j-1):
+                R[k,k] = 1.0
+            else :
+                R[k,k] = param.c[i,j]
+        # non-diagonal terms
+        if(cp != 0):
+            sd = np.sin(param.dcp[cp])
+            cd = np.cos(param.dcp[cp])
+            faseCP = complex(cd,sd)
+        else:
+            faseCP = complex(1.0,0.0)
+
+        R[i-1,j-1] = param.s[i,j]*faseCP.conjugate()
+        R[j-1,i-1] = -param.s[i,j]*faseCP
+        return R
+
+def calcU(param):
+        """ Defining the mixing matrix parametrization.
+        @type   param   :   PhysicsConstants
+        @param  param   :   set of physical parameters to be used.
+
+        @rtype          :   None
+        @return         :   Sets mixing matrix.
+        """
+        if(param.numneu == 2):
+	        return self.R(1,2,0,param)
+        elif(param.numneu == 3):
+            return np.dot(R(2,3,0,param),np.dot(R(1,3,1,param),R(1,2,0,param)))
+        elif(param.numneu == 4):
+            return np.dot(R(3,4,0,param),np.dot(R(2,4,2,param),np.dot(R(1,4,0,param),np.dot(R(2,3,0,param),np.dot(R(1,3,1,param),R(1,2,0,param))))))
+        elif(param.numneu == 5):
+            return np.dot(R(4,5,0,param),np.dot(R(3,5,0,param),np.dot(R(2,5,0,param),np.dot(R(1,5,3,param),np.dot(R(3,4,0,param),np.dot(R(2,4,0,param),np.dot(R(1,4,2,param),np.dot(R(2,3,0,param),np.dot(R(1,3,1,param),R(1,2,0,param))))))))))
+        elif(param.numneu == 6):
+            # 3+3 twin-sterile-neutrino model
+            return np.dot(R(3,6,0,param),np.dot(R(2,5,0,param),np.dot(R(1,4,0,param),np.dot(R(2,3,0,param),np.dot(R(1,3,1,param),R(1,2,0,param))))))
+        else:
+            print "Sorry, too many neutrinos. Not yet implemented! =(."
+            quit()
+
+        # antineutrino case
+        if param.neutype == "antineutrino" :
+            return self.U.conjugate()
+
+
+
+def massM2(param):
+    """ Mass term in the neutrino mass basis.
+
+    @type   param   :   PhysicsConstants
+    @param  param   :   set of physical parameters to be used.
+
+    @rtype          :   numpy array
+    @return         :   mass matrix in mass basis.
+    """
+    M2 = np.zeros([param.numneu,param.numneu],complex)
+    for k in np.arange(1,param.numneu,1):
+        M2[k,k] = param.dmsq[k+1]
+    return M2
+
+def flavorM2(param):
+    """ Mass term in the neutrino flavor basis.
+
+    @type   param   :   PhysicsConstants
+    @param  param   :   set of physical parameters to be used.
+
+    @rtype          :   numpy array
+    @return         :   mass matrix in flavor basis.
+    """
+    U = calcU(param)
+    return np.dot(U,np.dot(massM2(param),U.conjugate().T))
+
+class LVP(PhysicsConstants):
+    def __init__(self):
+        super(LVP,self).__init__()
+        self.th12 = 0.0
+        self.th13 = 0.0
+        self.th23 = 0.0
+        self.delta1 = 0.0
+        self.deltaCP = 0.0
+        super(LVP,self).Refresh()
+
+        # LVS
+        self.LVS21 = 0.0	                #
+        self.LVS31 = 0.0                    #
+        self.LVS41 = 0.0                    #
+        self.LVS51 = 0.0                    #
+        self.LVS61 = 0.0                    #
+        # SQUARED MASS DIFFERENCE MATRIX
+        self.LVS = np.zeros([self.numneumax+2],float)
+        self.LVS[2] = self.LVS21
+        self.LVS[3] = self.LVS31
+        self.LVS[4] = self.LVS41
+        self.LVS[5] = self.LVS51
+        self.LVS[6] = self.LVS61
+
+    def Refresh(self):
+        super(LVP,self).Refresh()
+        LVS = self.LVS
+        LVS[2] = self.LVS21
+        LVS[3] = self.LVS31
+        LVS[4] = self.LVS41
+        LVS[5] = self.LVS51
+        LVS[6] = self.LVS61
+
+def DiagonalMatrixLV(param):
+    DD = np.zeros([param.numneu,param.numneu],complex)
+    for k in np.arange(1,param.numneu,1):
+        DD[k,k] = param.LVS[k+1]
+    return DD
+
+def LVTerm(LVparam):
+    U = calcU(LVparam)
+    return np.dot(U,np.dot(DiagonalMatrixLV(LVparam),U.conjugate().T))
diff --git a/plot_llh/PhysConst.py b/plot_llh/PhysConst.py
new file mode 100644
index 0000000..89a0be0
--- /dev/null
+++ b/plot_llh/PhysConst.py
@@ -0,0 +1,390 @@
+"""
+Author  : C.A. Arguelles
+Date    : 10/MAY/2011
+
+Contains Physics constants and global variables.
+
+Log :
+- Modified on 23/ABR/2012 by C.Arguelles
+    + Changed the definition of PhysicsConstants to
+    include an __init__ to separate the class global
+    properties from its instances.
+"""
+
+# python standard modules
+import numpy as np
+
+class PhysicsConstants(object):
+    
+    def __init__(self):
+        ## PHYSICS CONSTANTS
+        #===========================================================================
+        # NAME
+        #===========================================================================
+        
+        self.name = "STD"                    # Default values
+        self.linestyle = "solid"             # Default linestyle in plots
+        self.markerstyle = "*"               # Default marker style
+        self.colorstyle = "red"              # Default color style
+        self.savefilename = "output.dat"     # Default color style
+        
+        #===============================================================================
+        # ## MATH
+        #===============================================================================
+        self.PI=3.14159265		            # Pi
+        self.PIby2=1.5707963268	            # Pi/2
+        self.sqr2=1.4142135624	            # Sqrt[2]
+        self.ln2 = np.log(2.0)
+        
+        #===============================================================================
+        # ## EARTH 
+        #===============================================================================
+        self.EARTHRADIUS = 6371.0	        # [km] Earth radius
+        #===============================================================================
+        # ## SUN 
+        #===============================================================================
+        self.SUNRADIUS = 109*self.EARTHRADIUS     # [km] Sun radius 
+        
+        #===============================================================================
+        # # PHYSICAL CONSTANTS
+        #===============================================================================
+        self.GF = 1.16639e-23	            # [eV^-2] Fermi Constant 
+        self.Na = 6.0221415e+23 		        # [mol cm^-3] Avogadro Number
+        self.sw_sq = 0.2312                  # [dimensionless] sin(th_weinberg) ^2
+        self.G  = 6.67300e-11                # [m^3 kg^-1 s^-2]
+        self.alpha = 1.0/137.0               # [dimensionless] fine-structure constant 
+        
+        #===============================================================================
+        # ## UNIT CONVERSION FACTORS
+        #===============================================================================
+        # Energy
+        self.TeV = 1.0e12                    # [eV/TeV]
+        self.GeV = 1.0e9                     # [eV/GeV]
+        self.MeV = 1.0e6                     # [eV/MeV]
+        self.keV = 1.0e3                     # [eV/keV]
+        self.Joule = 1/1.60225e-19           # [eV/J]
+        # Mass
+        self.kg = 5.62e35                    # [eV/kg]
+        self.gr = 1e-3*self.kg               # [eV/g] 
+        # Time
+        self.sec = 1.523e15                  # [eV^-1/s]
+        self.hour = 3600.0*self.sec          # [eV^-1/h]
+        self.day = 24.0*self.hour            # [eV^-1/d]
+        self.year = 365.0*self.day           # [eV^-1/yr]
+        self.yearstosec = self.sec/self.year # [s/yr]
+        # Distance
+        self.meter = 5.076e6                 # [eV^-1/m]
+        self.cm = 1.0e-2*self.meter          # [eV^-1/cm]
+        self.km = 1.0e3*self.meter           # [eV^-1/km]
+        self.fermi = 1.0e-15*self.meter      # [eV^-1/fm]
+        self.angstrom = 1.0e-10*self.meter   # [eV^-1/A]
+        self.AU = 149.60e9*self.meter        # [eV^-1/AU]
+        self.parsec = 3.08568025e16*self.meter# [eV^-1/parsec]
+        # Integrated Luminocity # review
+        self.picobarn = 1.0e-36*self.cm**2   # [eV^-2/pb]
+        self.femtobarn = 1.0e-39*self.cm**2  # [eV^-2/fb]
+        # Presure
+        self.Pascal = self.Joule/self.meter**3 # [eV^4/Pa]
+        self.hPascal = 100.0*self.Pascal     # [eV^4/hPa]
+        self.atm = 101325.0*self.Pascal      # [eV^4/atm]
+        self.psi = 6893.0*self.Pascal        # [eV^4/psi]
+        # Temperature
+        self.kelvin = 1/1.1604505e4          # [eV/K]
+        # Angle
+        self.degree = self.PI/180.0          # [rad/degree]
+        # magnetic field
+        self.T = 0.000692445                 # [eV^2/T]
+        
+        # old notation
+        self.cm3toev3 = 7.68351405e-15       # cm^3-> ev^3
+        self.KmtoEv =5.0677288532e+9         # km -> eV
+        self.yearstosec = 31536.0e3          # years -> sec
+        
+        #===============================================================================
+        # ## NEUTRINO OSCILLATION PARAMETERS ##
+        #===============================================================================
+        
+        self.numneu = 3                      # number of neutrinos
+        self.numneumax = 6                   # maximum neutrino number
+        self.neutype = 'neutrino'
+        #neutype = 'antineutrino'
+        
+        # values updated according to 1209.3023 Table 1 FreeFluxes + RSBL
+        
+        # MIXING ANGLES
+       
+        self.th12 = 0.579639
+        self.th13 = 0.150098
+        self.th23 = self.PIby2/2.0
+        self.th14 = 0.0
+        self.th24 = 0.0
+        self.th34 = 0.0
+        self.th15 = 0.0
+        self.th25 = 0.0
+        self.th35 = 0.0
+        self.th45 = 0.0
+        self.th16 = 0.0
+        self.th26 = 0.0
+        self.th36 = 0.0
+        self.th46 = 0.0
+        self.th56 = 0.0
+        
+        # mixing angles matrix array
+        self.th = np.zeros([self.numneumax+1,self.numneumax+1],float)
+        self.th[1,2] = self.th12
+        self.th[1,3] = self.th13
+        self.th[2,3] = self.th23
+        self.th[1,4] = self.th14
+        self.th[2,4] = self.th24
+        self.th[3,4] = self.th34
+        self.th[1,5] = self.th15
+        self.th[2,5] = self.th25
+        self.th[3,5] = self.th35
+        self.th[4,5] = self.th45
+        self.th[1,6] = self.th16
+        self.th[2,6] = self.th26
+        self.th[3,6] = self.th36
+        self.th[4,6] = self.th46
+        self.th[5,6] = self.th56
+        
+        self.s12 = np.sin(self.th12)
+        self.c12 = np.cos(self.th12)
+        self.s13 = np.sin(self.th13)
+        self.c13 = np.cos(self.th13)
+        self.s23 = np.sin(self.th23)
+        self.c23 = np.cos(self.th23)
+        self.s14 = np.sin(self.th14)
+        self.c14 = np.cos(self.th14)
+        self.s24 = np.sin(self.th24)
+        self.c24 = np.cos(self.th24)
+        self.s34 = np.sin(self.th34)
+        self.c34 = np.cos(self.th34)
+        self.s15 = np.sin(self.th15)
+        self.c15 = np.cos(self.th15)
+        self.s25 = np.sin(self.th25)
+        self.c25 = np.cos(self.th25)
+        self.s35 = np.sin(self.th35)
+        self.c35 = np.cos(self.th35)
+        self.s45 = np.sin(self.th45)
+        self.c45 = np.cos(self.th45)
+        self.s16 = np.sin(self.th16)
+        self.c16 = np.cos(self.th16)
+        self.s26 = np.sin(self.th26)
+        self.c26 = np.cos(self.th26)
+        self.s36 = np.sin(self.th36)
+        self.c36 = np.cos(self.th36)
+        self.s46 = np.sin(self.th46)
+        self.c46 = np.cos(self.th46)
+        self.s56 = np.sin(self.th56)
+        self.c56 = np.cos(self.th56)    
+        
+        # cos(th_ij) matrix array
+        self.c = np.zeros([self.numneumax+1,self.numneumax+1],float)
+        self.c[1,2] = self.c12
+        self.c[1,3] = self.c13
+        self.c[1,4] = self.c14
+        self.c[2,3] = self.c23
+        self.c[2,4] = self.c24
+        self.c[3,4] = self.c34
+        self.c[1,5] = self.c15
+        self.c[2,5] = self.c25
+        self.c[3,5] = self.c35
+        self.c[4,5] = self.c45
+        self.c[1,6] = self.c16
+        self.c[2,6] = self.c26
+        self.c[3,6] = self.c36
+        self.c[4,6] = self.c46
+        self.c[5,6] = self.c56    
+        
+        # sin(th_ij) matrix array
+        self.s = np.zeros([self.numneumax+1,self.numneumax+1],float)
+        self.s[1,2] = self.s12
+        self.s[1,3] = self.s13
+        self.s[1,4] = self.s14
+        self.s[2,3] = self.s23
+        self.s[2,4] = self.s24
+        self.s[3,4] = self.s34
+        self.s[1,5] = self.s15
+        self.s[2,5] = self.s25
+        self.s[3,5] = self.s35
+        self.s[4,5] = self.s45
+        self.s[1,6] = self.s16
+        self.s[2,6] = self.s26
+        self.s[3,6] = self.s36
+        self.s[4,6] = self.s46
+        self.s[5,6] = self.s56    
+        
+        # CP PHASES
+        #self.delta21=3.3e-5
+        #self.delta32=3.1e-3
+        #self.delta31=self.delta32+self.delta21
+        #self.deltaCP=self.PIby2
+        
+        # CP Phases
+        self.deltaCP = 5.235987
+        self.delta1 = self.deltaCP
+        self.delta2 = 0.0
+        self.delta3 = 0.0
+        
+        # d-CP phases
+        self.dcp = np.zeros([self.numneumax-2+1],complex)
+        self.dcp[0] = 1.0
+        self.dcp[1] = self.delta1
+        self.dcp[2] = self.delta2
+        self.dcp[3] = self.delta3
+        
+        # SQUARED MASS DIFFERENCE
+        self.dm21sq = 7.50e-5	             # [eV^2]
+        self.dm31sq = 2.47e-3	             # [eV^2]
+        self.dm32sq = -2.43e-3               # [eV^2]
+        # STERILE 
+        self.dm41sq = 0.0                    # [eV^2]
+        self.dm51sq = 0.0                    # [eV^2]
+        self.dm61sq = 0.0                    # [eV^2]
+        # SQUARED MASS DIFFERENCE MATRIX
+        self.dmsq = np.zeros([self.numneumax+2],float)
+        self.dmsq[2] = self.dm21sq
+        self.dmsq[3] = self.dm31sq
+        self.dmsq[4] = self.dm41sq
+        self.dmsq[5] = self.dm51sq
+        self.dmsq[6] = self.dm61sq    
+        
+        self.dm2 = np.zeros([self.numneumax+1,self.numneumax+1],float)
+        self.dm2[1,2] = self.dm21sq
+        self.dm2[1,3] = self.dm31sq
+        self.dm2[2,3] = self.dm32sq
+        self.dm2[1,4] = self.dm41sq
+        self.dm2[1,5] = self.dm51sq
+        self.dm2[1,6] = self.dm61sq    
+        
+        # MIXING MATRIX
+        self.U = None
+        
+        #===============================================================================
+        # # PARTICLE MASSES
+        #===============================================================================
+        self.muon_mass = 0.10565 	     # [GeV]
+        self.neutron_mass = 0.939565         # [GeV]
+        self.proton_mass = 0.938272          # [GeV]
+        self.electron_mass = 0.510998910e-3  # [GeV]
+        
+        self.atomic_mass_unit = 1.660538e-24 # [g]
+        
+        ## names
+        self.electron = 0
+        self.muon = 1
+        self.tau = 2
+        self.sterile1 = 3
+        self.sterile2 = 4
+        self.sterile3 = 5
+    
+    #===============================================================================
+    # REFRESH
+    #===============================================================================
+    
+    def Refresh(self):
+        # Refresh angles
+        self.s12 = np.sin(self.th12)
+        self.c12 = np.cos(self.th12)
+        self.s13 = np.sin(self.th13)
+        self.c13 = np.cos(self.th13)
+        self.s23 = np.sin(self.th23)
+        self.c23 = np.cos(self.th23)
+        self.s14 = np.sin(self.th14)
+        self.c14 = np.cos(self.th14)
+        self.s24 = np.sin(self.th24)
+        self.c24 = np.cos(self.th24)
+        self.s34 = np.sin(self.th34)
+        self.c34 = np.cos(self.th34)
+        self.s15 = np.sin(self.th15)
+        self.c15 = np.cos(self.th15)
+        self.s25 = np.sin(self.th25)
+        self.c25 = np.cos(self.th25)
+        self.s35 = np.sin(self.th35)
+        self.c35 = np.cos(self.th35)
+        self.s45 = np.sin(self.th45)
+        self.c45 = np.cos(self.th45)
+        self.s16 = np.sin(self.th16)
+        self.c16 = np.cos(self.th16)
+        self.s26 = np.sin(self.th26)
+        self.c26 = np.cos(self.th26)
+        self.s36 = np.sin(self.th36)
+        self.c36 = np.cos(self.th36)
+        self.s46 = np.sin(self.th46)
+        self.c46 = np.cos(self.th46)
+        self.s56 = np.sin(self.th56)
+        self.c56 = np.cos(self.th56)                
+        
+        th = self.th
+        th[1,2] = self.th12
+        th[1,3] = self.th13
+        th[2,3] = self.th23
+        th[1,4] = self.th14
+        th[2,4] = self.th24
+        th[3,4] = self.th34
+        th[1,5] = self.th15
+        th[2,5] = self.th25
+        th[3,5] = self.th35
+        th[4,5] = self.th45
+        th[1,6] = self.th16
+        th[2,6] = self.th26
+        th[3,6] = self.th36
+        th[4,6] = self.th46
+        th[5,6] = self.th56        
+        # Refresh cos(th_ij)
+        c = self.c
+        c[1,2] = self.c12
+        c[1,3] = self.c13
+        c[1,4] = self.c14
+        c[2,3] = self.c23
+        c[2,4] = self.c24
+        c[3,4] = self.c34
+        c[1,5] = self.c15
+        c[2,5] = self.c25
+        c[3,5] = self.c35
+        c[4,5] = self.c45
+        c[1,6] = self.c16
+        c[2,6] = self.c26
+        c[3,6] = self.c36
+        c[4,6] = self.c46
+        c[5,6] = self.c56        
+        # Refresh sin(th_ij)
+        s = self.s
+        self.s[1,2] = self.s12
+        self.s[1,3] = self.s13
+        self.s[1,4] = self.s14
+        self.s[2,3] = self.s23
+        self.s[2,4] = self.s24
+        self.s[3,4] = self.s34
+        self.s[1,5] = self.s15
+        self.s[2,5] = self.s25
+        self.s[3,5] = self.s35
+        self.s[4,5] = self.s45
+        self.s[1,6] = self.s16
+        self.s[2,6] = self.s26
+        self.s[3,6] = self.s36
+        self.s[4,6] = self.s46
+        self.s[5,6] = self.s56                
+        # Refresh CP-Phases
+        dcp = self.dcp
+        dcp[0] = 1.0
+        dcp[1] = self.delta1
+        dcp[2] = self.delta2
+        dcp[3] = self.delta3
+        #dcp[4] = self.delta2        
+        # Refresh Square mass differences
+        dmsq = self.dmsq
+        dmsq[2] = self.dm21sq
+        dmsq[3] = self.dm31sq
+        dmsq[4] = self.dm41sq
+        dmsq[5] = self.dm51sq
+        dmsq[6] = self.dm61sq        
+        
+        dm2 = self.dm2
+        dm2[1,2] = self.dm21sq
+        dm2[1,3] = self.dm31sq
+        dm2[2,3] = self.dm32sq
+        dm2[1,4] = self.dm41sq
+        dm2[1,5] = self.dm51sq
+        dm2[1,6] = self.dm61sq
+        
diff --git a/plot_llh/make_plots.py b/plot_llh/make_plots.py
new file mode 100644
index 0000000..9216396
--- /dev/null
+++ b/plot_llh/make_plots.py
@@ -0,0 +1,117 @@
+#!/usr/bin/python
+
+import sys
+sys.path.append('/Users/Teps/Theory')
+#import header as h
+#sys.path.append('/Users/Teps/Theory/HESE')
+#import anarchy_header as ah
+sys.path.append('/Users/Teps/Theory/HESE/Carlos')
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+from matplotlib import rc, rcParams
+import MinimalTools as MT
+import PhysConst as PC
+import LVCPT as LVCPT
+import numpy as np
+
+import sys,os
+
+rc('text', usetex=True)
+rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
+cols = ['#29A2C6','#FF6D31','#FFCB18','#73B66B','#EF597B', '#333333']
+
+font = {'family' : 'serif',
+        'weight' : 'bold',
+        'size'   : 18}
+
+if len(sys.argv)< 2:
+    print sys.argv
+    print "Usage: make_plots.py input_filepath."
+    exit(1)
+
+#colors_schemes = ["Greens","Reds","Blues","PuRd","summer"]
+colors_schemes = ["Greens","Reds","Blues","spring","summer"]
+#colors_schemes = ["Greens","Reds","Blues","cool","summer"]
+#colors_schemes = ["Greens","Reds","Blues","PuRd","summer"]
+#colors_schemes = ["Blues","Greens","Reds","PuRd","summer"]
+#colors_schemes = ["Greys","Greens","Reds","PuRd","summer"]
+#colors_schemes = ["Greys","Greys","Greys","Greys","summer"]
+#colors_schemes = ["PuRd","summer"]
+output_format = "pdf"
+
+# if True then will plot all files in the same figure
+use_same_canvas = True 
+figure = None
+
+for i in range(len(sys.argv)-1):
+    infile = sys.argv[i+1]
+    print "Load data: " + str(infile)
+    if infile[-3:] == 'txt':
+        flavor_list = np.genfromtxt(infile)
+    else:
+        flavor_list = np.load(infile)
+    if len(flavor_list[~np.isfinite(flavor_list)]) != 0:
+        fl = []
+        for x in flavor_list:
+            if np.sum(~np.isfinite(x)) == 0:
+                fl.append(x.tolist())
+        flavor_list = np.array(fl)
+    print flavor_list
+    print "Done loading data"
+
+    if not use_same_canvas :
+        filename = "triangle_plot_"+os.path.splitext(sys.argv[i+1])[0]
+    else :
+        filename = "triangle_plot"
+
+    # plots scale and diviions
+    scale = 8
+    divisions = 40
+
+    print "Begin making plot ..."
+    triangle_collection = LVCPT.GenerateTriangles(scale,divisions*2)
+    map(lambda f : LVCPT.AddPointToTriangleCollection(f,triangle_collection, scale = scale, divisions = divisions*2),flavor_list)
+
+    if use_same_canvas:
+        figure = LVCPT.MakeFlavorTriangle(flavor_list, divisions = 5, save_file=True, 
+                       filename = filename + "_hist", icolor = "g", icolormap = colors_schemes[i],
+                       triangle_collection=triangle_collection, figure = figure, alpha = 0.8,
+                       initial_flavor_ratio = [0,1,0], subdivisions = True, color_scale = "lin",
+                       output_format = output_format, inner_line_color ="silver",add_default_text = False,
+                       plot_color_bar =True)
+
+    else:
+        figure = LVCPT.MakeFlavorTriangle(flavor_list, divisions = 5, save_file=True,
+                       filename = filename + "_hist", icolor = "g", icolormap = colors_schemes[i],
+                       triangle_collection=triangle_collection, alpha = 0.8,
+                       initial_flavor_ratio = [0,1,0], subdivisions = True, color_scale = "lin",
+                       output_format = output_format, inner_line_color = "silver",add_default_text = False,
+                       plot_color_bar =True)
+
+    print "Done making plot: " + filename + "_hist."+output_format
+
+ax = figure.get_axes()[0]
+#ax.plot(6.5-0.35,5.5+0.3+0.35,"o", color = "grey", ms  = 10, markeredgewidth = 0.1, alpha = 0.9)
+#ax.text(6.7-0.35,5.44+0.3+0.35,r"$(1-x:x:0)$", fontsize = 16)
+#ax.axvline(x = 7.9)
+fsz = 32
+ax.plot(6.5-0.1,5.5+0.3+0.35+0.2+0.2,"o", color = "gold", ms  = 25, markeredgewidth = 0.1, alpha = 0.9)
+ax.text(6.7-0.1,5.44+0.3+0.35+0.2+0.2,r"$(1:2:0)$", fontsize = fsz)
+
+ax.plot(6.5-0.1,5.5+0.35+0.2,"o", color = "#2B653E", ms  = 25, markeredgewidth = 0.1, alpha = 0.9)
+ax.text(6.7-0.1,5.44+0.35+0.2,r"$(1:0:0)$", fontsize = fsz)
+
+ax.plot(6.5-0.1,5.5-0.3+0.35-0.2+0.2,"o", color = "#CA323D", ms  = 25, markeredgewidth = 0.1, alpha = 0.9)
+ax.text(6.7-0.1,5.44-0.3+0.35-0.2+0.2,r"$(0:1:0)$", fontsize = fsz)
+
+ax.plot(6.5-0.1,5.5-0.3+0.35-0.3-0.4+0.2,"o", color = "#2D4676", ms  = 25, markeredgewidth = 0.1, alpha = 0.9)
+ax.text(6.7-0.1,5.44-0.3+0.35-0.3-0.4+0.2,r"$(0:0:1)$", fontsize = fsz)
+
+plt.savefig("./plots/"+filename+"."+output_format, dpi = 600, bbox_inches='tight')
+
+exit(1)
+
+##os.system("cd plots")
+##os.system("gs triangle_plot_hist.eps")
+##os.system("cd ..")
+
diff --git a/plot_llh/mcmc_scan.py b/plot_llh/mcmc_scan.py
new file mode 100644
index 0000000..455d5dd
--- /dev/null
+++ b/plot_llh/mcmc_scan.py
@@ -0,0 +1,174 @@
+#! /usr/bin/env python
+"""
+Sample points from a gaussian likelihood
+"""
+
+from __future__ import absolute_import, division
+
+import sys
+
+import argparse
+import multiprocessing
+
+import numpy as np
+from scipy.stats import multivariate_normal
+
+import emcee
+import tqdm
+
+MEASURED_FR = [1, 1, 1]
+SIGMA = 0.001
+
+
+def angles_to_fr(angles):
+    sphi4, c2psi = 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 triangle_llh(theta):
+    """-Log likelihood function for a given theta."""
+    fr = angles_to_fr(theta)
+    fr_bf = MEASURED_FR
+    cov_fr = np.identity(3) * SIGMA
+    return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr))
+
+
+def lnprior(theta):
+    """Priors on theta."""
+    sphi4, c2psi = theta
+
+    # Flavour ratio bounds
+    if 0. <= sphi4 <= 1.0 and -1.0 <= c2psi <= 1.0:
+        pass
+    else: return -np.inf
+
+    return 0.
+
+
+def lnprob(theta):
+    """Prob function for mcmc."""
+    lp = lnprior(theta)
+    if not np.isfinite(lp):
+        return -np.inf
+    return lp + triangle_llh(theta)
+
+
+def parse_args():
+    """Parse command line arguments"""
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument(
+        '--measured-ratio', type=int, nargs=3, default=[1, 1, 1],
+        help='Set the central value for the measured flavour ratio at IceCube'
+    )
+    parser.add_argument(
+        '--sigma-ratio', type=float, default=0.01,
+        help='Set the 1 sigma for the measured flavour ratio'
+    )
+    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'
+    )
+    parser.add_argument(
+        '--seed', type=int, default=99,
+        help='Set the random seed value'
+    )
+    parser.add_argument(
+        '--outfile', type=str, default='./untitled',
+        help='Path to output chains'
+    )
+    args = parser.parse_args()
+    return args
+
+
+def main():
+    args = parse_args()
+
+    np.random.seed(args.seed)
+
+    global MEASURED_FR
+    global SIGMA
+
+    MEASURED_FR = np.array(args.measured_ratio) / float(np.sum(args.measured_ratio))
+    SIGMA = args.sigma_ratio
+
+    print 'MEASURED_FR = {0}'.format(MEASURED_FR)
+    print 'SIGMA = {0}'.format(SIGMA)
+
+    ndim = 2
+    nwalkers = args.nwalkers
+    ntemps = 1
+    burnin = args.burnin
+    betas = np.array([1e0, 1e-1, 1e-2, 1e-3, 1e-4])
+    p0_base = [0.5, 0.5]
+    p0_std = [.2, 0.2]
+
+    print 'p0_base', p0_base
+    print 'p0_std', p0_std
+    p0 = np.random.normal(p0_base, p0_std, size=[ntemps, nwalkers, ndim])
+    print map(lnprior, p0[0])
+
+    # threads = multiprocessing.cpu_count()
+    threads = 1
+    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"
+
+    nsteps = args.nsteps
+
+    print "Running"
+    for _ in tqdm.tqdm(sampler.sample(pos, iterations=nsteps), total=nsteps):
+        pass
+    print "Finished"
+
+    outfile = args.outfile+'_{0:03d}_{1:03d}_{2:03d}_{3:.1E}'.format(
+            int(MEASURED_FR[0]*100), int(MEASURED_FR[1]*100), int(MEASURED_FR[2]*100), SIGMA
+    )
+
+    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 FOR FILE ' + outfile
+    print 'outfile = ', outfile
+
+    fr_samples = np.array(map(angles_to_fr, samples))
+    np.save(outfile+'.npy', fr_samples)
+
+    print "DONE!"
+
+
+main.__doc__ = __doc__
+
+
+if __name__ == '__main__':
+    main()
+