Feature/differential sensitivity by chiarabellenghi · Pull Request #253 · icecube/skyllh

chiarabellenghi · 2026-03-12T15:50:41Z

Feature to address issue #249

This PR implements the option to inject a neutrino signal in a user-defined energy range.

The energy_range is now a property of the public data signal generator and of the analysis object.

chiarabellenghi · 2026-03-12T17:36:24Z

The PR should be working, but it needs some testing. Besides me, I assigned @adesai90 so that he can check that the code is now producing the differential sensitivity that he's willing to compute.

Maybe we could use Figure 5 of this paper to cross-check that the values we get are approximately correct, although that paper used only 7 years of data, whereas we have 10 (partially reprocessed) years in the public data release.

chiarabellenghi · 2026-03-24T18:00:17Z

Thinking about this a bit more... The energy range should be a property of the signal generator that we can change at run time without re-instantiating the analysis object.

The public data implementation would naturally support an architecture in which the energy_range is treated as a kwarg that is passed at runtime. However, the function that calculates the factor to convert a number of signal events into a flux is a property of the analysis object, which doesn't use the signal generator.

What we should do:

The general SkyLLH architecture should treat the energy_range as a property of the signal generator only
Then its value could be updated at run time when generating trials or simply simulating pseudo-data. The change would then be propagated to the signal_generator.energy_range property.
The conversion factor to go from ns to flux should be calculated by an analysis method. This analysis method should take energy_range as an optional argument that defaults to None. When None, the energy_range stored in the analysis's signal generator is used.

Note that the signal generator is created when the create_analysis method is called. So the user does not have to do much apart from creating the analysis and calling the methods they want.

…l_sensitivity

chiarabellenghi · 2026-03-27T17:34:45Z

@adesai90, the conversion from ns to flux should be fixed now!

There was some code restructuring involved. So here's a summary of how it works now:

You pass the energy range to the create_analysis method when you create the analysis:

ana = create_analysis(
    cfg=cfg,
    datasets=datasets,
    source=source,
    refplflux_gamma=2.0,
    energy_range=(1e2, 1e3))

Internally, this sets the signal generator's energy_range property. The signal generator uses the flux model specified when initializing the analysis. So if you initialized a power law with spectral index 2.0, that's what the signal generator uses to generate signal events. Then you can change the energy range by assigning a new value to the ana.energy_range property. Here's an example:

N = 100

rss = RandomStateService(0)
n_sig_erange, n_events_list_erange, events_list_erange = ana.generate_signal_events(rss, N)

rss = RandomStateService(0)
ana.energy_range = (1e2, 1e9)
n_sig, n_events_list, events_list = ana.generate_signal_events(rss, N)

which correctly generates the signal events:

So there is no sig_kwargs anymore. You have instead an analysis property that you can set at runtime. It will update the signal generator.

The calculate_fluxmodel_scaling_factor() method of the analysis now returns an dimentionless scaling factor. It does not take the model parameters as input anymore because it is bound to the flux model and energy range that are defined for the analysis.
But if you update the energy_range property, it will also be updated correctly. You just can't change the flux parameters (gamma in your case) anymore. For example:

ana.energy_range = (1e2, 1e9)
print(ana.calculate_fluxmodel_scaling_factor())
ana.energy_range = (1e2, 1e3)
print(ana.calculate_fluxmodel_scaling_factor())
ana.energy_range = (1e5, 1e6)
print(ana.calculate_fluxmodel_scaling_factor())

returns now different values:
5.627399490261528e-17
4.937587814014621e-21
2.8880126503503215e-18

I also added 2 utility methods that are just wrappers around ana.calculate_fluxmodel_scaling_factor(), so they also correctly use the energy_range property:

ana.mu2flux(mu) returns a flux normalization at the chosen refplflux_E0 (1TeV by default)
ana.flux2mu(flux_norm) returns the mean number of signal events for the given flux normalization value at the chosen refplflux_E0

…l_sensitivity

adesai90 · 2026-04-07T17:09:03Z

also be updated correctly. You just can't change the flux parameters (gamma in your case) anymore. Fo

Thanks, Just as an update. I tested it with 8 years of data and dec=0 at index of 2.0 to get the following result

The flux was calculated by
ana.calculate_fluxmodel_scaling_factor()*n_{sensitivity}
and then multiplied by (energy_centers/1e3)**2 where energy_centers/1e3 is the mean energy in the energy bin in TeV. If I did this I should get the flux in 1e3* TeV/cm2/s but it matches the TeV/cm2/s results from the plot below. Is the conversion right?

