add test scripts for Golem LV and NSI

12 files changed, 365 insertions, 111 deletions

diff --git a/fr.py b/fr.py
index c5dffaf..2de310c 100755
--- a/fr.py
+++ b/fr.py
@@ -21,7 +21,7 @@ from utils import gf as gf_utils
 from utils import likelihood as llh_utils
 from utils import mcmc as mcmc_utils
 from utils import misc as misc_utils
-from utils.enums import Likelihood, ParamTag
+from utils.enums import Likelihood, MCMCSeedType, ParamTag
 from utils.fr import MASS_EIGENVALUES, normalise_fr
 from utils.misc import Param, ParamSet
 from utils.plot import plot_argparse, chainer_plot
@@ -219,10 +219,14 @@ def main():
 
         ndim = len(mcmc_paramset)
         ntemps = 1
-        # p0 = mcmc_utils.gaussian_seed(
-        p0 = mcmc_utils.flat_seed(
-            mcmc_paramset, ntemps=ntemps, nwalkers=args.nwalkers
-        )
+        if args.mcmc_seed_type == MCMCSeedType.UNIFORM:
+            p0 = mcmc_utils.flat_seed(
+                mcmc_paramset, ntemps=ntemps, nwalkers=args.nwalkers
+            )
+        elif args.mcmc_seed_type == MCMCSeedType.GAUSSIAN:
+            p0 = mcmc_utils.gaussian_seed(
+                mcmc_paramset, ntemps=ntemps, nwalkers=args.nwalkers
+            )
 
         samples = mcmc_utils.mcmc(
             p0           = p0,
diff --git a/submitter/make_dag.py b/submitter/make_dag.py
index daddb65..33d05b3 100644
--- a/submitter/make_dag.py
+++ b/submitter/make_dag.py
@@ -36,6 +36,7 @@ nsteps   = 100
 nwalkers = 200
 seed     = 24
 threads  = 1
+mcmc_seed_type = 'uniform'
 
 # FR
 dimension    = [3]
@@ -46,6 +47,9 @@ sigma_ratio  = ['0.01']
 scale        = "1E-20 1E-30"
 scale_region = "1E10"
 
+# Likelihood
+likelihood = 'golemfit'
+
 # Nuisance
 astroDeltaGamma = 2.
 astroNorm       = 1.
@@ -54,11 +58,8 @@ muonNorm        = 1.
 promptNorm      = 0.
 
 # GolemFit
-aft    = 'hesespl'
-ast    = 'baseline'
-axs    = 'nom'
-data   = 'real'
-priors = 'uniform'
+ast  = 'p2_0'
+data = 'real'
 
 # Plot
 plot_angles   = 'True'
@@ -110,14 +111,13 @@ with open(outfile, 'w') as f:
                     f.write('VARS\tjob{0}\tmuonNorm="{1}"\n'.format(job_number, muonNorm))
                     f.write('VARS\tjob{0}\tpromptNorm="{1}"\n'.format(job_number, promptNorm))
                     f.write('VARS\tjob{0}\tdata="{1}"\n'.format(job_number, data))
-                    f.write('VARS\tjob{0}\tpriors="{1}"\n'.format(job_number, priors))
-                    f.write('VARS\tjob{0}\taft="{1}"\n'.format(job_number, aft))
                     f.write('VARS\tjob{0}\tast="{1}"\n'.format(job_number, ast))
-                    f.write('VARS\tjob{0}\taxs="{1}"\n'.format(job_number, axs))
                     f.write('VARS\tjob{0}\tplot_angles="{1}"\n'.format(job_number, plot_angles))
                     f.write('VARS\tjob{0}\tplot_elements="{1}"\n'.format(job_number, plot_elements))
                     f.write('VARS\tjob{0}\tseed="{1}"\n'.format(job_number, seed))
                     f.write('VARS\tjob{0}\tthreads="{1}"\n'.format(job_number, threads))
+                    f.write('VARS\tjob{0}\tlikelihood="{1}"\n'.format(job_number, likelihood))
+                    f.write('VARS\tjob{0}\tmcmc_seed_type="{1}"\n'.format(job_number, mcmc_seed_type))
                     job_number += 1
 
                 for frs in full_scan_mfr:
@@ -151,12 +151,11 @@ with open(outfile, 'w') as f:
                     f.write('VARS\tjob{0}\tmuonNorm="{1}"\n'.format(job_number, muonNorm))
                     f.write('VARS\tjob{0}\tpromptNorm="{1}"\n'.format(job_number, promptNorm))
                     f.write('VARS\tjob{0}\tdata="{1}"\n'.format(job_number, data))
-                    f.write('VARS\tjob{0}\tpriors="{1}"\n'.format(job_number, priors))
-                    f.write('VARS\tjob{0}\taft="{1}"\n'.format(job_number, aft))
                     f.write('VARS\tjob{0}\tast="{1}"\n'.format(job_number, ast))
-                    f.write('VARS\tjob{0}\taxs="{1}"\n'.format(job_number, axs))
                     f.write('VARS\tjob{0}\tplot_angles="{1}"\n'.format(job_number, plot_angles))
                     f.write('VARS\tjob{0}\tplot_elements="{1}"\n'.format(job_number, plot_elements))
                     f.write('VARS\tjob{0}\tseed="{1}"\n'.format(job_number, seed))
                     f.write('VARS\tjob{0}\tthreads="{1}"\n'.format(job_number, threads))
+                    f.write('VARS\tjob{0}\tlikelihood="{1}"\n'.format(job_number, likelihood))
+                    f.write('VARS\tjob{0}\tmcmc_seed_type="{1}"\n'.format(job_number, mcmc_seed_type))
                     job_number += 1
diff --git a/submitter/submit.sub b/submitter/submit.sub
index b57c64c..e932673 100644
--- a/submitter/submit.sub
+++ b/submitter/submit.sub
@@ -1,5 +1,5 @@
 Executable = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/fr.py
-Arguments = "--aft $(aft) --ast $(ast) --astroDeltaGamma $(astroDeltaGamma) --astroNorm $(astroNorm) --axs $(axs) --burnin $(burnin) --convNorm $(convNorm) --data $(data) --dimension $(dimension) --energy $(energy) --fix-mixing $(fix_mixing) --fix-scale $(fix_scale) --fix-source-ratio $(fix_source_ratio) --likelihood $(likelihood) --measured-ratio $(mr0) $(mr1) $(mr2) --muonNorm $(muonNorm) --no-bsm $(no_bsm) --nsteps $(nsteps) --nwalkers $(nwalkers) --outfile $(outfile) --plot-angles $(plot_angles) --plot-elements $(plot_elements) --priors $(priors) --promptNorm $(promptNorm) --run-mcmc $(run_mcmc) --scale $(scale) --scale-region $(scale_region) --seed $(seed) --sigma-ratio $(sigma_ratio) --source-ratio $(sr0) $(sr1) $(sr2) --threads $(threads)"
+Arguments = "--ast $(ast) --astroDeltaGamma $(astroDeltaGamma) --astroNorm $(astroNorm) --burnin $(burnin) --convNorm $(convNorm) --data $(data) --dimension $(dimension) --energy $(energy) --fix-mixing $(fix_mixing) --fix-scale $(fix_scale) --fix-source-ratio $(fix_source_ratio) --likelihood $(likelihood) --measured-ratio $(mr0) $(mr1) $(mr2) --muonNorm $(muonNorm) --no-bsm $(no_bsm) --nsteps $(nsteps) --nwalkers $(nwalkers) --outfile $(outfile) --plot-angles $(plot_angles) --plot-elements $(plot_elements) --promptNorm $(promptNorm) --run-mcmc $(run_mcmc) --scale $(scale) --scale-region $(scale_region) --seed $(seed) --sigma-ratio $(sigma_ratio) --source-ratio $(sr0) $(sr1) $(sr2) --threads $(threads) --likelihood $(likelihood) --mcmc-seed-type $(mcmc_seed_type)"
 
 # All logs will go to a single file 
 log    = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).log
