Blazar-neutrino chance coincidence simulation

We demonstrate how to run a simulation of a blazar population and diffuse neutrino flux that have no underlying connection. The simulation makes use of the configuration files provided with the package.

from popsynth.utils.logging import quiet_mode
from icecube_tools.detector.effective_area import EffectiveArea
from icecube_tools.detector.energy_resolution import EnergyResolution
from icecube_tools.detector.angular_resolution import AngularResolution

from nu_coincidence.blazar_nu.coincidence import BlazarNuCoincidenceSim
from nu_coincidence.blazar_nu.coincidence import BlazarNuCoincidenceResults
from nu_coincidence.utils.package_data import get_available_config

quiet_mode()  # silence popsynth logs/output

Firstly, we make sure that we have the necessary files for running an IceCube simulation by using the icecube_tools package

my_aeff = EffectiveArea.from_dataset("20181018")
my_aeff = EffectiveArea.from_dataset("20131121")
my_eres = EnergyResolution.from_dataset("20150820")
my_angres = AngularResolution.from_dataset("20181018")

20181018_All-sky_point-source_IceCube_da...: 100%|██████████| 11730609/11730609 [00:01<00:00, 9349396.29it/s]
20131121_Search_for_contained_neutrino_e...: 100%|██████████| 18812495/18812495 [00:01<00:00, 11761408.14it/s]
20150820_Astrophysical_muon_neutrino_flu...: 100%|██████████| 43711022/43711022 [00:09<00:00, 4593086.26it/s]

The inputs to BlazarNuCoincidenceSim are specified by the .yml config files provided with the nu_coincidence package. We can check what is available using the helper function get_available_config().

get_available_config()

['nu_connected_hese.yml',
 'fsrq_high.yml',
 'bllac_rbcu.yml',
 'nu_diffuse_tracks.yml',
 'fsrq_ref.yml',
 'nu_diffuse_ehe_hard.yml',
 'nu_diffuse_hese_soft.yml',
 'bllac_high.yml',
 'nu_connected_tracks.yml',
 'bllac_low.yml',
 'bllac_ref.yml',
 'fsrq_gpsel.yml',
 'fsrq_nobcu.yml',
 'nu_diffuse_ehe_soft.yml',
 'bllac_gpsel.yml',
 'fsrq_low.yml',
 'bllac_nobcu.yml',
 'nu_diffuse_hese.yml',
 'nu_diffuse_hese_hard.yml',
 'nu_connected_ehe.yml',
 'fsrq_rbcu.yml',
 'nu_diffuse_ehe.yml']

For this example, we can use the reference blazar model (bllac_ref.yml and fsrq_ref.yml) together with the diffuse IceCube HESE and EHE alerts models (nu_diffuse_hese.yml and nu_diffuse_ehe.yml).

sim = BlazarNuCoincidenceSim(
    file_name="output/test_coincidence_sim.h5",  # output file name
    N=2,  # number of independent simulations to run
    bllac_config="bllac_ref.yml",
    fsrq_config="fsrq_ref.yml",
    nu_hese_config="nu_diffuse_hese.yml",
    nu_ehe_config="nu_diffuse_ehe.yml",
    seed=42,
)

sim.run(parallel=False)  # define number of parallel jobs to run

The results are organised to store relevant information on the number of different coincidences, divided into the results.bllac and results.fsrq populations. We study spatial, variable and flaring coincidences. In just 2 simulations, we expect to see spatial and variable coincidences, but flaring coincidences are unlikely. In the event that flaring coincidences occur, several other population/neutrino properties are also stored for further analysis.

results = BlazarNuCoincidenceResults.load(["output/test_coincidence_sim.h5"])

results.bllac

OrderedDict([('n_spatial', array([ 5., 15.])),
             ('n_spatial_astro', array([1., 5.])),
             ('n_variable', array([0., 1.])),
             ('n_variable_astro', array([0., 0.])),
             ('n_flaring', array([0., 0.])),
             ('n_flaring_astro', array([0., 0.])),
             ('matched_flare_amplitudes', array([], dtype=float64)),
             ('matched_nu_ras', []),
             ('matched_nu_decs', []),
             ('matched_nu_ang_errs', []),
             ('matched_nu_times', []),
             ('pop_ras', []),
             ('pop_decs', []),
             ('pop_fluxes', [])])

As we ran 2 simulations, we have a coincidence count for each case.

results.bllac["n_spatial"]

array([ 5., 15.])

results.fsrq

OrderedDict([('n_spatial', array([2., 8.])),
             ('n_spatial_astro', array([0., 4.])),
             ('n_variable', array([0., 2.])),
             ('n_variable_astro', array([0., 1.])),
             ('n_flaring', array([0., 0.])),
             ('n_flaring_astro', array([0., 0.])),
             ('matched_flare_amplitudes', array([], dtype=float64)),
             ('matched_nu_ras', []),
             ('matched_nu_decs', []),
             ('matched_nu_ang_errs', []),
             ('matched_nu_times', []),
             ('pop_ras', []),
             ('pop_decs', []),
             ('pop_fluxes', [])])