cleanup benchmark analysis files

flarestack/analyses/benchmarks/benchmark_point_source_sensitivity_11yr_NT_KDE+gamma.py

@@ -0,0 +1,195 @@
+import numpy as np
+from flarestack import ResultsHandler, MinimisationHandler
+from flarestack.data.icecube import nt_v005_p01
+from flarestack.shared import plot_output_dir, flux_to_k
+from flarestack.utils.prepare_catalogue import ps_catalogue_name
+from flarestack.icecube_utils.reference_sensitivity import (
+    reference_sensitivity,
+    reference_discovery_potential,
+)
+import matplotlib.pyplot as plt
+from flarestack import analyse, wait_cluster
+import logging
+import sys, pickle
+
+logging.basicConfig(level=logging.INFO)
+
+# Initialise Injectors/LLHs
+
+injection_energy = {
+    "energy_pdf_name": "power_law",
+    "gamma": 2.0,
+}
+
+injection_time = {
+    "time_pdf_name": "steady",
+}
+
+injection_spatial = {
+    "spatial_pdf_name": "northern_tracks_kde",
+    "spatial_pdf_data": "/lustre/fs22/group/icecube/data_mirror/northern_tracks/version-005-p01/KDE_PDFs_v007/IC86_pass2/sig_E_psi_photospline_v006_4D.fits",
+}
+
+inj_kwargs = {
+    "injection_energy_pdf": injection_energy,
+    "injection_sig_time_pdf": injection_time,
+    "injection_spatial_pdf": injection_spatial,
+}
+
+
+llh_time = injection_time
+llh_spatial = injection_spatial
+llh_energy = injection_energy
+
+llh_kwargs = {
+    "llh_name": "standard_kde_enabled",
+    "llh_energy_pdf": llh_energy,
+    "llh_spatial_pdf": llh_spatial,
+    "llh_sig_time_pdf": llh_time,
+    "llh_bkg_time_pdf": {"time_pdf_name": "steady"},
+    "negative_ns_bool": True,
+}
+
+name = "analyses/benchmarks/ps_sens_10yr_NT_KDE+gamma_v2"
+
+# sindecs = np.linspace(0.90, -0.90, 3)
+sindecs = np.linspace(0.95, -0.05, 11)
+# sindecs = np.linspace(0.5, -0.5, 3)
+#
+analyses = []
+
+cluster = True
+
+job_ids = []
+
+for sindec in sindecs:
+    cat_path = ps_catalogue_name(sindec)
+
+    subname = name + "/sindec=" + "{0:.2f}".format(sindec) + "/"
+
+    scale = (flux_to_k(reference_sensitivity(sindec)) * 2)[0]
+
+    mh_dict = {
+        "name": subname,
+        "mh_name": "fixed_weights",
+        "dataset": nt_v005_p01,
+        "catalogue": cat_path,
+        "inj_dict": inj_kwargs,
+        "llh_dict": llh_kwargs,
+        "scale": scale,
+        "n_trials": 2000,
+        "n_steps": 15,
+    }
+
+    if len(sys.argv) > 1 and sys.argv[1] == "--run-jobs":
+        job_id = analyse(
+            mh_dict,
+            cluster=cluster,
+            n_cpu=1 if cluster else 16,
+            h_cpu="23:59:59",
+            ram_per_core="8.0G",
+            remove_old_logs=False,
+            trials_per_task=50,
+        )
+        job_ids.append(job_id)
+
+    analyses.append(mh_dict)
+
+if len(sys.argv) > 1 and sys.argv[1] == "--run-jobs":
+    wait_cluster(job_ids)
+
+
+sens = []
+sens_err = []
+disc_pots = []
+
+for rh_dict in analyses:
+    rh = ResultsHandler(rh_dict)
+    sens.append(rh.sensitivity)
+    sens_err.append(rh.sensitivity_err)
+    disc_pots.append(rh.disc_potential)
+
+with open(plot_output_dir(name) + "/PS10yrKDE+gamma.pkl", "wb") as pf:
+    pickle.dump(
+        {
+            "sin_dec": sindecs,
+            "sensitivity": list(sens),
+            "sens_error": [list(x) for x in list(sens_err)],
+            "disc_potential": list(disc_pots),
+        },
+        pf,
+    )
+
+sens_err = np.array(sens_err).T
+
+# sens = reference_sensitivity(sindecs, sample="7yr")
+# disc_pots = reference_discovery_potential(sindecs, sample="7yr")
+# sens_err = 0.1*sens
+
+plot_range = np.linspace(-0.99, 0.99, 1000)
+
+plt.figure()
+ax1 = plt.subplot2grid((4, 1), (0, 0), colspan=3, rowspan=3)
+ax1.plot(
+    sindecs,
+    reference_sensitivity(sindecs, sample="10yr"),
+    color="blue",
+    label=r"10-year Point Source analysis",
+)
+# ax1.plot(sindecs, reference_sensitivity(sindecs, sample="7yr"), color="green",
+#          label=r"7-year Point Source analysis")
+
+# ax1.plot(sindecs, sens, color='orange', label="Flarestack")
+
+ax1.errorbar(
+    sindecs, sens, yerr=sens_err, color="orange", label="Sensitivity", marker="o"
+)
+
+ax1.plot(
+    sindecs,
+    reference_discovery_potential(sindecs, sample="10yr"),
+    color="blue",
+    linestyle="--",
+)
+
+ax1.plot(
+    sindecs, disc_pots, color="orange", linestyle="--", label="Discovery Potential"
+)
+
+ax1.set_xlim(xmin=-1.0, xmax=1.0)
+# ax1.set_ylim(ymin=1.e-13, ymax=1.e-10)
+ax1.grid(True, which="both")
+ax1.semilogy(nonpositive="clip")
+ax1.set_ylabel(r"Flux Strength [ GeV$^{-1}$ cm$^{-2}$ s$^{-1}$ ]", fontsize=12)
+
+ax2 = plt.subplot2grid((4, 1), (3, 0), colspan=3, rowspan=1, sharex=ax1)
+
+sens_ratios = np.array(sens) / reference_sensitivity(sindecs, sample="10yr")
+sens_ratio_errs = sens_err / reference_sensitivity(sindecs, sample="10yr")
+
+disc_ratios = np.array(disc_pots) / reference_discovery_potential(
+    sindecs, sample="10yr"
+)
+
+ax2.errorbar(sindecs, sens_ratios, yerr=sens_ratio_errs, color="red", marker="o")
+ax2.scatter(sindecs, disc_ratios, color="k")
+ax2.plot(sindecs, disc_ratios, color="k", linestyle="--")
+ax2.set_ylabel(r"ratio", fontsize=12)
+ax2.set_xlabel(r"sin($\delta$)", fontsize=12)
+
+plt.suptitle("Point Source Sensitivity (10 year)")
+#
+ax1.set_xlim(xmin=-1.0, xmax=1.0)
+
+xticklabels = ax1.get_xticklabels()
+plt.setp(xticklabels, visible=False)
+plt.subplots_adjust(hspace=0.001)
+
+ax1.legend(loc="upper right", fancybox=True, framealpha=1.0)
+
+savefile = plot_output_dir(name) + "/PS10yrKDE+gamma.pdf"
+
+logging.info(f"Saving to {savefile}")
+
+plt.savefig(savefile)
+plt.close()

