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

flarestack/analyses/agn_cores/README.md

@@ -0,0 +1,87 @@
+# Search for neutrinos from AGN cores
+
+**Name:** AGN Cores Stacking Analysis
+
+**Analyser:** Federica Bradascio (federica.bradascio@desy.de)
+
+**Unblinding Date:** May, 2020
+
+**Flarestack Release:** v2.2.0 (Titan)
+
+**Datasets:** diffuse_8_year
+
+**Wiki Page:** https://wiki.icecube.wisc.edu/index.php/Neutrinos_from_AGN_cores
+
+**Results:** https://wiki.icecube.wisc.edu/index.php?title=Radio_Galaxies_Stacking_Analysis/Results&action=edit&redlink=1
+
+This analysis uses three AGN samples: Radio-selected AGN, IR-selected AGN
+and LLAGN. It performs a stacking analysis using the "fitting weights"
+method and the X-ray flux as source weight.
+
+## STEPS TO REPRODUCE AGN CORES ANALYSIS
+
+1. Create the 3 AGN samples:
+    - The `scripts_for_unblinding/create_agn_samples/create_*agn_north_catalogue.py` must be run, one for each AGN sample
+
+2. Reproduce sensitivity/dp study for each sub-sample of N brightest sources:
+    - The three `scripts scripts_for_unblinding/*agn_analysis.py` must be run.
+    In all cases, the fit results will be obtained.
+    - To reproduce sensitivity and DP plots, the `scripts_for_unblinding/make_plots_all_samples.ipynb`
+    notebook must be used.
+
+3. (OPTIONAL) Energy range study:
+    - To reproduce the sensitivity curves as a function of the minimum energy of the injected neutrinos,
+    `scripts_for_unblinding/calculate_energy_range_lower.py` (lower energy bound) and
+    `scripts_for_unblinding/calculate_energy_range_upper.py` (upper energy bound) must be run.
+    - To reproduce energy range plots, the `scripts_for_unblinding/make_energy_range_plots.ipynb`
+    notebook must be used.
+
+4. Unbind results in terms of integral significance/upper limits:
+    - Run the `scripts_for_unblinding/unblind_integral_*agn.py` for each catalogue.
+
+5. Differential sensitivity + differential significance/upper limits (unblinding):
+    - Run `scripts_for_unblinding/differential_sensitivity/diff_sens_*agn_analysis.py`
+    - For UNBLINDING, run `scripts_for_unblinding/differential_sensitivity/unblind_differential_*agn.py`
+
+All scripts are run on the DESY (Zeuthen) clusters.
+All plots can be found in: `/path/to/scratch/flarestack_data/output/plots/analyses/agn_cores/`
+
+
+## DESCRITPION OF SCRIPTS RELEVANT FOR UNBLINDING
+
+- `scripts_for_unblinding/create_agn_samples/create_radio_selected_agn_north_catalogue.py`,
+`scripts_for_unblinding/create_agn_samples/create_ir_selected_agn_north_catalogue.py`,
+`scripts_for_unblinding/create_agn_samples/create_llagn_north_catalogue.py`:
+    takes the original catalogue FITS files in /catalogues and produces the the .npy catalogue files.
+
+- `scripts_for_unblinding/radio_selected_agn_analysis.py`, `scripts_for_unblinding/ir_selected_agn_analysis.py`,
+`scripts_for_unblinding/llagn_analysis.py`:
+    calculate the sensitivity/dp for 3 AGN samples, and for 10 sub-selections of X-ray brightest sources
+    of each AGN sample.
+
+- `scripts_for_unblinding/unblind_integral_radio_selected_agn.py`,
+`scripts_for_unblinding/unblind_integral_ir_selected_agn.py`,
+`scripts_for_unblinding/unblind_integral_llagn.py`:
+    scripts to unblind the 3 AGN catalogues as integral significance/upper limits.
+
+- `scripts_for_unblinding/make_plots_all_samples.ipynb`: notebook to create sensitivity/dp plots.
+
+- `scripts_for_unblinding/differential_sensitivity/diff_sens_radio_selected_agn_analysis.py`,
+`scripts_for_unblinding/differential_sensitivity/diff_sens_ir_selected_agn_analysis.py`,
+`scripts_for_unblinding/differential_sensitivity/diff_sens_llagn_analysis.py`:
+    calculate the sensitivity/dp for the 3 AGN samples (total number of sources)
+    in energy decades between 100 GeV and 10 PeV.
+
+- `scripts_for_unblinding/differential_sensitivity/unblind_differential_radio_selected_agn.py`,
+`scripts_for_unblinding/differential_sensitivity/unblind_differential_ir_selected_agn.py`,
+`scripts_for_unblinding/differential_sensitivity/unblind_differential_llagn.py`:
+    scripts to unblind the 3 AGN catalogues as differential significance/upper limits.
+
+- `scripts_for_unblinding/energy_range/calculate_energy_range_lower.py`,
+`scripts_for_unblinding/energy_range/calculate_energy_range_upper.py`:
+    calculate energy range relevant for this analysis, using 100 radio-selected AGN sources.
+
+- `scripts_for_unblinding/make_energy_range_plots.ipynb`:
+    notebook to create sensitivity/dp plots as a function of min/max energy bound.
+
+**All other scripts are either not relevant for the unblinding or deprecated.**

