From 932a8691e16eb904e3eec61daae08d72c2039f10 Mon Sep 17 00:00:00 2001 From: shivesh Date: Sun, 6 May 2018 21:22:47 -0500 Subject: Sun May 6 21:22:46 CDT 2018 --- test/test_LV.py | 96 ++++++++++++++++++++++++++++++++++++++++++++ test/test_NSI.py | 106 +++++++++++++++++++++++++++++++++++++++++++++++++ test/test_all_gf.py | 79 ++++++++++++++++++++++++++++++++++++ test/test_gf.py | 21 ++++++++-- test/test_gf_simple.py | 40 ++++++++++++++++++- 5 files changed, 337 insertions(+), 5 deletions(-) create mode 100644 test/test_LV.py create mode 100644 test/test_NSI.py create mode 100644 test/test_all_gf.py (limited to 'test') diff --git a/test/test_LV.py b/test/test_LV.py new file mode 100644 index 0000000..72d0a9c --- /dev/null +++ b/test/test_LV.py @@ -0,0 +1,96 @@ +#!/usr/bin/env python + +from __future__ import absolute_import, division + +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib import rc + +import GolemFitPy as gf + +rc('text', usetex=True) +rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) + +dp = gf.DataPaths() +steer = gf.SteeringParams() + +fig = plt.figure(figsize=[12, 8]) +ax = fig.add_subplot(111) +ax.set_xscale('log') +ax.set_yscale('log') + +npp = gf.NewPhysicsParams() +npp.type = gf.NewPhysicsType.None +# npp.n_lv = 1 +# npp.lambda_1 = 1.e100 +# npp.lambda_2 = 1.e100 + +golem = gf.GolemFit(dp, steer, npp) + +binning = golem.GetEnergyBinsMC() +ax.set_xlim(binning[0], binning[-1]) +# ax.set_ylim(binning[0], binning[-1]) + +fit_params = gf.FitParameters(gf.sampleTag.HESE) +golem.SetupAsimov(fit_params) + +exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, + drawstyle='steps-pre', label='NULL', linestyle='--') + +print 'NULL min_llh', golem.MinLLH().likelihood +print 'NULL expectation', exp +print + +npp = gf.NewPhysicsParams() +npp.type = gf.NewPhysicsType.LorentzViolation +npp.n_lv = 1 +npp.lambda_1 = 1.e20 +npp.lambda_2 = 1.e20 + +golem.SetNewPhysicsParams(npp) + +exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, + drawstyle='steps-pre', label='1e-20 LV', linestyle='--') + +print '1e20 LV min_llh', golem.MinLLH().likelihood +print '1e20 LV expectation', exp + +npp = gf.NewPhysicsParams() +npp.type = gf.NewPhysicsType.LorentzViolation +npp.n_lv = 1 +npp.lambda_1 = 1.e10 +npp.lambda_2 = 1.e10 + +golem.SetNewPhysicsParams(npp) + +exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, + drawstyle='steps-pre', label='1e-10 LV', linestyle='--') + +print '1e10 LV min_llh', golem.MinLLH().likelihood +print '1e10 LV expectation', exp + +npp = gf.NewPhysicsParams() +npp.type = gf.NewPhysicsType.LorentzViolation +npp.n_lv = 1 +npp.lambda_1 = 1.e-20 +npp.lambda_2 = 1.e-20 + +golem.SetNewPhysicsParams(npp) + +ax.tick_params(axis='x', labelsize=12) +ax.tick_params(axis='y', labelsize=12) +ax.set_xlabel(r'Deposited energy / GeV') +ax.set_ylabel(r'Events') +for xmaj in ax.xaxis.get_majorticklocs(): + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.7, linewidth=1) +for ymaj in ax.yaxis.get_majorticklocs(): + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.7, linewidth=1) + +legend = ax.legend(prop=dict(size=12)) +fig.savefig('test.png', bbox_inches='tight', dpi=250) + diff --git a/test/test_NSI.py b/test/test_NSI.py new file mode 100644 index 0000000..d144420 --- /dev/null +++ b/test/test_NSI.py @@ -0,0 +1,106 @@ +#!/usr/bin/env python + +from __future__ import absolute_import, division + +import numpy as np +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib import rc + +import GolemFitPy as gf + +rc('text', usetex=True) +rc('font', **{'family':'serif', 'serif':['Computer Modern'], 'size':18}) + +dp = gf.DataPaths() +steer = gf.SteeringParams() +npp = gf.NewPhysicsParams() + +steer.quiet = False +steer.fastmode = True + +golem = gf.GolemFit(dp, steer, npp) + +fit_params = gf.FitParameters(gf.sampleTag.HESE) +golem.SetupAsimov(fit_params) + +fig = plt.figure(figsize=[6, 5]) +ax = fig.add_subplot(111) +ax.set_xscale('log') +ax.set_yscale('log') + +binning = golem.GetEnergyBinsMC() +ax.set_xlim(binning[0], binning[-1]) +# ax.set_ylim(binning[0], binning[-1]) + +print 'NULL min_llh', golem.MinLLH().likelihood + +# exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +# ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, +# drawstyle='steps-pre', label='NULL', linestyle='--') + +# print 'NULL expectation', exp +print + +npp.type = gf.NewPhysicsType.NonStandardInteraction +npp.epsilon_mutau = 0.1 +golem.SetNewPhysicsParams(npp) + +print '0.1 mutau min_llh', golem.MinLLH().likelihood + +# exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +# ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, +# drawstyle='steps-pre', label='0.1 mutau', linestyle='--') + +# print '0.1 mutau expectation', exp +print + +np.epsilon_mutau = 0.2 +golem.SetNewPhysicsParams(npp) + +print '0.2 mutau min_llh', golem.MinLLH().likelihood + +# exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +# ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, +# drawstyle='steps-pre', label='0.2 mutau', linestyle='--') + +# print '0.2 mutau expectation', exp +print + +np.epsilon_mutau = 0.3 +golem.SetNewPhysicsParams(npp) + +print '0.3 mutau min_llh', golem.MinLLH().likelihood + +# exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +# ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, +# drawstyle='steps-pre', label='0.3 mutau', linestyle='--') + +# print '0.3 mutau expectation', exp +print + +np.epsilon_mutau = 0.4 +golem.SetNewPhysicsParams(npp) + +print '0.4 mutau min_llh', golem.MinLLH().likelihood + +# exp = np.sum(golem.GetExpectation(fit_params), axis=(0, 1, 2, 3)) +# ax.step(binning, np.concatenate([[exp[0]], exp]), alpha=1, +# drawstyle='steps-pre', label='0.4 mutau', linestyle='--') + +# print '0.4 mutau expectation', exp +print + +ax.tick_params(axis='x', labelsize=12) +ax.tick_params(axis='y', labelsize=12) +ax.set_xlabel(r'Deposited energy / GeV') +ax.set_ylabel(r'Events') +for xmaj in ax.xaxis.get_majorticklocs(): + ax.axvline(x=xmaj, ls=':', color='gray', alpha=0.7, linewidth=1) +for ymaj in ax.yaxis.get_majorticklocs(): + ax.axhline(y=ymaj, ls=':', color='gray', alpha=0.7, linewidth=1) + +legend = ax.legend(prop=dict(size=12)) +fig.savefig('test_NSI.png', bbox_inches='tight', dpi=250) + diff --git a/test/test_all_gf.py b/test/test_all_gf.py new file mode 100644 index 0000000..04b0980 --- /dev/null +++ b/test/test_all_gf.py @@ -0,0 +1,79 @@ +import numpy as np +import matplotlib as mpl +mpl.use('Agg') +import matplotlib.pyplot as plt + +import GolemFitPy as gf + +FASTMODE = True +PARAMETERS = [ + # 'astroFlavorAngle1', 'astroFlavorAngle2', + 'convNorm', + # 'promptNorm', 'muonNorm', 'astroNorm' +] +DEFAULTS = [ + # 4/9., 0., 1., 0., 1., 6.9 + 1. +] +RANGES = [ + # (0, 1), (-1, 1), (0.01, 10), (0., 30), (0.01, 10), (0.01, 30) + (0.01, 10) +] +BINS = 50 + +def steering_params(): + steering_categ = 'p2_0' + params = gf.SteeringParams(gf.sampleTag.HESE) + if FASTMODE: + params.fastmode = True + params.quiet = False + params.simToLoad= steering_categ.lower() + return params + +def set_up_as(fitter, params): + print 'Injecting the model', params + asimov_params = gf.FitParameters(gf.sampleTag.HESE) + for x in params.iterkeys(): + asimov_params.__setattr__(x, float(params[x])) + fitter.SetupAsimov(asimov_params) + priors = gf.Priors() + priors.convNormWidth = 9e9 + fitter.SetFitPriors(priors) + +def setup_fitter(asimov_paramset): + datapaths = gf.DataPaths() + sparams = steering_params() + 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 x in params.iterkeys(): + fitparams.__setattr__(x, float(params[x])) + llh = -fitter.EvalLLH(fitparams) + return llh + +for ip, param in enumerate(PARAMETERS): + asimov_paramset = {param: DEFAULTS[ip]} + print 'injecting', asimov_paramset + fitter = setup_fitter(asimov_paramset) + binning = np.linspace(RANGES[ip][0], RANGES[ip][1], BINS) + llhs = [] + for b in binning: + test_paramset = {param: b} + print 'testing', test_paramset + llh = get_llh(fitter, test_paramset) + print 'llh', llh + llhs.append(llh) + plt.plot(binning, llhs) + plt.axvline(x=DEFAULTS[ip]) + plt.xlabel(param) + plt.ylabel('LLH') + outfile = 'llh_profile_noprior_' + if FASTMODE: + plt.savefig(outfile + 'fastmode_{0}.png'.format(param)) + else: + plt.savefig(outfile + '{0}.png'.format(param)) + plt.clf() diff --git a/test/test_gf.py b/test/test_gf.py index 6f8f463..4b77924 100644 --- a/test/test_gf.py +++ b/test/test_gf.py @@ -1,9 +1,14 @@ +import numpy.ma as ma + import GolemFitPy as gf +FASTMODE = True + def steering_params(): steering_categ = 'p2_0' params = gf.SteeringParams(gf.sampleTag.HESE) - params.fastmode = True + if FASTMODE: + params.fastmode = True params.quiet = False params.simToLoad= steering_categ.lower() return params @@ -50,7 +55,7 @@ test_paramset = {'astroFlavorAngle1': 4/9., 'astroFlavorAngle2': 0} print 'testing', test_paramset print 'llh', get_llh(fitter, test_paramset) -test_paramset = {'astroFlavorAngle1': 4/9, 'astroFlavorAngle2': 0} +test_paramset = {'astroFlavorAngle1': 4/9., 'astroFlavorAngle2': 0} print 'testing', test_paramset print 'llh', get_llh(fitter, test_paramset) @@ -68,9 +73,17 @@ for an1 in angle1_binning: llhs = [] for an2 in angle2_binning: test_paramset = {'astroFlavorAngle1': an1, 'astroFlavorAngle2': an2} - llhs.append(get_llh(fitter, test_paramset)) + llh = get_llh(fitter, test_paramset) + if abs(llh) > 9e9: + llhs.append(np.nan) + else: + llhs.append(llh) + llhs = ma.masked_invalid(llhs) plt.plot(angle2_binning, llhs, label='astroFlavorAngle1 = {0}'.format(an1)) plt.xlabel('astroFlavorAngle2') plt.ylabel('LLH') plt.legend() -plt.savefig('llh_profile_fastmode.png'.format(an1)) +if FASTMODE: + plt.savefig('llh_profile_fastmode.png'.format(an1)) +else: + plt.savefig('llh_profile.png'.format(an1)) diff --git a/test/test_gf_simple.py b/test/test_gf_simple.py index 51460c1..3a6ebb5 100644 --- a/test/test_gf_simple.py +++ b/test/test_gf_simple.py @@ -1,11 +1,20 @@ +import numpy as np +import matplotlib as mpl +mpl.use('Agg') +import matplotlib.pyplot as plt + import GolemFitPy as gf +FASTMODE = True + dp = gf.DataPaths() npp = gf.NewPhysicsParams() sp = gf.SteeringParams(gf.sampleTag.HESE) sp.quiet = False -# sp.fastmode = True +if FASTMODE: + sp.fastmode = True +sp.frequentist = True golem = gf.GolemFit(dp, sp, npp) @@ -20,4 +29,33 @@ fp_sh.astroFlavorAngle1 = 0.36 fp_sh.astroFlavorAngle2 = -0.57 print 'Eval fp = {0}'.format(golem.EvalLLH(fp)) + +# energy_centers = golem.GetEnergyBinsMC()[:-1]+ np.diff(golem.GetEnergyBinsMC())/2. + +# plt.hist(energy_centers,bins=golem.GetEnergyBinsMC(), +# weights=np.sum(golem.GetExpectation(fp),axis=(0,1,2,3)), +# histtype="step", lw = 2, label='injected') + +# data_energy_dist = np.sum(golem.GetDataDistribution(),axis=(0,1,2,3)) +# energy_centers=golem.GetEnergyBinsData()[:-1]+ np.diff(golem.GetEnergyBinsData())/2. +# plt.errorbar(energy_centers,data_energy_dist,yerr = np.sqrt(data_energy_dist),fmt='o') + print 'Eval fp_sh = {0}'.format(golem.EvalLLH(fp_sh)) + +# plt.hist(energy_centers,bins=golem.GetEnergyBinsMC(), +# weights=np.sum(golem.GetExpectation(fp_sh),axis=(0,1,2,3)), +# histtype="step", lw = 2, label='test') + +# data_energy_dist = np.sum(golem.GetDataDistribution(),axis=(0,1,2,3)) +# energy_centers=golem.GetEnergyBinsData()[:-1]+ np.diff(golem.GetEnergyBinsData())/2. +# plt.errorbar(energy_centers,data_energy_dist,yerr = np.sqrt(data_energy_dist),fmt='o') + +# plt.loglog(nonposy="clip") +# plt.xlabel(r"Deposited energy/GeV") +# plt.ylabel(r"Events") + +# outname = 'Expectation' +# if FASTMODE: +# plt.savefig(outname + 'fastmode.png') +# else: +# plt.savefig(outname + '.png') -- cgit v1.2.3