This can be compared to Figure 4 from https://arxiv.org/pdf/1811.07979 (attached here)

chiarabellenghi · 2026-04-08T07:15:30Z

Thanks, Just as an update. I tested it with 8 years of data and dec=0 at index of 2.0 to get the following result The flux was calculated by ana.calculate_fluxmodel_scaling_factor()*n_{sensitivity} and then multiplied by (energy_centers/1e3)**2 where energy_centers/1e3 is the mean energy in the energy bin in TeV. If I did this I should get the flux in 1e3* TeV/cm2/s but it matches the TeV/cm2/s results from the plot below. Is the conversion right?

This can be compared to Figure 4 from https://arxiv.org/pdf/1811.07979 (attached here)

I think the issue here is that you should not multiply your flux norm at 1TeV by (energy_centers/1e3)**2. Or if you do you should also multiply by the inverse.

By multiplying the flux norm by that factor, you are introducing an energy dependence on a $E^2 \times (E/E_0)^{-2}$, which should be flat in energy.

If you wanna transform your flux norm from 1/GeV/cm^2/s/sr to TeV/cm^2/s, you can simply multiply your flux norm by 1e3. This automatically also corresponds to a flux norm at 1TeV times E^2, because you are implicitly also multiplying by (1TeV/1TeV)^2. And because your spectrum is now flat, the flux norm is the same also in the tested energy bin.
(Note, the integral over the solid angle for a point source is just the integral of a delta-function, so that one is just another factor of 1.

Then, hopefully, your differential sensitivity will really look similar to Rene's!

adesai90 · 2026-04-10T16:12:36Z

-Update from Chiara for a case where gamma = 2.0 is injected, but the likelihood is free.

chiarabellenghi · 2026-04-10T16:33:39Z

With the last fix, I can reproduce Rene's sensitivity fairly well. Only the first bin (1TeV) seems better in this analysis. But I think this makes sense because Rene had a prior on the spectral index, which was forcing the likelihood to fit something around 2.13. And fitting a flat spectrum to a lot of low-energy events with no high-energy counterparts is... difficult 😅

I think this is ready for review now

Copilot

Pull request overview

This PR adds support for injecting / evaluating signal within a user-defined true-neutrino energy range (GeV), primarily targeting differential sensitivity workflows (issue #249), and refactors flux normalization conversion to use a signal-generator-provided scaling factor.

Changes:

Introduces a shared energy_range capability and propagates it through analysis + public-data signal generation components.
Adds energy-range correction factors for public-data signal generators and enforces consistent configuration across datasets when computing flux scaling.
Refactors sensitivity/discovery utilities to use calculate_fluxmodel_scaling_factor() as a per-event scaling factor and apply * mean_ns.

Reviewed changes

Copilot reviewed 9 out of 9 changed files in this pull request and generated 6 comments.

Show a summary per file

File	Description
`tests/core/test_signal_generator.py`	Adds tests for dataset-to-dataset `energy_range` consistency and kwarg policy around `MultiDatasetSignalGenerator`.
`skyllh/core/signal_generation.py`	Adds `HasEnergyRange` capability interface and clarifies energy-range units (GeV) for signal generation methods.
`skyllh/core/analysis.py`	Adds `Analysis.energy_range` property + propagation to energy-range-capable components; updates flux scaling conversion API and adds `mu2flux/flux2mu`.
`skyllh/core/signal_generator.py`	Adds `fluxmodel_scaling_factor()` to multi-dataset generators and integrates optional per-dataset/source correction factors.
`skyllh/i3/signal_generation.py`	Documentation clarification that `energy_range` is in true neutrino energy (GeV).
`skyllh/i3/utils/analysis.py`	Updates sensitivity/discovery curve utilities to compute flux scaling via `calculate_fluxmodel_scaling_factor() * mean_ns`.
`skyllh/analyses/i3/publicdata_ps/time_integrated_ps.py`	Adds `energy_range` parameter to analysis factory and applies it to the constructed analysis.
`skyllh/analyses/i3/publicdata_ps/time_integrated_ps_function_energy_spectrum.py`	Updates flux/ns conversion helpers to use the new scaling-factor API.
`skyllh/analyses/i3/publicdata_ps/signal_generator.py`	Adds `energy_range` property to the public-data dataset signal generator; rebuilds source-dependent structures and computes energy-range correction factors.

Comments suppressed due to low confidence (1)

skyllh/core/signal_generator.py:389

MultiDatasetSignalGenerator.generate_signal_events no longer accepts the src_detsigyield_weights_service keyword that is part of the SignalGenerator.generate_signal_events interface (it used to be swallowed by **kwargs). This makes the subclass signature incompatible with the base interface and breaks callers that pass that kwarg. Consider adding src_detsigyield_weights_service=None to the signature (even if unused) or otherwise maintaining signature compatibility while still rejecting unsupported kwargs like energy_range.

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

Copilot · 2026-04-10T16:39:48Z

+        # ref_N is the expected number of detected signal events for the reference flux normalization.
+        ref_N = np.sum(a_jk_eff)
+
+        if per_source:
+            a_k = np.sum(a_jk_eff, axis=0)
+            return a_k / (ref_N**2)
+
+        scaling_factor = 1.0 / ref_N
+


ref_N = np.sum(a_jk_eff) can become 0 (e.g. if an energy_range cut yields zero expected events), which will cause a division by zero in both the global and per-source branches. Please add a guard for ref_N <= 0 and raise a clear ValueError indicating that the configured signal hypothesis/energy_range produces zero detectable signal.

Fixed in c1d3407

Copilot · 2026-04-10T16:39:49Z

+        logger.warning(
+            'The fluxmodel_scaling_factor method of the MCMultiDatasetSignalGenerator class is currently only '
+            'implemented for the case of a single source or for multiple sources with the same flux model.'
+            'Multiple sources withdifferent flux models are currently not correctly handled.'
+        )


Spelling/formatting in the warning message: "withdifferent" is missing a space, and there should be a space between sentences after "flux model.".

Fixed in aa7ba5a

Copilot · 2026-04-10T16:39:49Z

+            # Convert the energy boundaries to the closest SM bin edges.
+            idx0 = self.sm.get_log10_true_e_idx(np.log10(r[0]))
+            idx1 = self.sm.get_log10_true_e_idx(np.log10(r[1]))
+            r_log10 = (self.sm.true_e_bin_edges[idx0], self.sm.true_e_bin_edges[idx1])
+
+            if r_log10[0] >= r_log10[1]:
+                raise ValueError(
+                    'The first element of the energy_range sequence must be strictly smaller than the second element!'
+                )


The bin-alignment of energy_range appears to compute the upper edge incorrectly: get_log10_true_e_idx returns a bin index, but true_e_bin_edges[idx1] is the lower edge of that bin. This makes ranges collapse when both bounds fall in the same bin (raising erroneously) and prevents selecting up to the last bin edge. Consider using true_e_bin_edges[idx1 + 1] (with clipping) for the upper boundary so the aligned range fully covers the requested E_max.

implemented in ad090bf

It makes sense to me as we also do this:

if self._energy_range_log10 is None: return (self.sm.true_e_bin_edges[0], self.sm.true_e_bin_edges[-1])

@chiarabellenghi please take a look.

I guess neither implementation passes in case the energy_range is outside true_e_bin_edges as the get_log10_true_e_idx raises an error instead of clipping indices

and also if passed energy_range exactly matches the SM binning this will take an additional bin, annoying 😓 . Now I don't know what's better

implemented in ad090bf

It makes sense to me as we also do this:

if self._energy_range_log10 is None: return (self.sm.true_e_bin_edges[0], self.sm.true_e_bin_edges[-1])

@chiarabellenghi please take a look.

I guess neither implementation passes in case the energy_range is outside true_e_bin_edges as the get_log10_true_e_idx raises an error instead of clipping indices

I think we can simply make this error clearer by saying what the full range supported by the SM is

and also if passed energy_range exactly matches the SM binning this will take an additional bin, annoying 😓 . Now I don't know what's better

Fixed in a683414. Now the sm helper function handles the two cases separately using an additional input flag.

Also added a dedicated unittest.

Copilot · 2026-04-10T16:39:49Z

+    def _calculate_flux_weight(self, src_idx, log_e_min, log_e_max):
+        """Calculates an unnormalized detectable flux weight for one source
+        and one true-energy range.
+        """
+        src = self._shg_mgr.source_list[src_idx]
+        fluxmodel = self.shg_mgr.get_fluxmodel_by_src_idx(src_idx=src_idx)
+
+        effA = PDAeff(
+            pathfilenames=self.ds.get_abs_pathfilename_list(self.ds.get_aux_data_definition('eff_area_datafile'))
+        )
+


_calculate_flux_weight constructs a new PDAeff instance for every call, and _calculate_energy_range_correction_factors calls it twice per source. Since PDAeff.__init__ loads the effective area file(s), this can become a significant performance bottleneck. Consider caching a single PDAeff instance (or reusing already loaded effective-area data) and reusing it across sources/ranges.

ah, yes. This is a good point

fixed in 18aa819

Copilot · 2026-04-10T16:39:49Z

            )

            mean_ns_arr[sin_dec_idx, iter_idx] = mean_ns
            mean_ns_err_arr[sin_dec_idx, iter_idx] = mean_ns_err
-            flux_scaling_arr[sin_dec_idx, iter_idx] = ana.calculate_fluxmodel_scaling_factor(
-                mean_ns=mean_ns, fitparam_values=np.array(fitparam_values)
-            )
+            flux_scaling_arr[sin_dec_idx, iter_idx] = ana.calculate_fluxmodel_scaling_factor() * mean_ns



fitparam_values is still a required parameter of estimate_ps_sin_dec_sensitivity_curve, but after switching to ana.calculate_fluxmodel_scaling_factor() it is no longer used (and the docstring still claims it is used for the conversion). Either remove/deprecate the parameter or update the docstring/signature to reflect the new behavior to avoid confusion.

Removed in dae2393

Yep... and this entire library has to be cleaned up at some point...

…testcase

Copilot

Pull request overview

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

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

Copilot · 2026-04-13T09:57:54Z

    )