diff --git a/test/LV.out b/test/LV.out
new file mode 100644
index 0000000..e9f2b7e
--- /dev/null
+++ b/test/LV.out
@@ -0,0 +1,74 @@
+test_LV.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 2u, std::allocator<double> > already registered; second conversion method ignored.
+  import GolemFitPy as gf
+test_LV.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 5u, std::allocator<double> > already registered; second conversion method ignored.
+  import GolemFitPy as gf
+reset_steering: 1
+reset_data: 1
+reset_npp: 1
+GolemFit constructor: checking paths
+Constructing DiffuseWeightMaker
+Loading data
+Loading HESE data.
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/data/HESEData.txt
+Loaded HESE 80 events.
+Loading MC
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 now main simulation size is 420660
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 now main simulation size is 616803
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 now main simulation size is 740041
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 now main simulation size is 740161
+Loaded 740161 events in main simulation set
+Loading Flux weighter
+Loading nusquids objects
+Weighting MC
+Initializing simulation weights
+Applying selfveto correction.
+HESE reshuffle has not been performed.
+Doing HESE reshuffle
+Making data hist
+Making sim hist
+Constructing likelihood problem with default settings
+Setting up Asimov set with parameters
+Conventional normalization 1
+Prompt normalization 1
+Atmospheric muon normalization 1
+Astro component overall normalization 8
+Astro E component 0.333333
+Astro Mu component 0.333333
+Astro Tau component 0.333333
+Astroparticle balance 1
+Astro gamma 2.5
+Astro cutoff 10
+CR gamma 0
+Conv pi/k ratio 1
+Conv particle balance 1
+Conv zenith correction 0
+Dark normalization 0
+DOM efficiency 0.99
+Astro component overall normalization second 0
+Astro gamma second 2
+Holeice forward 0
+Anisotropy scale1
+
+Remaking data hist
+Reconstrucing likelihood problem
+NULL min_llh 835.848062048
+NULL expectation [11.96478761 11.90169848 10.78840154  9.20158936  7.23605538  5.59952977
+  4.39850777  3.35653018  2.58491686  1.92834114  1.46813946  1.11581282
+  0.8492433   0.65154743  0.46727371  0.36211812  0.32718588  0.69746002
+  0.29475116  0.05150809]
+
+Applying new flavor physics to MC weights.
+Entering ApplyNewPhysics
+Entering ConstructNuSQuIDSObjectsForLV.
+Rotated Unitary Transform
+In for loop
+Entering inner for loop
+Leaving inner for loop
diff --git a/test/NSI.out b/test/NSI.out
new file mode 100644
index 0000000..657a9f3
--- /dev/null
+++ b/test/NSI.out
@@ -0,0 +1,75 @@
+test_NSI.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 2u, std::allocator<double> > already registered; second conversion method ignored.
+  import GolemFitPy as gf
+test_NSI.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 5u, std::allocator<double> > already registered; second conversion method ignored.
+  import GolemFitPy as gf
+reset_steering: 1
+reset_data: 1
+reset_npp: 1
+GolemFit constructor: checking paths
+Constructing DiffuseWeightMaker
+Loading data
+Loading HESE data.
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/data/HESEData.txt
+Loaded HESE 80 events.
+Loading MC
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 now main simulation size is 420660
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 now main simulation size is 616803
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 now main simulation size is 740041
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 now main simulation size is 740161
+Loaded 740161 events in main simulation set
+Loading Flux weighter
+Loading nusquids objects
+Weighting MC
+Initializing simulation weights
+Applying selfveto correction.
+HESE reshuffle has not been performed.
+Doing HESE reshuffle
+Making data hist
+Making sim hist
+Constructing likelihood problem with default settings
+Setting up Asimov set with parameters
+Conventional normalization 1
+Prompt normalization 1
+Atmospheric muon normalization 1
+Astro component overall normalization 8
+Astro E component 0.333333
+Astro Mu component 0.333333
+Astro Tau component 0.333333
+Astroparticle balance 1
+Astro gamma 2.5
+Astro cutoff 10
+CR gamma 0
+Conv pi/k ratio 1
+Conv particle balance 1
+Conv zenith correction 0
+Dark normalization 0
+DOM efficiency 0.99
+Astro component overall normalization second 0
+Astro gamma second 2
+Holeice forward 0
+Anisotropy scale1
+
+Remaking data hist
+Reconstrucing likelihood problem
+NULL min_llh 835.848062048
+NULL expectation [11.96478761 11.90169848 10.78840154  9.20158936  7.23605538  5.59952977
+  4.39850777  3.35653018  2.58491686  1.92834114  1.46813946  1.11581282
+  0.8492433   0.65154743  0.46727371  0.36211812  0.32718588  0.69746002
+  0.29475116  0.05150809]
+
+Applying new flavor physics to MC weights.
+Applying New Physics in normal mode to NonStandardInteraction
+HESE reshuffle has been performed.
+0.1 mutau min_llh 836.352921667
+0.1 mutau expectation [16.63858821 16.04604012 14.28034835 11.80228198  9.06453127  6.8815312
+  5.32704912  4.02384379  3.07832186  2.28036523  1.73329938  1.30341591
+  0.99239868  0.75237435  0.53312845  0.4093463   0.35723267  0.71433954
+  0.30173276  0.05276177]
diff --git a/test/test.png b/test/test.png
new file mode 100644
index 0000000..8aac9e3
--- /dev/null
+++ b/test/test.png
diff --git a/test/test_LV.py b/test/test_LV.py
new file mode 100644
index 0000000..96a1863
--- /dev/null
+++ b/test/test_LV.py
@@ -0,0 +1,95 @@
+#!/usr/bin/env python
+
+from __future__ import absolute_import, division
+
+import numpy as np
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+from matplotlib import rc
+
+import GolemFitPy as gf
+
+rc('text', usetex=True)
+rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
+
+dp = gf.DataPaths()
+steer = gf.SteeringParams()
+
+fig = plt.figure(figsize=[12, 8])
+ax = fig.add_subplot(111)
+ax.set_xscale('log')
+ax.set_yscale('log')
+
+npp = gf.NewPhysicsParams()
+npp.type = gf.NewPhysicsType.None
+# npp.n_lv = 1
+# npp.lambda_1 = 1.e100
+# npp.lambda_2 = 1.e100
+
+golem = gf.GolemFit(dp, steer, npp)
+
+binning = golem.GetEnergyBinsMC()
+ax.set_xlim(binning[0], binning[-1])
+# ax.set_ylim(binning[0], binning[-1])
+
+fit_params = gf.FitParameters(gf.sampleTag.HESE)
+golem.SetupAsimov(fit_params)
+
+exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
+ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
+        drawstyle='steps-pre', label='NULL', linestyle='--')
+
+print 'NULL min_llh', golem.MinLLH().likelihood
+print 'NULL expectation', exp
+print
+
+npp = gf.NewPhysicsParams()
+npp.type = gf.NewPhysicsType.LorentzViolation
+npp.n_lv = 1
+npp.lambda_1 = 1.e20
+npp.lambda_2 = 1.e20
+
+golem.SetNewPhysicsParams(npp)
+
+exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
+ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
+        drawstyle='steps-pre', label='1e-20 LV', linestyle='--')
+
+print '1e20 LV min_llh', golem.MinLLH().likelihood
+print '1e20 LV expectation', exp
+
+npp = gf.NewPhysicsParams()
+npp.type = gf.NewPhysicsType.LorentzViolation
+npp.n_lv = 1
+npp.lambda_1 = 1.e10
+npp.lambda_2 = 1.e10
+
+golem.SetNewPhysicsParams(npp)
+
+exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
+ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
+        drawstyle='steps-pre', label='1e-10 LV', linestyle='--')
+
+print '1e10 LV min_llh', golem.MinLLH().likelihood
+print '1e10 LV expectation', exp
+
+npp = gf.NewPhysicsParams()
+npp.type = gf.NewPhysicsType.LorentzViolation
+npp.n_lv = 1
+npp.lambda_1 = 1.e-20
+npp.lambda_2 = 1.e-20
+
+golem.SetNewPhysicsParams(npp)
+
+ax.tick_params(axis='x', labelsize=12)
+ax.tick_params(axis='y', labelsize=12)
+ax.set_xlabel(r'Deposited energy / GeV')
+ax.set_ylabel(r'Events')
+for xmaj in ax.xaxis.get_majorticklocs():
+    ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.7, linewidth=1)
+for ymaj in ax.yaxis.get_majorticklocs():
+    ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.7, linewidth=1)
+
+legend = ax.legend(prop=dict(size=12))
+fig.savefig('test.png', bbox_inches='tight', dpi=250)
diff --git a/test/test_NSI.png b/test/test_NSI.png
new file mode 100644
index 0000000..2eb4345
--- /dev/null
+++ b/test/test_NSI.png
diff --git a/test/test_NSI.py b/test/test_NSI.py
new file mode 100644
index 0000000..617c353
--- /dev/null
+++ b/test/test_NSI.py
@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+
+from __future__ import absolute_import, division
+
+import numpy as np
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+from matplotlib import rc
+
+import GolemFitPy as gf
+
+rc('text', usetex=True)
+rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18})
+
+dp = gf.DataPaths()
+steer = gf.SteeringParams()
+
+fig = plt.figure(figsize=[12, 8])
+ax = fig.add_subplot(111)
+ax.set_xscale('log')
+ax.set_yscale('log')
+
+npp = gf.NewPhysicsParams()
+npp.type = gf.NewPhysicsType.None
+# npp.epsilon_mutau = 0
+# npp.epsilon_prime = 0
+
+golem = gf.GolemFit(dp, steer, npp)
+
+binning = golem.GetEnergyBinsMC()
+ax.set_xlim(binning[0], binning[-1])
+# ax.set_ylim(binning[0], binning[-1])
+
+fit_params = gf.FitParameters(gf.sampleTag.HESE)
+golem.SetupAsimov(fit_params)
+
+exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
+ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
+        drawstyle='steps-pre', label='NULL', linestyle='--')
+
+print 'NULL min_llh', golem.MinLLH().likelihood
+print 'NULL expectation', exp
+print
+
+npp = gf.NewPhysicsParams()
+npp.type = gf.NewPhysicsType.NonStandardInteraction
+npp.epsilon_mutau = 0.1
+# npp.epsilon_prime = 0
+
+golem.SetNewPhysicsParams(npp)
+
+exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
+ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
+        drawstyle='steps-pre', label='0.1 mutau', linestyle='--')
+
+print '0.1 mutau min_llh', golem.MinLLH().likelihood
+print '0.1 mutau expectation', exp
+
+# npp = gf.NewPhysicsParams()
+# npp.epsilon_mutau = 0
+# # npp.epsilon_prime = 0
+
+# golem.SetNewPhysicsParams(npp)
+
+# exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3))
+# ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1,
+#         drawstyle='steps-pre', label='1e-10 LV', linestyle='--')
+
+# print '1e10 LV min_llh', golem.MinLLH().likelihood
+# print '1e10 LV expectation', exp
+
+# npp = gf.NewPhysicsParams()
+# npp.epsilon_mutau = 0
+# # npp.epsilon_prime = 0
+
+# golem.SetNewPhysicsParams(npp)
+
+ax.tick_params(axis='x', labelsize=12)
+ax.tick_params(axis='y', labelsize=12)
+ax.set_xlabel(r'Deposited energy / GeV')
+ax.set_ylabel(r'Events')
+for xmaj in ax.xaxis.get_majorticklocs():
+    ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.7, linewidth=1)
+for ymaj in ax.yaxis.get_majorticklocs():
+    ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.7, linewidth=1)
+
+legend = ax.legend(prop=dict(size=12))
+fig.savefig('test_NSI.png', bbox_inches='tight', dpi=250)
diff --git a/utils/enums.py b/utils/enums.py
index bce3da2..45f164d 100644
--- a/utils/enums.py
+++ b/utils/enums.py
@@ -16,19 +16,6 @@ class DataType(Enum):
     ASMIMOV = 3
 
 
