add test scripts for Golem LV and NSI

6 files changed, 333 insertions, 0 deletions

diff --git a/test/LV.out b/test/LV.out
new file mode 100644
index 0000000..e9f2b7e
--- /dev/null
+++ b/test/LV.out
@@ -0,0 +1,74 @@
+test_LV.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 2u, std::allocator<double> > already registered; second conversion method ignored.
+  import GolemFitPy as gf
+test_LV.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 5u, std::allocator<double> > already registered; second conversion method ignored.
+  import GolemFitPy as gf
+reset_steering: 1
+reset_data: 1
+reset_npp: 1
+GolemFit constructor: checking paths
+Constructing DiffuseWeightMaker
+Loading data
+Loading HESE data.
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/data/HESEData.txt
+Loaded HESE 80 events.
+Loading MC
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 now main simulation size is 420660
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 now main simulation size is 616803
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 now main simulation size is 740041
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 now main simulation size is 740161
+Loaded 740161 events in main simulation set
+Loading Flux weighter
+Loading nusquids objects
+Weighting MC
+Initializing simulation weights
+Applying selfveto correction.
+HESE reshuffle has not been performed.
+Doing HESE reshuffle
+Making data hist
+Making sim hist
+Constructing likelihood problem with default settings
+Setting up Asimov set with parameters
+Conventional normalization 1
+Prompt normalization 1
+Atmospheric muon normalization 1
+Astro component overall normalization 8
+Astro E component 0.333333
+Astro Mu component 0.333333
+Astro Tau component 0.333333
+Astroparticle balance 1
+Astro gamma 2.5
+Astro cutoff 10
+CR gamma 0
+Conv pi/k ratio 1
+Conv particle balance 1
+Conv zenith correction 0
+Dark normalization 0
+DOM efficiency 0.99
+Astro component overall normalization second 0
+Astro gamma second 2
+Holeice forward 0
+Anisotropy scale1
+
+Remaking data hist
+Reconstrucing likelihood problem
+NULL min_llh 835.848062048
+NULL expectation [11.96478761 11.90169848 10.78840154  9.20158936  7.23605538  5.59952977
+  4.39850777  3.35653018  2.58491686  1.92834114  1.46813946  1.11581282
+  0.8492433   0.65154743  0.46727371  0.36211812  0.32718588  0.69746002
+  0.29475116  0.05150809]
+
+Applying new flavor physics to MC weights.
+Entering ApplyNewPhysics
+Entering ConstructNuSQuIDSObjectsForLV.
+Rotated Unitary Transform
+In for loop
+Entering inner for loop
+Leaving inner for loop
diff --git a/test/NSI.out b/test/NSI.out
new file mode 100644
index 0000000..657a9f3
--- /dev/null
+++ b/test/NSI.out
@@ -0,0 +1,75 @@
+test_NSI.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 2u, std::allocator<double> > already registered; second conversion method ignored.
+  import GolemFitPy as gf
+test_NSI.py:11: RuntimeWarning: to-Python converter for nusquids::marray<double, 5u, std::allocator<double> > already registered; second conversion method ignored.
+  import GolemFitPy as gf
+reset_steering: 1
+reset_data: 1
+reset_npp: 1
+GolemFit constructor: checking paths
+Constructing DiffuseWeightMaker
+Loading data
+Loading HESE data.
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/data/HESEData.txt
+Loaded HESE 80 events.
+Loading MC
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nue_p2_0.h5 now main simulation size is 420660
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/numu_p2_0.h5 now main simulation size is 616803
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/nutau_p2_0.h5 now main simulation size is 740041
+Preparing to read /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 into main simulation
+Reading a file from path /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5
+Finish reading /data/user/smandalia/GolemTools/sources/GolemFit/monte_carlo//HESE/muongun.h5 now main simulation size is 740161
+Loaded 740161 events in main simulation set
+Loading Flux weighter
+Loading nusquids objects
+Weighting MC
+Initializing simulation weights
+Applying selfveto correction.
+HESE reshuffle has not been performed.
+Doing HESE reshuffle
+Making data hist
+Making sim hist
+Constructing likelihood problem with default settings
+Setting up Asimov set with parameters
+Conventional normalization 1
+Prompt normalization 1
+Atmospheric muon normalization 1
+Astro component overall normalization 8
+Astro E component 0.333333
+Astro Mu component 0.333333
+Astro Tau component 0.333333
+Astroparticle balance 1
+Astro gamma 2.5
+Astro cutoff 10
+CR gamma 0
+Conv pi/k ratio 1
+Conv particle balance 1
+Conv zenith correction 0
+Dark normalization 0
+DOM efficiency 0.99
+Astro component overall normalization second 0
+Astro gamma second 2
+Holeice forward 0
+Anisotropy scale1
+
+Remaking data hist
+Reconstrucing likelihood problem
+NULL min_llh 835.848062048
+NULL expectation [11.96478761 11.90169848 10.78840154  9.20158936  7.23605538  5.59952977
+  4.39850777  3.35653018  2.58491686  1.92834114  1.46813946  1.11581282
+  0.8492433   0.65154743  0.46727371  0.36211812  0.32718588  0.69746002
+  0.29475116  0.05150809]
+
+Applying new flavor physics to MC weights.
+Applying New Physics in normal mode to NonStandardInteraction
+HESE reshuffle has been performed.
+0.1 mutau min_llh 836.352921667
+0.1 mutau expectation [16.63858821 16.04604012 14.28034835 11.80228198  9.06453127  6.8815312
+  5.32704912  4.02384379  3.07832186  2.28036523  1.73329938  1.30341591
+  0.99239868  0.75237435  0.53312845  0.4093463   0.35723267  0.71433954
+  0.30173276  0.05276177]
diff --git a/test/test.png b/test/test.png
new file mode 100644
index 0000000..8aac9e3
--- /dev/null
+++ b/test/test.png
diff --git a/test/test_LV.py b/test/test_LV.py
new file mode 100644
index 0000000..96a1863
--- /dev/null
+++ b/test/test_LV.py
@@ -0,0 +1,95 @@
+#!/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.png b/test/test_NSI.png
new file mode 100644
index 0000000..2eb4345
--- /dev/null
+++ b/test/test_NSI.png
diff --git a/test/test_NSI.py b/test/test_NSI.py
new file mode 100644
index 0000000..617c353
--- /dev/null
+++ b/test/test_NSI.py
@@ -0,0 +1,89 @@
+#!/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.epsilon_mutau = 0
+# npp.epsilon_prime = 0
+
+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.NonStandardInteraction
+npp.epsilon_mutau = 0.1
+# npp.epsilon_prime = 0
+
+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='0.1 mutau', linestyle='--')
+
+print '0.1 mutau min_llh', golem.MinLLH().likelihood
+print '0.1 mutau expectation', exp
+
+# npp = gf.NewPhysicsParams()
+# npp.epsilon_mutau = 0
+# # npp.epsilon_prime = 0
+
+# 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.epsilon_mutau = 0
+# # npp.epsilon_prime = 0
+
+# 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_NSI.png', bbox_inches='tight', dpi=250)