+    # Set the energy range for signal generation if it is explicitly given.
+    ana.energy_range = energy_range


This comment says the energy range is set "if it is explicitly given", but the code assigns ana.energy_range = energy_range unconditionally (including when it is None). Consider guarding this with if energy_range is not None: (or update the comment), to avoid marking the analysis energy_range as explicitly set and to avoid unnecessary component rebuilds.

Suggested change

ana.energy_range = energy_range

if energy_range is not None:

ana.energy_range = energy_range

Fixed in 0ece010 by updating the comment

Copilot · 2026-04-13T09:57:54Z

+
+        with self.assertRaises(TypeError):
+            sig_gen.generate_signal_events(rss=rss, mean=10, poisson=False, energy_range=(1e2, 1e3))
+


This test enforces that MultiDatasetSignalGenerator.generate_signal_events rejects an energy_range kwarg. The PR description currently states that energy_range can be passed via sig_kwargs to trial/sensitivity helpers; that no longer matches the behavior. Please update the PR description (and/or user-facing docs) to reflect that energy_range must be set on the analysis/signal generator, not passed as sig_kwargs.

Copilot · 2026-04-13T09:57:55Z

+import os
+import unittest
+
+import numpy as np
+


These new public-data tests are placed under tests/publicdata_ps/, but the repository test runner (tests/run.sh) only discovers tests/core and tests/i3, so they won’t run in CI as-is. Either move them under a discovered package (e.g. tests/i3/...) or update the test runner to include this directory; if enabled in CI, also consider skipping when dataset downloads aren’t available (Config defaults download_from_origin=True).