-class FitCateg(Enum):
-    HESESPL  = 1
-    HESEDPL  = 2
-    ZPSPL    = 3
-    ZPDPL    = 4
-    NUNUBAR2 = 5
-
-
-class FitFlagCateg(Enum):
-    DEFAULT = 1
-    XS      = 2
-
-
 class FitPriorsCateg(Enum):
     DEFAULT = 1
     XS      = 2
@@ -49,10 +36,9 @@ class ParamTag(Enum):
     NONE      = 6
 
 
-class Priors(Enum):
-    UNIFORM    = 1
-    LOGUNIFORM = 2
-    JEFFREYS   = 3
+class MCMCSeedType(Enum):
+    UNIFORM  = 1
+    GAUSSIAN = 2
 
 
 class SteeringCateg(Enum):
@@ -69,10 +55,3 @@ class SteeringCateg(Enum):
     LONGLIFE              = 10
     DPL                   = 11
 
-class XSCateg(Enum):
-    HALF    = 1
-    NOM     = 2
-    TWICE   = 3
-    TWICE01 = 4
-    TWICE02 = 5
-
diff --git a/utils/gf.py b/utils/gf.py
index eaa9450..3fb063b 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -15,7 +15,7 @@ from functools import partial
 
 import GolemFitPy as gf
 
