diff options
| -rwxr-xr-x | fr.py | 40 | ||||
| -rw-r--r-- | plot_llh/LVCPT.py | 26 | ||||
| -rw-r--r-- | plot_llh/make_plots.py | 4 | ||||
| -rw-r--r-- | submitter/make_dag.py | 27 | ||||
| -rw-r--r-- | submitter/submit.sub | 2 | ||||
| -rw-r--r-- | utils/enums.py | 10 | ||||
| -rw-r--r-- | utils/fr.py | 1 | ||||
| -rw-r--r-- | utils/likelihood.py | 66 | ||||
| -rw-r--r-- | utils/plot.py | 49 |
9 files changed, 162 insertions, 63 deletions
@@ -21,9 +21,9 @@ 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, MCMCSeedType, ParamTag +from utils.enums import EnergyDependance, Likelihood, MCMCSeedType, ParamTag from utils.fr import MASS_EIGENVALUES, normalise_fr -from utils.misc import Param, ParamSet +from utils.misc import enum_parse, Param, ParamSet from utils.plot import plot_argparse, chainer_plot @@ -108,12 +108,25 @@ def process_args(args): if args.fix_source_ratio: args.source_ratio = normalise_fr(args.source_ratio) - if not args.fix_scale: - args.scale = np.power( - 10, np.round(np.log10(MASS_EIGENVALUES[1]/args.energy)) - \ - np.log10(args.energy**(args.dimension-3)) + if args.energy_dependance is EnergyDependance.SPECTRAL: + args.binning = np.logspace( + np.log10(args.binning[0]), np.log10(args.binning[1]), args.binning[2]+1 ) - """Estimate the scale of when to expect the BSM term to contribute.""" + + if not args.fix_scale: + if args.energy_dependance is EnergyDependance.MONO: + args.scale = np.power( + 10, np.round(np.log10(MASS_EIGENVALUES[1]/args.energy)) - \ + np.log10(args.energy**(args.dimension-3)) + ) + elif args.energy_dependance is EnergyDependance.SPECTRAL: + args.scale = np.power( + 10, np.round( + np.log10(MASS_EIGENVALUES[1]/np.power(10, np.average(np.log10(args.binning)))) \ + - np.log10(np.power(10, np.average(np.log10(args.binning)))**(args.dimension-3)) + ) + ) + """Estimate the scale of when to expect the BSM term to contribute.""" def parse_args(): @@ -135,6 +148,19 @@ def parse_args(): help='Fix the source flavour ratio' ) parser.add_argument( + '--energy-dependance', default='spectral', type=partial(enum_parse, c=EnergyDependance), + choices=EnergyDependance, + help='Type of energy dependance to use in the BSM fit' + ) + parser.add_argument( + '--spectral-index', default=-2, type=int, + help='Spectral index for spectral energy dependance' + ) + parser.add_argument( + '--binning', default=[4, 7, 50], type=int, nargs=3, + help='Binning for spectral energy dependance' + ) + parser.add_argument( '--source-ratio', type=int, nargs=3, default=[2, 1, 0], help='Set the source flavour ratio for the case when you want to fix it' ) diff --git a/plot_llh/LVCPT.py b/plot_llh/LVCPT.py index f5b02e6..e55c3b4 100644 --- a/plot_llh/LVCPT.py +++ b/plot_llh/LVCPT.py @@ -127,11 +127,11 @@ def MakeFlavorTriangle(list_of_flavor_ratios, scale = 8, # horizontalalignment="center",fontsize = 50, zorder = 10, rotation = -60) ax.text((p+s0*0.5+s0*0.025)[0], (p+s0*0.5-np.array([0,0.15*scale]))[1],r"$f_{e,\oplus}$", - horizontalalignment="center",fontsize = 50, zorder = 10) - ax.text((p+s1*0.5-0.2*s0)[0], (p+s1*0.5+0.1*s0)[1]+scale*0.1,r"$f_{\tau, \oplus}$", horizontalalignment="center",fontsize = 50, zorder = 10, rotation = 60.) + ax.text((p+s1*0.5-0.2*s0)[0], (p+s1*0.5+0.1*s0)[1]+scale*0.1,r"$f_{\tau, \oplus}$", + horizontalalignment="center",fontsize = 50, zorder = 10, rotation = -60.) ax.text((p+s1*0.5 + 0.7*s0)[0], (p+s1*0.5 + 0.6*s0)[1]+0.05*scale,r"$f_{\mu, \oplus}$", - horizontalalignment="center",fontsize = 50, zorder = 10, rotation = -60) + horizontalalignment="center",fontsize = 50, zorder = 10) # construct triangle grid fsl = 35 @@ -155,6 +155,11 @@ def MakeFlavorTriangle(list_of_flavor_ratios, scale = 8, plt.plot(*PointToList(p+s0*subsize*float(divisions-(i-0.5)), p+s1-s1*subsize*float(i-0.5)), color = inner_line_color, lw = slw, ls = "dotted", zorder = -1) plt.plot(*PointToList(p+s1*subsize*float(divisions-(i-0.5)), p+s1+s2*subsize*float(i-0.5)), color = inner_line_color, lw = slw, ls = "dotted", zorder = -1) + # chi2 + if (triangle_collection != None): + total_points = float(sum([ triangle.number_of_points for triangle in triangle_collection])) + max_points = float(max([ triangle.number_of_points for triangle in triangle_collection])) + delta_llh = [-2*(np.log10(triangle.number_of_points) - np.log10(max_points)) for triangle in triangle_collection] # plot triangle collection if (triangle_collection != None): @@ -162,14 +167,18 @@ def MakeFlavorTriangle(list_of_flavor_ratios, scale = 8, total_points = float(sum([ triangle.number_of_points for triangle in triangle_collection])) max_points = float(max([ triangle.number_of_points for triangle in triangle_collection])) color_map = plt.get_cmap(icolormap) - for triangle in triangle_collection: + for i_triangle, triangle in enumerate(triangle_collection): if triangle.number_of_points > 0: xx,yy = zip(*triangle.coordinates) if color_scale == "lin": plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(0.75), alpha = np.sqrt(float(triangle.number_of_points)/max_points)) #plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(float(triangle.number_of_points)/max_points), alpha = alpha) elif color_scale == "log": - plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(0.75), alpha = (np.log10(float(triangle.number_of_points))/np.log10(max_points))) + # plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(0.75), alpha = (np.log10(float(triangle.number_of_points))/np.log10(max_points))) + if delta_llh[i_triangle] > 2.30: + plt.fill(xx,yy,lw = 0., zorder = -0.8, color = 'blue', alpha = np.sqrt(float(triangle.number_of_points)/max_points)) + else: + plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(0.75), alpha = np.sqrt(float(triangle.number_of_points)/max_points)) #plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(0.7), alpha = (np.log10(float(triangle.number_of_points))/np.log10(max_points))**(2./3.)) #plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(np.log10(float(triangle.number_of_points))/np.log10(max_points)), alpha = alpha) #plt.fill(xx,yy,lw = 0., zorder = -0.8, color = color_map(np.log10(float(triangle.number_of_points)))) @@ -186,15 +195,18 @@ def MakeFlavorTriangle(list_of_flavor_ratios, scale = 8, cbaxes = fig.add_axes([left,bottom,width,height]) if color_scale == "lin": norm = mpl.colors.Normalize(vmin = 0., vmax = max_points) + label_format='%.0e' elif color_scale == "log": - norm = mpl.colors.Normalize(vmin = 0., vmax = 1.0) + # norm = mpl.colors.Normalize(vmin = 0., vmax = 1.0) + norm = mpl.colors.LogNorm(vmin = 1., vmax = max_points) + label_format=None else : raise NameError('Error. Love CA.') mpl.rcParams.update({'font.size': 10}) triangle_colorbar = mpl.colorbar.ColorbarBase(cbaxes, cmap = color_map, norm = norm, orientation = "horizontal", spacing = "proportional", # ) - format ='%.0e') + format =label_format) cbaxes.set_xlabel("Likelihood", fontsize = 12) # plot flavor ratio points diff --git a/plot_llh/make_plots.py b/plot_llh/make_plots.py index 9216396..8cf8318 100644 --- a/plot_llh/make_plots.py +++ b/plot_llh/make_plots.py @@ -76,7 +76,7 @@ for i in range(len(sys.argv)-1): figure = LVCPT.MakeFlavorTriangle(flavor_list, divisions = 5, save_file=True, filename = filename + "_hist", icolor = "g", icolormap = colors_schemes[i], triangle_collection=triangle_collection, figure = figure, alpha = 0.8, - initial_flavor_ratio = [0,1,0], subdivisions = True, color_scale = "lin", + initial_flavor_ratio = [0,1,0], subdivisions = True, color_scale = "log", output_format = output_format, inner_line_color ="silver",add_default_text = False, plot_color_bar =True) @@ -84,7 +84,7 @@ for i in range(len(sys.argv)-1): figure = LVCPT.MakeFlavorTriangle(flavor_list, divisions = 5, save_file=True, filename = filename + "_hist", icolor = "g", icolormap = colors_schemes[i], triangle_collection=triangle_collection, alpha = 0.8, - initial_flavor_ratio = [0,1,0], subdivisions = True, color_scale = "lin", + initial_flavor_ratio = [0,1,0], subdivisions = True, color_scale = "log", output_format = output_format, inner_line_color = "silver",add_default_text = False, plot_color_bar =True) diff --git a/submitter/make_dag.py b/submitter/make_dag.py index 106a1b7..216ca72 100644 --- a/submitter/make_dag.py +++ b/submitter/make_dag.py @@ -39,13 +39,16 @@ threads = 12 mcmc_seed_type = 'uniform' # FR -dimension = [3, 6] -energy = [1e6] -likelihood = 'golemfit' -no_bsm = 'False' -sigma_ratio = ['0.01'] -scale = "1E-20 1E-30" -scale_region = "1E10" +dimension = [3, 6] +energy = [1e6] +likelihood = 'golemfit' +no_bsm = 'False' +sigma_ratio = ['0.01'] +scale = "1E-20 1E-30" +scale_region = "1E10" +energy_dependance = 'spectral' +spectral_index = -2 +binning = [4, 7, 51] # Likelihood likelihood = 'golemfit' @@ -118,6 +121,11 @@ with open(outfile, 'w') as f: 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)) + f.write('VARS\tjob{0}\tenergy_dependance="{1}"\n'.format(job_number, energy_dependance)) + f.write('VARS\tjob{0}\tspectral_index="{1}"\n'.format(job_number, spectral_index)) + f.write('VARS\tjob{0}\tbinning_0="{1}"\n'.format(job_number, binning[0])) + 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])) job_number += 1 for frs in full_scan_mfr: @@ -158,4 +166,9 @@ with open(outfile, 'w') as f: 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)) + f.write('VARS\tjob{0}\tenergy_dependance="{1}"\n'.format(job_number, energy_dependance)) + f.write('VARS\tjob{0}\tspectral_index="{1}"\n'.format(job_number, spectral_index)) + f.write('VARS\tjob{0}\tbinning_0="{1}"\n'.format(job_number, binning[0])) + 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])) job_number += 1 diff --git a/submitter/submit.sub b/submitter/submit.sub index 51a97f6..b47ac5c 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)" +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)" # 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/utils/enums.py b/utils/enums.py index a5875a8..4637eeb 100644 --- a/utils/enums.py +++ b/utils/enums.py @@ -16,6 +16,11 @@ class DataType(Enum): ASMIMOV = 3 +class EnergyDependance(Enum): + MONO = 1 + SPECTRAL = 2 + + class FitPriorsCateg(Enum): DEFAULT = 1 XS = 2 @@ -42,6 +47,5 @@ class MCMCSeedType(Enum): class SteeringCateg(Enum): - P2_0 = 1 - P2_1 = 2 - + P2_0 = 1 + P2_1 = 2 diff --git a/utils/fr.py b/utils/fr.py index ddcb5d2..4242e97 100644 --- a/utils/fr.py +++ b/utils/fr.py @@ -344,7 +344,6 @@ def u_to_fr(source_fr, matrix): array([ 0.33740075, 0.33176584, 0.33083341]) """ - # TODO(shivesh): energy dependence composition = np.einsum( 'ai, bi, a -> b', abs(matrix)**2, abs(matrix)**2, source_fr ) diff --git a/utils/likelihood.py b/utils/likelihood.py index 599ff8b..f590200 100644 --- a/utils/likelihood.py +++ b/utils/likelihood.py @@ -59,26 +59,64 @@ def triangle_llh(theta, args, asimov_paramset, mcmc_paramset, fitter): for param in mcmc_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.fix_source_ratio: - fr1, fr2, fr3 = args.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, SPECTRAL_INDEX) + for fr in args.source_ratio] + ).T else: - fr1, fr2, fr3 = fr_utils.angles_to_fr( - mcmc_paramset.from_tag(ParamTag.SRCANGLES, values=True) - ) + if args.energy_dependance is EnergyDependance.MONO: + source_flux = fr_utils.angles_to_fr( + mcmc_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 angles_to_fr(theta[-2:])] + ).T + bsm_angles = mcmc_paramset.from_tag( [ParamTag.SCALE, ParamTag.MMANGLES], values=True ) - u = fr_utils.params_to_BSMu( - theta = bsm_angles, - dim = args.dimension, - energy = args.energy, - no_bsm = args.no_bsm, - fix_mixing = args.fix_mixing, - fix_scale = args.fix_scale, - scale = args.scale - ) - fr = fr_utils.u_to_fr((fr1, fr2, fr3), u) + if args.energy_dependance is EnergyDependance.MONO: + u = fr_utils.params_to_BSMu( + theta = bsm_angles, + dim = args.dimension, + energy = args.energy, + no_bsm = args.no_bsm, + fix_mixing = args.fix_mixing, + 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 = args.energy, + no_bsm = args.no_bsm, + fix_mixing = args.fix_mixing, + 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) + for idx, param in enumerate(hypo_paramset.from_tag(ParamTag.BESTFIT)): param.value = fr[idx] diff --git a/utils/plot.py b/utils/plot.py index 968d10f..87a6f5c 100644 --- a/utils/plot.py +++ b/utils/plot.py @@ -21,7 +21,7 @@ from getdist import plots from getdist import mcsamples from utils import misc as misc_utils -from utils.enums import ParamTag +from utils.enums import Likelihood, ParamTag from utils.fr import angles_to_u, angles_to_fr rc('text', usetex=False) @@ -105,40 +105,47 @@ def plot_Tchain(Tchain, axes_labels, ranges): def gen_figtext(args): """Generate the figure text.""" + t = '' mr1, mr2, mr3 = args.measured_ratio if args.fix_source_ratio: sr1, sr2, sr3 = args.source_ratio if args.fix_scale: - return 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \ + t += 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \ 'observed flavour ratio = [{3:.2f}, {4:.2f}, ' \ - '{5:.2f}]\nSigma = {6:.3f}\nDimension = {7}\nEnergy = ' \ - '{8} GeV\nScale = {9}'.format( - sr1, sr2, sr3, mr1, mr2, mr3, args.sigma_ratio, - args.dimension, int(args.energy), args.scale + '{5:.2f}]\nDimension = {7}\nScale = {9}'.format( + sr1, sr2, sr3, mr1, mr2, mr3, args.dimension, + int(args.energy), args.scale ) else: - return 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \ + t += 'Source flavour ratio = [{0:.2f}, {1:.2f}, {2:.2f}]\nIC ' \ 'observed flavour ratio = [{3:.2f}, {4:.2f}, ' \ - '{5:.2f}]\nSigma = {6:.3f}\nDimension = {7}\nEnergy = {8} ' \ - 'GeV'.format( - sr1, sr2, sr3, mr1, mr2, mr3, args.sigma_ratio, - args.dimension, int(args.energy) + '{5:.2f}]\nDimension = {7}'.format( + sr1, sr2, sr3, mr1, mr2, mr3, args.dimension, + int(args.energy) ) else: if args.fix_scale: - return 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \ - '{2:.2f}]\nSigma = {3:.3f}\nDimension = {4}\nEnergy = {5} ' \ - 'GeV\nScale = {6}'.format( - mr1, mr2, mr3, args.sigma_ratio, args.dimension, - int(args.energy), args.scale + t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \ + '{2:.2f}]\nDimension = {4}\nScale = {6}'.format( + mr1, mr2, mr3, args.dimension, int(args.energy), + args.scale ) else: - return 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \ - '{2:.2f}]\nSigma = {3:.3f}\nDimension = {4}\nEnergy = {5} ' \ - 'GeV'.format( - mr1, mr2, mr3, args.sigma_ratio, args.dimension, - int(args.energy) + t += 'IC observed flavour ratio = [{0:.2f}, {1:.2f}, ' \ + '{2:.2f}]\nDimension = {4}'.format( + mr1, mr2, mr3, args.dimension, int(args.energy) ) + if args.likelihood is Likelihood.GAUSSIAN: + t += '\nSigma = {0:.3f}'.format(args.sigma_ratio) + if args.energy_dependance is EnergyDependance.SPECTRAL: + t += '\nSpectral Index = {0}\nBinning = [{7}, {8}] GeV - {9} bins'.format( + int(args.spectral_index), int(10**args.binning[0]), + int(10**args.binning[1]), int(args.binning[2]) + ) + elif args.energy_dependance is EnergyDependance.MONO: + t += '\nEnergy = {0} GeV'.format(args.energy) + return t + def chainer_plot(infile, outfile, outformat, args, mcmc_paramset): """Make the triangle plot.""" |