add planck scale in plot

5 files changed, 52 insertions, 16 deletions

diff --git a/plot_sens.py b/plot_sens.py
index 893e7a0..610ee95 100755
--- a/plot_sens.py
+++ b/plot_sens.py
@@ -231,7 +231,8 @@ def main():
                 infile += '/gaussian/'
             if args.likelihood is Likelihood.GAUSSIAN:
                 infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_'))
-            infile += '/DIM{0}/fix_ifr/{1}/{2}/{3}/fr_stat'.format(
+            # infile += '/DIM{0}/fix_ifr/{1}/{2}/{3}/fr_stat'.format(
+            infile += '/DIM{0}/fix_ifr/seed2/{1}/{2}/{3}/fr_stat'.format(
             # infile += '/DIM{0}/fix_ifr/100TeV/{1}/{2}/{3}/fr_stat'.format(
                 dim, *map(misc_utils.parse_enum, [args.stat_method, args.run_method, args.data])
             ) + misc_utils.gen_identifier(argsc)
@@ -278,7 +279,8 @@ def main():
                 base_infile += '/gaussian/'
             if args.likelihood is Likelihood.GAUSSIAN:
                 base_infile += '{0}/'.format(str(args.sigma_ratio).replace('.', '_'))
-            base_infile += '/DIM{0}/fix_ifr'.format(dim)
+            # base_infile += '/DIM{0}/fix_ifr'.format(dim)
+            base_infile += '/DIM{0}/fix_ifr/seed2'.format(dim)
             # base_infile += '/DIM{0}/fix_ifr/100TeV'.format(dim)
 
             for isrc, src in enumerate(args.source_ratios):
@@ -345,14 +347,14 @@ def main():
             plot_utils.plot_sens_full(
                 data      = data,
                 outfile   = baseoutfile + '/FULL',
-                outformat = ['png'],
+                outformat = ['png', 'pdf'],
                 args      = args,
             )
         elif args.run_method in fixed_angle_categ:
             plot_utils.plot_sens_fixed_angle_pretty(
                 data      = data,
                 outfile   = baseoutfile + '/fixed_angle_pretty',
-                outformat = ['png'],
+                outformat = ['png', 'pdf'],
                 args      = args,
             )
             # plot_utils.plot_sens_fixed_angle(
@@ -365,7 +367,7 @@ def main():
             plot_utils.plot_sens_corr_angle(
                 data      = data,
                 outfile   = baseoutfile + '/CORR_ANGLE',
-                outformat = ['png'],
+                outformat = ['png', 'pdf'],
                 args      = args,
             )
 
diff --git a/sens.py b/sens.py
index 5b79d0d..36e4333 100755
--- a/sens.py
+++ b/sens.py
@@ -56,7 +56,7 @@ def define_nuisance():
         Param(name='promptNorm',      value=0.,  seed=[0. , 6. ], ranges=[0. , 50.], std=0.05, tag=tag),
         Param(name='muonNorm',        value=1.,  seed=[0.1, 2. ], ranges=[0. , 50.], std=0.1,  tag=tag),
         Param(name='astroNorm',       value=6.9, seed=[0.1, 10.], ranges=[0. , 50.], std=0.1,  tag=tag),
-        Param(name='astroDeltaGamma', value=2.5, seed=[1. , 3. ], ranges=[-5., 5. ], std=0.1,  tag=tag)
+        Param(name='astroDeltaGamma', value=2.5, seed=[2.4, 3. ], ranges=[-5., 5. ], std=0.1,  tag=tag)
     ])
     return ParamSet(nuisance)
 
@@ -288,6 +288,9 @@ def main():
                         print '## Evidence = {0}'.format(stat)
                     elif args.stat_method is StatCateg.FREQUENTIST:
                         def fn(x):
+                            # Force prior ranges to be inside "seed"
+                            for el in x:
+                                if el < 0 or el > 1: return np.inf
                             pranges = sens_paramset.seeds
                             for i, name in enumerate(sens_paramset.names):
                                 llh_paramset[name].value = \
@@ -316,8 +319,8 @@ def main():
                                 if s > stat:
                                     stat = s
                         except AssertionError:
-                            print 'Failed run, continuing'
-                            # raise
+                            # print 'Failed run, continuing'
+                            raise
                             continue
                         print '=== final llh', stat
                     if args.run_method is SensitivityCateg.FULL:
diff --git a/submitter/sens_dag.py b/submitter/sens_dag.py
index 6734278..1dc0ed8 100644
--- a/submitter/sens_dag.py
+++ b/submitter/sens_dag.py
@@ -30,7 +30,7 @@ GLOBAL_PARAMS.update(dict(
     run_method    = 'fixed_angle', # full, fixed_angle, corr_angle
     stat_method   = 'bayesian',
     sens_bins     = 20,
-    seed          = 'None'
+    seed          = None
 ))
 
 # MultiNest
@@ -41,9 +41,9 @@ GLOBAL_PARAMS.update(dict(
 ))
 
 # FR