It will be fixed in #280

tomaskontrimas

Looks good to me! 🥳

chiarabellenghi added 4 commits March 11, 2026 19:26

add energy range for signal injeciton

5ab07eb

energy_range as keyword argument for the signal generator

a7fd685

log energy_range

dde0a35

cache input energy_range state

4b7459c

chiarabellenghi requested a review from tomaskontrimas March 12, 2026 15:50

chiarabellenghi self-assigned this Mar 12, 2026

tomaskontrimas reviewed Mar 12, 2026

View reviewed changes

Comment thread skyllh/analyses/i3/publicdata_ps/signal_generator.py Outdated

Comment thread skyllh/analyses/i3/publicdata_ps/signal_generator.py Outdated

fix property definition

8d1b753

chiarabellenghi assigned adesai90 and unassigned chiarabellenghi and adesai90 Mar 12, 2026

chiarabellenghi linked an issue Mar 20, 2026 that may be closed by this pull request

Addition of Differential Sensitivity/DP #249

Closed

chiarabellenghi added 5 commits March 25, 2026 16:31

Merge remote-tracking branch 'origin/master' into feature/differentia…

acf787b

…l_sensitivity

formatting

db2abab

WIP fix flux scaling factor

9ecf3a8

formatting

f8c6103

ns-to-flux conversion with energy range

