GOLEMFIT takes in Haar measure params, fix. Add Bayes Factor calculation

5 files changed, 37 insertions, 8 deletions

diff --git a/utils/enums.py b/utils/enums.py
index 4637eeb..90ca17d 100644
--- a/utils/enums.py
+++ b/utils/enums.py
@@ -30,6 +30,7 @@ class Likelihood(Enum):
     FLAT     = 1
     GAUSSIAN = 2
     GOLEMFIT = 3
+    GF_FREQ  = 4
 
 
 class ParamTag(Enum):
diff --git a/utils/fr.py b/utils/fr.py
index 7f9d855..8f71f78 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -194,6 +194,23 @@ def normalise_fr(fr):
     return np.array(fr) / float(np.sum(fr))
 
 
+def fr_to_angles(ratios):
+    """Convert from flavour ratio into the angular projection of the flavour
+    ratios.
+
+    Parameters
+    ----------
+    TODO(shivesh)
+    """
+    fr0, fr1, fr2 = normalise_fr(ratios)
+    sphi4 = (fr2 - 1.0)**2
+    if (fr2 - 1.0) == 0:
+        c2psi = 0
+    else:
+        c2psi = (fr1*2.0 + fr2 - 1.0) * (fr2 - 1.0)
+    return sphi4, c2psi
+
+
 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)"""
 
diff --git a/utils/gf.py b/utils/gf.py
index ebc8538..e575094 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -44,14 +44,21 @@ def set_up_as(fitter, params):
 
 def get_llh(fitter, params):
     fitparams = gf.FitParameters(gf.sampleTag.HESE)
-    # print params
     for parm in params:
         fitparams.__setattr__(parm.name, parm.value)
     llh = -fitter.EvalLLH(fitparams)
-    # print '=== llh = {0}'.format(llh)
     return llh
 
 
+def get_llh_freq(fitter, params):
+    print 'setting to {0}'.format(params)
+    fitparams = gf.FitParameters(gf.sampleTag.HESE)
+    for parm in params:
+        fitparams.__setattr__(parm.name, parm.value)
+    fitter.SetFitParametersSeed(fitparams)
+    return -fitter.MinLLH().likelihood
+
+
 def data_distributions(fitter):
     hdat = fitter.GetDataDistribution()
     binedges = np.asarray([fitter.GetZenithBinsData(), fitter.GetEnergyBinsData()])
diff --git a/utils/likelihood.py b/utils/likelihood.py
index fe4301d..c1208ab 100644
--- a/utils/likelihood.py
+++ b/utils/likelihood.py
@@ -101,7 +101,7 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter):
     elif args.energy_dependance is EnergyDependance.SPECTRAL:
         mf_perbin = []
         for i_sf, sf_perbin in enumerate(source_flux):
-            u                     = fr_utils.params_to_BSMu(
+            u = fr_utils.params_to_BSMu(
                 theta             = bsm_angles,
                 dim               = args.dimension,
                 energy            = args.energy,
@@ -119,10 +119,9 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter):
             intergrated_measured_flux
         fr = averaged_measured_flux / np.sum(averaged_measured_flux)
 
+    flavour_angles = fr_utils.fr_to_angles(fr)
     for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)):
-        param.value = fr[idx]
-
-    # print 'hypo_paramset', hypo_paramset
+        param.value = flavour_angles[idx]
 
     if args.likelihood is Likelihood.FLAT:
         return 1.
@@ -131,6 +130,9 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter):
         return gaussian_llh(fr, fr_bf, args.sigma_ratio)
     elif args.likelihood is Likelihood.GOLEMFIT:
         return gf_utils.get_llh(fitter, hypo_paramset)
+    elif args.likelihood is Likelihood.GF_FREQ:
+        return gf_utils.get_llh_freq(fitter, hypo_paramset)
+
 
 def ln_prob(theta, args, fitter, asimov_paramset, mcmc_paramset):
     lp = lnprior(theta, paramset=mcmc_paramset)
diff --git a/utils/plot.py b/utils/plot.py
index c70ffec..4c5ff1c 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -171,6 +171,7 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
     ranges = mcmc_paramset.ranges
 
     if args.plot_angles:
+        print "Making triangle plots"
         Tchain = raw
         g = plot_Tchain(Tchain, axes_labels, ranges)
 
@@ -182,8 +183,8 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
                     itv = interval(Tchain[:,i_ax_1], percentile=90.)
                     for l in itv:
                         ax_2.axvline(l, color='gray', ls='--')
-                        ax_2.set_title(r'${0:.2f}_{{-{1:.2f}}}^{{+{2:.2f}}}$'.format(
-                            itv[1], itv[0], itv[2]
+                        ax_2.set_title(r'${0:.2f}_{{{1:.2f}}}^{{+{2:.2f}}}$'.format(
+                            itv[1], itv[0]-itv[1], itv[2]-itv[1]
                         ), fontsize=10)
 
         # if not args.fix_mixing:
@@ -198,6 +199,7 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset):
             g.export(outfile+'_angles.'+of)
 
     if args.plot_elements:
+        print "Making triangle plots"
         if args.fix_mixing_almost:
             raise NotImplementedError
         nu_index = mcmc_paramset.from_tag(ParamTag.NUISANCE, index=True)