flarestack/data/icecube/northern_tracks/nt_v005_p00.py

@@ -1,32 +1,60 @@
-from flarestack.data.icecube.ic_season import IceCubeDataset, get_dataset_dir
+from flarestack.data.icecube.ic_season import IceCubeDataset, icecube_dataset_dir
 from flarestack.data.icecube.northern_tracks import (
     NTSeasonNewStyle,
     get_diffuse_binning,
 )
 
-icecube_dataset_dir = get_dataset_dir()
+nt_data_dir = icecube_dataset_dir / "northern_tracks/version-005-p00"
 
-nt_data_dir = icecube_dataset_dir + "northern_tracks/version-005-p00/"
+sample_name = "northern_tracks_v005_p00"
 
-nt_v005_p00 = IceCubeDataset()
+dataset_name = "icecube." + sample_name
 
-sample_name = "northern_tracks_v005_p00"
+nt_v005_p00 = IceCubeDataset(name=dataset_name)
+
+IC86_start_year = 2011
+IC86_stop_year = 2019
+IC86_timerange = range(IC86_start_year, IC86_stop_year + 1)
+
+seasons = [f"IC86_{yr}" for yr in IC86_timerange]
+
+# ==================================
+# Add individual years as subseasons
+# ==================================
 
 
 def generate_diffuse_season(name):
     season = NTSeasonNewStyle(
         season_name=name,
         sample_name=sample_name,
-        exp_path=nt_data_dir + f"{name}_exp.npy",
-        mc_path=nt_data_dir + "IC86_pass2_MC.npy",
-        grl_path=nt_data_dir + f"GRL/{name}_exp.npy",
+        exp_path=nt_data_dir / f"{name}_exp.npy",
+        mc_path=nt_data_dir / "IC86_pass2_MC.npy",
+        grl_path=nt_data_dir / f"GRL/{name}_exp.npy",
         sin_dec_bins=get_diffuse_binning(name)[0],
         log_e_bins=get_diffuse_binning(name)[1],
     )
-    nt_v005_p00.add_season(season)
-
+    return season
 
-seasons = [f"IC86_201{i}" for i in [1, 2, 3, 4, 5, 6, 7, 8, 9]]
 
 for season in seasons:
-    generate_diffuse_season(season)
+    subseason = generate_diffuse_season(season)
+    nt_v005_p00.add_subseason(subseason)
+
+
+# ==================================
+# Add combo season
+# ==================================
+
+name = "IC86_1-9"
+
+combo_season = NTSeasonNewStyle(
+    season_name=name,
+    sample_name=sample_name,
+    exp_path=[nt_data_dir / f"IC86_{yr}_exp.npy" for yr in IC86_timerange],
+    mc_path=nt_data_dir / "IC86_pass2_MC.npy",
+    grl_path=[nt_data_dir / f"GRL/IC86_{yr}_exp.npy" for yr in IC86_timerange],
+    sin_dec_bins=get_diffuse_binning(name)[0],
+    log_e_bins=get_diffuse_binning(name)[1],
+)
+
+nt_v005_p00.add_season(combo_season)