diff options
| -rwxr-xr-x | fr.py | 40 | ||||
| -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 |
7 files changed, 141 insertions, 54 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/submitter/make_dag.py b/submitter/make_dag.py index 63a94d7..f5a06dd 100644 --- a/submitter/make_dag.py +++ b/submitter/make_dag.py @@ -39,13 +39,16 @@ threads = 12 mcmc_seed_type = 'uniform' # FR -dimension = [3] -energy = [1e6] -likelihood = 'gaussian' -no_bsm = 'False' -sigma_ratio = ['0.01'] -scale = "1E-20 1E-30" -scale_region = "1E10" +dimension = [3] +energy = [1e6] +likelihood = 'gaussian' +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 fe58eef..f7c4e7e 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 5c9af51..9bba79a 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 a24c69c..a659b78 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) @@ -68,40 +68,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.""" |