From 5bcc51c59700968f8dad4359a133ab6a64f85f01 Mon Sep 17 00:00:00 2001 From: shivesh Date: Mon, 26 Nov 2018 20:43:08 -0600 Subject: begin contour.py and rename some util files --- utils/fr.py | 7 +- utils/likelihood.py | 212 ---------------------------------------------------- utils/llh.py | 212 ++++++++++++++++++++++++++++++++++++++++++++++++++++ utils/mn.py | 109 +++++++++++++++++++++++++++ utils/multinest.py | 109 --------------------------- utils/param.py | 68 +++++++++-------- 6 files changed, 362 insertions(+), 355 deletions(-) delete mode 100644 utils/likelihood.py create mode 100644 utils/llh.py create mode 100644 utils/mn.py delete mode 100644 utils/multinest.py (limited to 'utils') diff --git a/utils/fr.py b/utils/fr.py index d899d42..7c05f85 100644 --- a/utils/fr.py +++ b/utils/fr.py @@ -347,8 +347,11 @@ def fr_to_angles(ratios): cphi2 = fr2 sphi2 = (1.0 - cphi2) - spsi2 = fr1 / sphi2 - cpsi2 = fr0 / sphi2 + if sphi2 == 0.: + return (0., 0.) + else: + spsi2 = fr1 / sphi2 + cpsi2 = fr0 / sphi2 sphi4 = sphi2**2 c2psi = COS(ACOS(SQRT(cpsi2))*2) diff --git a/utils/likelihood.py b/utils/likelihood.py deleted file mode 100644 index f4404c4..0000000 --- a/utils/likelihood.py +++ /dev/null @@ -1,212 +0,0 @@ -# author : S. Mandalia -# s.p.mandalia@qmul.ac.uk -# -# date : April 04, 2018 - -""" -Likelihood functions for the BSM flavour ratio analysis -""" - -from __future__ import absolute_import, division - -from functools import partial - -import numpy as np -from scipy.stats import multivariate_normal, rv_continuous - -from utils import fr as fr_utils -from utils import gf as gf_utils -from utils.enums import EnergyDependance, Likelihood, ParamTag, PriorsCateg -from utils.misc import enum_parse, gen_identifier, parse_bool - - -class Gaussian(rv_continuous): - """Gaussian for one dimension.""" - def _pdf(self, x, mu, sig): - return (1./np.sqrt(2*np.pi*sig**2))*np.exp(-((x-mu)**2)/(2*sig**2)) - - -def multi_gaussian(fr, fr_bf, sigma, offset=-320): - """Multivariate gaussian likelihood.""" - cov_fr = np.identity(3) * sigma - return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr)) + offset - - -def likelihood_argparse(parser): - parser.add_argument( - '--likelihood', default='gaussian', type=partial(enum_parse, c=Likelihood), - choices=Likelihood, help='likelihood contour' - ) - parser.add_argument( - '--sigma-ratio', type=float, default=0.01, - help='Set the 1 sigma for the measured flavour ratio for a gaussian LLH' - ) - parser.add_argument( - '--save-measured-fr', type=parse_bool, default='False', - help='Output the measured flavour ratios' - ) - parser.add_argument( - '--output-measured-fr', type=str, default='./frs', - help='Output of the measured flavour ratios' - ) - - -def lnprior(theta, paramset): - """Priors on theta.""" - if len(theta) != len(paramset): - raise AssertionError( - 'Length of MCMC scan is not the same as the input ' - 'params\ntheta={0}\nparamset={1}'.format(theta, paramset) - ) - for idx, param in enumerate(paramset): - param.value = theta[idx] - ranges = paramset.ranges - for value, range in zip(theta, ranges): - if range[0] <= value <= range[1]: - pass - else: return -np.inf - - prior = 0 - for param in paramset: - if param.prior is PriorsCateg.GAUSSIAN: - prior += Gaussian().logpdf( - param.nominal_value, param.value, param.std - ) - elif param.prior is PriorsCateg.HALFGAUSS: - prior += Gaussian().logpdf( - param.nominal_value, param.value, param.std - ) + Gaussian().logcdf(1, param.value, param.std) - return prior - - -def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): - """Log likelihood function for a given theta.""" - if len(theta) != len(llh_paramset): - raise AssertionError( - 'Length of MCMC scan is not the same as the input ' - 'params\ntheta={0}\nparamset]{1}'.format(theta, llh_paramset) - ) - for idx, param in enumerate(llh_paramset): - param.value = theta[idx] - hypo_paramset = asimov_paramset - for param in llh_paramset.from_tag(ParamTag.NUISANCE): - hypo_paramset[param.name].value = param.value - - if args.energy_dependance is EnergyDependance.SPECTRAL: - bin_centers = np.sqrt(args.binning[:-1]*args.binning[1:]) - bin_width = np.abs(np.diff(args.binning)) - if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ] \ - and args.fold_index: - args.spectral_index = -hypo_paramset['astroDeltaGamma'].value - - if args.fix_source_ratio: - if args.energy_dependance is EnergyDependance.MONO: - source_flux = args.source_ratio - elif args.energy_dependance is EnergyDependance.SPECTRAL: - source_flux = np.array( - [fr * np.power(bin_centers, args.spectral_index) - for fr in args.source_ratio] - ).T - else: - if args.energy_dependance is EnergyDependance.MONO: - source_flux = fr_utils.angles_to_fr( - llh_paramset.from_tag(ParamTag.SRCANGLES, values=True) - ) - elif args.energy_dependance is EnergyDependance.SPECTRAL: - source_flux = np.array( - [fr * np.power(bin_centers, args.spectral_index) - for fr in fr_utils.angles_to_fr(theta[-2:])] - ).T - - bsm_angles = llh_paramset.from_tag( - [ParamTag.SCALE, ParamTag.MMANGLES], values=True - ) - - m_eig_names = ['m21_2', 'm3x_2'] - ma_names = ['s_12_2', 'c_13_4', 's_23_2', 'dcp'] - - if set(m_eig_names+ma_names).issubset(set(llh_paramset.names)): - mass_eigenvalues = [x.value for x in llh_paramset if x.name in m_eig_names] - sm_u = fr_utils.angles_to_u( - [x.value for x in llh_paramset if x.name in ma_names] - ) - else: - mass_eigenvalues = fr_utils.MASS_EIGENVALUES - sm_u = fr_utils.NUFIT_U - - if args.no_bsm: - fr = fr_utils.u_to_fr(source_flux, np.array(sm_u, dtype=np.complex256)) - elif args.energy_dependance is EnergyDependance.MONO: - u = fr_utils.params_to_BSMu( - theta = bsm_angles, - dim = args.dimension, - energy = args.energy, - mass_eigenvalues = mass_eigenvalues, - sm_u = sm_u, - no_bsm = args.no_bsm, - fix_mixing = args.fix_mixing, - fix_mixing_almost = args.fix_mixing_almost, - fix_scale = args.fix_scale, - scale = args.scale - ) - fr = fr_utils.u_to_fr(source_flux, u) - elif args.energy_dependance is EnergyDependance.SPECTRAL: - mf_perbin = [] - for i_sf, sf_perbin in enumerate(source_flux): - u = fr_utils.params_to_BSMu( - theta = bsm_angles, - dim = args.dimension, - energy = bin_centers[i_sf], - mass_eigenvalues = mass_eigenvalues, - sm_u = sm_u, - no_bsm = args.no_bsm, - fix_mixing = args.fix_mixing, - fix_mixing_almost = args.fix_mixing_almost, - fix_scale = args.fix_scale, - scale = args.scale - ) - fr = fr_utils.u_to_fr(sf_perbin, u) - mf_perbin.append(fr) - measured_flux = np.array(mf_perbin).T - intergrated_measured_flux = np.sum(measured_flux * bin_width, axis=1) - averaged_measured_flux = (1./(args.binning[-1] - args.binning[0])) * \ - intergrated_measured_flux - fr = averaged_measured_flux / np.sum(averaged_measured_flux) - - flavour_angles = fr_utils.fr_to_angles(fr) - # print 'flavour_angles', map(float, flavour_angles) - for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)): - param.value = flavour_angles[idx] - - if args.likelihood is Likelihood.FLAT: - llh = 1. - elif args.likelihood is Likelihood.GAUSSIAN: - fr_bf = args.measured_ratio - llh = multi_gaussian(fr, fr_bf, args.sigma_ratio) - elif args.likelihood is Likelihood.GOLEMFIT: - llh = gf_utils.get_llh(fitter, hypo_paramset) - elif args.likelihood is Likelihood.GF_FREQ: - llh = gf_utils.get_llh_freq(fitter, hypo_paramset) - - if args.save_measured_fr and args.burnin is False: - n = gen_identifier(args) + '.txt' - with open(args.output_measured_fr + n, 'a') as f: - f.write(r'{0:.3f} {1:.3f} {2:.3f}'.format( - float(fr[0]), float(fr[1]), float(fr[2]) - )) - for p in llh_paramset: - f.write(r' {0:.3f}'.format(p.value)) - f.write(' LLH = {0:.3f}'.format(llh)) - f.write('\n') - - return llh - - -def ln_prob(theta, args, asimov_paramset, llh_paramset, fitter): - lp = lnprior(theta, paramset=llh_paramset) - if not np.isfinite(lp): - return -np.inf - return lp + triangle_llh( - theta, args=args, asimov_paramset=asimov_paramset, - llh_paramset=llh_paramset, fitter=fitter - ) diff --git a/utils/llh.py b/utils/llh.py new file mode 100644 index 0000000..f4404c4 --- /dev/null +++ b/utils/llh.py @@ -0,0 +1,212 @@ +# author : S. Mandalia +# s.p.mandalia@qmul.ac.uk +# +# date : April 04, 2018 + +""" +Likelihood functions for the BSM flavour ratio analysis +""" + +from __future__ import absolute_import, division + +from functools import partial + +import numpy as np +from scipy.stats import multivariate_normal, rv_continuous + +from utils import fr as fr_utils +from utils import gf as gf_utils +from utils.enums import EnergyDependance, Likelihood, ParamTag, PriorsCateg +from utils.misc import enum_parse, gen_identifier, parse_bool + + +class Gaussian(rv_continuous): + """Gaussian for one dimension.""" + def _pdf(self, x, mu, sig): + return (1./np.sqrt(2*np.pi*sig**2))*np.exp(-((x-mu)**2)/(2*sig**2)) + + +def multi_gaussian(fr, fr_bf, sigma, offset=-320): + """Multivariate gaussian likelihood.""" + cov_fr = np.identity(3) * sigma + return np.log(multivariate_normal.pdf(fr, mean=fr_bf, cov=cov_fr)) + offset + + +def likelihood_argparse(parser): + parser.add_argument( + '--likelihood', default='gaussian', type=partial(enum_parse, c=Likelihood), + choices=Likelihood, help='likelihood contour' + ) + parser.add_argument( + '--sigma-ratio', type=float, default=0.01, + help='Set the 1 sigma for the measured flavour ratio for a gaussian LLH' + ) + parser.add_argument( + '--save-measured-fr', type=parse_bool, default='False', + help='Output the measured flavour ratios' + ) + parser.add_argument( + '--output-measured-fr', type=str, default='./frs', + help='Output of the measured flavour ratios' + ) + + +def lnprior(theta, paramset): + """Priors on theta.""" + if len(theta) != len(paramset): + raise AssertionError( + 'Length of MCMC scan is not the same as the input ' + 'params\ntheta={0}\nparamset={1}'.format(theta, paramset) + ) + for idx, param in enumerate(paramset): + param.value = theta[idx] + ranges = paramset.ranges + for value, range in zip(theta, ranges): + if range[0] <= value <= range[1]: + pass + else: return -np.inf + + prior = 0 + for param in paramset: + if param.prior is PriorsCateg.GAUSSIAN: + prior += Gaussian().logpdf( + param.nominal_value, param.value, param.std + ) + elif param.prior is PriorsCateg.HALFGAUSS: + prior += Gaussian().logpdf( + param.nominal_value, param.value, param.std + ) + Gaussian().logcdf(1, param.value, param.std) + return prior + + +def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter): + """Log likelihood function for a given theta.""" + if len(theta) != len(llh_paramset): + raise AssertionError( + 'Length of MCMC scan is not the same as the input ' + 'params\ntheta={0}\nparamset]{1}'.format(theta, llh_paramset) + ) + for idx, param in enumerate(llh_paramset): + param.value = theta[idx] + hypo_paramset = asimov_paramset + for param in llh_paramset.from_tag(ParamTag.NUISANCE): + hypo_paramset[param.name].value = param.value + + if args.energy_dependance is EnergyDependance.SPECTRAL: + bin_centers = np.sqrt(args.binning[:-1]*args.binning[1:]) + bin_width = np.abs(np.diff(args.binning)) + if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ] \ + and args.fold_index: + args.spectral_index = -hypo_paramset['astroDeltaGamma'].value + + if args.fix_source_ratio: + if args.energy_dependance is EnergyDependance.MONO: + source_flux = args.source_ratio + elif args.energy_dependance is EnergyDependance.SPECTRAL: + source_flux = np.array( + [fr * np.power(bin_centers, args.spectral_index) + for fr in args.source_ratio] + ).T + else: + if args.energy_dependance is EnergyDependance.MONO: + source_flux = fr_utils.angles_to_fr( + llh_paramset.from_tag(ParamTag.SRCANGLES, values=True) + ) + elif args.energy_dependance is EnergyDependance.SPECTRAL: + source_flux = np.array( + [fr * np.power(bin_centers, args.spectral_index) + for fr in fr_utils.angles_to_fr(theta[-2:])] + ).T + + bsm_angles = llh_paramset.from_tag( + [ParamTag.SCALE, ParamTag.MMANGLES], values=True + ) + + m_eig_names = ['m21_2', 'm3x_2'] + ma_names = ['s_12_2', 'c_13_4', 's_23_2', 'dcp'] + + if set(m_eig_names+ma_names).issubset(set(llh_paramset.names)): + mass_eigenvalues = [x.value for x in llh_paramset if x.name in m_eig_names] + sm_u = fr_utils.angles_to_u( + [x.value for x in llh_paramset if x.name in ma_names] + ) + else: + mass_eigenvalues = fr_utils.MASS_EIGENVALUES + sm_u = fr_utils.NUFIT_U + + if args.no_bsm: + fr = fr_utils.u_to_fr(source_flux, np.array(sm_u, dtype=np.complex256)) + elif args.energy_dependance is EnergyDependance.MONO: + u = fr_utils.params_to_BSMu( + theta = bsm_angles, + dim = args.dimension, + energy = args.energy, + mass_eigenvalues = mass_eigenvalues, + sm_u = sm_u, + no_bsm = args.no_bsm, + fix_mixing = args.fix_mixing, + fix_mixing_almost = args.fix_mixing_almost, + fix_scale = args.fix_scale, + scale = args.scale + ) + fr = fr_utils.u_to_fr(source_flux, u) + elif args.energy_dependance is EnergyDependance.SPECTRAL: + mf_perbin = [] + for i_sf, sf_perbin in enumerate(source_flux): + u = fr_utils.params_to_BSMu( + theta = bsm_angles, + dim = args.dimension, + energy = bin_centers[i_sf], + mass_eigenvalues = mass_eigenvalues, + sm_u = sm_u, + no_bsm = args.no_bsm, + fix_mixing = args.fix_mixing, + fix_mixing_almost = args.fix_mixing_almost, + fix_scale = args.fix_scale, + scale = args.scale + ) + fr = fr_utils.u_to_fr(sf_perbin, u) + mf_perbin.append(fr) + measured_flux = np.array(mf_perbin).T + intergrated_measured_flux = np.sum(measured_flux * bin_width, axis=1) + averaged_measured_flux = (1./(args.binning[-1] - args.binning[0])) * \ + intergrated_measured_flux + fr = averaged_measured_flux / np.sum(averaged_measured_flux) + + flavour_angles = fr_utils.fr_to_angles(fr) + # print 'flavour_angles', map(float, flavour_angles) + for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)): + param.value = flavour_angles[idx] + + if args.likelihood is Likelihood.FLAT: + llh = 1. + elif args.likelihood is Likelihood.GAUSSIAN: + fr_bf = args.measured_ratio + llh = multi_gaussian(fr, fr_bf, args.sigma_ratio) + elif args.likelihood is Likelihood.GOLEMFIT: + llh = gf_utils.get_llh(fitter, hypo_paramset) + elif args.likelihood is Likelihood.GF_FREQ: + llh = gf_utils.get_llh_freq(fitter, hypo_paramset) + + if args.save_measured_fr and args.burnin is False: + n = gen_identifier(args) + '.txt' + with open(args.output_measured_fr + n, 'a') as f: + f.write(r'{0:.3f} {1:.3f} {2:.3f}'.format( + float(fr[0]), float(fr[1]), float(fr[2]) + )) + for p in llh_paramset: + f.write(r' {0:.3f}'.format(p.value)) + f.write(' LLH = {0:.3f}'.format(llh)) + f.write('\n') + + return llh + + +def ln_prob(theta, args, asimov_paramset, llh_paramset, fitter): + lp = lnprior(theta, paramset=llh_paramset) + if not np.isfinite(lp): + return -np.inf + return lp + triangle_llh( + theta, args=args, asimov_paramset=asimov_paramset, + llh_paramset=llh_paramset, fitter=fitter + ) diff --git a/utils/mn.py b/utils/mn.py new file mode 100644 index 0000000..48ccc26 --- /dev/null +++ b/utils/mn.py @@ -0,0 +1,109 @@ +# author : S. Mandalia +# s.p.mandalia@qmul.ac.uk +# +# date : April 19, 2018 + +""" +Useful functions to use MultiNest for the BSM flavour ratio analysis +""" + +from __future__ import absolute_import, division + +from functools import partial + +import numpy as np + +from pymultinest import analyse, run + +from utils import llh as llh_utils +from utils.misc import gen_identifier, make_dir, solve_ratio + + +def CubePrior(cube, ndim, n_params): + pass + + +def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset, + args, fitter): + if ndim != len(mn_paramset): + raise AssertionError( + 'Length of MultiNest scan paramset is not the same as the input ' + 'params\ncube={0}\nmn_paramset]{1}'.format(cube, mn_paramset) + ) + pranges = mn_paramset.seeds + for i in xrange(ndim): + mn_paramset[i].value = (pranges[i][1]-pranges[i][0])*cube[i] + pranges[i][0] + for pm in mn_paramset.names: + llh_paramset[pm].value = mn_paramset[pm].value + 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 + return llh + + +def mn_argparse(parser): + parser.add_argument( + '--mn-live-points', type=int, default=3000, + help='Number of live points for MultiNest evaluations' + ) + parser.add_argument( + '--mn-tolerance', type=float, default=0.01, + help='Tolerance for MultiNest' + ) + parser.add_argument( + '--mn-output', type=str, default='./mnrun/', + help='Folder to store MultiNest evaluations' + ) + + +def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter, + identifier='mn'): + """Run the MultiNest algorithm to calculate the evidence.""" + n_params = len(mn_paramset) + + for n in mn_paramset.names: + llh_paramset[n].value = mn_paramset[n].value + + lnProbEval = partial( + lnProb, + mn_paramset = mn_paramset, + llh_paramset = llh_paramset, + asimov_paramset = asimov_paramset, + args = args, + fitter = fitter + ) + + # prefix = './mnrun/DIM{0}/{1}_{2}_{3:>010}_'.format( + # args.dimension, args.likelihood, gen_identifier(args), np.random.randint(0, 2**30) + # ) + llh = '{0}'.format(args.likelihood).split('.')[1] + data = '{0}'.format(args.data).split('.')[1] + sr1, sr2, sr3 = solve_ratio(args.source_ratio) + prefix = './mnrun/DIM{0}/{1}/{2}/s{3}{4}{5}/{6}'.format( + args.dimension, data, llh, sr1, sr2, sr3, identifier + ) + make_dir(prefix) + print 'Running evidence calculation for {0}'.format(prefix) + run( + LogLikelihood = lnProbEval, + Prior = CubePrior, + n_dims = n_params, + n_live_points = args.mn_live_points, + evidence_tolerance = args.mn_tolerance, + outputfiles_basename = prefix, + importance_nested_sampling = True, + resume = False, + verbose = True + ) + + analyser = analyse.Analyzer( + outputfiles_basename=prefix, n_params=n_params + ) + return analyser.get_stats()['global evidence'] diff --git a/utils/multinest.py b/utils/multinest.py deleted file mode 100644 index fdd87cd..0000000 --- a/utils/multinest.py +++ /dev/null @@ -1,109 +0,0 @@ -# author : S. Mandalia -# s.p.mandalia@qmul.ac.uk -# -# date : April 19, 2018 - -""" -Useful functions to use MultiNest for the BSM flavour ratio analysis -""" - -from __future__ import absolute_import, division - -from functools import partial - -import numpy as np - -from pymultinest import analyse, run - -from utils import likelihood -from utils.misc import gen_identifier, make_dir, solve_ratio - - -def CubePrior(cube, ndim, n_params): - pass - - -def lnProb(cube, ndim, n_params, mn_paramset, llh_paramset, asimov_paramset, - args, fitter): - if ndim != len(mn_paramset): - raise AssertionError( - 'Length of MultiNest scan paramset is not the same as the input ' - 'params\ncube={0}\nmn_paramset]{1}'.format(cube, mn_paramset) - ) - pranges = mn_paramset.seeds - for i in xrange(ndim): - mn_paramset[i].value = (pranges[i][1]-pranges[i][0])*cube[i] + pranges[i][0] - for pm in mn_paramset.names: - llh_paramset[pm].value = mn_paramset[pm].value - theta = llh_paramset.values - # print 'llh_paramset', llh_paramset - llh = likelihood.ln_prob( - theta=theta, - args=args, - asimov_paramset=asimov_paramset, - llh_paramset=llh_paramset, - fitter=fitter - ) - # print 'llh', llh - return llh - - -def mn_argparse(parser): - parser.add_argument( - '--mn-live-points', type=int, default=3000, - help='Number of live points for MultiNest evaluations' - ) - parser.add_argument( - '--mn-tolerance', type=float, default=0.01, - help='Tolerance for MultiNest' - ) - parser.add_argument( - '--mn-output', type=str, default='./mnrun/', - help='Folder to store MultiNest evaluations' - ) - - -def mn_evidence(mn_paramset, llh_paramset, asimov_paramset, args, fitter, - identifier='mn'): - """Run the MultiNest algorithm to calculate the evidence.""" - n_params = len(mn_paramset) - - for n in mn_paramset.names: - llh_paramset[n].value = mn_paramset[n].value - - lnProbEval = partial( - lnProb, - mn_paramset = mn_paramset, - llh_paramset = llh_paramset, - asimov_paramset = asimov_paramset, - args = args, - fitter = fitter - ) - - # prefix = './mnrun/DIM{0}/{1}_{2}_{3:>010}_'.format( - # args.dimension, args.likelihood, gen_identifier(args), np.random.randint(0, 2**30) - # ) - llh = '{0}'.format(args.likelihood).split('.')[1] - data = '{0}'.format(args.data).split('.')[1] - sr1, sr2, sr3 = solve_ratio(args.source_ratio) - prefix = './mnrun/DIM{0}/{1}/{2}/s{3}{4}{5}/{6}'.format( - args.dimension, data, llh, sr1, sr2, sr3, identifier - ) - make_dir(prefix) - print 'Running evidence calculation for {0}'.format(prefix) - run( - LogLikelihood = lnProbEval, - Prior = CubePrior, - n_dims = n_params, - n_live_points = args.mn_live_points, - evidence_tolerance = args.mn_tolerance, - outputfiles_basename = prefix, - importance_nested_sampling = True, - resume = False, - verbose = True - ) - - analyser = analyse.Analyzer( - outputfiles_basename=prefix, n_params=n_params - ) - return analyser.get_stats()['global evidence'] diff --git a/utils/param.py b/utils/param.py index 62f59b1..8c17541 100644 --- a/utils/param.py +++ b/utils/param.py @@ -230,43 +230,47 @@ def get_paramsets(args, nuisance_paramset): for parm in llh_paramset: parm.value = args.__getattribute__(parm.name) tag = ParamTag.BESTFIT - flavour_angles = fr_to_angles(args.measured_ratio) + try: + flavour_angles = fr_to_angles(args.measured_ratio) + except: + flavour_angles = fr_to_angles(args.injected_ratio) asimov_paramset.extend([ Param(name='astroFlavorAngle1', value=flavour_angles[0], ranges=[0., 1.], std=0.2, tag=tag), Param(name='astroFlavorAngle2', value=flavour_angles[1], ranges=[-1., 1.], std=0.2, tag=tag), ]) asimov_paramset = ParamSet(asimov_paramset) - if args.fix_mixing is MixingScenario.NONE and not args.fix_mixing_almost: - tag = ParamTag.MMANGLES - llh_paramset.extend([ - Param(name='np_s_12^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag), - Param(name='np_c_13^4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag), - Param(name='np_s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^2', tag=tag), - Param(name='np_dcp', value=np.pi, ranges=[0., 2*np.pi], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag) - ]) - if args.fix_mixing_almost: - tag = ParamTag.MMANGLES - llh_paramset.extend([ - Param(name='np_s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag) - ]) - if not args.fix_scale: - tag = ParamTag.SCALE - if hasattr(args, 'dimension'): - llh_paramset.append( - Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, - tex=r'{\rm log}_{10}\left (\Lambda^{-1}'+get_units(args.dimension)+r'\right )', tag=tag) - ) - elif hasattr(args, 'dimensions'): - llh_paramset.append( - Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, - tex=r'{\rm log}_{10}\left (\Lambda^{-1} / GeV^{-d+4}\right )', tag=tag) - ) - if not args.fix_source_ratio: - tag = ParamTag.SRCANGLES - llh_paramset.extend([ - Param(name='s_phi4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'sin^4(\phi)', tag=tag), - Param(name='c_2psi', value=0.5, ranges=[0., 1.], std=0.2, tex=r'cos(2\psi)', tag=tag) - ]) + if hasattr(args, 'fix_source_ratio'): + if args.fix_mixing is MixingScenario.NONE and not args.fix_mixing_almost: + tag = ParamTag.MMANGLES + llh_paramset.extend([ + Param(name='np_s_12^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{12}^2', tag=tag), + Param(name='np_c_13^4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{c}_{13}^4', tag=tag), + Param(name='np_s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^2', tag=tag), + Param(name='np_dcp', value=np.pi, ranges=[0., 2*np.pi], std=0.2, tex=r'\tilde{\delta_{CP}}', tag=tag) + ]) + if args.fix_mixing_almost: + tag = ParamTag.MMANGLES + llh_paramset.extend([ + Param(name='np_s_23^2', value=0.5, ranges=[0., 1.], std=0.2, tex=r'\tilde{s}_{23}^4', tag=tag) + ]) + if not args.fix_scale: + tag = ParamTag.SCALE + if hasattr(args, 'dimension'): + llh_paramset.append( + Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, + tex=r'{\rm log}_{10}\left (\Lambda^{-1}'+get_units(args.dimension)+r'\right )', tag=tag) + ) + elif hasattr(args, 'dimensions'): + llh_paramset.append( + Param(name='logLam', value=np.log10(args.scale), ranges=np.log10(args.scale_region), std=3, + tex=r'{\rm log}_{10}\left (\Lambda^{-1} / GeV^{-d+4}\right )', tag=tag) + ) + if not args.fix_source_ratio: + tag = ParamTag.SRCANGLES + llh_paramset.extend([ + Param(name='s_phi4', value=0.5, ranges=[0., 1.], std=0.2, tex=r'sin^4(\phi)', tag=tag), + Param(name='c_2psi', value=0.5, ranges=[0., 1.], std=0.2, tex=r'cos(2\psi)', tag=tag) + ]) llh_paramset = ParamSet(llh_paramset) return asimov_paramset, llh_paramset -- cgit v1.2.3