diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rwxr-xr-x | fr.py | 64 | ||||
| -rw-r--r-- | mnrun/plot.py | 2 | ||||
| -rw-r--r-- | mnrun/test.png | bin | 64972 -> 0 bytes | |||
| -rwxr-xr-x | sens.py | 48 | ||||
| -rwxr-xr-x | submitter/clean.sh | 5 | ||||
| -rw-r--r-- | submitter/make_dag.py | 55 | ||||
| -rw-r--r-- | submitter/submit.sub | 7 | ||||
| -rw-r--r-- | utils/gf.py | 9 | ||||
| -rw-r--r-- | utils/plot.py | 113 |
10 files changed, 225 insertions, 79 deletions
@@ -7,3 +7,4 @@ plots/ dagman_FR.submit* submitter/logs/ mnrun_* +*.png @@ -21,10 +21,10 @@ 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 import plot as plot_utils from utils.enums import EnergyDependance, Likelihood, MCMCSeedType, ParamTag from utils.fr import MASS_EIGENVALUES, normalise_fr, fr_to_angles from utils.misc import enum_parse, Param, ParamSet -from utils.plot import plot_argparse, chainer_plot def define_nuisance(): @@ -89,7 +89,8 @@ def get_paramsets(args): lL_range = (logLam-scale_region, logLam+scale_region) tag = ParamTag.SCALE mcmc_paramset.append( - Param(name='logLam', value=logLam, ranges=lL_range, std=3, tex=r'{\rm log}_{10}\Lambda', tag=tag) + Param(name='logLam', value=logLam, ranges=lL_range, std=3, + tex=r'{\rm log}_{10}\Lambda'+plot_utils.get_units(args.dimension), tag=tag) ) if not args.fix_source_ratio: tag = ParamTag.SRCANGLES @@ -98,7 +99,6 @@ def get_paramsets(args): Param(name='c_2psi', value=0.5, ranges=[0., 1.], std=0.2, tex=r'cos(2\psi)', tag=tag) ]) mcmc_paramset = ParamSet(mcmc_paramset) - # TODO(shivesh): unify return asimov_paramset, mcmc_paramset @@ -110,9 +110,9 @@ def process_args(args): raise NotImplementedError( '--fix-mixing and --fix-mixing-almost cannot be used together' ) - if args.bayes_factor and args.fix_scale: + if args.run_bayes_factor and args.fix_scale: raise NotImplementedError( - '--bayes-factor and --fix-scale cannot be used together' + '--run-bayes-factor and --fix-scale cannot be used together' ) args.measured_ratio = normalise_fr(args.measured_ratio) @@ -172,14 +172,18 @@ def parse_args(): help='Binning for spectral energy dependance' ) parser.add_argument( - '--bayes-factor', type=misc_utils.parse_bool, default='False', + '--run-bayes-factor', type=misc_utils.parse_bool, default='False', help='Make the bayes factor plot for the scale' ) parser.add_argument( - '--bayes-bins', type=int, default=100, + '--bayes-bins', type=int, default=10, help='Number of bins for the Bayes factor plot' ) parser.add_argument( + '--bayes-live-points', type=int, default=400, + help='Number of live points for MultiNest evaluations' + ) + parser.add_argument( '--bayes-output', type=str, default='./mnrun/', help='Folder to store MultiNest evaluations' ) @@ -233,9 +237,9 @@ def parse_args(): ) llh_utils.likelihood_argparse(parser) mcmc_utils.mcmc_argparse(parser) - nuisance_argparse(parser) gf_utils.gf_argparse(parser) - plot_argparse(parser) + plot_utils.plot_argparse(parser) + nuisance_argparse(parser) return parser.parse_args() @@ -252,12 +256,7 @@ def main(): if args.run_mcmc: if args.likelihood is Likelihood.GOLEMFIT: - datapaths = gf.DataPaths() - sparams = gf_utils.steering_params(args) - npp = gf.NewPhysicsParams() - fitter = gf.GolemFit(datapaths, sparams, npp) - gf_utils.set_up_as(fitter, asimov_paramset) - # fitter.WriteCompact() + fitter = gf_utils.setup_fitter(args, asimov_paramset) else: fitter = None @@ -289,7 +288,7 @@ def main(): ) mcmc_utils.save_chains(samples, outfile) - chainer_plot( + plot_utils.chainer_plot( infile = outfile+'.npy', outfile = outfile[:5]+outfile[5:].replace('data', 'plots'), outformat = ['pdf'], @@ -297,18 +296,14 @@ def main(): mcmc_paramset = mcmc_paramset ) - if args.bayes_factor: + if args.run_bayes_factor: # TODO(shivesh) import pymultinest from pymultinest.solve import solve from pymultinest.watch import ProgressPrinter if not args.run_mcmc and args.likelihood is Likelihood.GOLEMFIT: - datapaths = gf.DataPaths() - sparams = gf_utils.steering_params(args) - npp = gf.NewPhysicsParams() - fitter = gf.GolemFit(datapaths, sparams, npp) - gf_utils.set_up_as(fitter, asimov_paramset) + fitter = gf_utils.setup_fitter(args, asimov_paramset) else: fitter = None @@ -318,9 +313,8 @@ def main(): ) p = mcmc_paramset.from_tag(ParamTag.SCALE, invert=True) - prior_ranges = p.ranges n_params = len(p) - hypo_paramset = asimov_paramset + prior_ranges = p.ranges def CubePrior(cube, ndim, nparams): # default are uniform priors @@ -333,7 +327,7 @@ def main(): for i in range(ndim): prange = prior_ranges[i][1] - prior_ranges[i][0] theta[i] = prange*cube[i] + prior_ranges[i][0] - theta_ = np.concatenate([theta, [s]]) + theta_ = np.array(theta.tolist() + [s]) return llh_utils.triangle_llh( theta=theta_, args=args, @@ -342,25 +336,31 @@ def main(): fitter=fitter ) - prefix = 'mnrun_{0:.0E}'.format(np.power(10, s)) + prefix = 'mnrun_{0:.0E}'.format(np.power(10, s)) + '_' + outfile[2:] print 'begin running evidence calculation for {0}'.format(prefix) result = pymultinest.run( - LogLikelihood=lnProb, Prior=CubePrior, n_dims=n_params, - outputfiles_basename=prefix#, verbose=True + LogLikelihood=lnProb, + Prior=CubePrior, + n_dims=n_params, + n_live_points=args.bayes_live_points, + outputfiles_basename=prefix, + # resume=False ) - print 'end running evidence calculation for {0}'.format(prefix) - print 'begin analysis for {0}'.format(prefix) analyzer = pymultinest.Analyzer(outputfiles_basename=prefix, n_params=n_params) a_lnZ = analyzer.get_stats()['global evidence'] - print 'end analysis for {0}'.format(prefix) print 'Evidence = {0}'.format(a_lnZ) arr.append([s, a_lnZ]) - out = args.bayes_output+'fr_evidence_'+outfile[2:]+'.npy' + out = args.bayes_output+'fr_evidence_'+outfile[2:]+'.npy' misc_utils.make_dir(out) np.save(out, np.array(arr)) + b_out = args.bayes_output+'fr_evidence_'+outfile[2:] + plot_utils.bayes_factor_plot( + infile=b_out+'.npy', outfile=b_out, outformat=['png'], args=args + ) + print "DONE!" diff --git a/mnrun/plot.py b/mnrun/plot.py index e010257..f7e750b 100644 --- a/mnrun/plot.py +++ b/mnrun/plot.py @@ -7,7 +7,7 @@ from matplotlib import rc rc('text', usetex=False) rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) -arr = np.load('fr_evidence_test_050_050_000_0100_sfr_033_066_000_DIM3_single_scale.npy') +arr = np.load('fr_evidence_test_033_033_033_0010_sfr_033_066_000_DIM6_single_scale.npy') print arr ev = arr.T[1] diff --git a/mnrun/test.png b/mnrun/test.png Binary files differdeleted file mode 100644 index 3f68e17..0000000 --- a/mnrun/test.png +++ /dev/null @@ -23,12 +23,13 @@ from utils import gf as gf_utils from utils import likelihood as llh_utils from utils import misc as misc_utils from utils.enums import Likelihood, ParamTag +from utils.plot import plot_BSM_angles_limit rc('text', usetex=False) rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) -RUN = True +RUN = False z = 0+1E-6 @@ -77,6 +78,8 @@ def main(): fr.process_args(args) misc_utils.print_args(args) + bins = 10 + asimov_paramset, mcmc_paramset = fr.get_paramsets(args) outfile = misc_utils.gen_outfile_name(args) print '== {0:<25} = {1}'.format('outfile', outfile) @@ -85,9 +88,10 @@ def main(): mcmc_paramset = mcmc_paramset.from_tag(ParamTag.NUISANCE, invert=True) sc_range = mcmc_paramset.from_tag(ParamTag.SCALE)[0].ranges - scan_scales = np.linspace(sc_range[0], sc_range[1], 100) + scan_scales = np.linspace(sc_range[0], sc_range[1], bins+1) print 'scan_scales', scan_scales + outfile = './sens' if RUN: datapaths = gf.DataPaths() sparams = gf_utils.steering_params(args) @@ -96,10 +100,10 @@ def main(): fitter.SetFitParametersFlag(fit_flags(default_flags)) gf_utils.set_up_as(fitter, asimov_paramset) - x = [] - y = [] + data = np.zeros((len(scenarios), bins+1, 2)) mm_angles = mcmc_paramset.from_tag(ParamTag.MMANGLES) for idx, scen in enumerate(scenarios): + arr = [] scales = [] llhs = [] for yidx, an in enumerate(mm_angles): @@ -114,28 +118,20 @@ def main(): print 'llh', llh scales.append(sc) llhs.append(llh) - min_llh = np.min(llhs) - delta_llh = 2*(np.array(llhs) - min_llh) - print 'scales', scales - print 'delta_llh', delta_llh - - n_arr = [] - for i, d in enumerate(delta_llh): - # 90% for 1 DOF - if d < 2.71: - x.append(idx) - y.append(scales[i]) - - np.save('sens.npy', np.array([x, y])) - else: - x, y = np.load('sens.npy') - - plt.plot(x, y, linewidth=0., marker='.') - plt.xticks(range(len(scenarios)), xticks) - plt.xlim(-1, len(scenarios)) - plt.ylim(sc_range[0], sc_range[1]) - plt.ylabel(r'${\rm log}_{10}\Lambda / GeV$') - plt.savefig('./sens.png', bbox_inches='tight', dpi=150) + + for i, d in enumerate(llhs): + data[idx][i][0] = scales[i] + data[idx][i][1] = d + + np.save(outfile + '.npy', data) + + plot_BSM_angles_limit( + infile=outfile+'.npy', + outfile=outfile, + xticks=xticks, + outformat=['png'], + args=args + ) main.__doc__ = __doc__ diff --git a/submitter/clean.sh b/submitter/clean.sh new file mode 100755 index 0000000..edf0eea --- /dev/null +++ b/submitter/clean.sh @@ -0,0 +1,5 @@ +#!/bin/bash + +rm -f dagman_FR.submit.* +rm -f logs/* +rm -f metaouts/* diff --git a/submitter/make_dag.py b/submitter/make_dag.py index 735d213..641e00e 100644 --- a/submitter/make_dag.py +++ b/submitter/make_dag.py @@ -16,30 +16,30 @@ full_scan_mfr = [ ] 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) + # (1, 1, 0, 1, 2, 0), + # (1, 2, 0, 1, 2, 0), + # (1, 1, 1, 1, 0, 0), + # (1, 1, 0, 1, 0, 0), + # (1, 0, 0, 1, 0, 0), + # (1, 1, 1, 0, 1, 0), + # (1, 1, 0, 0, 1, 0), + # (0, 1, 0, 0, 1, 0), + # (1, 2, 0, 0, 1, 0), + # (1, 1, 1, 0, 0, 1), ] # MCMC -run_mcmc = 'True' +run_mcmc = 'False' burnin = 500 nsteps = 2000 nwalkers = 60 seed = 24 -threads = 4 +threads = 12 mcmc_seed_type = 'uniform' # FR -dimension = [4, 5, 7, 8] +dimension = [3, 6] energy = [1e6] likelihood = 'golemfit' no_bsm = 'False' @@ -66,9 +66,15 @@ promptNorm = 0. ast = 'p2_0' data = 'real' +# Bayes Factor +run_bayes_factor = 'True' +bayes_bins = 10 +bayes_live_points = 200 + # Plot -plot_angles = 'True' +plot_angles = 'False' plot_elements = 'False' +plot_bayes = 'True' outfile = 'dagman_FR.submit' golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit' @@ -86,11 +92,15 @@ with open(outfile, 'w') as f: elif energy_dependance == 'spectral': outchain_head = '/data/user/smandalia/flavour_ratio/data/{0}/DIM{1}/SI_{2}'.format(likelihood, dim, spectral_index) + bayes_output = 'None' for sig in sigma_ratio: print 'sigma', sig for frs in fix_sfr_mfr: print frs - outchains = outchain_head + '/fix_ifr/{0}/mcmc_chain'.format(str(sig).replace('.', '_')) + outchains = outchain_head + '/fix_ifr/{0}/'.format(str(sig).replace('.', '_')) + if run_bayes_factor == 'True': + bayes_output = outchains + '/bayes_factor/' + outchains += 'mcmc_chain' 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])) @@ -132,11 +142,19 @@ with open(outfile, 'w') as f: f.write('VARS\tjob{0}\tbinning_1="{1}"\n'.format(job_number, binning[1])) f.write('VARS\tjob{0}\tbinning_2="{1}"\n'.format(job_number, binning[2])) f.write('VARS\tjob{0}\tfix_mixing_almost="{1}"\n'.format(job_number, fix_mixing_almost)) + f.write('VARS\tjob{0}\trun_bayes_factor="{1}"\n'.format(job_number, run_bayes_factor)) + f.write('VARS\tjob{0}\tbayes_bins="{1}"\n'.format(job_number, bayes_bins)) + f.write('VARS\tjob{0}\tbayes_output="{1}"\n'.format(job_number, bayes_output)) + f.write('VARS\tjob{0}\tbayes_live_points="{1}"\n'.format(job_number, bayes_live_points)) + f.write('VARS\tjob{0}\tplot_bayes="{1}"\n'.format(job_number, plot_bayes)) job_number += 1 # for frs in full_scan_mfr: # print frs - # outchains = outchain_head + '/full_scan/{0}/mcmc_chain'.format(str(sig).replace('.', '_')) + # outchains = outchain_head + '/full_scan/{0}'.format(str(sig).replace('.', '_')) + # if run_bayes_factor == 'True': + # bayes_output = outchains + '/bayes_factor/' + # outchains += 'mcmc_chain' # 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])) @@ -178,4 +196,9 @@ with open(outfile, 'w') as f: # f.write('VARS\tjob{0}\tbinning_1="{1}"\n'.format(job_number, binning[1])) # f.write('VARS\tjob{0}\tbinning_2="{1}"\n'.format(job_number, binning[2])) # f.write('VARS\tjob{0}\tfix_mixing_almost="{1}"\n'.format(job_number, fix_mixing_almost)) + # f.write('VARS\tjob{0}\trun_bayes_factor="{1}"\n'.format(job_number, run_bayes_factor)) + # f.write('VARS\tjob{0}\tbayes_bins="{1}"\n'.format(job_number, bayes_bins)) + # f.write('VARS\tjob{0}\tbayes_output="{1}"\n'.format(job_number, bayes_output)) + # f.write('VARS\tjob{0}\tbayes_live_points="{1}"\n'.format(job_number, bayes_live_points)) + # f.write('VARS\tjob{0}\tplot_bayes="{1}"\n'.format(job_number, plot_bayes)) # job_number += 1 diff --git a/submitter/submit.sub b/submitter/submit.sub index b984c89..e563924 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 = "--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) --energy-dependance $(energy_dependance) --spectral-index $(spectral_index) --binning $(binning_0) $(binning_1) $(binning_2) --fix-mixing-almost $(fix_mixing_almost)" +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) --energy-dependance $(energy_dependance) --spectral-index $(spectral_index) --binning $(binning_0) $(binning_1) $(binning_2) --fix-mixing-almost $(fix_mixing_almost) --run-bayes-factor $(run_bayes_factor) --bayes-bins $(bayes_bins) --bayes-output $(bayes_output) --bayes-live-points $(bayes_live_points) --plot-bayes $(plot_bayes)" # All logs will go to a single file log = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/logs/job_$(Cluster).log @@ -14,9 +14,10 @@ getenv = True # but do not try to copy outputs back (see: https://htcondor-wiki.cs.wisc.edu/index.cgi/tktview?tn=3081) # +TransferOutput="" +transfer_output_files = /data/user/smandalia/GolemTools/sources/GolemFit/scripts/flavour_ratio/submitter/metaouts/ -request_memory = 5GB -request_cpus = 4 +request_memory = 30GB +request_cpus = 12 Universe = vanilla Notification = never diff --git a/utils/gf.py b/utils/gf.py index e575094..3f770eb 100644 --- a/utils/gf.py +++ b/utils/gf.py @@ -42,6 +42,15 @@ def set_up_as(fitter, params): fitter.SetupAsimov(asimov_params) +def setup_fitter(args, asimov_paramset): + datapaths = gf.DataPaths() + sparams = steering_params(args) + npp = gf.NewPhysicsParams() + fitter = gf.GolemFit(datapaths, sparams, npp) + set_up_as(fitter, asimov_paramset) + return fitter + + def get_llh(fitter, params): fitparams = gf.FitParameters(gf.sampleTag.HESE) for parm in params: diff --git a/utils/plot.py b/utils/plot.py index 4c5ff1c..4180eb3 100644 --- a/utils/plot.py +++ b/utils/plot.py @@ -15,6 +15,8 @@ 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 import getdist from getdist import plots @@ -32,6 +34,15 @@ def centers(x): return (x[:-1]+x[1:])*0.5 +def get_units(dimension): + if dimension == 3: return r' / GeV' + if dimension == 4: return r'' + if dimension == 5: return r' / GeV^{-1}' + if dimension == 6: return r' / GeV^{-2}' + if dimension == 7: return r' / GeV^{-3}' + if dimension == 8: return r' / GeV^{-4}' + + def calc_nbins(x): n = (np.max(x) - np.min(x)) / (2 * len(x)**(-1./3) * (np.percentile(x, 75) - np.percentile(x, 25))) return np.floor(n) @@ -84,6 +95,10 @@ def plot_argparse(parser): '--plot-elements', type=misc_utils.parse_bool, default='False', help='Plot MCMC triangle in the mixing elements space' ) + parser.add_argument( + '--plot-bayes', type=misc_utils.parse_bool, default='False', + help='Plot Bayes factor' + ) def flat_angles_to_u(x): @@ -139,7 +154,7 @@ def gen_figtext(args): mr1, mr2, mr3, args.dimension, int(args.energy), args.scale ) - else: + else: t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \ '{2:.2f}]\nDimension = {3}'.format( mr1, mr2, mr3, args.dimension, int(args.energy) @@ -242,3 +257,99 @@ def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): mpl.pyplot.figtext(0.5, 0.7, fig_text, fontsize=15) for of in outformat: g.export(outfile+'_elements.'+of) + + +def bayes_factor_plot(infile, outfile, outformat, args): + """Make Bayes factor plot.""" + if not args.plot_bayes: return + print "Making Bayes Factor plot" + fig_text = gen_figtext(args) + + raw = np.load(infile) + scales, evidences = raw.T + null = evidences[0] + + reduced_ev = -(evidences - null) + + fig = plt.figure(figsize=(7, 5)) + ax = fig.add_subplot(111) + + ax.set_xlabel(r'{\rm log}_{10} \Lambda ' + get_units(args.dimension)) + ax.set_ylabel(r'Bayes Factor') + + ax.plot(scales, reduced_ev) + + for ymaj in ax.yaxis.get_majorticklocs(): + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.7, linewidth=1) + for xmaj in ax.xaxis.get_majorticklocs(): + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.7, linewidth=1) + + at = AnchoredText( + fig_text, prop=dict(size=7), frameon=True, loc=2 + ) + at.patch.set_boxstyle("round,pad=0.,rounding_size=0.5") + ax.add_artist(at) + for of in outformat: + fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150) + + +def myround(x, base=5, up=False, down=False): + if up == down and up is True: assert 0 + if up: return int(base * np.round(float(x)/base-0.5)) + elif down: return int(base * np.round(float(x)/base+0.5)) + else: int(base * np.round(float(x)/base)) + + +def plot_BSM_angles_limit(infile, outfile, xticks, outformat, args): + """Make BSM angles vs scale limit plot.""" + print "Making BSM angles limit plot.""" + fig_text = gen_figtext(args) + + raw = np.load(infile) + sc_ranges = ( + myround(np.min(raw[0][:,0]), down=True), + myround(np.max(raw[0][:,0]), up=True) + ) + + proc = [] + for idx, theta in enumerate(raw): + scale, llh = theta.T + delta_llh = -2 * (llh - np.max(llh)) + # 90% CL for 1 dof + al = scale[delta_llh > 2.71] + proc.append((idx+1, al[0])) + + limits = np.array(proc) + + fig = plt.figure(figsize=(7, 5)) + ax = fig.add_subplot(111) + + ax.set_xticklabels(['']+xticks+['']) + ax.set_xlim(0, len(xticks)+1) + ax.set_ylim(sc_ranges[0], sc_ranges[-1]) + ax.set_xlabel(r'BSM angle') + ylabel = r'${\rm log}_{10} \Lambda' + get_units(args.dimension) + r'$' + ax.set_ylabel(ylabel) + + for l in limits: + line = plt.Line2D( + (l[0]-0.1, l[0]+0.1), (l[1], l[1]), lw=3, color='r' + ) + ax.add_line(line) + # ax.arrow( + # l[0], l[1], 0, -1.5, head_width=0.05, head_length=0.2, fc='r', + # ec='r', lw=2 + # ) + ax.annotate( + s='', xy=l, xytext=(l[0], l[1]-1.5), + arrowprops={'arrowstyle': '<-', 'lw': 1.5, 'color':'r'} + ) + + for ymaj in ax.yaxis.get_majorticklocs(): + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.4, linewidth=1) + for xmaj in ax.xaxis.get_majorticklocs(): + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.4, linewidth=1) + + for of in outformat: + fig.savefig(outfile+'.'+of, bbox_inches='tight', dpi=150) + |
