Implement CutoffPowerLaw Fluxes with public data

[ssclafani949]

This pull request allows for the use of the cutoff power law function with public data. The major fix is the implementation of the get_integral method in the flux_model. This uses the default implementation. It could be sped up by implementing an analytic form.

The second update is to edit the create_analysis script so that if you give a cutoff value it will set the analysis up using the CutoffPowerLawEnergyFluxProfile, otherwise it will default to the power law implementation.

[ssclafani949]

Seems like there is an error in the integration when the cutoff is too low. Not sure where it is breaking yet. Full traceback:

Cutoff
/home/ssclafani/SNR_public_data_release/./skyllh/skyllh/analyses/i3/publicdata_ps/signalpdf.py:345: RuntimeWarning: divide by zero encountered in divide
  my_fluxmodel.energy_profile.get_integral(
/home/ssclafani/SNR_public_data_release/./skyllh/skyllh/analyses/i3/publicdata_ps/signalpdf.py:345: RuntimeWarning: invalid value encountered in divide
  my_fluxmodel.energy_profile.get_integral(
Traceback (most recent call last):
  File "/home/ssclafani/SNR_public_data_release/SNR_Example.py", line 83, in <module>
    mu, sens_flux = get_sens(args.ra, args.dec, args.gamma, args.cutoff, args.nsigma)
  File "/home/ssclafani/SNR_public_data_release/SNR_Example.py", line 38, in get_sens
    ana = create_analysis(cfg=cfg, datasets=datasets, 
  File "/home/ssclafani/SNR_public_data_release/./skyllh/skyllh/analyses/i3/publicdata_ps/time_integrated_ps.py", line 407, in create_analysis
    energy_sigpdfset = PDSignalEnergyPDFSet(
  File "/home/ssclafani/SNR_public_data_release/./skyllh/skyllh/analyses/i3/publicdata_ps/signalpdf.py", line 415, in __init__
    pdf_list = parallelize(
  File "/home/ssclafani/SNR_public_data_release/./skyllh/skyllh/core/multiproc.py", line 254, in parallelize
    result_list = master_wrapper(
  File "/home/ssclafani/SNR_public_data_release/./skyllh/skyllh/core/multiproc.py", line 227, in master_wrapper
    result_list.append(func(*args, **kwargs))
  File "/home/ssclafani/SNR_public_data_release/./skyllh/skyllh/analyses/i3/publicdata_ps/signalpdf.py", line 398, in create_energy_pdf
    [
  File "/home/ssclafani/SNR_public_data_release/./skyllh/skyllh/analyses/i3/publicdata_ps/signalpdf.py", line 399, in <listcomp>
    create_reco_e_pdf_for_true_e(i, true_e_idx)
  File "/home/ssclafani/SNR_public_data_release/./skyllh/skyllh/analyses/i3/publicdata_ps/signalpdf.py", line 392, in create_reco_e_pdf_for_true_e
    spline = FctSpline1D(p, log10_reco_e_binedges)
  File "/home/ssclafani/SNR_public_data_release/./skyllh/skyllh/analyses/i3/publicdata_ps/utils.py", line 49, in __init__
    self.spl_f = interpolate.PchipInterpolator(
  File "/home/ssclafani/.local/lib/python3.10/site-packages/scipy/interpolate/_cubic.py", line 232, in __init__
    x, _, y, axis, _ = prepare_input(x, y, axis)
  File "/home/ssclafani/.local/lib/python3.10/site-packages/scipy/interpolate/_cubic.py", line 54, in prepare_input
    raise ValueError("`y` must contain only finite values.")
ValueError: `y` must contain only finite values.

It occurs on the public data with the source:
ra = 350.8
dec = 58.8
gamma = 2.219
cutoff = 2.86

    cfg = Config()

    from skyllh.datasets.i3.PublicData_10y_ps import create_dataset_collection

    dsc = create_dataset_collection(
        cfg=cfg,
        base_path='/home/ssclafani/publicdata/')

    datasets = dsc['IC40', 'IC59', 'IC79', 'IC86_I', 'IC86_II-VII']

    from skyllh.analyses.i3.publicdata_ps.time_integrated_ps import create_analysis
    from skyllh.core.source_model import PointLikeSource

    source = PointLikeSource(ra=np.deg2rad(ra), dec=np.deg2rad(dec))
    ana = create_analysis(cfg=cfg, datasets=datasets, 
                    source=source, gamma_min = 1, gamma_max = 4, 
                    refplflux_gamma=gamma, refplflux_Ec=cutoff*1000)

[chiarabellenghi]

Thanks Steve. I don't think this is related to the integration, but to the values used to create the interpolated signal energy PDF. The PChipInterpolator is quite peaky with the input values it can take. I have run into this issue already in the past, and it might already be fixed in the new data release branch. I'll have a look today

[chiarabellenghi]

Good news. I have found that the issue is related to the quad integration failing for very small values of the integral

skyllh/skyllh/core/flux_model.py

Lines 816 to 817 in 14b5485

	for (i, (E1_i, E2_i)) in enumerate(zip(E1, E2)):
	integral[i] = quad(self, E1_i, E2_i, full_output=True)[0]

Possibly easily solved by replacing that integration with something like numpy.trapz on log energies (using log energies should also improve numerical stability since we're moving across 7 orders of magnitude)

[chiarabellenghi]

This looks good after the flux integral fix and can be merged.

[ssclafani949]

Thanks Chiara. Looks good from me too. Should I merge now or wait for someone else to do so?

[chiarabellenghi]

Go ahead and merge @ssclafani949 :)