Point source likelihood

icecube_tools also provides an interface to the likelihoods often used in point source searches of the neutrino data (see this paper by Braun et al.).

\[\mathcal{L} = \prod_{i=1}^N \Bigg[ \frac{n_s}{N} \mathcal{S}(\theta_i, E_i, \gamma) + (1-\frac{n_s}{N}) \mathcal{B}(\theta_i, E_i) \Bigg],\]

where \(N\) is the total number of detected neutrino events, \(n_s\) is the expected number of source events, \(\theta\) is the neutrino direction, \(E\) is the reconstructed neutrino energy and \(\gamma\) is the source spectral index.

The point source likelihood is a mixture model with two components: one representing possible astrophysical neutrino sources, \(\mathcal{S}(\theta, E)\), and the other known background, \(\mathcal{B}(\theta, E)\). Each component has terms depending on the directional or spatial source–neutrino relationship and also the energy of the neutrinos, as higher energy neutrinos are more likely to come from astrophysical sources. Depending on the search, the energy dependence may be omitted. Also, there may be a temporal dependence added, but this is not yet implemented in icecube_tools.

Here we implement a simple likelihood and apply it to some simulated data. There are several likelihoods available, and more information can be found in the API documentation.

[1]:
import numpy as np
from matplotlib import pyplot as plt
import h5py

from icecube_tools.point_source_likelihood.spatial_likelihood import (
    EventDependentSpatialGaussianLikelihood,
)
from icecube_tools.point_source_likelihood.energy_likelihood import (
    MarginalisedEnergyLikelihood2021,
    read_input_from_file,
)
from icecube_tools.point_source_likelihood.point_source_likelihood import (
    PointSourceLikelihood,
    TimeDependentPointSourceLikelihood,
)

Spatial likelihood

We can start with the spatial/directional term. Let’s use the energy dependent spatial likelihood. It is build from a Gaussian with an event-wise uncertainty sampled from the IRF data. The background case will simply be an isotropic distribution on the sphere.

[2]:
angular_resolution = 1  # deg
spatial_likelihood = EventDependentSpatialGaussianLikelihood(angular_resolution)

We show the likelihood profile for a single event with an assumed uncertainty of 1 degree.

[3]:
source_coord = (np.pi, np.deg2rad(30))
test_coords = [(np.pi + _, source_coord[1]) for _ in np.linspace(-0.1, 0.1, 100)]

fig, ax = plt.subplots()
ax.plot(
    np.rad2deg([tc[0] for tc in test_coords]),
    [spatial_likelihood(1.0, tc, source_coord) for tc in test_coords],
)
ax.axvline(
    np.rad2deg(source_coord[0]), color="k", linestyle=":", label="Source location"
)
ax.set_xlabel("RA [deg]")
ax.set_ylabel("Spatial likelihood")
ax.legend()
[3]:
<matplotlib.legend.Legend at 0x7f17d4067dc0>
../_images/notebooks_point_source_likelihood_6_1.png

Energy likelihood

Now let’s think about the energy-dependent term. The way this is handled is to marginalise over the true neutrino energies, to directly connect the reconstructed neutrino energies to the spectral index of a simple power-law source model.

Doing this properly requires a knowledge of the relationship between the true and reconstructed energies as well as the details of the power law model. The most straightforward way to implement this is to simulate the a large number of events using the Simulator and build a likelihood using the output of this simulation and MarginalisedEnergyLikelihoodFromSim. We do exactly this with pre-computed lists of events, to be found in the data subdirectory: sim_output_{index}.h5. These were simulated using point sources with spectral index index at 45 degrees declination. The likelihood is restricted to a small band of declination around the assumed source. Using the same declination for our test source, this is fine. For different source declinations further simulations would be needed to account for the declination dependence of the detector acceptance.

[4]:
index_list = [1.5, 1.7, 1.9, 2.1, 2.3, 2.5, 2.7, 2.9, 3.1, 3.3, 3.5, 3.7, 2.0]
energy_likelihood = MarginalisedEnergyLikelihood2021(
    index_list, "data", "p_IC86_II", np.deg2rad(30)
)
data/p_IC86_II_index_1.5.h5
data/p_IC86_II_index_1.7.h5
data/p_IC86_II_index_1.9.h5
data/p_IC86_II_index_2.0.h5
data/p_IC86_II_index_2.1.h5
data/p_IC86_II_index_2.3.h5
data/p_IC86_II_index_2.5.h5
data/p_IC86_II_index_2.7.h5
data/p_IC86_II_index_2.9.h5
data/p_IC86_II_index_3.1.h5
data/p_IC86_II_index_3.3.h5
data/p_IC86_II_index_3.5.h5
data/p_IC86_II_index_3.7.h5
[5]:
energy_likelihood.likelihood["1.5"].likelihood.shape
[5]:
(49,)
[6]:
test_energies = np.geomspace(100, 1e7)  # GeV
test_indices = [2.0, 2.5, 3.1, 3.5]

