fix bug in energy spectrum

10 files changed, 62 insertions, 73 deletions

diff --git a/plot_sens.py b/plot_sens.py
index 2e3b1c6..9aa5089 100755
--- a/plot_sens.py
+++ b/plot_sens.py
@@ -145,10 +145,6 @@ def parse_args(args=None):
         help='Make sensitivity plots'
     )
     parser.add_argument(
-        '--data', type=misc_utils.parse_bool, default='False',
-        help='Use HESE7 data'
-    )
-    parser.add_argument(
         '--plot-statistic', type=misc_utils.parse_bool, default='False',
         help='Plot MultiNest evidence or LLH value'
     )
@@ -231,13 +227,11 @@ def main():
                 infile += '/golemfit/'
             elif args.likelihood is Likelihood.GAUSSIAN:
                 infile += '/gaussian/'
-            # infile += '/DIM{0}/fix_ifr/HESESim'.format(dim)
-            # infile += '/DIM{0}/fix_ifr/'.format(dim)
-            infile += '/DIM{0}/fix_ifr/data'.format(dim)
             if args.likelihood is Likelihood.GAUSSIAN:
                 infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_'))
-            infile += 'fr_stat/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \
-                + misc_utils.gen_identifier(argsc)
+            infile += '/{0}/{1}/{2}/fr_stat'.format(
+                *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data])
+            ) + misc_utils.gen_identifier(argsc)
             print '== {0:<25} = {1}'.format('infile', infile)
 
             if args.split_jobs:
@@ -279,17 +273,14 @@ def main():
             argsc.dimension = dim
             _, scale_region = fr_utils.estimate_scale(argsc)
             argsc.scale_region = scale_region
-            # infile = base_infile + '/DIM{0}/fix_ifrHESESim/'.format(dim)
-            # infile = base_infile + '/DIM{0}/fix_ifr/'.format(dim)
-            infile = base_infile + '/DIM{0}/fix_ifr/data'.format(dim)
             if args.likelihood is Likelihood.GAUSSIAN:
                 infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_'))
-            infile += 'fr_stat/'
 
             for isrc, src in enumerate(args.source_ratios):
                 argsc.source_ratio = src
-                finfile = infile +'/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \
-                    + misc_utils.gen_identifier(argsc)
+                finfile = infile +'/{0}/{1}/{2}/fr_stat'.format(
+                    *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data])
+                ) + misc_utils.gen_identifier(argsc)
                 basename = os.path.dirname(finfile) + '/statrun/'
                 baseoutfile = basename[:5]+basename[5:].replace('data', 'plots')
                 if args.data:
@@ -343,11 +334,9 @@ def main():
         print 'Plotting sensitivities'
 
         basename = args.infile[:5]+args.infile[5:].replace('data', 'plots')
-        if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]:
-            baseoutfile = basename + 'golemfit/{1}'.format(args.likelihood, args.stat_method)
-            if args.data: baseoutfile += '/data'
-        elif args.likelihood is Likelihood.GAUSSIAN:
-            baseoutfile = basename + 'gaussian/{1}'.format(args.likelihood, args.stat_method)
+        baseoutfile = basename + '/{1}/{2}/{3}'.format(
+            *map(misc_utils.parse_enum, [args.likelihood, args.stat_method, args.data])
+        )
 
         if args.run_method is SensitivityCateg.FULL:
             plot_utils.plot_sens_full(
diff --git a/sens.py b/sens.py
index 6145f4a..581a3ca 100755
--- a/sens.py
+++ b/sens.py
@@ -207,8 +207,9 @@ def main():
         eval_dim = args.sens_bins
     else: eval_dim = 1
 
-    out = args.outfile+'/{0}/{1}/fr_stat'.format(args.stat_method, args.run_method) \
-        + misc_utils.gen_identifier(args)
+    out = args.outfile+'/{0}/{1}/{2}/fr_stat'.format(
+        *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data])
+    ) + misc_utils.gen_identifier(args)
     if args.sens_run:
         if args.likelihood in [Likelihood.GOLEMFIT, Likelihood.GF_FREQ]:
             fitter = gf_utils.setup_fitter(args, asimov_paramset)
