Thu 8 Nov 15:02:27 CST 2018

4 files changed, 45 insertions, 10 deletions

diff --git a/utils/fr.py b/utils/fr.py
index a2c2356..81ba3b2 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -248,7 +248,10 @@ 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 args.energy_dependance is EnergyDependance.MONO:
+    if args.scale is not None:
+        scale = args.scale
+        scale_region = (scale/args.scale_region, scale*args.scale_region)
+    elif 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))
@@ -287,7 +290,7 @@ def fr_argparse(parser):
         help='Type of energy dependance to use in the BSM fit'
     )
     parser.add_argument(
-        '--spectral-index', default=-2, type=int,
+        '--spectral-index', default=-2, type=float,
         help='Spectral index for spectral energy dependance'
     )
     parser.add_argument(
@@ -316,7 +319,7 @@ def fr_argparse(parser):
         help='Fix the new physics scale'
     )
     parser.add_argument(
-        '--scale', type=float, default=1e-30,
+        '--scale', type=float, default=None,
         help='Set the new physics scale'
     )
     parser.add_argument(
@@ -415,6 +418,9 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
             '--fix-mixing and --fix-mixing-almost cannot be used together'
         )
 
+    if not isinstance(theta, list):
+        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:
@@ -429,7 +435,11 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
         sc2 = np.log10(scale)
     else:
         s12_2, c13_4, s23_2, dcp, sc2 = theta
-    sc2 = np.power(10., sc2)
+
+    if len(theta) != 0:
+        sc2 = np.power(10., sc2)
+    else:
+        sc2 = np.log10(scale)
     sc1 = sc2 / 100.
 
     mass_matrix = np.array(
@@ -515,8 +525,16 @@ def u_to_fr(source_fr, matrix):
         array([ 0.33740075,  0.33176584,  0.33083341])
 
     """
-    composition = np.einsum(
-        'ai, bi, a -> b', np.abs(matrix)**2, np.abs(matrix)**2, source_fr,
-    )
+    try:
+        composition = np.einsum(
+            'ai, bi, a -> b', np.abs(matrix)**2, np.abs(matrix)**2, source_fr,
+        )
+    except:
+        matrix = np.array(matrix, dtype=np.complex256)
+        composition = np.einsum(
+            'ai, bi, a -> b', np.abs(matrix)**2, np.abs(matrix)**2, source_fr,
+        )
+        pass
+
     ratio = composition / np.sum(source_fr)
     return ratio
diff --git a/utils/likelihood.py b/utils/likelihood.py
index 4580b67..93f3aea 100644
--- a/utils/likelihood.py
+++ b/utils/likelihood.py
@@ -17,7 +17,7 @@ from scipy.stats import multivariate_normal, rv_continuous
 from utils import fr as fr_utils
 from utils import gf as gf_utils
 from utils.enums import EnergyDependance, Likelihood, ParamTag, PriorsCateg
-from utils.misc import enum_parse
+from utils.misc import enum_parse, gen_identifier, parse_bool
 
 
 class Gaussian(rv_continuous):
@@ -41,6 +41,14 @@ def likelihood_argparse(parser):
         '--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'
+    )
 
 
 def lnprior(theta, paramset):
@@ -165,6 +173,14 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
             intergrated_measured_flux
         fr = averaged_measured_flux / np.sum(averaged_measured_flux)
 
+    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} {3:.1f}'.format(
+                fr[0], fr[1], fr[2], llh_paramset['logLam'].value
+            ))
+            f.write('\n')
+
     flavour_angles = fr_utils.fr_to_angles(fr)
     # print 'flavour_angles', map(float, flavour_angles)
     for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)):
diff --git a/utils/mcmc.py b/utils/mcmc.py
index fa044ea..e5bd8da 100644
--- a/utils/mcmc.py
+++ b/utils/mcmc.py
@@ -20,7 +20,7 @@ from utils.enums import MCMCSeedType
 from utils.misc import enum_parse, make_dir, parse_bool
 
 
-def mcmc(p0, ln_prob, ndim, nwalkers, burnin, nsteps, threads=1):
+def mcmc(p0, ln_prob, ndim, nwalkers, burnin, nsteps, args, threads=1):
     """Run the MCMC."""
     sampler = emcee.EnsembleSampler(
         nwalkers, ndim, ln_prob, threads=threads
@@ -31,6 +31,7 @@ def mcmc(p0, ln_prob, ndim, nwalkers, burnin, nsteps, threads=1):
         pos, prob, state = result
     sampler.reset()
     print "Finished burn-in"
+    args.burnin = False
 
     print "Running"
     for _ in tqdm.tqdm(sampler.sample(pos, iterations=nsteps), total=nsteps):
diff --git a/utils/param.py b/utils/param.py
index fe0a0a0..62f59b1 100644
--- a/utils/param.py
+++ b/utils/param.py
@@ -237,7 +237,7 @@ def get_paramsets(args, nuisance_paramset):
     ])
     asimov_paramset = ParamSet(asimov_paramset)
 
-    if args.fix_mixing is not MixingScenario.NONE and not args.fix_mixing_almost:
+    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),