-from utils.enums import *
+from utils.enums import DataType, SteeringCateg
 from utils.misc import enum_parse, thread_factors
 
 
@@ -70,60 +70,6 @@ def data_distributions(fitter):
     return hdat, binedges
 
 
-def fit_flags(fitflag_categ):
-    flags = gf.FitParametersFlag()
-    if fitflag_categ is FitFlagCateg.xs:
-        # False means it's not fixed in minimization
-        flags.NeutrinoAntineutrinoRatio = False
-    return flags
-
-
-def fit_params(fit_categ):
-    params = gf.FitParameters()
-    params.astroNorm = 7.5
-    params.astroDeltaGamma = 0.9
-    if fit_categ is FitCateg.hesespl:
-        params.astroNormSec = 0
-    elif fit_categ is FitCateg.hesedpl:
-        params.astroNormSec=7.0
-    elif fit_categ is FitCateg.zpspl:
-        # zero prompt, single powerlaw
-        params.promptNorm = 0
-        params.astroNormSec = 0
-    elif fit_categ is FitCateg.zpdpl:
-        # zero prompt, double powerlaw
-        params.promptNorm = 0
-        params.astroNormSec=7.0
-    elif fit_categ is FitCateg.nunubar2:
-        params.NeutrinoAntineutrinoRatio = 2
-
-
-def fit_priors(fitpriors_categ):
-    priors = gf.Priors()
-    if fitpriors_categ == FitPriorsCateg.xs:
-        priors.promptNormCenter = 1
-        priors.promptNormWidth = 3
-        priors.astroDeltaGammaCenter = 0
-        priors.astroDeltaGammaWidth = 1
-    return priors
-
-
-def gen_steering_params(steering_categ, quiet=False):
-    params = gf.SteeringParams()
-    if quiet: params.quiet = True
-    params.fastmode = False
-    params.do_HESE_reshuffle = False
-    params.numc_tag = steering_categ.name
-    params.baseline_astro_spectral_index = -2.
-    if steering_categ is SteeringCateg.LONGLIFE:
-        params.years = [999]
-        params.numc_tag = 'std_half1'
-    if steering_categ is SteeringCateg.DPL:
-        params.diffuse_fit_type = gf.DiffuseFitType.DoublePowerLaw
-        params.numc_tag = 'std_half1'
-    return params
-
-
 def gf_argparse(parser):
     parser.add_argument(
         '--data', default='real', type=partial(enum_parse, c=DataType),
@@ -134,18 +80,4 @@ def gf_argparse(parser):
         choices=SteeringCateg,
         help='use asimov/fake dataset with specific steering'
     )