flarestack/analyses/agn_cores/scripts_for_unblinding/create_agn_samples/create_ir_selected_agn_north_catalogue.py

@@ -0,0 +1,76 @@
+"""
+Script to create catalogue entries for radio-selected AGN sample.
+
+The catalogue is given by the positional cross-match between 2RXS and NVSS,
+and XMMSL2 and NVSS. From the two subsamples 3LAC blazars have been removed.
+Also double X-Ray sources and double radio sources have been removed.
+"""
+from flarestack.analyses.agn_cores.shared_agncores import raw_cat_dir, agn_catalogue_name, agn_cores_output_dir
+from shared_agncores import create_random_src, plot_catalogue
+from flarestack.utils.prepare_catalogue import cat_dtype
+import astropy.io.fits as pyfits
+from astropy.table import Table
+import numpy as np
+import pandas as pd
+import os
+
+def select_nrandom_sources(cat, n_random=100):
+    # select n_random random sources
+    # import random
+    # list_of_random_srcs = random.sample(raw_cat, 100)
+    # raw_cat = list_of_random_srcs
+
+    df = cat.to_pandas()
+    df_random = df.sample(n=n_random)
+    cat_new = Table.from_pandas(df_random)
+    print (cat_new)
+    return cat_new
+
+def select_n_brightest_srcs(cat, nr_srcs):
+    """
+    Select the first nr_srcs brightest sources
+
+    :param cat: original catalogue of sources
+    :param nr_srcs: number of sources to select
+    :return: catalogue after selection
+    """
+    print ("Selecting", nr_srcs, "brightest sources.Length after cuts:", len(raw_cat))
+    return cat[-nr_srcs:]
+
+'''Open original (complete) catalogue'''
+raw_cat = pyfits.open(raw_cat_dir+'RadioLoudAGN_IRSelected_2rxsXmmsl2AllWise_no3LACbl_June2019_small.fits')
+raw_cat = Table(raw_cat[1].data, masked=True)
+print ("Catalogue length:", len(raw_cat))
+
+raw_cat = raw_cat[raw_cat['DEC_DEG']>-5]   # Select Northen sky sources only
+raw_cat = raw_cat.group_by('XRay_FLUX')  # order catalog by flux
+print ("Catalogue length after cut:", len(raw_cat))
+
+new_cat = np.empty(len(raw_cat), dtype=cat_dtype)
+new_cat["ra_rad"] = np.deg2rad(raw_cat["RA_DEG"])  # rosat RA in radians    #np.deg2rad(random_ra)
+new_cat["dec_rad"] = np.deg2rad(raw_cat["DEC_DEG"])  # rosat DEC in radians  #np.deg2rad(random_dec)
+new_cat["distance_mpc"] = np.ones(len(raw_cat))
+new_cat["ref_time_mjd"] = np.ones(len(raw_cat))
+new_cat["start_time_mjd"] = np.ones(len(raw_cat))
+new_cat["end_time_mjd"] = np.ones(len(raw_cat))
+
+new_cat["base_weight"] = raw_cat["XRay_FLUX"] * 1e13
+new_cat["injection_weight_modifier"] = np.ones(len(raw_cat))
+
+src_name = []
+for src, vv10 in enumerate(raw_cat['2RXS_ID']):
+    # if (vv10!='N/A'):
+    #     src_name.append(vv10)
+    if (raw_cat['2RXS_ID'][src] != 'N/A'):
+        src_name.append(raw_cat['2RXS_ID'][src])
+    elif (raw_cat['XMMSL2_ID'][src] != 'N/A'):
+        src_name.append(raw_cat['XMMSL2_ID'][src])
+    else:
+        print ("No valid name found for source nr ", src)
+        break
+new_cat["source_name"] = src_name
+print (len(new_cat))
+save_path = agn_catalogue_name("radioloud", "irselected_north")
+
+np.save(save_path, new_cat)
+print ("Saving to", save_path)