diff --git a/submitter/mcmc_dag.py b/submitter/mcmc_dag.py
index 7a0cde9..e40c043 100644
--- a/submitter/mcmc_dag.py
+++ b/submitter/mcmc_dag.py
@@ -38,7 +38,7 @@ dimension         = [3]
 # dimension         = [4, 5, 7, 8]
 GLOBAL_PARAMS.update(dict(
     threads           = 1,
-    binning           = '6e4 1e7 5',
+    binning           = '6e4 1e7 20',
     no_bsm            = 'False',
     scale_region      = "1E10",
     energy_dependance = 'spectral',
@@ -57,7 +57,7 @@ GLOBAL_PARAMS.update(dict(
 # GolemFit
 GLOBAL_PARAMS.update(dict(
     ast  = 'p2_0',
-    data = 'real'
+    data = 'asimov'
 ))
 
 # Plot
@@ -82,6 +82,7 @@ with open(outfile, 'w') as f:
             outchains = outchain_head + '/fix_ifr/'
             if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
                 outchains += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
+            outchains += '{0}/'.format(GLOBAL_PARAMS['data'].lower())
             outchains += 'mcmc_chain'
             f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
             f.write('VARS\tjob{0}\tdimension="{1}"\n'.format(job_number, dim))
diff --git a/submitter/mcmc_submit.sub b/submitter/mcmc_submit.sub
index 998932e..fa36250 100644
--- a/submitter/mcmc_submit.sub
+++ b/submitter/mcmc_submit.sub
@@ -32,9 +32,10 @@ Notification    = never
 
 # +TransferOutput=""
 
++NATIVE_OS = True
 # Requirements = IS_GLIDEIN && HAS_CVMFS_icecube_opensciencegrid_org && (OpSysAndVer =?= "CentOS6" || OpSysAndVer =?= "RedHat6" || OpSysAndVer =?= "SL6")
 # Requirements = IS_GLIDEIN
-Requirements = (OpSysMajorVer =?= 6)
+# Requirements = (OpSysMajorVer =?= 6)
 
 # GO! 
 queue
diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py
index 666a0dd..5e00dfd 100644
--- a/submitter/sens_dag.py
+++ b/submitter/sens_dag.py
@@ -35,27 +35,25 @@ GLOBAL_PARAMS.update(dict(
 
 # MultiNest
 GLOBAL_PARAMS.update(dict(
-    mn_live_points = 300,
+    mn_live_points = 1000,
     mn_tolerance   = 0.1,
     mn_output      = './mnrun'
 ))
 
 # FR
-# dimension         = [3]
-dimension         = [3, 6]
-# dimension         = [4, 5, 7, 8]
-# dimension         = [3, 4, 5, 6, 7, 8]
+dimension         = [3]
+# dimension         = [3, 6]
 GLOBAL_PARAMS.update(dict(
     threads           = 1,
-    # binning           = '6e4 1e7 5',
-    binning           = '1e5 1e7 5',
+    binning           = '6e4 1e7 20',
+    # binning           = '1e5 1e7 5',
     no_bsm            = 'False',
     scale_region      = "1E10",
     energy_dependance = 'spectral',
     spectral_index    = -2,
     fix_mixing        = 'False',
     fix_mixing_almost = 'False',
-    fold_index        = 'False'
+    fold_index        = 'True'
 ))
 
 # Likelihood
@@ -67,7 +65,7 @@ GLOBAL_PARAMS.update(dict(
 # GolemFit
 GLOBAL_PARAMS.update(dict(
     ast  = 'p2_0',
-    data = 'real'
+    data = 'asimov'
 ))
 
 # Plot
@@ -75,8 +73,9 @@ GLOBAL_PARAMS.update(dict(
     plot_statistic = 'True'
 ))
 
-outfile = 'dagman_FR_SENS_{0}_{1}_{2}_data_100TeV.submit'.format(
-    GLOBAL_PARAMS['stat_method'], GLOBAL_PARAMS['run_method'], GLOBAL_PARAMS['likelihood']
+outfile = 'dagman_FR_SENS_{0}_{1}_{2}_{3}.submit'.format(
+    GLOBAL_PARAMS['stat_method'], GLOBAL_PARAMS['run_method'],
+    GLOBAL_PARAMS['likelihood'], GLOBAL_PARAMS['data']
 )
 golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
 condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub'
@@ -97,14 +96,9 @@ with open(outfile, 'w') as f:
         )
         for frs in fix_sfr_mfr:
             print 'frs', frs
-            # output = outchain_head + '/fix_ifr/'
-            # output = outchain_head + '/fix_ifr/HESESim'
-            # output = outchain_head + '/fix_ifr/sim'
-            # output = outchain_head + '/fix_ifr/data'
-            output = outchain_head + '/fix_ifr/data/100TeV/'
+            output = outchain_head + '/fix_ifr/'
             if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
                 output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
-            output += 'fr_stat'
             for r in xrange(sens_runs):
                 print 'run', r
                 f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
@@ -130,7 +124,6 @@ with open(outfile, 'w') as f:
             output = outchain_head + '/full/'
             if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
                 output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
-            output += 'fr_stat'
             for r in xrange(sens_runs):
                 print 'run', r
                 f.write('JOB\tjob{0}\t{1}\n'.format(job_number, condor_script))
diff --git a/utils/enums.py b/utils/enums.py
index b80b712..c40d681 100644
--- a/utils/enums.py
+++ b/utils/enums.py
@@ -11,9 +11,8 @@ from enum import Enum
 
 
 class DataType(Enum):
-    REAL    = 1
-    FAKE    = 2
-    ASMIMOV = 3
+    REAL   = 1
+    ASIMOV = 2
 
 
 class EnergyDependance(Enum):
diff --git a/utils/fr.py b/utils/fr.py
index 64e2706..639927f 100644
--- a/utils/fr.py
+++ b/utils/fr.py
@@ -140,16 +140,16 @@ def cardano_eqn(ham):
         )
 
     a = -np.trace(ham)
-    b = (0.5) * ((np.trace(ham))**2 - np.trace(np.dot(ham, ham)))
+    b = DTYPE(1)/2 * ((np.trace(ham))**DTYPE(2) - np.trace(np.dot(ham, ham)))
     c = -np.linalg.det(ham)
 
-    Q = (1/9.) * (a**2 - 3*b)
-    R = (1/54.) * (2*a**3 - 9*a*b + 27*c)
-    theta = np.arccos(R / np.sqrt(Q**3))
+    Q = (DTYPE(1)/9) * (a**DTYPE(2) - DTYPE(3)*b)
+    R = (DTYPE(1)/54) * (DTYPE(2)*a**DTYPE(3) - DTYPE(9)*a*b + DTYPE(27)*c)
+    theta = np.arccos(R / np.sqrt(Q**DTYPE(3)))
 
-    E1 = -2 * np.sqrt(Q) * np.cos(theta/3.) - (1/3.)*a
-    E2 = -2 * np.sqrt(Q) * np.cos((theta - 2.*PI)/3.) - (1/3.)*a
-    E3 = -2 * np.sqrt(Q) * np.cos((theta + 2.*PI)/3.) - (1/3.)*a
+    E1 = -DTYPE(2) * np.sqrt(Q) * np.cos(theta/DTYPE(3)) - (DTYPE(1)/3)*a
+    E2 = -DTYPE(2) * np.sqrt(Q) * np.cos((theta - DTYPE(2)*PI)/DTYPE(3)) - (DTYPE(1)/3)*a
+    E3 = -DTYPE(2) * np.sqrt(Q) * np.cos((theta + DTYPE(2)*PI)/DTYPE(3)) - (DTYPE(1)/3)*a
 
     A1 = ham[1][2] * (ham[0][0] - E1) - ham[1][0]*ham[0][2]
     A2 = ham[1][2] * (ham[0][0] - E2) - ham[1][0]*ham[0][2]
@@ -383,7 +383,7 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
     mass_matrix = np.array(
         [[0, 0, 0], [0, mass_eigenvalues[0], 0], [0, 0, mass_eigenvalues[1]]]
     )
-    sm_ham = (1./(2*energy))*np.dot(sm_u, np.dot(mass_matrix, sm_u.conj().T))
+    sm_ham = (1./(2*energy))*np.matmul(sm_u, np.matmul(mass_matrix, sm_u.conj().T))
     if no_bsm:
         eg_vector = cardano_eqn(sm_ham)
     else:
@@ -391,23 +391,16 @@ def params_to_BSMu(theta, dim, energy, mass_eigenvalues=MASS_EIGENVALUES,
         scale_matrix = np.array(
             [[0, 0, 0], [0, sc1, 0], [0, 0, sc2]]
         )
-        bsm_term = (energy**(dim-3)) * np.dot(new_physics_u, np.dot(scale_matrix, new_physics_u.conj().T))
-
+        bsm_term = (energy**(dim-3)) * np.matmul(new_physics_u, np.matmul(scale_matrix, new_physics_u.conj().T))
         bsm_ham = sm_ham + bsm_term
         eg_vector = cardano_eqn(bsm_ham)
 
     if check_uni:
-        tu = test_unitarity(eg_vector)
-        if not np.abs(np.trace(tu) - 3.) < epsilon or \
-           not np.abs(np.sum(tu) - 3.) < epsilon:
-            raise AssertionError(
-                'Matrix is not unitary!\neg_vector\n{0}\ntest '
-                'u\n{1}'.format(eg_vector, tu)
-            )
+        test_unitarity(eg_vector, rse=True, epsilon=epsilon)
     return eg_vector
 
 
-def test_unitarity(x, prnt=False):
+def test_unitarity(x, prnt=False, rse=False, epsilon=None):
     """Test the unitarity of a matrix.
 
     Parameters
@@ -418,6 +411,9 @@ def test_unitarity(x, prnt=False):
     prnt : bool
         Print the result
 
+    rse : bool
+        Raise Assertion if matrix is not unitary
+
     Returns
     ----------
     numpy ndarray
@@ -435,6 +431,13 @@ def test_unitarity(x, prnt=False):
     f = np.abs(np.dot(x, x.conj().T), dtype=DTYPE)
     if prnt:
         print 'Unitarity test:\n{0}'.format(f)
+    if rse:
+        if not np.abs(np.trace(f) - 3.) < epsilon or \
+           not np.abs(np.sum(f) - 3.) < epsilon:
+            raise AssertionError(
+                'Matrix is not unitary!\nx\n{0}\ntest '
+                'u\n{1}'.format(x, f)
+            )
     return f
 
 
diff --git a/utils/gf.py b/utils/gf.py
index 7ded152..cc093e1 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -55,7 +55,6 @@ def fit_flags(llh_paramset):
 
 def steering_params(args):
     steering_categ = args.ast
-    # params = gf.SteeringParams(gf.sampleTag.HESE)
     params = gf.SteeringParams(gf.sampleTag.MagicTau)
     params.quiet = False
     params.fastmode = True
@@ -66,14 +65,13 @@ def steering_params(args):
     # For Tianlu
     # params.years = [999]
 
-    params.minFitEnergy = 1.0e5 # GeV
+    # params.minFitEnergy = 1.0e5 # GeV
 
     return params
 
 
 def set_up_as(fitter, params):
     print 'Injecting the model', params
-    # asimov_params = gf.FitParameters(gf.sampleTag.HESE)
     asimov_params = gf.FitParameters(gf.sampleTag.MagicTau)
     for parm in params:
         asimov_params.__setattr__(parm.name, float(parm.value))
@@ -85,13 +83,14 @@ def setup_fitter(args, asimov_paramset):
     sparams = steering_params(args)
     npp = gf.NewPhysicsParams()
     fitter = gf.GolemFit(datapaths, sparams, npp)
-    # comment to use data
-    # set_up_as(fitter, asimov_paramset)
+    if args.data is DataType.ASIMOV:
+        set_up_as(fitter, asimov_paramset)
+    elif args.data is DataType.REAL:
+        print 'Using MagicTau DATA'
     return fitter
 
 
 def get_llh(fitter, params):
-    # fitparams = gf.FitParameters(gf.sampleTag.HESE)
     fitparams = gf.FitParameters(gf.sampleTag.MagicTau)
     for parm in params:
         fitparams.__setattr__(parm.name, float(parm.value))
@@ -101,7 +100,6 @@ def get_llh(fitter, params):
 
 def get_llh_freq(fitter, params):
     print 'setting to {0}'.format(params)
-    # fitparams = gf.FitParameters(gf.sampleTag.HESE)
     fitparams = gf.FitParameters(gf.sampleTag.MagicTau)
     for parm in params:
         fitparams.__setattr__(parm.name, float(parm.value))
@@ -118,7 +116,7 @@ def data_distributions(fitter):
 
 def gf_argparse(parser):
     parser.add_argument(
-        '--data', default='real', type=partial(enum_parse, c=DataType),
+        '--data', default='asimov', type=partial(enum_parse, c=DataType),
         choices=DataType, help='select datatype'
     )
     parser.add_argument(
diff --git a/utils/likelihood.py b/utils/likelihood.py
index ebc3539..ed598c8 100644
--- a/utils/likelihood.py
+++ b/utils/likelihood.py
@@ -146,7 +146,7 @@ def triangle_llh(theta, args, asimov_paramset, llh_paramset, fitter):
             u = fr_utils.params_to_BSMu(
                 theta             = bsm_angles,
                 dim               = args.dimension,
-                energy            = args.energy,
+                energy            = bin_centers[i_sf],
                 mass_eigenvalues  = mass_eigenvalues,
                 sm_u              = sm_u,
                 no_bsm            = args.no_bsm,
diff --git a/utils/misc.py b/utils/misc.py
index 5483675..7e40d13 100644
--- a/utils/misc.py
+++ b/utils/misc.py
@@ -106,6 +106,10 @@ def parse_bool(s):
         raise ValueError
 
 
+def parse_enum(e):
+    return '{0}'.format(e).split('.')[1].lower()
+
+
 def print_args(args):
     """Print the input arguments.