Feature/differential sensitivity

[chiarabellenghi]

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.
It can also be passed to the methods that generate trials or estimate sensitivity as sig_kwargs

[chiarabellenghi]

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]

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.

[chiarabellenghi]

@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