flarestack/analyses/agn_cores/scripts_for_unblinding/create_agn_samples/create_llagn_north_catalogue.py

@@ -0,0 +1,69 @@
+"""
+Script to create catalogue entries for LLAGN sample.
+
+The catalogue is given by the positional cross-match between 2RXS and AllWISE,
+and removing the 3LAC blazars. A Seyferntess PDF is assigned and only sources with
+Seyfertness larger than 0.5 are selected in the final sample.
+"""
+from flarestack.analyses.agn_cores.shared_agncores import raw_cat_dir, agn_catalogue_name, agn_cores_output_dir
+from shared_agncores import create_random_src, plot_catalogue
+from flarestack.utils.prepare_catalogue import cat_dtype
+import astropy.io.fits as pyfits
+from astropy.table import Table
+import numpy as np
+import pandas as pd
+import os
+
+def select_nrandom_sources(cat, n_random=100):
+    df = cat.to_pandas()
+    df_random = df.sample(n=n_random)
+    cat_new = Table.from_pandas(df_random)
+    print (cat_new)
+    return cat_new
+
+def select_n_brightest_srcs(cat, nr_srcs):
+    """
+    Select the first nr_srcs brightest sources
+
+    :param cat: original catalogue of sources
+    :param nr_srcs: number of sources to select
+    :return: catalogue after selection
+    """
+    print ("Selecting", nr_srcs, "brightest sources.Length after cuts:", len(raw_cat))
+    return cat[-nr_srcs:]
+
+'''Open original (complete) catalogue'''
+raw_cat = pyfits.open(raw_cat_dir+'LLAGN_2rxs2AllWiseSayfertness_no3LACbl_April2020_small.fits')
+raw_cat = Table(raw_cat[1].data, masked=True)
+print ("Catalogue length:", len(raw_cat))
+
+raw_cat = raw_cat[raw_cat['DEC_DEG']>-5]   # Select Northen sky sources only
+raw_cat = raw_cat.group_by('XRay_FLUX')  # order catalog by flux
+print ("Catalogue length after cut:", len(raw_cat))
+
+new_cat = np.empty(len(raw_cat), dtype=cat_dtype)
+new_cat["ra_rad"] = np.deg2rad(raw_cat["RA_DEG"])  # rosat RA in radians    #np.deg2rad(random_ra)
+new_cat["dec_rad"] = np.deg2rad(raw_cat["DEC_DEG"])  # rosat DEC in radians  #np.deg2rad(random_dec)
+new_cat["distance_mpc"] = np.ones(len(raw_cat))
+new_cat["ref_time_mjd"] = np.ones(len(raw_cat))
+new_cat["start_time_mjd"] = np.ones(len(raw_cat))
+new_cat["end_time_mjd"] = np.ones(len(raw_cat))
+new_cat["base_weight"] = raw_cat["XRay_FLUX"] * 1e13
+new_cat["injection_weight_modifier"] = np.ones(len(raw_cat))
+
+src_name = []
+for src, vv10 in enumerate(raw_cat['2RXS_ID']):
+    # if (vv10!='N/A'):
+    #     src_name.append(vv10)
+    if (raw_cat['2RXS_ID'][src] != 'N/A'):
+        src_name.append(raw_cat['2RXS_ID'][src])
+    elif (raw_cat['XMMSL2_ID'][src] != 'N/A'):
+        src_name.append(raw_cat['XMMSL2_ID'][src])
+    else:
+        print ("No valid name found for source nr ", src)
+        break
+new_cat["source_name"] = src_name
+
+save_path = agn_catalogue_name("lowluminosity", "irselected_north")
+
+np.save(save_path, new_cat)

