Merge branch 'master' of github.com:IceCubeOpenSource/flarestack

flarestack/analyses/ccsn/__init__.py

@@ -34,13 +34,11 @@ def plot_dec_sens(dec, sens, fname):
     plt.savefig(fname)
 
 
-sn_colors = {'IIp': 'gold', 'IIn': 'red', 'Ibc': 'blue'}
+sn_colors = {'IIp': 'orange', 'IIn': 'red', 'Ibc': 'blue'}
 
 
 def get_sn_color(type):
-    color = sn_colors.get(type, None)
-    if not color and type == 'IIP':
-        sn_colors.get('IIp', None)
-    if not color:
-        raise Exception('SN type {0} not in dictionary!'.format(type))
-    return color
+    if type == 'IIP':
+        return sn_colors['IIp']
+    else:
+        return sn_colors[type]

flarestack/analyses/ccsn/necker_2019/README.md

@@ -1,7 +1,11 @@
-Name: Core Collapse Supernova Analysis
-Analyser: Jannis Necker (updated) Alex Stasik (original)
-Wiki Page: https://wiki.icecube.wisc.edu/index.php/Point_Source_Supernova_Stacking_Analysis
-Datasets: ps_tracks (v002_p01)
+**Name**: Core Collapse Supernova Analysis
+
+**Analyser**: Jannis Necker (updated) Alex Stasik (original)
+
+**Wiki Page**: https://wiki.icecube.wisc.edu/index.php/Point_Source_Supernova_Stacking_Analysis/Update
+
+**Datasets**: ps_tracks (v002_p01)
+
 
 This analysis was originally performed by Alexander Stasik, using three CCSNe
 catalogues with multiple time PDFs. Both fixed weight analyses, as
@@ -11,4 +15,27 @@ with a nearby sample containing 70% of the weight being analysed
 separately from the remaining 30%. Individual analyses on specific CCSNe were
 also performed. 
 
-More info provided at some point.
+The original supernova catalogue contained some errors. I (Jannis Necker) went through the nearby
+sample and and checked each supernova and collected references. A link to the catalogue and 
+references can be found on the Wiki. To also try to reproduce the sensitivities assuming the decay model, 
+the corresponding PDF was implemented.
+
+
+**REPRODUCICIBILITY GUIDE**
+
+For the sensitivity calculation run `box_sensitivity.py` and `decay_sensitivity.py` respectively.
+
+If the .npy catalogue files are not present under 
+`/path/to/flaresatck/flaresatck/analyses/ccsn/necker2019/catalogues`
+download the updated and flarestack readable .csv files from the link in the Wiki and store 
+them under
+`/path/to/flaresatck/flaresatck/analyses/ccsn/necker2019/catalogues/raw_necker/<sn_type>.csv`.
+Then run `build_catalogues_from_raw.py`
+
+The plots illustrating the difference to Alex Stasik's original results, 
+the reproduced ones with _Flaresatck_ and the ones from this analysis 
+are made with `make_ratio_plots.py`. Note that before that, the analysis scripts under 
+`/path/to/flaresatck/flaresatck/analyses/ccsn/stasik2017/` have to be run.
+
+All plots can be seen under 
+/path/to/scratch/flarestack_data/output/plots/analyses/ccsn/stasik2107/fit_weights/

flarestack/analyses/ccsn/necker_2019/__init__.py

@@ -1,46 +0,0 @@
-# import os
-# from flarestack.analyses.ccsn.stasik_2017.shared_ccsn import sn_catalogue_name
-#
-# ccsn_dir = os.path.abspath(os.path.dirname(__file__))
-# ccsn_cat_dir = ccsn_dir + "/catalogues/"
-# raw_cat_dir = ccsn_cat_dir + "raw/"
-#
-# # sn_cats = ["IIn", "IIp", "Ibc"]
-# sn_cats = ["IIn"]
-#
-# sn_times = [100., 300., 1000.]
-# sn_times = [300.]
-#
-# def updated_sn_catalogue_name(sn_type, nearby=True):
-#     sn_name = sn_type + "_"
-#
-#     if nearby:
-#         sn_name += "nearby.npy"
-#     else:
-#         sn_name += "distant.npy"
-#
-#     return ccsn_cat_dir + sn_name
-#
-# def sn_time_pdfs(sn_type):
-#
-#     time_pdfs = []
-#
-#     for i in sn_times:
-#         time_pdfs.append(
-#             {
-#                 "time_pdf_name": "box",
-#                 "pre_window": 0,
-#                 "post_window": i
-#             }
-#         )
-#
-#     if sn_type == "Ibc":
-#         time_pdfs.append(
-#             {
-#                 "time_pdf_name": "box",
-#                 "pre_window": 20,
-#                 "post_window": 0
-#             }
-#         )
-#
-#     return time_pdfs

flarestack/analyses/ccsn/necker_2019/box_sensitivity.py

