%load_ext autoreload
%autoreload 2

import numpy as np

from skyllh.core.config import Config

cfg = Config()

PDFs

The most basic building block of a likelihood function is a probability density function (PDF). SkyLLH provides the skyllh.core.pdf module, which defines base classes for PDFs. The abstract base class of all PDFs is PDF. It defines the abstract get_pd() method. In the derived PDF class this method needs to return the probability density values for each event. In case of a signal PDF it will have to return those values also for all the sources.

Time PDFs

SkyLLH provides the very generic time PDF abstract base class TimePDF, which takes a Livetime instance and a TimeFluxProfile instance to calculate the probability density values of the events.

For signal and background time PDFs the SignalTimePDF and BackgroundTimePDF classes exist, respectively.

As an example we are creating a signal time PDF with a gaussian time profile:

from skyllh.core.flux_model import GaussianTimeFluxProfile
from skyllh.core.livetime import Livetime
from skyllh.core.signalpdf import SignalTimePDF

livetime = Livetime(np.array([[0.0, 3.0], [6.0, 10.0]], dtype=np.float64))
gaussian_tfp = GaussianTimeFluxProfile(t0=7, sigma_t=0.2, cfg=cfg)
pdf = SignalTimePDF(pmm=None, livetime=livetime, time_flux_profile=gaussian_tfp, cfg=cfg)

The gaussian time flux profile is restricted to the following time range:

print(f'Time range: [{gaussian_tfp.t_start}, {gaussian_tfp.t_stop}]')

Time range: [5.513231124460065, 8.486768875539935]

In order to evaluate the time PDF, we need to create a TrialDataManager instance and the structured array holding the local source parameter names and values.

from skyllh.core.detsigyield import NullDetSigYieldBuilder
from skyllh.core.flux_model import NullFluxModel
from skyllh.core.parameters import ParameterModelMapper
from skyllh.core.source_hypo_grouping import (
    SourceHypoGroup,
    SourceHypoGroupManager,
)
from skyllh.core.source_model import SourceModel
from skyllh.core.storage import DataFieldRecordArray
from skyllh.core.trialdata import TrialDataManager

sources = [SourceModel()]

shg_mgr = SourceHypoGroupManager(
    SourceHypoGroup(
        sources=sources, fluxmodel=NullFluxModel(cfg=cfg), detsigyield_builders=NullDetSigYieldBuilder(cfg=cfg)
    )
)

pmm = ParameterModelMapper(models=sources)
src_params_recarray = pmm.create_src_params_recarray([])

events = DataFieldRecordArray({'time': np.array([0.0, 1.1, 5.6, 8.1, 9.5])})

tdm = TrialDataManager()
tdm.initialize_trial(shg_mgr=shg_mgr, pmm=pmm, events=events)

Via the get_pd() method we can now get the PDF values:

(pd, grads) = pdf.get_pd(tdm=tdm, params_recarray=src_params_recarray)
print(pd)

[0.00000000e+00 0.00000000e+00 0.00000000e+00 5.38488156e-07
 0.00000000e+00]

As can be seen, only the 4th event time (8.1) has a non-zero PDF value, because all other values are either outside the detector on-time or are not covered by the time flux profile.