From ccffb521195eb5f41471e166e1ba8f695740bcb3 Mon Sep 17 00:00:00 2001 From: shivesh Date: Fri, 6 Apr 2018 17:21:57 -0500 Subject: add test scripts for Golem LV and NSI --- fr.py | 14 +++++--- submitter/make_dag.py | 21 ++++++----- submitter/submit.sub | 2 +- test/LV.out | 74 +++++++++++++++++++++++++++++++++++++++ test/NSI.out | 75 +++++++++++++++++++++++++++++++++++++++ test/test.png | Bin 0 -> 104565 bytes test/test_LV.py | 95 ++++++++++++++++++++++++++++++++++++++++++++++++++ test/test_NSI.png | Bin 0 -> 107718 bytes test/test_NSI.py | 89 ++++++++++++++++++++++++++++++++++++++++++++++ utils/enums.py | 27 ++------------ utils/gf.py | 70 +------------------------------------ utils/mcmc.py | 9 ++++- 12 files changed, 365 insertions(+), 111 deletions(-) create mode 100644 test/LV.out create mode 100644 test/NSI.out create mode 100644 test/test.png create mode 100644 test/test_LV.py create mode 100644 test/test_NSI.png create mode 100644 test/test_NSI.py 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 > already registered; second conversion method ignored. + import GolemFitPy as gf +test_LV.py:11: RuntimeWarning: to-Python converter for nusquids::marray > 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 > already registered; second conversion method ignored. + import GolemFitPy as gf +test_NSI.py:11: RuntimeWarning: to-Python converter for nusquids::marray > 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 Binary files /dev/null and b/test/test.png differ 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 Binary files /dev/null and b/test/test_NSI.png differ 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): -- cgit v1.2.3