diff options
| -rwxr-xr-x | fr.py | 10 | ||||
| -rwxr-xr-x | plot_sens.py | 30 | ||||
| -rwxr-xr-x | sens.py | 17 | ||||
| -rw-r--r-- | submitter/mcmc_dag.py | 12 | ||||
| -rw-r--r-- | submitter/sens_dag.py | 22 | ||||
| -rw-r--r-- | utils/likelihood.py | 3 |
6 files changed, 68 insertions, 26 deletions
@@ -30,11 +30,13 @@ def define_nuisance(): """Define the nuisance parameters.""" tag = ParamTag.SM_ANGLES g_prior = PriorsCateg.GAUSSIAN + hg_prior = PriorsCateg.HALFGAUSS + e = 1e-9 nuisance = [ - Param(name='s_12_2', value=0.307, seed=[0.29, 0.31], ranges=[0., 1.], std=0.013, tex=r's_{12}^2', prior=g_prior, tag=tag), - Param(name='c_13_4', value=1-(0.02206)**2, seed=[0.998, 1.0], ranges=[0., 1.], std=0.00147, tex=r'c_{13}^4', prior=g_prior, tag=tag), - Param(name='s_23_2', value=0.538, seed=[0.46, 0.61], ranges=[0., 1.], std=0.069, tex=r's_{23}^2', prior=g_prior, tag=tag), - Param(name='dcp', value=4.08404, seed=[0, 2*np.pi], ranges=[0., 2*np.pi], std=2.0, tex=r'\delta_{CP}', tag=tag), + Param(name='s_12_2', value=0.307, seed=[0.26, 0.35], ranges=[0., 1.], std=0.013, tex=r's_{12}^2', prior=g_prior, tag=tag), + Param(name='c_13_4', value=1-(0.02206)**2, seed=[0.995, 1-e], ranges=[0., 1.], std=0.00147, tex=r'c_{13}^4', prior=hg_prior, tag=tag), + Param(name='s_23_2', value=0.538, seed=[0.31, 0.75], ranges=[0., 1.], std=0.069, tex=r's_{23}^2', prior=g_prior, tag=tag), + Param(name='dcp', value=4.08404, seed=[0+e, 2*np.pi-e], ranges=[0., 2*np.pi], std=2.0, tex=r'\delta_{CP}', tag=tag), Param( name='m21_2', value=7.40E-23, seed=[7.2E-23, 7.6E-23], ranges=[6.80E-23, 8.02E-23], std=2.1E-24, tex=r'\Delta m_{21}^2{\rm GeV}^{-2}', prior=g_prior, tag=tag diff --git a/plot_sens.py b/plot_sens.py new file mode 100755 index 0000000..e51d55d --- /dev/null +++ b/plot_sens.py @@ -0,0 +1,30 @@ +#! /usr/bin/env python +# author : S. Mandalia +# s.p.mandalia@qmul.ac.uk +# +# date : April 28, 2018 + +""" +HESE BSM flavour ratio analysis plotting script +""" + +from __future__ import absolute_import, division + +import numpy as np +import matplotlib as mpl +mpl.use('Agg') +from matplotlib import rc +from matplotlib import pyplot as plt +from matplotlib.offsetbox import AnchoredText + +rc('text', usetex=False) +rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) + +FRS = [ + (1, 1, 1, 1, 2, 0), + (1, 1, 1, 0, 1, 0), + # (1, 1, 1, 1, 0, 0), + # (1, 1, 1, 0, 0, 1), +] + +DIMENSION = [3, 6] @@ -84,10 +84,6 @@ def process_args(args): '--sens-run and --fix-scale cannot be used together' ) - if args.sens_eval_bin is not None and args.plot_statistic: - print 'Cannot make plot when specific scale bin is chosen' - args.plot_statistic = False - args.measured_ratio = fr_utils.normalise_fr(args.measured_ratio) if args.fix_source_ratio: args.source_ratio = fr_utils.normalise_fr(args.source_ratio) @@ -96,6 +92,10 @@ def process_args(args): args.binning = np.logspace( np.log10(args.binning[0]), np.log10(args.binning[1]), args.binning[2]+1 ) + if args.likelihood is Likelihood.GOLEMFIT: + print 'GolemFit selected with spectral index energy dependance, ' \ + 'will attempt to use the astroDeltaGamma systematic to fold ' \ + 'in the spectral index.' if not args.fix_scale: args.scale = fr_utils.estimate_scale(args) @@ -106,6 +106,10 @@ def process_args(args): else: args.sens_eval_bin = int(args.sens_eval_bin) + if args.sens_eval_bin is not None and args.plot_statistic: + print 'Cannot make plot when specific scale bin is chosen' + args.plot_statistic = False + if args.stat_method is StatCateg.FREQUENTIST and \ args.likelihood is Likelihood.GOLEMFIT: args.likelihood = Likelihood.GF_FREQ @@ -285,12 +289,12 @@ def main(): llh_paramset[name].value = \ (pranges[i][1]-pranges[i][0])*x[i] + pranges[i][0] theta = llh_paramset.values - print 'llh_paramset', llh_paramset llh = llh_utils.ln_prob( theta=theta, args=args, asimov_paramset=asimov_paramset, llh_paramset=llh_paramset, fitter=fitter ) - print 'llh', llh + # print 'llh_paramset', llh_paramset + # print 'llh', llh return -llh n_params = len(sens_paramset) @@ -321,6 +325,7 @@ def main(): np.save(out+'.npy', statistic_arr) if args.plot_statistic: + print 'Plotting statistic' if args.sens_run: raw = statistic_arr else: raw = np.load(out+'.npy') data = ma.masked_invalid(raw, 0) diff --git a/submitter/mcmc_dag.py b/submitter/mcmc_dag.py index 2866887..827ab9e 100644 --- a/submitter/mcmc_dag.py +++ b/submitter/mcmc_dag.py @@ -9,9 +9,9 @@ full_scan_mfr = [ fix_sfr_mfr = [ (1, 1, 1, 1, 2, 0), - # (1, 1, 1, 1, 0, 0), - # (1, 1, 1, 0, 1, 0), - # (1, 1, 1, 0, 0, 1), + (1, 1, 1, 1, 0, 0), + (1, 1, 1, 0, 1, 0), + (1, 1, 1, 0, 0, 1), # (1, 1, 0, 1, 2, 0), # (1, 1, 0, 1, 0, 0), # (1, 1, 0, 0, 1, 0), @@ -48,7 +48,7 @@ GLOBAL_PARAMS.update(dict( # Likelihood GLOBAL_PARAMS.update(dict( - likelihood = 'gaussian', + likelihood = 'golemfit', sigma_ratio = '0.01' )) @@ -72,8 +72,8 @@ with open(outfile, 'w') as f: job_number = 1 for dim in dimension: print 'dimension', dim - outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format( - GLOBAL_PARAMS['likelihood'], dim, GLOBAL_PARAMS['spectral_index'] + outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/'.format( + GLOBAL_PARAMS['likelihood'], dim ) for frs in fix_sfr_mfr: print 'frs', frs diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py index 63aaaba..17b063d 100644 --- a/submitter/sens_dag.py +++ b/submitter/sens_dag.py @@ -9,8 +9,8 @@ full_scan_mfr = [ fix_sfr_mfr = [ (1, 1, 1, 1, 2, 0), - # (1, 1, 1, 1, 0, 0), - # (1, 1, 1, 0, 1, 0), + (1, 1, 1, 1, 0, 0), + (1, 1, 1, 0, 1, 0), # (1, 1, 1, 0, 0, 1), # (1, 1, 0, 1, 2, 0), # (1, 1, 0, 1, 0, 0), @@ -27,7 +27,7 @@ GLOBAL_PARAMS = {} sens_eval_bin = 'all' # set to 'all' to run normally GLOBAL_PARAMS.update(dict( sens_run = 'True', - run_method = 'corr_angle', + run_method = 'fixed_angle', # full, fixed_angle, corr_angle stat_method = 'frequentist', sens_bins = 10, seed = 'None' @@ -55,7 +55,7 @@ GLOBAL_PARAMS.update(dict( # Likelihood GLOBAL_PARAMS.update(dict( - likelihood = 'gaussian', + likelihood = 'golemfit', sigma_ratio = '0.01' )) @@ -87,8 +87,8 @@ with open(outfile, 'w') as f: job_number = 1 for dim in dimension: print 'dimension', dim - outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format( - GLOBAL_PARAMS['likelihood'], dim, GLOBAL_PARAMS['spectral_index'] + outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}'.format( + GLOBAL_PARAMS['likelihood'], dim ) for frs in fix_sfr_mfr: print 'frs', frs @@ -107,7 +107,10 @@ with open(outfile, 'w') as f: 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}\tsens_eval_bin="{1}"\n'.format(job_number, r)) + if sens_eval_bin.lower() != 'all': + f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, r)) + else: + f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, 'all')) for key in GLOBAL_PARAMS.iterkeys(): f.write('VARS\tjob{0}\t{1}="{2}"\n'.format(job_number, key, GLOBAL_PARAMS[key])) f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, output)) @@ -130,7 +133,10 @@ with open(outfile, 'w') as f: 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}\tsens_eval_bin="{1}"\n'.format(job_number, r)) + if sens_eval_bin.lower() != 'all': + f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, r)) + else: + f.write('VARS\tjob{0}\tsens_eval_bin="{1}"\n'.format(job_number, 'all')) for key in GLOBAL_PARAMS.iterkeys(): f.write('VARS\tjob{0}\t{1}="{2}"\n'.format(job_number, key, GLOBAL_PARAMS[key])) f.write('VARS\tjob{0}\toutfile="{1}"\n'.format(job_number, output)) diff --git a/utils/likelihood.py b/utils/likelihood.py index 9629b65..70b54c9 100644 --- a/utils/likelihood.py +++ b/utils/likelihood.py @@ -89,8 +89,7 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): bin_width = np.abs(np.diff(args.binning)) if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]: if 'astroDeltaGamma' in hypo_paramset.names: - args.spectral_index = hypo_paramset['astroDeltaGamma'].value - print 'args.spectral_index', args.spectral_index + args.spectral_index = -hypo_paramset['astroDeltaGamma'].value if args.fix_source_ratio: if args.energy_dependance is EnergyDependance.MONO: |