@@ -43,10 +43,10 @@
 
 # gammas = [1.8, 1.9, 2.0, 2.1, 2.3, 2.5, 2.7]
 # gammas = [1.8, 2.0, 2.5]
-gammas = [2., 2.5]
+# gammas = [2., 2.5]
+gammas = [2.]
 
 # Base name
-
 mh_name = 'fit_weights'
 pdf_type = 'box'
 
@@ -125,8 +125,8 @@
                     np.sin(closest_src["dec_rad"]), gamma=gamma
                 ) * 40 * math.sqrt(float(len(catalogue)))) * 200.) / length
 
-                if cat == 'IIP' and pdf_type != 100:
-                    scale *= 0.4
+                if cat == 'IIP':# and pdf_type != 100:
+                    scale *= 0.2
                 if cat == 'Ibc':
                     scale *= 0.5
                 if cat == 'IIn' and gamma == 2.5 and pdf_time == 1000:

flarestack/analyses/ccsn/necker_2019/decay_sensitivity.py

@@ -37,14 +37,14 @@
     "energy_pdf_name": "power_law",
 }
 
-cluster = 50
+cluster = 500
 
 # Spectral indices to loop over
 
 # gammas = [1.8, 1.9, 2.0, 2.1, 2.3, 2.5, 2.7]
 # gammas = [1.8, 2.0, 2.5]
-gammas = [2., 2.5]
-
+# gammas = [2., 2.5]
+gammas = [2.]
 # Base name
 
 mh_name = 'fit_weights'
@@ -121,12 +121,14 @@
 
                 length = float(time_key)
 
-                scale = 0.5 * (flux_to_k(reference_sensitivity(
+                scale = 0.025 * (flux_to_k(reference_sensitivity(
                     np.sin(closest_src["dec_rad"]), gamma=gamma
                 ) * 40 * math.sqrt(float(len(catalogue)))) * 200.) / length
 
                 if (gamma == 2.5) and (cat == 'IIn'):
                     scale *= 1.8
+                if cat == 'IIn':
+                    scale *= 5
 
                 injection_energy = dict(llh_energy)
                 injection_energy["gamma"] = gamma
@@ -146,17 +148,17 @@
                     "inj_dict": inj_dict,
                     "llh_dict": llh_dict,
                     "scale": scale,
-                    "n_trials": 50/cluster if cluster else 1000,
+                    "n_trials": 500/cluster if cluster else 1000,
                     "n_steps": 10
                 }
                 # !!! number of total trial is ntrials * n_jobs !!!
 
                 job_id = None
-                job_id = analyse(mh_dict,
-                                 cluster=True if cluster else False,
-                                 n_cpu=1 if cluster else 32,
-                                 n_jobs=cluster,
-                                 h_cpu='00:59:59')
+                # job_id = analyse(mh_dict,
+                #                  cluster=True if cluster else False,
+                #                  n_cpu=1 if cluster else 32,
+                #                  n_jobs=cluster,
+                #                  h_cpu='00:59:59')
                 job_ids.append(job_id)
 
                 time_res[gamma] = mh_dict

flarestack/analyses/ccsn/stasik_2017/README.txt

@@ -1,5 +1,5 @@
 Name: Core Collapse Supernova Analysis
-Analyser: Alex Stasik (original) Robert Stein (reproduced) (robert.stein@desy.de)
+Analyser: Alex Stasik (original) Robert Stein (reproduced) (robert.stein@desy.de) Jannis Necker (reproduced) (jannis.necker@desy.de)
 Wiki Page: https://wiki.icecube.wisc.edu/index.php/Point_Source_Supernova_Stacking_Analysis
 Datasets: ps_tracks (v002_p01)
 
@@ -18,3 +18,20 @@ calculated using calculate_sensitivity.py, and the results were unblinded.
 As in the original analysis, no significant excess was found, and resulting
 upper limits were calculated accordingly.
 
+After the results of Alex Stasik and Flarestack were found to not agree, I (Jannis Necker), did a thorough investigation.
+I copied the original catalogue from the wiki page on order to be able to also reproduce
+the sensitivities for the decay model. A link to the catalogue file I used can be found on the Wiki of the Update
+(https://wiki.icecube.wisc.edu/index.php/Point_Source_Supernova_Stacking_Analysis/Update)
+The results of the sensitivity calculation and a comparison to Alex Stasik's results can be found on the wiki.
+
+REPRODUCICIBILITY GUIDE
+
+If the .npy catalogue files are not present, download the original catalogues as .csv from the file referenced in the
+update wiki. Copy them to /path/to/flarestack/analyses/ccsn/stasik2017/catalogues/raw and execute
+build_catalogues_from_raw.py in the same directory.
+
+For the sensitivity calculation run calculate_sensitivity_reproduce_stasik.py
+(calculate_sensitivity_reproduce_stasik_decay.py respectively for the decay model sensitivities)
+
+All plots can be found under
+/path/to/scratch/flarestack_data/output/plots/analyses/ccsn/stasik2107/fit_weights/