-    parser.add_argument(
-        '--aft', default='hesespl', type=partial(enum_parse, c=FitCateg),
-        choices=FitCateg,
-        help='use asimov/fake dataset with specific Fit'
-    )
-    parser.add_argument(
-        '--axs', default='nom', type=partial(enum_parse, c=XSCateg),
-        choices=XSCateg,
-        help='use asimov/fake dataset with xs scaling'
-    )
-    parser.add_argument(
-        '--priors', default='uniform', type=partial(enum_parse, c=Priors),
-        choices=Priors, help='Bayesian priors'
-    )
 
diff --git a/utils/mcmc.py b/utils/mcmc.py
index d2036da..0b78f1e 100644
--- a/utils/mcmc.py
+++ b/utils/mcmc.py
@@ -10,13 +10,15 @@ Useful functions to use an MCMC for the BSM flavour ratio analysis
 from __future__ import absolute_import, division
 
 import argparse
+from functools import partial
 
 import emcee
 import tqdm
 
 import numpy as np
 
-from utils.misc import make_dir, parse_bool
+from utils.enums import MCMCSeedType
+from utils.misc import enum_parse, make_dir, parse_bool
 
 
 def mcmc(p0, triangle_llh, lnprior, ndim, nwalkers, burnin, nsteps, ntemps=1, threads=1):
@@ -65,6 +67,11 @@ def mcmc_argparse(parser):
         '--nsteps', type=int, default=2000,
         help='Number of steps to run'
     )
+    parser.add_argument(
+        '--mcmc-seed-type', default='uniform',
+        type=partial(enum_parse, c=MCMCSeedType), choices=MCMCSeedType,
+        help='Type of distrbution to make the initial MCMC seed'
+    )
 
 
 def flat_seed(paramset, ntemps, nwalkers):