aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitignore6
-rw-r--r--__init__.py0
-rwxr-xr-xchainer_plot.py265
-rwxr-xr-xmcmc_scan.py563
-rw-r--r--nufit_u.txt14
-rw-r--r--plot_llh/LVCPT.py431
-rw-r--r--plot_llh/MinimalTools.py166
-rw-r--r--plot_llh/PhysConst.py390
-rw-r--r--plot_llh/make_plots.py117
-rw-r--r--plot_llh/mcmc_scan.py174
-rw-r--r--submitter/make_dag.py100
-rw-r--r--submitter/submit.sub37
12 files changed, 2263 insertions, 0 deletions
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..546aa78
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,6 @@
+*.npy
+nohup.out
+*.pyc
+*.nfs*
+plots/
+*.pdf
diff --git a/__init__.py b/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/__init__.py
diff --git a/chainer_plot.py b/chainer_plot.py
new file mode 100755
index 0000000..26dc164
--- /dev/null
+++ b/chainer_plot.py
@@ -0,0 +1,265 @@
+#! /usr/bin/env python
+"""
+From an MCMC chains file, make a triangle plot.
+"""
+
+from __future__ import absolute_import, division
+
+import sys, os
+import errno
+
+import numpy as np
+import matplotlib as mpl
+mpl.use('Agg')
+from matplotlib import rc, rcParams
+
+import getdist
+from getdist import plots
+from getdist import mcsamples
+
+import mcmc_scan
+
+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}
+
+
+def plot(infile, angles, outfile, measured_ratio, sigma_ratio, fix_sfr,
+ fix_mixing, fix_scale, source_ratio, scale, dimension, energy, scale_bounds):
+ """Make the triangle plot"""
+ if not angles:
+ if fix_mixing:
+ labels = [r'{\rm log}_{10}\Lambda', r'\phi_e', r'\phi_\mu', r'\phi_\tau']
+ elif fix_sfr:
+ if fix_scale:
+ labels = [r'\mid \tilde{U}_{e1} \mid', r'\mid \tilde{U}_{e2} \mid', r'\mid \tilde{U}_{e3} \mid', \
+ r'\mid \tilde{U}_{\mu1} \mid', r'\mid \tilde{U}_{\mu2} \mid', r'\mid \tilde{U}_{\mu3} \mid', \
+ r'\mid \tilde{U}_{\tau1} \mid', r'\mid \tilde{U}_{\tau2} \mid', r'\mid \tilde{U}_{\tau3} \mid']
+ else:
+ labels = [r'\mid \tilde{U}_{e1} \mid', r'\mid \tilde{U}_{e2} \mid', r'\mid \tilde{U}_{e3} \mid', \
+ r'\mid \tilde{U}_{\mu1} \mid', r'\mid \tilde{U}_{\mu2} \mid', r'\mid \tilde{U}_{\mu3} \mid', \
+ r'\mid \tilde{U}_{\tau1} \mid', r'\mid \tilde{U}_{\tau2} \mid', r'\mid \tilde{U}_{\tau3} \mid', \
+ r'{\rm log}_{10}(\Lambda)']
+ else:
+ if fix_scale:
+ labels = [r'\mid \tilde{U}_{e1} \mid', r'\mid \tilde{U}_{e2} \mid', r'\mid \tilde{U}_{e3} \mid', \
+ r'\mid \tilde{U}_{\mu1} \mid', r'\mid \tilde{U}_{\mu2} \mid', r'\mid \tilde{U}_{\mu3} \mid', \
+ r'\mid \tilde{U}_{\tau1} \mid', r'\mid \tilde{U}_{\tau2} \mid', r'\mid \tilde{U}_{\tau3} \mid', \
+ r'{\rm log}_{10}\Lambda', r'\phi_e', r'\phi_\mu', r'\phi_\tau']
+ else:
+ labels = [r'\mid \tilde{U}_{e1} \mid', r'\mid \tilde{U}_{e2} \mid', r'\mid \tilde{U}_{e3} \mid', \
+ r'\mid \tilde{U}_{\mu1} \mid', r'\mid \tilde{U}_{\mu2} \mid', r'\mid \tilde{U}_{\mu3} \mid', \
+ r'\mid \tilde{U}_{\tau1} \mid', r'\mid \tilde{U}_{\tau2} \mid', r'\mid \tilde{U}_{\tau3} \mid', \
+ r'\phi_e', r'\phi_\mu', r'\phi_\tau']
+ else:
+ if fix_sfr:
+ if fix_mixing:
+ assert 0
+ labels=[r'\tilde{s}_{12}^2', r'{\rm log}_{10}\Lambda']
+ elif fix_scale:
+ labels=[r'\tilde{s}_{12}^2', r'\tilde{c}_{13}^4',
+ r'\tilde{s}_{23}^2', r'\tilde{\delta_{CP}}']
+ else:
+ labels=[r'\tilde{s}_{12}^2', r'\tilde{c}_{13}^4',
+ r'\tilde{s}_{23}^2', r'\tilde{\delta_{CP}}',
+ r'{\rm log}_{10}\Lambda']
+ else:
+ if fix_mixing:
+ labels=[r'{\rm log}_{10}\Lambda', r'sin^4(\phi)', r'cos(2\psi)']
+ elif fix_scale:
+ labels=[r'\tilde{s}_{12}^2', r'\tilde{c}_{13}^4',
+ r'\tilde{s}_{23}^2', r'\tilde{\delta_{CP}}',
+ r'sin^4(\phi)', r'cos(2\psi)']
+ else:
+ labels=[r'\tilde{s}_{12}^2', r'\tilde{c}_{13}^4',
+ r'\tilde{s}_{23}^2', r'\tilde{\delta_{CP}}',
+ r'{\rm log}_{10}\Lambda', r'sin^4(\phi)', r'cos(2\psi)']
+ print labels
+
+ if not fix_scale:
+ s2 = np.log10(scale_bounds)
+
+ if not angles:
+ if fix_mixing:
+ ranges = [s2, (0, 1), (0, 1), (0, 1)]
+ elif fix_sfr:
+ if fix_scale:
+ ranges = [(0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1)]
+ else:
+ ranges = [(0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), s2]
+ else:
+ if fix_scale:
+ ranges = [(0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1)]
+ else:
+ ranges = [(0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), (0, 1), s2, (0, 1), (0, 1), (0, 1)]
+ else:
+ if fix_sfr:
+ if fix_mixing:
+ ranges = [(0, 1), s2]
+ elif fix_scale:
+ ranges = [(0, 1), (0, 1), (0, 1), (0, 2*np.pi)]
+ else:
+ ranges = [(0, 1), (0, 1), (0, 1), (0, 2*np.pi), s2]
+ else:
+ if fix_mixing:
+ ranges = [s2, (0, 1), (-1, 1)]
+ elif fix_scale:
+ ranges = [(0, 1), (0, 1), (0, 1), (0, 2*np.pi), (0, 1), (-1, 1)]
+ else:
+ ranges = [(0, 1), (0, 1), (0, 1), (0, 2*np.pi), s2, (0, 1), (-1, 1)]
+
+ def flat_angles_to_u(x):
+ return abs(mcmc_scan.angles_to_u(x)).astype(np.float32).flatten().tolist()
+
+ raw = np.load(infile)
+ print 'raw.shape', raw.shape
+ if not angles:
+ if fix_mixing:
+ fr_elements = np.array(map(mcmc_scan.angles_to_fr, raw[:,-2:]))
+ sc_elements = raw[:,:-2]
+ Tchain = np.column_stack([sc_elements, fr_elements])
+ elif fix_sfr:
+ if fix_scale:
+ Tchain = np.array(map(flat_angles_to_u, raw))
+ else:
+ sc_elements = raw[:,-1:]
+ m_elements = np.array(map(flat_angles_to_u, raw[:,:-1]))
+ Tchain = np.column_stack([m_elements, sc_elements])
+ else:
+ if fix_scale:
+ fr_elements = np.array(map(mcmc_scan.angles_to_fr, raw[:,-2:]))
+ m_elements = np.array(map(flat_angles_to_u, raw[:,:-2]))
+ Tchain = np.column_stack([m_elements, fr_elements])
+ else:
+ fr_elements = np.array(map(mcmc_scan.angles_to_fr, raw[:,-2:]))
+ sc_elements = raw[:,-3:-2]
+ m_elements = np.array(map(flat_angles_to_u, raw[:,:-3]))
+ Tchain = np.column_stack([m_elements, sc_elements, fr_elements])
+ else:
+ Tchain = raw
+
+ if fix_sfr:
+ if fix_scale:
+ label = 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC observed flavour ratio = [{3:.2f}, {4:.2f}, {5:.2f}]\nSigma = {6:.3f}\nDimension = {7}\nEnergy = {8} GeV\nScale = {9}'.format(
+ source_ratio[0], source_ratio[1], source_ratio[2],
+ measured_ratio[0], measured_ratio[1], measured_ratio[2], sigma_ratio,
+ dimension, int(energy), scale
+ )
+ else:
+ label = 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC observed flavour ratio = [{3:.2f}, {4:.2f}, {5:.2f}]\nSigma = {6:.3f}\nDimension = {7}\nEnergy = {8} GeV'.format(
+ source_ratio[0], source_ratio[1], source_ratio[2],
+ measured_ratio[0], measured_ratio[1], measured_ratio[2], sigma_ratio,
+ dimension, int(energy)
+ )
+ else:
+ if fix_scale:
+ label = 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nSigma = {3:.3f}\nDimension = {4}\nEnergy = {5} GeV\nScale = {6}'.format(
+ measured_ratio[0], measured_ratio[1], measured_ratio[2], sigma_ratio,
+ dimension, int(energy), scale
+ )
+ else:
+ label = 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nSigma = {3:.3f}\nDimension = {4}\nEnergy = {5} GeV'.format(
+ measured_ratio[0], measured_ratio[1], measured_ratio[2], sigma_ratio,
+ dimension, int(energy)
+ )
+
+ Tsample = mcsamples.MCSamples(
+ samples=Tchain, labels=labels, ranges=ranges
+ )
+
+ Tsample.updateSettings({'contours': [0.90, 0.99]})
+ Tsample.num_bins_2D=500
+ Tsample.fine_bins_2D=500
+ Tsample.smooth_scale_2D=0.03
+
+ g = plots.getSubplotPlotter()
+ g.settings.num_plot_contours = 2
+ g.settings.axes_fontsize = 10
+ g.settings.figure_legend_frame = False
+ g.triangle_plot(
+ [Tsample], filled=True,
+ )
+ if fix_mixing and fix_sfr:
+ mpl.pyplot.figtext(0.4, 0.7, label, fontsize=4)
+ else:
+ mpl.pyplot.figtext(0.5, 0.7, label, fontsize=15)
+ print 'outfile = {0}'.format(outfile)
+ 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
+ g.export(outfile)
+
+
+def parse_args():
+ """Parse command line arguments"""
+ parser = argparse.ArgumentParser(description=__doc__)
+ parser.add_argument(
+ '--infile', type=str, required=True,
+ help='Path to MCMC chains'
+ )
+ parser.add_argument(
+ '--angles', default=False, action='store_true',
+ help='Plot in terms of mixing angles'
+ )
+ parser.add_argument(
+ '--outfile', type=str, default='./untitled.pdf',
+ help='Path to output plot'
+ )
+ parser.add_argument(
+ '--bestfit-ratio', type=int, nargs=3, required=False,
+ help='Set the bestfit flavour ratio'
+ )
+ parser.add_argument(
+ '--sigma-ratio', type=float, required=False,
+ help='Set the 1 sigma for the flavour ratio'
+ )
+ parser.add_argument(
+ '--fix-sfr', action='store_true',
+ help='Fix the source flavour ratio'
+ )
+ parser.add_argument(
+ '--fix-mixing', action='store_true',
+ help='Fix the new physics mixing values to a single term, s_12^2'
+ )
+ parser.add_argument(
+ '--source-ratio', type=int, nargs=3, default=[2, 1, 0],
+ help='Set the source flavour ratio for the case when you want to fix it'
+ )
+ parser.add_argument(
+ '--scale', type=float, required=False,
+ help='Fix the scale to this value'
+ )
+ parser.add_argument(
+ '--dimension', type=int, default=3, help='Dimension'
+ )
+ parser.add_argument(
+ '--energy', type=float, default=1000, help='Energy'
+ )
+ parser.add_argument(
+ '--scale-bounds', type=float, nargs=2,
+ help='Upper and lower limits to plot the new physics scale'
+ )
+ args = parser.parse_args()
+ return args
+
+
+def main():
+ args = vars(parse_args())
+ plot(**args)
+
+ print "DONE!"
+
+
+main.__doc__ = __doc__
+
+
+if __name__ == '__main__':
+ main()
diff --git a/mcmc_scan.py b/mcmc_scan.py
new file mode 100755
index 0000000..f79e2af
--- /dev/null
+++ b/mcmc_scan.py
@@ -0,0 +1,563 @@
+#! /usr/bin/env python
+"""
+From a gaussian likelihood run an MCMC scan to find the posteriors
+"""
+
+from __future__ import absolute_import, division
+
+import os, sys
+import errno
+
+import argparse
+import multiprocessing
+
+import numpy as np
+from scipy.stats import multivariate_normal
+from scipy import linalg
+
+import emcee
+import tqdm
+
+import chainer_plot
+
+
+RUN_SCAN = True
+
+FIX_MIXING = False
+FIX_SFR = True
+SOURCE_FR = [1, 2, 0]
+FIX_SCALE = True
+
+DIMENSION = 3
+ENERGY = 1000000 # GeV
+MEASURED_FR = [1, 1, 1]
+SIGMA = 0.001
+SCALE_RATIO = 100.
+MASS_EIGENVALUES = [7.40E-23, 2.515E-21]
+# MASS_EIGENVALUES = [7.40E-100, 2.515E-100]
+FLAT = False
+
+CHANGE_MASS = False
+if CHANGE_MASS:
+ SCALE = 1E-27
+else:
+ SCALE = np.power(10, np.round(np.log10(MASS_EIGENVALUES[1]/ENERGY)) - np.log10(ENERGY**(DIMENSION-3)))
+SCALE2_BOUNDS = (SCALE*1E-10, SCALE*1E10)
+
+
+def test_uni(x):
+ """Test the unitarity of a matrix."""
+ # print 'Unitarity test:\n{0}'.format(abs(np.dot(x, x.conj().T)))
+ return abs(np.dot(x, x.conj().T))
+
+
+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 angles_to_u(angles):
+ s12_2, c13_4, s23_2, dcp = 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)
+ """
+ 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
+
+
+# s_12^2, c_13^4, s_23^2, dcp
+NUFIT_U = angles_to_u((0.307, (1-0.02195)**2, 0.565, 3.97935))
+
+
+def params_to_BSMu(theta):
+ 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 / SCALE_RATIO
+
+ 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))
+
+ 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**(DIMENSION-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)
+ tu = test_uni(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 u_to_fr(initial_fr, matrix):
+ """Compute the observed flavour ratio assuming decoherence."""
+ # TODO(shivesh): energy dependence
+ composition = np.einsum(
+ 'ai, bi, a -> b', abs(matrix)**2, abs(matrix)**2, initial_fr
+ )
+ ratio = composition / np.sum(initial_fr)
+ return ratio
+
+
+def triangle_llh(theta):
+ """-Log likelihood function for a given theta."""
+ if FIX_SFR:
+ fr1, fr2, fr3 = SOURCE_FR
+ u = params_to_BSMu(theta)
+ else:
+ fr1, fr2, fr3 = angles_to_fr(theta[-2:])
+ u = params_to_BSMu(theta[:-2])
+
+ fr = u_to_fr((fr1, fr2, fr3), u)
+ fr_bf = MEASURED_FR
+ cov_fr = np.identity(3) * SIGMA
+ if FLAT:
+ return 10.
+ else:
+ return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr))
+
+
+def lnprior(theta):
+ """Priors on theta."""
+ if FIX_SFR:
+ if FIX_MIXING:
+ s12_2, sc2 = theta
+ elif FIX_SCALE:
+ s12_2, c13_4, s23_2, dcp = theta
+ else:
+ s12_2, c13_4, s23_2, dcp, sc2 = theta
+ else:
+ if FIX_MIXING:
+ sc2, sphi4, c2psi = theta
+ elif FIX_SCALE:
+ s12_2, c13_4, s23_2, dcp, sphi4, c2psi = theta
+ else:
+ s12_2, c13_4, s23_2, dcp, sc2, sphi4, c2psi = theta
+ if not FIX_SCALE:
+ sc2 = np.power(10., sc2)
+
+ if not FIX_SFR:
+ # Flavour ratio bounds
+ if 0. <= sphi4 <= 1.0 and -1.0 <= c2psi <= 1.0:
+ pass
+ else: return -np.inf
+
+ # Mixing angle bounds
+ if not FIX_MIXING:
+ if 0. <= s12_2 <= 1. and 0. <= c13_4 <= 1. and 0. <= s23_2 <= 1. \
+ and 0. <= dcp <= 2*np.pi:
+ pass
+ else: return -np.inf
+
+ # Scale bounds
+ if not FIX_SCALE:
+ if SCALE2_BOUNDS[0] <= sc2 <= SCALE2_BOUNDS[1]:
+ 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(
+ '--fix-source-ratio', type=str, default='False',
+ help='Fix the source flavour ratio'
+ )
+ parser.add_argument(
+ '--source-ratio', type=int, nargs=3, default=[2, 1, 0],
+ help='Set the source flavour ratio for the case when you want to fix it'
+ )
+ parser.add_argument(
+ '--fix-mixing', type=str, default='False',
+ help='Fix all mixing parameters except one'
+ )
+ parser.add_argument(
+ '--fix-scale', type=str, default='False',
+ help='Fix the new physics scale'
+ )
+ parser.add_argument(
+ '--scale', type=float, default=1e-30,
+ help='Set the new physics scale'
+ )
+ parser.add_argument(
+ '--dimension', type=int, default=3,
+ help='Set the new physics dimension to consider'
+ )
+ parser.add_argument(
+ '--energy', type=float, default=1000,
+ help='Set the energy scale'
+ )
+ parser.add_argument(
+ '--flat-llh', type=str, default='False',
+ help='Run with a flat likelihood'
+ )
+ 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 DIMENSION
+ global ENERGY
+ global MEASURED_FR
+ global SIGMA
+ global FLAT
+ global FIX_SFR
+ global SOURCE_FR
+ global FIX_MIXING
+ global FIX_SCALE
+ global SCALE
+ global SCALE2_BOUNDS
+
+ DIMENSION = args.dimension
+ ENERGY = args.energy
+
+ MEASURED_FR = np.array(args.measured_ratio) / float(np.sum(args.measured_ratio))
+ SIGMA = args.sigma_ratio
+
+ if args.flat_llh.lower() == 'true':
+ FLAT = True
+ elif args.flat_llh.lower() == 'false':
+ FLAT = False
+ else:
+ raise ValueError
+
+ if args.fix_source_ratio.lower() == 'true':
+ FIX_SFR = True
+ elif args.fix_source_ratio.lower() == 'false':
+ FIX_SFR = False
+ else:
+ raise ValueError
+
+ if args.fix_mixing.lower() == 'true':
+ FIX_MIXING = True
+ elif args.fix_mixing.lower() == 'false':
+ FIX_MIXING = False
+ else:
+ raise ValueError
+
+ if args.fix_scale.lower() == 'true':
+ FIX_SCALE = True
+ elif args.fix_scale.lower() == 'false':
+ FIX_SCALE = False
+ else:
+ raise ValueError
+
+ if FIX_SFR:
+ SOURCE_FR = np.array(args.source_ratio) / float(np.sum(args.source_ratio))
+
+ if FIX_SCALE:
+ SCALE = args.scale
+ else:
+ if CHANGE_MASS:
+ SCALE = 1E-27
+ else:
+ SCALE = np.power(10, np.round(np.log10(MASS_EIGENVALUES[1]/ENERGY)) - np.log10(ENERGY**(DIMENSION-3)))
+
+ if FIX_MIXING and FIX_SFR:
+ raise NotImplementedError('Fixed mixing and sfr not implemented')
+ if FIX_MIXING and FIX_SCALE:
+ raise NotImplementedError('Fixed mixing and scale not implemented')
+
+ SCALE2_BOUNDS = (SCALE*1E-10, SCALE*1E10)
+
+ print 'RUN_SCAN = {0}'.format(RUN_SCAN)
+ print 'MEASURED_FR = {0}'.format(MEASURED_FR)
+ print 'SIGMA = {0}'.format(SIGMA)
+ print 'FLAT = {0}'.format(FLAT)
+ print 'ENERGY = {0}'.format(ENERGY)
+ print 'DIMENSION = {0}'.format(DIMENSION)
+ print 'SCALE2_BOUNDS = {0}'.format(SCALE2_BOUNDS)
+ print 'FIX_SFR = {0}'.format(FIX_SFR)
+ if FIX_SFR:
+ print 'SOURCE_FR = {0}'.format(SOURCE_FR)
+ print 'FIX_MIXING = {0}'.format(FIX_MIXING)
+ print 'FIX_SCALE = {0}'.format(FIX_SCALE)
+ if FIX_SCALE:
+ print 'SCALE = {0}'.format(SCALE)
+
+ if FIX_SFR:
+ if FIX_MIXING:
+ ndim = 2
+ elif FIX_SCALE:
+ ndim = 4
+ else:
+ ndim = 5
+ else:
+ if FIX_MIXING:
+ ndim = 3
+ elif FIX_SCALE:
+ ndim = 6
+ else:
+ ndim = 7
+ nwalkers = args.nwalkers
+ ntemps = 1
+ burnin = args.burnin
+ betas = np.array([1e0, 1e-1, 1e-2, 1e-3, 1e-4])
+ s2 = np.average(np.log10(SCALE2_BOUNDS))
+ if FIX_SFR:
+ if FIX_MIXING:
+ p0_base = [0.5, s2]
+ p0_std = [0.2, 3]
+ elif FIX_SCALE:
+ p0_base = [0.5, 0.5, 0.5, np.pi]
+ p0_std = [0.2, 0.2, 0.2, 0.2]
+ else:
+ p0_base = [0.5, 0.5, 0.5, np.pi, s2]
+ p0_std = [0.2, 0.2, 0.2, 0.2, 3]
+ else:
+ if FIX_MIXING:
+ p0_base = [s2, 0.5, 0.0]
+ p0_std = [3, 0.2, 0.2]
+ elif FIX_SCALE:
+ p0_base = [0.5, 0.5, 0.5, np.pi, 0.5, 0.0]
+ p0_std = [0.2, 0.2, 0.2, 0.2, 0.2, 0.2]
+ else:
+ p0_base = [0.5, 0.5, 0.5, np.pi, s2, 0.5, 0.0]
+ p0_std = [0.2, 0.2, 0.2, 0.2, 3, 0.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])
+
+ if RUN_SCAN:
+ # threads = multiprocessing.cpu_count()
+ threads = 1
+ sampler = emcee.PTSampler(
+ ntemps, nwalkers, ndim, triangle_llh, lnprior, threads=threads
+ )
+
+ if RUN_SCAN:
+ 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"
+
+ if FIX_SFR:
+ if 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(MEASURED_FR[0]*100), int(MEASURED_FR[1]*100), int(MEASURED_FR[2]*100), int(SIGMA*1000),
+ int(SOURCE_FR[0]*100), int(SOURCE_FR[1]*100), int(SOURCE_FR[2]*100), DIMENSION
+ )
+ elif 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(MEASURED_FR[0]*100), int(MEASURED_FR[1]*100), int(MEASURED_FR[2]*100), int(SIGMA*1000),
+ int(SOURCE_FR[0]*100), int(SOURCE_FR[1]*100), int(SOURCE_FR[2]*100), DIMENSION, 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(MEASURED_FR[0]*100), int(MEASURED_FR[1]*100), int(MEASURED_FR[2]*100), int(SIGMA*1000),
+ int(SOURCE_FR[0]*100), int(SOURCE_FR[1]*100), int(SOURCE_FR[2]*100), DIMENSION
+ )
+ else:
+ if FIX_MIXING:
+ outfile = args.outfile+'_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fix_mixing'.format(
+ int(MEASURED_FR[0]*100), int(MEASURED_FR[1]*100), int(MEASURED_FR[2]*100),
+ int(SIGMA*1000), DIMENSION
+ )
+ elif FIX_SCALE:
+ outfile = args.outfile+'_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}_fixed_scale_{5}'.format(
+ int(MEASURED_FR[0]*100), int(MEASURED_FR[1]*100), int(MEASURED_FR[2]*100),
+ int(SIGMA*1000), DIMENSION, SCALE
+ )
+ else:
+ outfile = args.outfile+'_{0:03d}_{1:03d}_{2:03d}_{3:04d}_DIM{4}'.format(
+ int(MEASURED_FR[0]*100), int(MEASURED_FR[1]*100), int(MEASURED_FR[2]*100),
+ int(SIGMA*1000), DIMENSION
+ )
+
+ if RUN_SCAN:
+ 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
+
+ if FLAT:
+ outfile += '_flat'
+
+ if RUN_SCAN:
+ 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
+ np.save(outfile+'.npy', samples)
+
+ # print "Making triangle plots"
+ # chainer_plot.plot(
+ # infile=outfile+'.npy',
+ # angles=True,
+ # outfile=outfile[:5]+outfile[5:].replace('data', 'plots')+'_angles.pdf',
+ # measured_ratio=MEASURED_FR,
+ # sigma_ratio=SIGMA,
+ # fix_sfr=FIX_SFR,
+ # fix_mixing=FIX_MIXING,
+ # fix_scale=FIX_SCALE,
+ # source_ratio=SOURCE_FR,
+ # scale=SCALE,
+ # dimension=DIMENSION,
+ # energy=ENERGY,
+ # scale_bounds=SCALE2_BOUNDS,
+ # )
+ # chainer_plot.plot(
+ # infile=outfile+'.npy',
+ # angles=False,
+ # outfile=outfile[:5]+outfile[5:].replace('data', 'plots')+'.pdf',
+ # measured_ratio=MEASURED_FR,
+ # sigma_ratio=SIGMA,
+ # fix_sfr=FIX_SFR,
+ # fix_mixing=FIX_MIXING,
+ # fix_scale=FIX_SCALE,
+ # source_ratio=SOURCE_FR,
+ # scale=SCALE,
+ # dimension=DIMENSION,
+ # energy=ENERGY,
+ # scale_bounds=SCALE2_BOUNDS,
+ # )
+ print "DONE!"
+
+
+main.__doc__ = __doc__
+
+
+if __name__ == '__main__':
+ main()
+
diff --git a/nufit_u.txt b/nufit_u.txt
new file mode 100644
index 0000000..a9c9cb2
--- /dev/null
+++ b/nufit_u.txt
@@ -0,0 +1,14 @@
+
+[[ 0.73544986+0.j 0.48950332+0.j -0.31349347+0.34816928j]
+ [-0.16927408+0.2178619j 0.67961271+0.14500529j 0.66406513+0.j ]
+ [ 0.58860262+0.19116206j -0.51117311+0.12723432j 0.58268130+0.j ]]
+
+abs:
+[[ 0.82327921 0.54796108 0.14815532]
+ [ 0.31112912 0.59211521 0.74336952]
+ [ 0.47477365 0.5908792 0.65226662]]
+
+M
+[[0, 0, 0 ],
+ [0, 7.40E-23, 0 ],
+ [0, 0, 2.515E-21]]
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()
+
diff --git a/submitter/make_dag.py b/submitter/make_dag.py
new file mode 100644
index 0000000..10381f6
--- /dev/null
+++ b/submitter/make_dag.py
@@ -0,0 +1,100 @@
+#! /usr/bin/env python
+
+import numpy as np
+
+a_fr = (1, 2, 0)
+b_fr = (1, 0, 0)
+c_fr = (0, 1, 0)
+d_fr = (0, 0, 1)
+e_fr = (1, 1, 1)
+f_fr = (2, 1, 0)
+g_fr = (1, 1, 0)
+
+full_scan_mfr = [
+ (1, 1, 1), (1, 1, 0)
+]
+
+fix_sfr_mfr = [
+ (1, 1, 1, 1, 0, 0),
+ (1, 1, 1, 0, 1, 0),
+ (1, 1, 1, 0, 0, 1),
+ (1, 1, 1, 1, 2, 0),
+ (1, 1, 0, 0, 1, 0),
+ (1, 1, 0, 1, 2, 0),
+ (1, 1, 0, 1, 0, 0),
+ (1, 0, 0, 1, 0, 0),
+ (0, 1, 0, 0, 1, 0),
+ (1, 2, 0, 0, 1, 0),
+ (1, 2, 0, 1, 2, 0)
+]
+
+sigmas = ['0.1', '0.01']
+dimensions = [6]
+energy = [1e4, 1e6, 1e7]
+flat = False
+burnin = 1000
+nwalkers = 200
+nsteps = 10000
+scales = "1E-20 1E-30"
+
+outfile = 'dagman_FR.submit'
+condor_script = '/home/smandalia/Documents/flavour_ratio/submitter/submit.sub'
+
+with open(outfile, 'w') as f:
+ job_number = 1
+ for dim in dimensions:
+ print 'dimension', dim
+ for en in energy:
+ print 'energy {0:.0E}'.format(en)
+
+ outchain_head = '/data/user/smandalia/flavour_ratio/data/DIM{0}/{1:.0E}'.format(dim, en)
+
+ for sig in sigmas:
+ print 'sigma', sig
+ for frs in fix_sfr_mfr:
+ print frs
+ outchains = outchain_head + '/fix_ifr/{0}/mcmc_chain'.format(str(sig).replace('.', '_'))
+ f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
+ f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
+ f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
+ f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
+ f.write('VARS\tjob{0}\tsigma="{1}"\n'.format(job_number, sig))
+ f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, 'True'))
+ f.write('VARS\tjob{0}\tsr0="{1}"\n'.format(job_number, frs[3]))
+ f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, frs[4]))
+ f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, frs[5]))
+ f.write('VARS\tjob{0}\tfix_scale="{1}"\n'.format(job_number, 'False'))
+ f.write('VARS\tjob{0}\tscale="{1}"\n'.format(job_number, 0))
+ f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
+ f.write('VARS\tjob{0}\tenergy="{1}"\n'.format(job_number, en))
+ f.write('VARS\tjob{0}\tflat_llh="{1}"\n'.format(job_number, flat))
+ f.write('VARS\tjob{0}\tburnin="{1}"\n'.format(job_number, burnin))
+ f.write('VARS\tjob{0}\tnwalkers="{1}"\n'.format(job_number, nwalkers))
+ f.write('VARS\tjob{0}\tnsteps="{1}"\n'.format(job_number, nsteps))
+ f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outchains))
+ f.write('VARS\tjob{0}\tfix_mixing="{1}"\n'.format(job_number, 'False'))
+ job_number += 1
+
+ for frs in full_scan_mfr:
+ print frs
+ outchains = outchain_head + '/full_scan/{0}/mcmc_chain'.format(str(sig).replace('.', '_'))
+ f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
+ f.write('VARS\tjob{0}\tmr0="{1}"\n'.format(job_number, frs[0]))
+ f.write('VARS\tjob{0}\tmr1="{1}"\n'.format(job_number, frs[1]))
+ f.write('VARS\tjob{0}\tmr2="{1}"\n'.format(job_number, frs[2]))
+ f.write('VARS\tjob{0}\tsigma="{1}"\n'.format(job_number, sig))
+ f.write('VARS\tjob{0}\tfix_source_ratio="{1}"\n'.format(job_number, 'False'))
+ f.write('VARS\tjob{0}\tsr0="{1}"\n'.format(job_number, 0))
+ f.write('VARS\tjob{0}\tsr1="{1}"\n'.format(job_number, 0))
+ f.write('VARS\tjob{0}\tsr2="{1}"\n'.format(job_number, 0))
+ f.write('VARS\tjob{0}\tfix_scale="{1}"\n'.format(job_number, 'False'))
+ f.write('VARS\tjob{0}\tscale="{1}"\n'.format(job_number, 0))
+ f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
+ f.write('VARS\tjob{0}\tenergy="{1}"\n'.format(job_number, en))
+ f.write('VARS\tjob{0}\tflat_llh="{1}"\n'.format(job_number, flat))
+ f.write('VARS\tjob{0}\tburnin="{1}"\n'.format(job_number, burnin))
+ f.write('VARS\tjob{0}\tnwalkers="{1}"\n'.format(job_number, nwalkers))
+ f.write('VARS\tjob{0}\tnsteps="{1}"\n'.format(job_number, nsteps))
+ f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, outchains))
+ f.write('VARS\tjob{0}\tfix_mixing="{1}"\n'.format(job_number, 'False'))
+ job_number += 1
diff --git a/submitter/submit.sub b/submitter/submit.sub
new file mode 100644
index 0000000..9e21fc5
--- /dev/null
+++ b/submitter/submit.sub
@@ -0,0 +1,37 @@
+Executable = /home/smandalia/Documents/flavour_ratio/mcmc_scan.py
+Arguments = --measured-ratio $(mr0) $(mr1) $(mr2) --sigma-ratio $(sigma) --fix-source-ratio $(fix_source_ratio) --source-ratio $(sr0) $(sr1) $(sr2) --fix-scale $(fix_scale) --scale $(scale) --dimension $(dimension) --energy $(energy) --flat-llh $(flat_llh) --burnin $(burnin) --nwalkers $(nwalkers) --nsteps $(nsteps) --seed 24 --outfile $(outfile) --fix-mixing $(fix_mixing)
+
+# All logs will go to a single file
+log = /home/smandalia/Documents/flavour_ratio/submitter/logs/job_$(Cluster).log
+output = /home/smandalia/Documents/flavour_ratio/submitter/logs/job_$(Cluster).out
+error = /home/smandalia/Documents/flavour_ratio/submitter/logs/job_$(Cluster).err
+
+getenv = True
+# environment = "X509_USER_PROXY=x509up_u14830"
+
+# Stage user cert to the node (Gridftp-Users is already on CVMFS)
+# transfer_input_files = /tmp/x509up_u14830
+
+# but do not try to copy outputs back (see: https://htcondor-wiki.cs.wisc.edu/index.cgi/tktview?tn=3081)
+# +TransferOutput=""
+
+request_memory = 1GB
+request_cpus = 1
+
+Universe = vanilla
+Notification = never
+
+# run on both SL5 and 6
+# +WantRHEL6 = True
+# +WantSLC6 = False
+
+# # run on OSG
+# +WantGlidein = True
+
+# +TransferOutput=""
+
+# Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6")
+# Requirements = IS_GLIDEIN
+
+# GO!
+queue
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-17 11:17:12 +0000