fig, ax = plt.subplots()
for index in test_indices:
    ax.plot(
        test_energies,
        [energy_likelihood(e, index) for e in test_energies],
        label=f"index {index:.1f}",
    )
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("E_reco [GeV]")
ax.set_ylabel("Energy likelihood")
ax.legend()
[6]:
<matplotlib.legend.Legend at 0x7f17c9ee0f70>
../_images/notebooks_point_source_likelihood_11_1.png

Point source likelihood

Now we can bring together the spatial and energy terms to build a full PointSourceLikelihood. First, let’s load some data from the simulation notebook to allow us to demonstrate.

[7]:
data = {}
with h5py.File("data/p_IC86_I_test_sim.h5", "r") as f:
    for key in f:
        if "source_0" not in key and "source_1" not in key:
            data[key] = f[key][()]

Let’s see how many source events are included in the simulated sample.

[8]:
np.where(data["source_label"] == 1)[0].size
[8]:
36

Now lets put our likelihood structure and data in together, along with a proposed source location:

[9]:
likelihood = PointSourceLikelihood(
    spatial_likelihood,
    energy_likelihood,
    data["ra"],
    data["dec"],
    data["reco_energy"],
    source_coord,
)

The likelihood will automatically select a declination band around the proposed source location. Because of the Gaussian spatial likelihood, neutrinos far from the source will have negligible contribution. We can control the width of this band with the optional argument band_width_factor. Let’s see how many events ended up in the band, compared to the total number:

[10]:
likelihood.N
[10]:
55
[11]:
likelihood.Ntot
[11]:
249

We also note that the background likelihood is implemented automatically, for more information on the options here, check out the API docs. This is just a function of energy, with a constant factor to account for the isotropic directional likelihood.

[12]:
fig, ax = plt.subplots()
ax.plot(test_energies, [likelihood._background_likelihood(e) for e in test_energies])
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("E_reco [GeV]")
ax.set_ylabel("Background likelihood")
[12]:
Text(0, 0.5, 'Background likelihood')
../_images/notebooks_point_source_likelihood_23_1.png

Time dependent point source likelihood

By repeating the exercise with the time dependent likelihood we can analyse data spanning multiple data periods. We define the source and its coordinates as above, pass the appropriate event_files, an index_list containing the indices of simulated spectra used to build the energy likelihoods and a directory in which these files are stored.

[21]:
source_coords = (np.pi, np.deg2rad(30))
# index_list = list(np.arange(1.5, 4.25, 0.25))
event_files = ["data/p_IC86_I_test_sim.h5", "data/p_IC86_II_test_sim.h5"]
tllh = TimeDependentPointSourceLikelihood(
    source_coords, ["IC86_I", "IC86_II"], event_files, index_list, "data"
)
IC86_I
data/p_IC86_I_index_1.5.h5
data/p_IC86_I_index_1.7.h5
data/p_IC86_I_index_1.9.h5
data/p_IC86_I_index_2.0.h5
data/p_IC86_I_index_2.1.h5
data/p_IC86_I_index_2.3.h5
data/p_IC86_I_index_2.5.h5
data/p_IC86_I_index_2.7.h5
data/p_IC86_I_index_2.9.h5
data/p_IC86_I_index_3.1.h5
data/p_IC86_I_index_3.3.h5
data/p_IC86_I_index_3.5.h5
data/p_IC86_I_index_3.7.h5
IC86_II
data/p_IC86_II_index_1.5.h5
data/p_IC86_II_index_1.7.h5
data/p_IC86_II_index_1.9.h5
data/p_IC86_II_index_2.0.h5
data/p_IC86_II_index_2.1.h5
data/p_IC86_II_index_2.3.h5
data/p_IC86_II_index_2.5.h5
data/p_IC86_II_index_2.7.h5
data/p_IC86_II_index_2.9.h5
data/p_IC86_II_index_3.1.h5
data/p_IC86_II_index_3.3.h5
data/p_IC86_II_index_3.5.h5
data/p_IC86_II_index_3.7.h5
[22]:
m = tllh._minimize()
[23]:
m
[23]:
Migrad
FCN = -503.5 Nfcn = 365
EDM = 6.59e-05 (Goal: 0.0001) time = 3.4 sec
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 n0 36 5 0 134
1 n1 41 6 0 158
2 index 2.5000 0.0025 1.55 3.65
n0 n1 index
n0 26.7 0.0185 (0.001) 7.85e-07
n1 0.0185 (0.001) 30.7 2.27e-06
index 7.85e-07 2.27e-06 6.21e-06
[ ]: