Fixing the Spectral Index

In this tutorial, we will show how to initialise an analysis instance with a fixed spectral index.

First, we initialize the logging, config, datasets, and source objects. Details on these can be found in the Fitting a steady point-source with the public 14-year IceCube track data and Dataset collections tutorials.

import numpy as np

import skyllh
from skyllh.core.config import Config
from skyllh.core.logging import setup_logging
from skyllh.core.source_model import PointLikeSource

cfg = Config()
logger = setup_logging(cfg=cfg, name='fixed_spectral_index')
datasets = skyllh.create_datasets('IceTracks-DR2', cfg=cfg)

# Location of NGC 1068
src_ra = 40.67  # degrees
src_dec = -0.01  # degrees

source = PointLikeSource(ra=np.radians(src_ra), dec=np.radians(src_dec))

2026-05-13 19:15:18,993 MainProcess skyllh.datasets.datasets INFO: Loaded 4 dataset(s) from sample "IceTracks-DR2": IC40, IC59, IC79, IC86_I-XI

Then we initialize the analysis object as already shown in the Fitting a steady point-source with the public 14-year IceCube track data tutorial, but we fix the seed, minimum and maximum values defining the spectral index parameter. This effectively fixes the spectral index to the desired value.

Let’s say we want to fix the spectral index to 2.0. Then we need to pass gamma_seed, gamma_min, and gamma_max equal to 2.0 to the create_analysis method.

Additionally, one would also set the refplflux_gamma parameter to the same value of 2.0, as that parameter sets the power-law spectral index that is used for generating pseudo signal events.

from skyllh.analyses.i3.publicdata_ps.time_integrated_ps import create_analysis

ana = create_analysis(
    cfg=cfg,
    datasets=datasets,
    source=source,
    refplflux_gamma=2.0,
    gamma_seed=2.0,
    gamma_min=2.0,
    gamma_max=2.0,
)

100%|██████████| 33/33 [00:10<00:00,  3.12it/s]
100%|██████████| 33/33 [00:09<00:00,  3.61it/s]
100%|██████████| 33/33 [00:07<00:00,  4.67it/s]
100%|██████████| 33/33 [00:15<00:00,  2.14it/s]
100%|██████████| 4/4 [04:45<00:00, 71.42s/it]
100%|██████████| 136/136 [00:00<00:00, 2999.89it/s]

This now initialises the analysis instance with a fixed power law spectral index \(\gamma = 2.0\).

..todo: add a cell where signal events are injected and fit and then plot the distribution of TS, ns, and gamma. This shows that the outcome is what one expects.

..todo: show how the fit of NGC 1068 worsens if one fix the gamma to 2.0, because NGC 1068 is a soft source!

Consequently, 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 the energy_range property is updated, the flux normalisation will also be updated correctly. You just can’t change the flux parameters (gamma in this case) anymore.

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())

3.0837040854223075e-17
2.004341561536863e-15
1.041381138284864e-16