Add utility functions to get distribution of event counts in reconstructed energy given a flux model and a declination by chiarabellenghi · Pull Request #322 · icecube/skyllh

chiarabellenghi · 2026-06-01T15:12:37Z

This PR addresses issue #301.

It has been tested using the PowerLawEnergyFluxProfile flux model and a custom callable that implements the same power law. Passing one of the two gives the same results.

This is an initial implementation. For instance, it could be extended in the future to allow integration over a declination band.

Copilot

Pull request overview

This PR adds helper utilities (requested in issue #301) to compute the expected distribution of reconstructed-energy event counts for IceCube public-data point-source datasets, given a flux model and a fixed declination.

Changes:

Add utilities to fold a flux model through effective area + smearing matrix to obtain expected reco-energy counts per dataset and summed over seasons.
Make get_bin_indices_from_lower_and_upper_binedges accept scalar inputs by normalizing values with np.atleast_1d.
Add a generic math_function_str property for FunctionEnergyFluxProfile.

Reviewed changes

Copilot reviewed 3 out of 3 changed files in this pull request and generated 12 comments.

File	Description
`skyllh/core/flux_model.py`	Adds `math_function_str` for callable-based energy flux profiles.
`skyllh/core/binning.py`	Makes bin-index lookup robust to scalar `values` inputs.
`skyllh/analyses/i3/publicdata_ps/utils.py`	Introduces the new reco-energy expected-counts helper functions.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

chiarabellenghi · 2026-06-01T15:20:45Z

+    # Build a uniform 0.1-step grid in log10(E/GeV) covering all valid reco-E bins.
+    reco_e_edges = np.arange(0.0, 9.0 + 0.05, 0.1)
+    reco_e_centers = get_bincenters_from_binedges(reco_e_edges)
+    out_bin_width = reco_e_edges[1] - reco_e_edges[0]


I don't wanna do this because then combining counts from different datasets is non-trivial. And that's because different smearing matrices have different reco energy bins.

chiarabellenghi · 2026-06-02T09:41:45Z

+    # Resolve energy range limits (in GeV).
+    E_min = energy_range[0] if energy_range is not None else 10.0 ** true_e_edges[0]
+    E_max = energy_range[1] if energy_range is not None else 10.0 ** true_e_edges[-1]


This is now checked in the public function. See 72441d8

+        log10_E_centre = 0.5 * (log10_E1 + log10_E2)
+        aeff_e_idx = get_bin_indices_from_lower_and_upper_binedges(aeff_e_bins_lo, aeff_e_bins_hi, log10_E_centre)[0]
+        A_eff = aeff_dec[aeff_e_idx]  # cm^2


chiarabellenghi · 2026-06-02T09:42:08Z

+        # Marginalise smearing matrix over psi and ang_err axes:
+        p_reco = sm.histogram[i, sm_dec_idx, :, :, :].sum(axis=(-1, -2))  # (n_reco_e,)
+        if p_reco.sum() == 0.0:
+            continue
+
+        # Convert probability to probability density in log10(E_reco) space.
+        bw = reco_e_hi[i, sm_dec_idx, :] - reco_e_lo[i, sm_dec_idx, :]  # bin widths
+        pdf = np.where(bw > 0, p_reco / bw, 0.0)  # (n_reco_e,)
+
+        # Original reco-energy bin edges for this (true_e_i, dec) slice.
+        # Note: this might break for zero-binwidth bins, but they should never happen for the reco energies.
+        original_edges = np.append(reco_e_lo[i, sm_dec_idx, :], reco_e_hi[i, sm_dec_idx, -1])
+
+        # 1D spline interpolation of the PDF onto the uniform output grid.
+        spl = FctSpline1D(pdf, original_edges)
+        counts_per_sec += weight * spl(reco_e_centers) * out_bin_width


Fixed in 72441d8

chiarabellenghi · 2026-06-02T09:42:25Z

+    if livetimes is not None and len(livetimes) != len(datasets):
+        raise ValueError('Length of livetimes must match the number of datasets.')
+
+    counts_total = 0
+
+    # Loop over datasets and sum contributions accounting for livetimes.
+    for i, ds in enumerate(datasets):
+        if livetimes is not None:
+            if livetimes[i] is None:
+                data = ds.load_data()
+                livetime = data.livetime
+            else:
+                livetime = livetimes[i]
+        else:
+            data = ds.load_data()
+            livetime = data.livetime
+
+        reco_e_edges, counts_per_sec = reco_energy_counts_per_second_per_ds(
+            ds, flux, dec, Phi0=Phi0, energy_range=energy_range
+        )
+        counts_total += counts_per_sec * livetime * 24 * 3600  # Convert days to seconds
+    return reco_e_edges, counts_total


Fixed in 72441d8

chiarabellenghi · 2026-06-02T09:42:47Z

+    # Resolve energy range limits (in GeV).
+    E_min = energy_range[0] if energy_range is not None else 10.0 ** true_e_edges[0]
+    E_max = energy_range[1] if energy_range is not None else 10.0 ** true_e_edges[-1]


Fixed in 72441d8

chiarabellenghi · 2026-06-02T08:00:58Z

+        log10_E_centre = 0.5 * (log10_E1 + log10_E2)
+        aeff_e_idx = get_bin_indices_from_lower_and_upper_binedges(aeff_e_bins_lo, aeff_e_bins_hi, log10_E_centre)[0]
+        A_eff = aeff_dec[aeff_e_idx]  # cm^2


This also does not happen because in true energy, the sm and aeff ranges are the same

chiarabellenghi · 2026-06-01T15:27:51Z

+        # Marginalise smearing matrix over psi and ang_err axes:
+        p_reco = sm.histogram[i, sm_dec_idx, :, :, :].sum(axis=(-1, -2))  # (n_reco_e,)
+        if p_reco.sum() == 0.0:
+            continue
+
+        # Convert probability to probability density in log10(E_reco) space.
+        bw = reco_e_hi[i, sm_dec_idx, :] - reco_e_lo[i, sm_dec_idx, :]  # bin widths
+        pdf = np.where(bw > 0, p_reco / bw, 0.0)  # (n_reco_e,)
+
+        # Original reco-energy bin edges for this (true_e_i, dec) slice.
+        # Note: this might break for zero-binwidth bins, but they should never happen for the reco energies.
+        original_edges = np.append(reco_e_lo[i, sm_dec_idx, :], reco_e_hi[i, sm_dec_idx, -1])
+
+        # 1D spline interpolation of the PDF onto the uniform output grid.
+        spl = FctSpline1D(pdf, original_edges)
+        counts_per_sec += weight * spl(reco_e_centers) * out_bin_width


As written in the comment note, I believe this never actually happens. But yes, it is a liability in principle. Just not sure I want to overcomplicate the code for something that currently can't happen...

chiarabellenghi · 2026-06-02T09:43:18Z

+    if livetimes is not None and len(livetimes) != len(datasets):
+        raise ValueError('Length of livetimes must match the number of datasets.')
+
+    counts_total = 0
+
+    # Loop over datasets and sum contributions accounting for livetimes.
+    for i, ds in enumerate(datasets):
+        if livetimes is not None:
+            if livetimes[i] is None:
+                data = ds.load_data()
+                livetime = data.livetime
+            else:
+                livetime = livetimes[i]
+        else:
+            data = ds.load_data()
+            livetime = data.livetime
+
+        reco_e_edges, counts_per_sec = reco_energy_counts_per_second_per_ds(
+            ds, flux, dec, Phi0=Phi0, energy_range=energy_range
+        )
+        counts_total += counts_per_sec * livetime * 24 * 3600  # Convert days to seconds
+    return reco_e_edges, counts_total


Fixed in 72441d8

chiarabellenghi · 2026-06-02T09:43:26Z

+    reco_e_edges : ndarray, shape (n_bins + 1,)
+        Bin edges of the output reconstructed log10(E/GeV) axis.
+    counts_per_sec : ndarray, shape (n_bins,)
+        Expected counts per second in each reco-energy bin.
+    """


Fixed in 72441d8

tomaskontrimas · 2026-06-01T19:33:01Z

What do you think about making reco_energy_counts_per_second_per_ds an internal function with e.g. _ prefix and then making sure that compute_expected_reco_energy_counts accepts one dataset too?

You could then load aeff/sm (and temporarily cache them) there to figure out the common bounds for reco_e_edges and then pass them to _ reco_energy_counts_per_second_per_ds. This way also the energy_range argument validation could be performed only once, as I think we should assume that internal function calls are correct within the framework and only validate user input

tomaskontrimas

Looks good to me! 🚀

Address #301

1d041d2

chiarabellenghi requested review from Copilot and tomaskontrimas June 1, 2026 15:12

chiarabellenghi added the enhancement New feature or request label Jun 1, 2026

Copilot started reviewing on behalf of chiarabellenghi June 1, 2026 15:12 View session

Copilot AI reviewed Jun 1, 2026

View reviewed changes

tomaskontrimas assigned chiarabellenghi Jun 1, 2026

tomaskontrimas linked an issue Jun 1, 2026 that may be closed by this pull request

Helper function to get number of observed events as a function of reco energy given a flux model #301

Closed

chiarabellenghi and others added 5 commits June 2, 2026 11:37

Address PR comments

72441d8

Fix docstring

436b8ce

avoid initializing empty list for smearing matrices

79a6321

Avoid initializing lists anywhere if possible

a39304c

Fix docstring

a02d2da

tomaskontrimas approved these changes Jun 2, 2026

View reviewed changes

chiarabellenghi merged commit 93f86a9 into master Jun 2, 2026
5 checks passed

chiarabellenghi deleted the reco_events_from_fluxmodel branch June 2, 2026 10:06

Uh oh!

Conversation

chiarabellenghi commented Jun 1, 2026

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

tomaskontrimas commented Jun 1, 2026

Uh oh!

tomaskontrimas left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants