Sun May 6 21:22:46 CDT 2018

5 files changed, 337 insertions, 5 deletions

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')