-# dimension         = [3]
-dimension         = [3, 6]
-# dimension         = [4, 5, 7, 8]
+dimension         = [3]
+# dimension         = [3, 6]
+# dimension         = [3, 4, 5, 6, 7, 8]
 GLOBAL_PARAMS.update(dict(
     threads           = 1,
     binning           = '6e4 1e7 20',
@@ -77,7 +77,7 @@ outfile = 'dagman_FR_SENS_{0}_{1}_{2}_{3}'.format(
     GLOBAL_PARAMS['stat_method'], GLOBAL_PARAMS['run_method'],
     GLOBAL_PARAMS['likelihood'], GLOBAL_PARAMS['data']
 )
-# outfile += '_100TeV'
+outfile += '_seed2'
 outfile += '.submit'
 golemfitsourcepath = os.environ['GOLEMSOURCEPATH'] + '/GolemFit'
 condor_script = golemfitsourcepath + '/scripts/flavour_ratio/submitter/sens_submit.sub'
@@ -101,7 +101,7 @@ with open(outfile, 'w') as f:
             output = outchain_head + '/fix_ifr/'
             if GLOBAL_PARAMS['likelihood'].lower() == 'gaussian':
                 output += '{0}/'.format(str(GLOBAL_PARAMS['sigma_ratio']).replace('.', '_'))
-            # output += '100TeV/'
+            output += 'seed2/'
             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/gf.py b/utils/gf.py
index 02231f0..a288622 100644
--- a/utils/gf.py
+++ b/utils/gf.py
@@ -95,7 +95,7 @@ def setup_fitter(args, asimov_paramset):
     if args.data is DataType.ASIMOV:
         setup_asimov(fitter, asimov_paramset)
     elif args.data is DataType.REALISATION:
-        seed = args.seed if args.seed is not None else 0
+        seed = args.seed if args.seed is not None else 1
         setup_realisation(fitter, asimov_paramset, seed)
     elif args.data is DataType.REAL:
         print 'Using MagicTau DATA'
diff --git a/utils/plot.py b/utils/plot.py
index 7a82f2b..ea3c852 100644
--- a/utils/plot.py
+++ b/utils/plot.py
@@ -408,6 +408,13 @@ def plot_sens_fixed_angle_pretty(data, outfile, outformat, args):
         7: (3.6E-41, 1.77E-32),
         8: (1.4E-45, 1.00E-34)
     }
+    PS = 8.203e-20 # GeV^{-1}
+    planck_scale = {
+        5: PS,
+        6: PS**2,
+        7: PS**3,
+        8: PS**4
+    }
 
     best_limits = {
         3: ('CMB polarization', 1E-43, 'brown'),
@@ -497,6 +504,23 @@ def plot_sens_fixed_angle_pretty(data, outfile, outformat, args):
                 bc_limit = best_limits[dim]
                 # ax.axvline(x=np.log10(bc_limit[1]), color=bc_limit[2], alpha=0.7, linewidth=1.5)
 
+                if dim in planck_scale:
+                    ps = np.log10(planck_scale[dim])
+                    if ps < xlims[0]:
+                        ax.annotate(
+                            s='', xy=(xlims[0], 1), xytext=(xlims[0]+1, 1),
+                            arrowprops={'arrowstyle': '->, head_length=0.2',
+                                        'lw': 1, 'color':'purple'}
+                        )
+                    elif ps > xlims[1]:
+                        ax.annotate(
+                            s='', xy=(xlims[1]-1, 1), xytext=(xlims[1], 1),
+                            arrowprops={'arrowstyle': '<-, head_length=0.2',
+                                        'lw': 1, 'color':'purple'}
+                        )
+                    else:
+                        ax.axvline(x=ps, color='purple', alpha=1., linewidth=1.5)
+
                 scales, statistic = ma.compress_rows(data[idim][isrc][ian]).T
                 tck, u = interpolate.splprep([scales, statistic], s=0)
                 scales, statistic = interpolate.splev(np.linspace(0, 1, 1000), tck)
@@ -557,6 +581,10 @@ def plot_sens_fixed_angle_pretty(data, outfile, outformat, args):
     legend_elements.append(
         Patch(facecolor=col, alpha=alpha+0.1, edgecolor='k', label='maximum reach')
     )
+    purple = [0.5019607843137255, 0.0, 0.5019607843137255]
+    legend_elements.append(
+        Patch(facecolor=purple+[0.7], edgecolor=purple+[1], label='Planck Scale')
+    )
     legend_elements.append(
         Patch(facecolor='none', hatch='\\\\', edgecolor='k', label='IceCube arXiv:1709.03434')
     )
@@ -570,12 +598,15 @@ def plot_sens_fixed_angle_pretty(data, outfile, outformat, args):
     plt.setp(legend.get_title(), fontsize='11')
     legend.get_frame().set_linestyle('-')
 
-    if show_data: ybound = 0.65
+    if show_data: ybound = 0.595
     else: ybound = 0.73
     if args.data is DataType.REAL and show_data:
         # fig.text(0.295, 0.684, 'IceCube Preliminary', color='red', fontsize=13,
         fig.text(0.278, ybound, r'\bf IceCube Preliminary', color='red', fontsize=13,
                  ha='center', va='center', zorder=11)
+    elif args.data is DataType.REALISATION:
+        fig.text(0.278, ybound-0.05, r'\bf IceCube Simulation', color='red', fontsize=13,
+                 ha='center', va='center', zorder=11)
     else:
         fig.text(0.278, ybound, r'\bf IceCube Simulation', color='red', fontsize=13,
                  ha='center', va='center', zorder=11)