Fitting a steady point-source with the public 14-year IceCube track data#

This tutorial shows how to fit a neutrino point-like steady source (time-integrated) using the IceCube public 14-year point-source data with SkyLLH.

We assume that you have installed the SkyLLH package in your environment by following the Installation guide.

[1]:

import numpy as np

Initializing a configuration dictionary#

First we have to create a config instance.

[2]:

from skyllh.core.config import Config

cfg = Config()

The Config object stores a dictionary with global settings, including computing, logging, and analysis settings. It is a required argument to several methods.

[3]:

cfg

[3]:

{'multiproc': {'ncpu': None},
 'logging': {'log_level': 'INFO',
  'log_format': '%(asctime)s %(processName)s %(name)s %(levelname)s: %(message)s',
  'enable_tracing': False},
 'project': {'working_directory': '.'},
 'repository': {'base_path': '/Users/tomas/.cache/skyllh',
  'download_from_origin': True},
 'units': {'internal': {'angle': Unit("rad"),
   'energy': Unit("GeV"),
   'length': Unit("cm"),
   'time': Unit("s")},
  'defaults': {'fluxes': {'angle': Unit("rad"),
    'energy': Unit("GeV"),
    'length': Unit("cm"),
    'time': Unit("s")}}},
 'datafields': {'run': 4,
  'ra': 4,
  'dec': 4,
  'ang_err': 4,
  'time': 4,
  'log_energy': 4,
  'true_ra': 8,
  'true_dec': 8,
  'true_energy': 8,
  'mcweight': 8},
 'caching': {'pdf': {'MultiDimGridPDF': False}}}

Multiprocessing#

The multiprocessing in skyllh can be enabled globally by updating the configuration with a number of CPU cores. It almost linearly speeds up creating an analysis and running trials, with a trade-off of a slightly higher memory usage:

[4]:

cfg.set_ncpu(6)

cfg

[4]:

{'multiproc': {'ncpu': 6},
 'logging': {'log_level': 'INFO',
  'log_format': '%(asctime)s %(processName)s %(name)s %(levelname)s: %(message)s',
  'enable_tracing': False},
 'project': {'working_directory': '.'},
 'repository': {'base_path': '/Users/tomas/.cache/skyllh',
  'download_from_origin': True},
 'units': {'internal': {'angle': Unit("rad"),
   'energy': Unit("GeV"),
   'length': Unit("cm"),
   'time': Unit("s")},
  'defaults': {'fluxes': {'angle': Unit("rad"),
    'energy': Unit("GeV"),
    'length': Unit("cm"),
    'time': Unit("s")}}},
 'datafields': {'run': 4,
  'ra': 4,
  'dec': 4,
  'ang_err': 4,
  'time': 4,
  'log_energy': 4,
  'true_ra': 8,
  'true_dec': 8,
  'true_energy': 8,
  'mcweight': 8},
 'caching': {'pdf': {'MultiDimGridPDF': False}}}

Setting up logging#

The user sets up the main logger. The logger name, logging format and level, output file, and other options can be set.

Default logging format and level are defined in the Config instance, but they can be overwritten by the arguments of logging.setup_logging.

[5]:

from skyllh.core.logging import setup_logging

logger = setup_logging(cfg=cfg, name='fitting_a_source', log_level='INFO')

Getting the datasets#

The datasets are automatically downloaded from dataverse.harvard.edu to a local cache directory (~/.cache/skyllh) on first run. To use a custom dataset location, set cfg['repository']['base_path'] to the desired path.

Note

IceTracks-DR2 is a pretty heavy dataset (~2GB) and it might take a while to initialize the analysis when it has to be downloaded for the first time. Be patient, and the next time it will be loaded from ~/.cache/skyllh :)

The public datasets are available through the top-level skyllh package:

[6]:

import skyllh

The collection of datasets is created with skyllh.create_datasets. With the default configuration, the data is downloaded automatically from dataverse.harvard.edu.

[7]:

help(skyllh.create_datasets)

Help on function create_datasets in module skyllh.datasets.datasets:

create_datasets(sample_name, cfg, names=None, base_path=None, sub_path_fmt=None)
    Creates a list of Dataset instances for a named data sample.

    Parameters
    ----------
    sample_name : str
        The name of the data sample. Available samples are the keys of
        ``skyllh.datasets.data_samples``.
    cfg : instance of Config
        The instance of Config holding the local configuration.
    names : sequence of str | None
        The dataset names to return. If None, the module's ``DATASET_NAMES``
        default is used (combined IC86 seasons where applicable).
    base_path : str | None
        The base path of the data files. If None, uses
        ``cfg['repository']['base_path']``.
    sub_path_fmt : str | None
        The sub path format override. If None, uses the module default.

    Returns
    -------
    datasets : list of Dataset

If no names are passed to skyllh.create_datasets, it loads all data sets that are part of the IceTracks-DR2 dataset collection.

A tutorial on how to access specific information about which datasets are included in, e.g., IceTracks-DR2 is provided in Dataset collections tutorial.

[8]:

datasets = skyllh.create_datasets('IceTracks-DR2', cfg=cfg)

2026-10-28 21:35:47,797 MainProcess skyllh.datasets.datasets INFO: Loaded 4 dataset(s) from sample "IceTracks-DR2": IC40, IC59, IC79, IC86_I-XI

By default, this imports all available data sets in IceTracks-DR2 (all 14 years of data).

Creating the analysis object#

The analysis for doing a point-source search as in, e.g., https://doi.org/10.1103/PhysRevLett.124.051103 is pre-defined. This analysis instance can be created via the create_analysis method of the time_integrated_ps module of the public data interface.

The function requires the config, the datasets, and a source object. It also takes a number of additional optional arguments for which we refer the reader to the docstring.

[9]:

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

help(create_analysis)

Help on function create_analysis in module skyllh.analyses.i3.publicdata_ps.time_integrated_ps:

create_analysis(
    cfg: Config,
    datasets: list[Dataset],
    source: PointLikeSource,
    refplflux_Phi0: float = 1,
    refplflux_E0: float = 1000.0,
    refplflux_gamma: float = 2.0,
    refplflux_Ec: float = inf,
    energy_range: tuple[float, float] | None = None,
    ns_seed: float = 10.0,
    ns_min: float = 0.0,
    ns_max: float = 1000.0,
    gamma_seed: float | None = 3.0,
    gamma_min: float = 1.0,
    gamma_max: float = 4.0,
    kde_smoothing: bool = False,
    minimizer_impl: str = 'LBFGS',
    minimizer_max_rep: int = 100,
    compress_data: bool = False,
    keep_data_fields: list[str] | None = None,
    evt_sel_delta_angle_deg: float = 10,
    construct_sig_generator: bool = True,
    tl: TimeLord | None = None,
    ppbar: ProgressBar | None = None,
    logger_name: str | None = None
) -> SingleSourceMultiDatasetLLHRatioAnalysis
    Creates the Analysis instance for this particular analysis.

    Parameters
    ----------
    cfg
        The instance of Config holding the local configuration.
    datasets
        The list of Dataset instances, which should be used in the
        analysis.
    source
        The PointLikeSource instance defining the point source position.
    refplflux_Phi0
        The flux normalization to use for the reference power law flux model.
    refplflux_E0
        The reference energy to use for the reference power law flux model.
    refplflux_gamma
        The spectral index to use for the reference power law flux model.
    refplflux_Ec
        The cutoff energy for the cutoff power law flux model.
    energy_range
        The energy range for signal generation. Both low and high energies are
        given in GeV. If set to ``None``, the entire energy range of the
        dataset is used.
    ns_seed
        Value to seed the minimizer with for the ns fit.
    ns_min
        Lower bound for ns fit.
    ns_max
        Upper bound for ns fit.
    gamma_seed
        Value to seed the minimizer with for the gamma fit. If set to None,
        the refplflux_gamma value will be set as gamma_seed.
    gamma_min
        Lower bound for gamma fit.
    gamma_max
        Upper bound for gamma fit.
    kde_smoothing
        Deprecated: use of ``kde_smoothing=True`` is deprecated and will be
        removed in a future version.
        Apply a KDE-based smoothing to the data-driven background pdf.
        Default
    minimizer_impl
        Minimizer implementation to be used. Supported options are ``"LBFGS"``
        (L-BFG-S minimizer used from the :mod:`scipy.optimize` module),
        or ``"minuit"`` (Minuit minimizer used by the :mod:`iminuit` module).
        Default: "LBFGS".
    minimizer_max_rep
        In case the minimization process did not converge at the first time
        this option specifies the maximum number of repetitions with
        different initials. Default is 100.
    compress_data
        Flag if the data should get converted from float64 into float32.
    keep_data_fields
        List of additional data field names that should get kept when loading
        the data.
    evt_sel_delta_angle_deg
        The delta angle in degrees for the event selection optimization methods.
    construct_sig_generator
        Flag if the signal generator should be constructed (``True``) or not
        (``False``).
    tl
        The TimeLord instance to use to time the creation of the analysis.
    ppbar
        The instance of ProgressBar for the optional parent progress bar.
    logger_name
        The name of the logger to be used. If set to ``None``, ``__name__`` will
        be used.

    Returns
    -------
    ana
        The Analysis instance for this analysis.

Defining the Source: Here we use the location of NGC 1068 as an example source.

Instantiating a PointLikeSource object requires the source location (R.A. and Dec.) in radians. Optionally, the source can have a name and a weight. The latter is in principle useful when performing a stacking analysis which is not supported for public data yet.

[10]:

from skyllh.core.source_model import PointLikeSource

help(PointLikeSource)

Help on class PointLikeSource in module skyllh.core.source_model:

class PointLikeSource(SourceModel, IsPointlike)
 |  PointLikeSource(
 |      ra: float,
 |      dec: float,
 |      name: str | None = None,
 |      weight: float | None = None,
 |      **kwargs
 |  )
 |
 |  The PointLikeSource class is a source model for a point-like source
 |  object in the sky at a given location (right-ascention and declination).
 |
 |  Method resolution order:
 |      PointLikeSource
 |      SourceModel
 |      skyllh.core.model.Model
 |      IsPointlike
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  __init__(
 |      self,
 |      ra: float,
 |      dec: float,
 |      name: str | None = None,
 |      weight: float | None = None,
 |      **kwargs
 |  )
 |      Creates a new PointLikeSource instance for defining a point-like
 |      source.
 |
 |      Parameters
 |      ----------
 |      ra
 |          The right-ascention coordinate of the source in radians.
 |      dec
 |          The declination coordinate of the source in radians.
 |      name
 |          The name of the source.
 |      weight
 |          The relative weight of the source w.r.t. other sources.
 |          If set to None, unity will be used.
 |
 |  __str__(self)
 |      Pretty string representation.
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from SourceModel:
 |
 |  classification
 |      The astronomical classification of the source.
 |
 |  weight
 |      The weight of the source.
 |
 |  ----------------------------------------------------------------------
 |  Readonly properties inherited from skyllh.core.model.Model:
 |
 |  id
 |      (read-only) The ID of the model. It's an integer generated with
 |      Python's `id` function. Hence, it's related to the memory address
 |      of the object.
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from skyllh.core.model.Model:
 |
 |  __dict__
 |      dictionary for instance variables
 |
 |  __weakref__
 |      list of weak references to the object
 |
 |  name
 |      The name of the model.
 |
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from IsPointlike:
 |
 |  dec
 |      The declination coordinate of the point-like source.
 |
 |  ra
 |      The right-ascention coordinate of the point-like source.


[11]:

src_ra = 40.67  # degrees
src_dec = -0.01  # degrees

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

If not otherwise specified through the optional arguments, the time_integrated_ps.create_analysis method assumes a signal hypothesis in which the defined source emits a powerlaw flux with spectral index \(\gamma=2.0\):

\[\phi(E_{\nu}) = \mathrm{refplflux\_Phi0} \cdot \left(\frac{E_{\nu}}{\mathrm{refplflux\_E0}}\right)^{-\mathrm{refplflux\_gamma}}\]
[12]:

ana = create_analysis(cfg=cfg, datasets=datasets, source=source)