flarestack/analyses/agn_cores/scripts_for_unblinding/create_agn_samples/create_radio_selected_agn_north_catalogue.py

@@ -0,0 +1,73 @@
+"""
+Script to create catalogue entries for IR-selected AGN sample.
+
+The catalogue is given by the positional cross-match between 2RXS and AllWISE,
+and XMMSL2 and AllWISE. From the two subsamples 3LAC blazars have been removed.
+Also double X-Ray sources have been removed.
+"""
+from flarestack.analyses.agn_cores.shared_agncores import raw_cat_dir, agn_catalogue_name, agn_cores_output_dir
+from shared_agncores import create_random_src, plot_catalogue
+from flarestack.utils.prepare_catalogue import cat_dtype
+import astropy.io.fits as pyfits
+from astropy.table import Table
+import numpy as np
+import pandas as pd
+import os
+
+def select_nrandom_sources(cat, n_random=100):
+    # select n_random random sources
+    # import random
+    # list_of_random_srcs = random.sample(raw_cat, 100)
+    # raw_cat = list_of_random_srcs
+
+    df = cat.to_pandas()
+    df_random = df.sample(n=n_random)
+    cat_new = Table.from_pandas(df_random)
+    print (cat_new)
+    return cat_new
+
+def select_n_brightest_srcs(cat, nr_srcs):
+    """
+    Select the first nr_srcs brightest sources
+
+    :param cat: original catalogue of sources
+    :param nr_srcs: number of sources to select
+    :return: catalogue after selection
+    """
+    print ("Selecting", nr_srcs, "brightest sources.Length after cuts:", len(raw_cat))
+    return cat[-nr_srcs:]
+
+'''Open original (complete) catalogue'''
+raw_cat = pyfits.open(raw_cat_dir+'RadioLoudAGN_RadioSelected_2rxsXmmsl2NvssAllWise_no3LACbl.fits')
+raw_cat = Table(raw_cat[1].data, masked=True)
+
+raw_cat = raw_cat[raw_cat['DEC']>-5]   # Select Northen sky sources only
+raw_cat = raw_cat.group_by('XRay_FLUX')  # order catalog by flux
+print ("Catalogue length:", len(raw_cat))
+
+new_cat = np.empty(len(raw_cat), dtype=cat_dtype)
+new_cat["ra_rad"] = np.deg2rad(raw_cat["RA"])  # NVSS RA in radians    #np.deg2rad(random_ra)
+new_cat["dec_rad"] = np.deg2rad(raw_cat["DEC"])  # NVSS DEC in radians  #np.deg2rad(random_dec)
+new_cat["distance_mpc"] = np.ones(len(raw_cat))
+new_cat["ref_time_mjd"] = np.ones(len(raw_cat))
+new_cat["start_time_mjd"] = np.ones(len(raw_cat))
+new_cat["end_time_mjd"] = np.ones(len(raw_cat))
+
+new_cat["base_weight"] = raw_cat["XRay_FLUX"] * 1e13
+new_cat["injection_weight_modifier"] = np.ones(len(raw_cat))
+
+src_name = []
+for src, vv10 in enumerate(raw_cat['2RXS_ID']):
+    # if (vv10!='N/A'):
+    #     src_name.append(vv10)
+    if (raw_cat['2RXS_ID'][src] != 'N/A'):
+        src_name.append(raw_cat['2RXS_ID'][src])
+    elif (raw_cat['XMMSL2_ID'][src] != 'N/A'):
+        src_name.append(raw_cat['XMMSL2_ID'][src])
+    else:
+        print ("No valid name found for source nr ", src)
+        break
+new_cat["source_name"] = src_name
+
+np.save(save_path, new_cat)
+print ("Saving to", save_path)