0203fde

Merge remote-tracking branch 'origin/master' into feature/differentia…

dfbe695

…l_sensitivity

chiarabellenghi added 2 commits April 8, 2026 08:53

Merge branch 'master' into feature/differential_sensitivity

885d00a

Merge branch 'master' into feature/differential_sensitivity

2af7ad0

chiarabellenghi added 2 commits April 8, 2026 09:36

Merge branch 'master' into feature/differential_sensitivity

04e5024

Merge branch 'master' into feature/differential_sensitivity

31b3383

adesai90 closed this Apr 10, 2026

adesai90 reopened this Apr 10, 2026

adesai90 marked this pull request as ready for review April 10, 2026 16:14

fix scaling factors for ns to flux conversion

b642753

Copilot AI review requested due to automatic review settings April 10, 2026 16:29

Copilot started reviewing on behalf of chiarabellenghi April 10, 2026 16:29 View session

chiarabellenghi requested a review from tomaskontrimas April 10, 2026 16:33

Copilot AI reviewed Apr 10, 2026

View reviewed changes

tomaskontrimas added 6 commits April 10, 2026 20:19

Remove no longer usedfitparam_values argument

dae2393

Fix grammar and typos

aa7ba5a

Add public data time_integrated_ps test with additional energy_range …

46536de

…testcase

Check that ref_N is strictly positive

c1d3407

Use upper true_e_bin_edges for upper energy range limit

ad090bf

Additional testcase for changing energy_range

443eaba

tomaskontrimas mentioned this pull request Apr 11, 2026

Vectorize _renorm_per_y_pairs method #298

Merged

chiarabellenghi added 2 commits April 13, 2026 10:06

Merge branch 'master' into feature/differential_sensitivity

7249636

Merge branch 'master' into feature/differential_sensitivity

9660df1

chiarabellenghi added the enhancement New feature or request label Apr 13, 2026

chiarabellenghi added 3 commits April 13, 2026 11:13

control energy_range bin alignment edge cases

a683414

cache effective area data

18aa819

Better error message

6fdc6ea

chiarabellenghi requested a review from Copilot April 13, 2026 09:52

Copilot started reviewing on behalf of chiarabellenghi April 13, 2026 09:52 View session

Copilot AI reviewed Apr 13, 2026

View reviewed changes

chiarabellenghi added 2 commits April 13, 2026 12:20

Add comment for forgetful future self

9861548

Fix create_analysis comment on energy_range

0ece010

tomaskontrimas approved these changes Apr 13, 2026

View reviewed changes

chiarabellenghi merged commit e244567 into master Apr 13, 2026
6 checks passed

chiarabellenghi deleted the feature/differential_sensitivity branch April 13, 2026 12:08

	ana.energy_range = energy_range
	if energy_range is not None:
	ana.energy_range = energy_range


		with self.assertRaises(TypeError):
		sig_gen.generate_signal_events(rss=rss, mean=10, poisson=False, energy_range=(1e2, 1e3))

Conversation

chiarabellenghi commented Mar 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Uh oh!

chiarabellenghi commented Mar 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

chiarabellenghi commented Mar 24, 2026

Uh oh!

chiarabellenghi commented Mar 27, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

adesai90 commented Apr 7, 2026

Uh oh!

chiarabellenghi commented Apr 8, 2026

Uh oh!

adesai90 commented Apr 10, 2026

Uh oh!

chiarabellenghi commented Apr 10, 2026

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

Uh oh!

Copilot AI Apr 10, 2026

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 10, 2026

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 10, 2026

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

tomaskontrimas Apr 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

chiarabellenghi Apr 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 10, 2026

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!

Copilot AI Apr 10, 2026

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!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

Uh oh!

Copilot AI Apr 13, 2026

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Copilot AI Apr 13, 2026

Choose a reason for hiding this comment

Uh oh!

chiarabellenghi commented Mar 12, 2026 •

edited

Loading

chiarabellenghi commented Mar 12, 2026 •

edited

Loading

chiarabellenghi commented Mar 27, 2026 •

edited

Loading

tomaskontrimas Apr 10, 2026 •

edited

Loading

chiarabellenghi Apr 13, 2026 •

edited

Loading