2026-10-28 21:35:48,696 MainProcess skyllh.analyses.i3.publicdata_ps.time_integrated_ps INFO: SourceHypoGroupManager
    Source Hypothesis Groups:
        0: SourceHypoGroup:
            sources (1):
                0: PointLikeSource: "4659955424": { ra=40.670 deg, dec=-0.010 deg }
            fluxmodel:
                1.000e+00 * (E / (1000 GeV))^-2 * 1 (GeV cm^2 s)^-1
            detector signal yield builders (1):
                PDSingleParamFluxPointLikeSourceI3DetSigYieldBuilder
            signal generation method:
                NoneType
2026-10-28 21:35:48,697 MainProcess skyllh.analyses.i3.publicdata_ps.time_integrated_ps INFO: ParameterModelMapper: 2 global parameters, 2 models (1 source)
    Parameters:
        ns [floating (0 <= 10 <= 1000)]
            in models:
            - IceCube: ns

        gamma [floating (1 <= 3 <= 4)]
            in models:
            - 4659955424: gamma


100%|██████████| 33/33 [00:00<00:00, 480.94it/s]
100%|██████████| 33/33 [00:00<00:00, 487.45it/s]
100%|██████████| 33/33 [00:00<00:00, 370.73it/s]
100%|██████████| 33/33 [00:00<00:00, 480.19it/s]
100%|██████████| 4/4 [00:24<00:00,  6.12s/it]
100%|██████████| 136/136 [00:00<00:00, 836.92it/s]

Understanding how the analysis works#

One usually wants to get best-fit values (maximum likelihood estimators) for the source parameters and a test-statistic, using methods like ana.do_trial or ana.unblind (see below). However, here we describe what happens under the hood when the analysis is run either on simulated or real experimental data.

Initializing a trial#

After the Analysis instance was created, the point-source analysis can be run. To do so the analysis is initialized with some trial data. For instance we could initialize the analysis with the experimental data but the same applies for simulated data.

The Analysis object provides the method initialize_trial to initialize a trial with data. It takes a list of DataFieldRecordArray instances holding the events. If we want to initialize a trial with the experimental data, we can get that list from the Analysis instance itself:

[13]:

events_list = [data.exp for data in ana.data_list]
ana.initialize_trial(events_list)

Maximizing the log-likelihood ratio function#

After initializing a trial, we can maximize the LLH ratio function using the ana.llhratio.maximize method. This method requires a RandomStateService instance in case the minimizer does not succeed and a new set of initial values for the fit parameters need to get generated. The method returns a 3-element tuple. The first element is the value of the LLH ratio function at its maximum. The second element is the set of fit parameters that maximize the LLH ratio. The third element is the status dictionary of the minimizer.

[14]:

from skyllh.core.random import RandomStateService

rss = RandomStateService(seed=1)

[15]:

(log_lambda_max, fitparam_values, status) = ana.llhratio.maximize(rss)

[16]:

print(f'log_lambda_max = {log_lambda_max}')
print(f'fitparam_values = {fitparam_values}')
print(f'status = {status}')

log_lambda_max = 14.383579243318888
fitparam_values = [79.54729206  3.20155111]
status = {'grad': array([-3.08035275e-06, -3.17262427e-04]), 'task': 'CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH', 'funcalls': 13, 'nit': 10, 'warnflag': 0, 'skyllh_minimizer_n_reps': 0, 'n_llhratio_func_calls': 13}

Calculating the test-statistic#

Using the maximum of the LLH ratio function and the fit parameter values at the maximum we can calculate the test-statistic using the ana.calculate_test_statistic method:

[17]:

TS = ana.calculate_test_statistic(log_lambda_max, fitparam_values)
print(f'TS = {TS:.3f}')

TS = 28.767

Unblinding the data#

Above, we perform the analysis following pedagogical steps: we initialize the analysis with a trial of the experimental data, maximize the log-likelihood ratio function for all given experimental data events, and calculate the test-statistic value.

However, the same calculations can be performed all together using the ana.unblind method of the analysis instance, which returns the test statistic and the best fit parameters directly.

Note that the ana.unblind method runs the analysis on the experimental data.

[18]:

help(ana.unblind)

Help on method unblind in module skyllh.core.analysis:

unblind(minimizer_rss: RandomStateService, tl: TimeLord | None = None) -> tuple[float, dict, dict] method of skyllh.core.analysis.SingleSourceMultiDatasetLLHRatioAnalysis instance
    Evaluates the unscrambled data, i.e. unblinds the data.

    Parameters
    ----------
    minimizer_rss
        The instance of RandomStateService that should be used by the
        minimizer to generate new random initial fit parameter values.
    tl
        The optional instance of TimeLord that should be used to time the
        maximization of the LLH ratio function.

    Returns
    -------
    TS
        The test-statistic value.
    global_params_dict
        The dictionary holding the global parameter names and their
        best fit values. It includes fixed and floating parameters.
    status
        The status dictionary with information about the performed
        minimization process of the negative of the log-likelihood ratio
        function.


[19]:

from skyllh.core.random import RandomStateService

rss = RandomStateService(seed=1)

ts, fitparam_values, status = ana.unblind(rss)

[20]:

print(f'TS = {ts:.3f}')
print(f'ns = {fitparam_values["ns"]:.2f}')
print(f'gamma = {fitparam_values["gamma"]:.2f}')
print(f'status = {status}')

TS = 28.767
ns = 79.55
gamma = 3.20
status = {'grad': array([-3.08035275e-06, -3.17262427e-04]), 'task': 'CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH', 'funcalls': 13, 'nit': 10, 'warnflag': 0, 'skyllh_minimizer_n_reps': 0, 'n_llhratio_func_calls': 13}

Converting the fit parameters into a flux normalization#

By default the analysis assuming a powerlaw signal spectrum is created with a flux normalization of 1 \(\text{GeV}^{-1}\text{cm}^{-2}\text{s}^{-1}\) at 1 TeV. Note that the differential flux of a point source has no solid-angle term, hence no \(\text{sr}^{-1}\) in its unit (see refplflux_Phi0 and refplflux_E0 arguments of the create_analysis method).

The analysis instance offers the method ana.calculate_fluxmodel_scaling_factor that calculates the scaling factor that, for a given spectral model, converts a mean number of signal events of 1 into a flux normalization refplflux_Phi0 at refplflux_E0.

[21]:

# The default value of `create_analysis()`
refplflux_gamma = 2.0

scaling_factor = ana.calculate_fluxmodel_scaling_factor()
print(f'1 event = {scaling_factor:.3e} (E/1000 GeV)^{{-{refplflux_gamma:.1f}}} 1/(GeV s cm^2)')

1 event = 3.083e-17 (E/1000 GeV)^{-2.0} 1/(GeV s cm^2)

By default, the conversion factor is returned for the spectral index with which the analysis was created, i.e. refplflux_gamma (2.0 by default).

The method optionally takes a fitparam_values argument where a specific spectral index can be specified. Note that this won’t change the internal refplflux_gamma value with which the analysis was initialized.

[22]:

refplflux_gamma = fitparam_values['gamma']

scaling_factor = ana.calculate_fluxmodel_scaling_factor(fitparam_values=[1, refplflux_gamma])
print(
    f'best-fit flux =  {fitparam_values["ns"] * scaling_factor:.3e} (E/1000 GeV)^{{-{refplflux_gamma:.1f}}} 1/(GeV s cm^2)'
)

best-fit flux =  2.845e-14 (E/1000 GeV)^{-3.2} 1/(GeV s cm^2)

Alternatively, the ana.mu2flux and ana.flux2mu methods return:

  • a flux normalization at refplflux_E0 given a certain mean number of signal events passed to it as an argument

  • a mean number of signal events corresponding to a flux normalization at refplflux_E0 passed to it as an argument

respectively.

They also take different spectral parameter values via the optional fitparam_values argument.

[23]:

phi0 = ana.mu2flux(10)
print(f'Flux normalization for 10 events = {phi0:.3e} 1/(GeV s cm^2)')

Flux normalization for 10 events = 3.083e-16 1/(GeV s cm^2)

[24]:

mean_ns = ana.flux2mu(5e-15)
print(f'Mean number of events for a flux normalization at 1000 GeV of 5e-15 1/(GeV s cm^2) = {mean_ns:.2f}')

Mean number of events for a flux normalization at 1000 GeV of 5e-15 1/(GeV s cm^2) = 162.18