API Reference

analysis

Analysis cache objects for csky.

This module provides caching of key objects required for virtually all analyses. These include background space parameterizations, signal acceptance parameterizations, and energy PDF ratio models.

Classes:

`Analysis`(repo, specs[, load_sig])	Analysis cache for all datasets (i.e. seasons) used in an analysis.
`SubAnalysis`(repo, spec[, src, mask_plane, ...])	Analysis cache for a single dataset (i.e. "season").

class csky.analysis.Analysis(repo, specs, load_sig=True, **kw)[source]

Analysis cache for all datasets (i.e. seasons) used in an analysis.

This class manages an ensemble of SubAnalysis instances so they can be treated as a unit when producing a csky.trial.TrialRunner in order to perform analysis tasks.

These should not be constructed directly; instead, see csky.conf.get_analysis().

Methods:

__init__(repo, specs[, load_sig])

__init__(repo, specs, load_sig=True, **kw)[source]

class csky.analysis.SubAnalysis(repo, spec, src=[], mask_plane=False, gp_width=0.08726646259971647, mask_radius=0.0, src_tomask=[], dir='', gammas=array([1., 1.0625, 1.125, 1.1875, 1.25, 1.3125, 1.375, 1.4375, 1.5, 1.5625, 1.625, 1.6875, 1.75, 1.8125, 1.875, 1.9375, 2., 2.0625, 2.125, 2.1875, 2.25, 2.3125, 2.375, 2.4375, 2.5, 2.5625, 2.625, 2.6875, 2.75, 2.8125, 2.875, 2.9375, 3., 3.0625, 3.125, 3.1875, 3.25, 3.3125, 3.375, 3.4375, 3.5, 3.5625, 3.625, 3.6875, 3.75, 3.8125, 3.875, 3.9375, 4. ]), energy_kw={}, acc_kw={}, space_bg_kw={}, min_sigma=0.003490658503988659, compress=True, load_sig=True, mmap=True, energy_pdf_ratio_model_cls=<class 'csky.pdf.EnergyPDFRatioModel'>, acc_param_cls=<class 'csky.pdf.SinDecAccParameterization'>, bg_space_param_cls=<class 'csky.pdf.BgSinDecParameterization'>)[source]

Analysis cache for a single dataset (i.e. “season”).

This class holds references to the background space parameterizations, signal acceptance parameterizations, and energy PDF ratio models for a single dataset. These can be saved to and loaded from disk. Parameterizations for custom signal fluxes can also be requested on-demand, in which case these will be cached during execution but not included in the on-disk representation.

These should not be constructed directly; instead, see csky.conf.get_analysis().

See csky.selections for pre-defined csky.selections.DataSpec options.

Methods:

__init__(repo, spec[, src, mask_plane, ...])

Construct a SubAnalysis.

__init__(repo, spec, src=[], mask_plane=False, gp_width=0.08726646259971647, mask_radius=0.0, src_tomask=[], dir='', gammas=array([1., 1.0625, 1.125, 1.1875, 1.25, 1.3125, 1.375, 1.4375, 1.5, 1.5625, 1.625, 1.6875, 1.75, 1.8125, 1.875, 1.9375, 2., 2.0625, 2.125, 2.1875, 2.25, 2.3125, 2.375, 2.4375, 2.5, 2.5625, 2.625, 2.6875, 2.75, 2.8125, 2.875, 2.9375, 3., 3.0625, 3.125, 3.1875, 3.25, 3.3125, 3.375, 3.4375, 3.5, 3.5625, 3.625, 3.6875, 3.75, 3.8125, 3.875, 3.9375, 4. ]), energy_kw={}, acc_kw={}, space_bg_kw={}, min_sigma=0.003490658503988659, compress=True, load_sig=True, mmap=True, energy_pdf_ratio_model_cls=<class 'csky.pdf.EnergyPDFRatioModel'>, acc_param_cls=<class 'csky.pdf.SinDecAccParameterization'>, bg_space_param_cls=<class 'csky.pdf.BgSinDecParameterization'>)[source]

Construct a SubAnalysis.

Parameters:

repo (csky.selections.Repository) – the repository for loading the data
spec (csky.selections.DataSpec) – the dataset specification
src (list of csky.utils.Sources) – if given, sources whose [t_start,t_stop] time windows must be excluded from background parameterizations
dir (str) – directory from which to load cached PDF parameterizations
gammas (array) – gammas to parameterize in energy PDFs and signal acceptance
energy_kw (mapping) – extra kwargs for the energy_pdf_ratio_model
acc_kw (mapping) – extra kwargs for the acc_param
space_bg_kw (mapping) – extra kwargs for the bg_space_param
min_sigma (float) – minimum sigma (angular uncertainty estimate) to allow
compress (bool) – whether to compress the dataset (drop unneeded columns and use float32 where possible)
load_sig (bool) – whether to load signal MC (set to False to save RAM in background trial jobs)
mmap (bool) – whether to use memmap disk reading (should reduce peak RAM usage)
energy_pdf_ratio_model_cls (type) – class to use for energy_pdf_ratio_model
acc_param_cls (type) – class to use for acc_param
bg_space_param_cls (type) – class to use for bg_space_param_cls
gp_width (float) – size of galactic plane gp_width to mask in radians.

bk

Bookkeeping helper functions.

Functions:

`get_all`(root_dir, pattern[, merge, ...])	Get a tree of nested dictionaries with files loaded as leaf items.
`get_best`(d, *values)	Descend into dict-tree d, obtaining the best-match subtree or item.

csky.bk.get_all(root_dir, pattern, merge=<function sum>, pre_convert=None, post_convert=None, log=True)[source]

Get a tree of nested dictionaries with files loaded as leaf items.

Parameters:

root_dir (str) – the root of the directory tree
pattern (str) – a glob pattern for the files to load
merge (func(list) -> object) – function that merges loaded objects

csky.bk.get_best(d, *values)[source]: Descend into dict-tree d, obtaining the best-match subtree or item.

conf

Configuration tools for csky.

Functions:

`describe`(o[, visited, d, path])	Describe a csky object and any more csky objects inside it.
`do_bind`(f, conf[, is_class])	Partially bind `f` by applying configuration `conf`.
`get_analysis`(repo, args, *kw)	Get an `csky.analysis.Analysis` instance.
`get_injs`(a, llh_model, conf[, do_inj, llh_conf])	Get `csky.inj.Injector` instances for TRUTH, background, and signal
`get_llh`(a, conf, **kw)	Get a `csky.llh.LLHModel`.
`get_multiflare_trial_runner`([conf])	Get a `csky.trial.MultiflareTrialRunner`.
`get_obj`(cls, conf[, subana, is_class])	Get object of type `cls`.
`get_sky_scan_trial_runner`([conf, inj_conf, ...])	Get a `csky.trial.SkyScanTrialRunner`.
`get_spatial_prior_trial_runner`([conf, ...])	Get a `csky.trial.SpatialPriorTrialRunner`.
`get_trial_runner`([conf, inj_conf])	Get a `csky.trial.TrialRunner`.
`overlay`(conf[, new_conf])	Override elements of `conf` with other arguments.

csky.conf.describe(o, visited=None, d=0, path='')[source]: Describe a csky object and any more csky objects inside it.

csky.conf.do_bind(f, conf, is_class=True)[source]: Partially bind f by applying configuration conf.

csky.conf.get_analysis(repo, *args, **kw)[source]

Get an csky.analysis.Analysis instance.

Parameters:

repo (csky.selections.Repository) – the repository for loading the data
*args – one or more data specifications, with each one optionally preceeded by a string indicating the desired dataset version (see csky.selections)
**kw – other keyword arguments passed to csky.analysis.Analysis constructor

csky.conf.get_injs(a, llh_model, conf, do_inj=True, llh_conf={}, **kw)[source]: Get csky.inj.Injector instances for TRUTH, background, and signal

csky.conf.get_llh(a, conf, **kw)[source]: Get a csky.llh.LLHModel.

csky.conf.get_multiflare_trial_runner(conf={}, **kw)[source]: Get a csky.trial.MultiflareTrialRunner.

csky.conf.get_obj(cls, conf, subana=None, is_class=True, **kw)[source]: Get object of type cls.

csky.conf.get_sky_scan_trial_runner(conf={}, inj_conf={}, multiflare=False, src_tr=False, **kw)[source]: Get a csky.trial.SkyScanTrialRunner.

csky.conf.get_spatial_prior_trial_runner(conf={}, inj_conf={}, multiflare=False, src_tr=False, **kw)[source]: Get a csky.trial.SpatialPriorTrialRunner.

csky.conf.get_trial_runner(conf={}, inj_conf={}, **kw)[source]: Get a csky.trial.TrialRunner.

csky.conf.overlay(conf, new_conf={})[source]: Override elements of conf with other arguments.

coord

Coordinate transformations.

Functions:

`delta_angle`(zenith1, azimuth1, zenith2, azimuth2)	Return the space angle between `fit1` and `fit2`.
`rotate_source_to_xaxis`(source_zenith, ...[, ...])	Rotate coordinates such that points at the source location move to the x axis.
`rotate_source_to_zaxis`(source_zenith, ...[, ...])	Rotate coordinates such that points at the source location move to the z axis.
`rotate_xaxis_to_source`(source_zenith, ...[, ...])	Rotate coordinates such that points on the xaxis move to a source location.
`rotate_xzplane_to_source`(axis_zenith, ...[, ...])	Rotate coordinates such that points in the x-z plane move to a source location.
`rotate_zaxis_to_source`(source_zenith, ...[, ...])	Rotate coordinates such that points on the zaxis move to a source location.
`switch_ra_azimuth`(phi_in, mjd)	Givin MJD, transform azimuth->RA or RA->azimuth.

csky.coord.delta_angle(zenith1, azimuth1, zenith2, azimuth2, cos=False, latlon=False)[source]

Return the space angle between fit1 and fit2.

Parameters:

zenith1 (float or ndarray) – zenith 1 or declination 1
azimuth1 (float or ndarray) – azimuth 1
zenith2 (float or ndarray) – zenith 2 or declination 2
azimuth2 (float or ndarray) – azimuth 2
cos (bool) – if True, return cos(delta_angle); otherwise return delta_angle
latlon (bool) – if True, zenith1 and zenith2 are interpreted as declinations rather than polar zeniths

Returns:

of delta_angle or cos(delta_angle) values.

Return type:

numpy.ndarray

This function seems redundant but is useful in its treatment of nontrivial numpy broadcasting use cases.

csky.coord.rotate_source_to_xaxis(source_zenith, source_azimuth, zenith, azimuth, latlon=False)[source]

Rotate coordinates such that points at the source location move to the x axis.

Parameters:

source_zenith (float or ndarray) – the destination source zenith
source_azimuth (float or ndarray) – the destination source azimuth
zenith (float or ndarray) – the zenith of the map point(s)
azimuth (float or ndarray) – the azimuth of the map point(s)
latlon (bool) – if True, source_zenith and zenith are interpreted as declinations rather than polar zeniths

Input coordinates are assumed to be standard spherical (either DC or CC).

csky.coord.rotate_source_to_zaxis(source_zenith, source_azimuth, zenith, azimuth, latlon=False)[source]

Rotate coordinates such that points at the source location move to the z axis.

Parameters:

source_zenith (float or ndarray) – the destination source zenith
source_azimuth (float or ndarray) – the destination source azimuth
zenith (float or ndarray) – the zenith of the map point(s)
azimuth (float or ndarray) – the azimuth of the map point(s)
latlon (bool) – if True, source_zenith and zenith are interpreted as declinations rather than polar zeniths

Input coordinates are assumed to be standard spherical (either DC or CC).

csky.coord.rotate_xaxis_to_source(source_zenith, source_azimuth, zenith, azimuth, latlon=False)[source]

Rotate coordinates such that points on the xaxis move to a source location.

Parameters:

source_zenith (float or ndarray) – the destination source zenith
source_azimuth (float or ndarray) – the destination source azimuth
zenith (float or ndarray) – the zenith of the map point(s)
azimuth (float or ndarray) – the azimuth of the map point(s)
latlong (bool) – if True, source_zenith and zenith are interpreted as declinations rather than polar zeniths

Input coordinates are assumed to be standard spherical(either DC or CC).

csky.coord.rotate_xzplane_to_source(axis_zenith, source_zenith, source_azimuth, zenith, azimuth, reverse=False, latlon=False)[source]

Rotate coordinates such that points in the x-z plane move to a source location.

Parameters:

axis_zenith (float) – the destination source zenith
source_zenith (float or ndarray) – the destination source zenith
source_azimuth (float or ndarray) – the destination source azimuth
zenith (float or ndarray) – the zenith of the map point(s)
azimuth (float or ndarray) – the azimuth of the map point(s)
reverse (bool) – if True, perform the inverse coordinate rotation
latlon (bool) – if True, source_zenith and zenith are interpreted as declinations rather than polar zeniths

Input coordinates are assumed to be standard spherical(in DC) or equatorial.

csky.coord.rotate_zaxis_to_source(source_zenith, source_azimuth, zenith, azimuth, latlon=True)[source]

Rotate coordinates such that points on the zaxis move to a source location.

Parameters:

source_zenith (float or ndarray) – the destination source zenith
source_azimuth (float or ndarray) – the destination source azimuth
zenith (float or ndarray) – the zenith of the map point(s)
azimuth (float or ndarray) – the azimuth of the map point(s)
latlon (bool) – if True, source_zenith and zenith are interpreted as declinations rather than polar zeniths

Input coordinates are assumed to be standard spherical (either DC or CC).

csky.coord.switch_ra_azimuth(phi_in, mjd)[source]: Givin MJD, transform azimuth->RA or RA->azimuth.

dists

Test statistic distributions models.

Classes:

`BinnedTSD`(values[, n_zero])	Histogram-based representation of a test statistic distribution.
`Chi2TSD`(values[, n_zero, threshold])	\(\chi^2\)-fit representation of a test statistic distribution.
`TSD`(values[, n_zero])	Trials-based representation of a test statistic distribution.

Functions:

ts_to_p(bg, key, ts[, fit])

Convert TS to p-values, interpolating over the background grid.

class csky.dists.BinnedTSD(values, n_zero=0, **kw)[source]

Bases: TSD

Histogram-based representation of a test statistic distribution.

This class tries to reproduce the behavior of TSD without storing information about every single trial, thereby reducing RAM and disk space requirements.

Methods:

`__init__`(values[, n_zero])	Construct a Chi2TSD.
`cdf`(x[, fit])	The cumulative distribution function, similar to `scipy.stats` distributions.
`get_hist`(**kw)	Get a `histlite.Hist` representation of the distribution.
`isf`(p[, fit])	The inverse survival function, similar to `scipy.stats` distributions.
`median`()	The median TS value.
`sf`(x[, fit])	The survival function, similar to `scipy.stats` distributions.
`sf_nsigma`(x, **kw)	The survival function, but returning a value in terms of number of sigmas.

__init__(values, n_zero=0, **kw)[source]

Construct a Chi2TSD.

Parameters:

values (array of float or utis.Arrays) – the test statistic values
n_zero (int) – number of trials not contained by values for which TS=0
kw (mapping) – additional arguments for scipy.stats.chi2.fit. By default, if loc and floc are not set, then floc=0 will be used; if scale and fscale are not set, then fscale=1 will be used.

cdf(x, fit=True)[source]: The cumulative distribution function, similar to scipy.stats distributions.

get_hist(**kw): Get a histlite.Hist representation of the distribution.

isf(p, fit=True)[source]: The inverse survival function, similar to scipy.stats distributions.

median()[source]: The median TS value.

sf(x, fit=True)[source]: The survival function, similar to scipy.stats distributions.

sf_nsigma(x, **kw): The survival function, but returning a value in terms of number of sigmas.

class csky.dists.Chi2TSD(values, n_zero=0, threshold=0, **kw)[source]

Bases: TSD

\(\chi^2\)-fit representation of a test statistic distribution.

Information about trials is retained, and a chi2 fit is applied to the region TS > 0.

Methods:

`__init__`(values[, n_zero, threshold])	Construct a Chi2TSD.
`cdf`(x[, fit])	The cumulative distribution function, similar to `scipy.stats` distributions.
`get_hist`(**kw)	Get a `histlite.Hist` representation of the distribution.
`isf`(p[, fit])	The inverse survival function, similar to `scipy.stats` distributions.
`median`([fit])	The median TS value.
`sf`(x[, fit])	The survival function, similar to `scipy.stats` distributions.
`sf_nsigma`(x, **kw)	The survival function, but returning a value in terms of number of sigmas.

Attributes:

chi2

The fitted scipy.stats.chi2 object.

__init__(values, n_zero=0, threshold=0, **kw)[source]

Construct a Chi2TSD.

Parameters:

values (array of float or utis.Arrays) – the test statistic values
n_zero (int) – number of trials not contained by values for which TS=0
kw (mapping) – additional arguments for scipy.stats.chi2.fit. By default, if loc and floc are not set, then floc=0 will be used; if scale and fscale are not set, then fscale=1 will be used.

cdf(x, fit=True)[source]: The cumulative distribution function, similar to scipy.stats distributions.

property chi2: The fitted scipy.stats.chi2 object.

get_hist(**kw): Get a histlite.Hist representation of the distribution.

isf(p, fit=True)[source]: The inverse survival function, similar to scipy.stats distributions.

median(fit=False)[source]: The median TS value.

sf(x, fit=True)[source]: The survival function, similar to scipy.stats distributions.

sf_nsigma(x, **kw): The survival function, but returning a value in terms of number of sigmas.

class csky.dists.TSD(values, n_zero=0, **kw)[source]

Bases: object

Trials-based representation of a test statistic distribution.

Several methods are in the spirit of scipy.stats distribution objects, but note that at least currently, the interface is not a perfect match.

Methods:

`__init__`(values[, n_zero])	Construct a Chi2TSD.
`cdf`(x)	The cumulative distribution function, similar to `scipy.stats` distributions.
`get_hist`(**kw)	Get a `histlite.Hist` representation of the distribution.
`isf`(p[, fit])	The inverse survival function, similar to `scipy.stats` distributions.
`median`()	The median TS value.
`sf`(x[, fit])	The survival function, similar to `scipy.stats` distributions.
`sf_nsigma`(x, **kw)	The survival function, but returning a value in terms of number of sigmas.

__init__(values, n_zero=0, **kw)[source]

Construct a Chi2TSD.

Parameters:

values (array of float or utis.Arrays) – the test statistic values
n_zero (int) – number of trials not contained by values for which TS=0
kw (mapping) – additional arguments for scipy.stats.chi2.fit. By default, if loc and floc are not set, then floc=0 will be used; if scale and fscale are not set, then fscale=1 will be used.

cdf(x)[source]: The cumulative distribution function, similar to scipy.stats distributions.

get_hist(**kw)[source]: Get a histlite.Hist representation of the distribution.

isf(p, fit=False)[source]: The inverse survival function, similar to scipy.stats distributions.

median()[source]: The median TS value.

sf(x, fit=False)[source]: The survival function, similar to scipy.stats distributions.

sf_nsigma(x, **kw)[source]: The survival function, but returning a value in terms of number of sigmas.

csky.dists.ts_to_p(bg, key, ts, fit=True)[source]

Convert TS to p-values, interpolating over the background grid.

Parameters:

bg (dict) – the backgound distributions
key (array of float) – the values that indicate where along the background parameterization each TS value originates from
fit (bool) – whether to request sf(ts_value, fit=True)

hyp

Spectral hypothesis models.

Classes:

`BinnedFlux`(bins_energy, flux[, ...])	Energy-binned spectrum.
`Flux`()	Abstract base class for fluxes.
`PowerLawFlux`(gamma[, norm, energy_range, ...])	Power law spectrum.
`SeyfertCoreCoronaFlux`(psp_table, ...[, ...])	Implements the Core-Corona Seyfert Galaxy neutrino flux model of A.
`SplineFluxModel`(norm, psp_table, ...)	Spline flux spectrum.

class csky.hyp.BinnedFlux(bins_energy, flux, max_energy_approx_gamma=1000000.0)[source]

Bases: Flux

Energy-binned spectrum.

Methods:

`__call__`(E)	Get dN/dE given E.
`__init__`(bins_energy, flux[, ...])	Construct a BinnedFlux.
`to_E2dNdE`(ns, acc_total[, E0, unit])	Get E0^2 * dN/dE for particular reference energy and units.
`to_dNdE`(ns, acc_total[, E0, unit])	Get dN/dE for particular reference energy and units.
`to_model_norm`(ns, acc_total)	Get the model normalization resulting in ns, given the total acceptance.
`to_ns`(flux, acc_total[, E0, unit, use_E2dNdE])	Get number of signal events corresponding to a particular flux.

Attributes:

`as_params`	Convert to str->value fit parameter kwargs.
`energy_range`	Get the energy range containing nonzero flux, (Emin,Emax), in GeV.

__call__(E)[source]

Get dN/dE given E.

Parameters:: E (float) – energy in GeV
Returns:: dN/dE in units 1/GeV/cm2/s
Return type:: float

__init__(bins_energy, flux, max_energy_approx_gamma=1000000.0)[source]

Construct a BinnedFlux.

Parameters:

bins_energy (ndarray of float) – energy bin edges(up to or including Emax of highest energy bin)
flux (ndarray of float) – average dN/dE, in 1/GeV/cm2/s, for each energy bin
max_energy_approx_gamma (float) – maximum energy(in GeV) to consider when computing approximate average spectral index(used by :meth:BinnedFlux.as_params)

property as_params: Convert to str->value fit parameter kwargs.

property energy_range: Get the energy range containing nonzero flux, (Emin,Emax), in GeV.

to_E2dNdE(ns, acc_total, E0=1, unit=1)

Get E0^2 * dN/dE for particular reference energy and units.

Just like Flux.to_dNdE(), except multiplied by E0^2.

to_dNdE(ns, acc_total, E0=1, unit=1)

Get dN/dE for particular reference energy and units.

Parameters:

ns (float) – number of signal events
acc_total (float) – total signal acceptance
E0 (float) – reference energy in units unit * GeV
unit (float) – energy unit for E0 as well as in denominator of return value

Returns:

dN/dE in units 1/(unit * GeV)/cm2/s

Return type:

float

For example, to get dN/dE at 100 TeV in units 1/TeV/cm2/s, use to_dNdE(ns, acc_total, E0=100, unit=1000).

Todo

Check units of acc_total

to_model_norm(ns, acc_total)

Get the model normalization resulting in ns, given the total acceptance.

Parameters:

ns (float) – number of signal events
acc_total (float) – total signal acceptance

Returns:

model normalization

Return type:

float

Todo

Check units of acc_total

This is just a semantically meaningful way of expressing ns / acc_total.

to_ns(flux, acc_total, E0=1, unit=1, use_E2dNdE=True)

Get number of signal events corresponding to a particular flux.

Parameters:

flux (float) – dN/dE or E0^2 * dN/dE(see argument E2dNdE)
acc_total (float) – total signal acceptance
E0 (float) – reference energy in units unit * GeV
unit (float) – energy unit for E0 as well as in denominator of return value
E2dNdE (bool) – if true(default), flux gives E0^2 * dN/dE. otherwise, flux gives dN/dE.

Returns:

number of signal events

Return type:

float

Todo

Check units of acc_total

class csky.hyp.Flux[source]

Bases: object

Abstract base class for fluxes.

This class provides basic operations for flux functions of the form dN/dE=f(E), where E is in GeV and dN/dE is in units 1/GeV/cm2/s. The main value-added from this class is counts <–> flux calculations, which are performed in terms of total signal acceptance.

Todo

Check and then document the exact units of “total signal acceptance”.

Methods:

`__call__`(E)	Get dN/dE given E.
`__init__`()
`to_E2dNdE`(ns, acc_total[, E0, unit])	Get E0^2 * dN/dE for particular reference energy and units.
`to_dNdE`(ns, acc_total[, E0, unit])	Get dN/dE for particular reference energy and units.
`to_model_norm`(ns, acc_total)	Get the model normalization resulting in ns, given the total acceptance.
`to_ns`(flux, acc_total[, E0, unit, use_E2dNdE])	Get number of signal events corresponding to a particular flux.

Attributes:

`as_params`	Convert to str->value fit parameter kwargs.
`energy_range`	Get the energy range containing nonzero flux, (Emin,Emax), in GeV.

abstract __call__(E)[source]

Get dN/dE given E.

Parameters:: E (float) – energy in GeV
Returns:: dN/dE in units 1/GeV/cm2/s
Return type:: float

__init__()

abstract property as_params: Convert to str->value fit parameter kwargs.

abstract property energy_range: Get the energy range containing nonzero flux, (Emin,Emax), in GeV.

to_E2dNdE(ns, acc_total, E0=1, unit=1)[source]

Get E0^2 * dN/dE for particular reference energy and units.

Just like Flux.to_dNdE(), except multiplied by E0^2.

to_dNdE(ns, acc_total, E0=1, unit=1)[source]

Get dN/dE for particular reference energy and units.

Parameters:

ns (float) – number of signal events
acc_total (float) – total signal acceptance
E0 (float) – reference energy in units unit * GeV
unit (float) – energy unit for E0 as well as in denominator of return value

Returns:

dN/dE in units 1/(unit * GeV)/cm2/s

Return type:

float

For example, to get dN/dE at 100 TeV in units 1/TeV/cm2/s, use to_dNdE(ns, acc_total, E0=100, unit=1000).

Todo

Check units of acc_total

to_model_norm(ns, acc_total)[source]

Get the model normalization resulting in ns, given the total acceptance.

Parameters:

ns (float) – number of signal events
acc_total (float) – total signal acceptance

Returns:

model normalization

Return type:

float

Todo

Check units of acc_total

This is just a semantically meaningful way of expressing ns / acc_total.

to_ns(flux, acc_total, E0=1, unit=1, use_E2dNdE=True)[source]

Get number of signal events corresponding to a particular flux.

Parameters:

flux (float) – dN/dE or E0^2 * dN/dE(see argument E2dNdE)
acc_total (float) – total signal acceptance
E0 (float) – reference energy in units unit * GeV
unit (float) – energy unit for E0 as well as in denominator of return value
E2dNdE (bool) – if true(default), flux gives E0^2 * dN/dE. otherwise, flux gives dN/dE.

Returns:

number of signal events

Return type:

float

Todo

Check units of acc_total

class csky.hyp.PowerLawFlux(gamma, norm=1, energy_range=(0, inf), energy_cutoff=inf)[source]

Bases: Flux

Power law spectrum.

Methods:

`__call__`(E)	Get dN/dE given E.
`__init__`(gamma[, norm, energy_range, ...])	param gamma: spectral index
`to_E2dNdE`(ns, acc_total[, E0, unit])	Get E0^2 * dN/dE for particular reference energy and units.
`to_dNdE`(ns, acc_total[, E0, unit])	Get dN/dE for particular reference energy and units.
`to_model_norm`(ns, acc_total)	Get the model normalization resulting in ns, given the total acceptance.
`to_ns`(flux, acc_total[, E0, unit, use_E2dNdE])	Get number of signal events corresponding to a particular flux.

Attributes:

`as_params`	Convert to str->value fit parameter kwargs.
`energy_range`	Get the energy range containing nonzero flux, (Emin,Emax), in GeV.

__call__(E)[source]

Get dN/dE given E.

Parameters:: E (float) – energy in GeV
Returns:: dN/dE in units 1/GeV/cm2/s
Return type:: float

__init__(gamma, norm=1, energy_range=(0, inf), energy_cutoff=inf)[source]

Parameters:

gamma (float) – spectral index
norm (float) – normalization factor(at 1 GeV)
energy_range (tuple of (float, float)) – hard energy bounds(flux is zero elsewhere)
energy_cutoff (float) – exponential cutoff energy(in GeV)

property as_params: Convert to str->value fit parameter kwargs.

property energy_range: Get the energy range containing nonzero flux, (Emin,Emax), in GeV.

to_E2dNdE(ns, acc_total, E0=1, unit=1)

Get E0^2 * dN/dE for particular reference energy and units.

Just like Flux.to_dNdE(), except multiplied by E0^2.

to_dNdE(ns, acc_total, E0=1, unit=1)

Get dN/dE for particular reference energy and units.

Parameters:

ns (float) – number of signal events
acc_total (float) – total signal acceptance
E0 (float) – reference energy in units unit * GeV
unit (float) – energy unit for E0 as well as in denominator of return value

Returns:

dN/dE in units 1/(unit * GeV)/cm2/s

Return type:

float

For example, to get dN/dE at 100 TeV in units 1/TeV/cm2/s, use to_dNdE(ns, acc_total, E0=100, unit=1000).

Todo

Check units of acc_total

to_model_norm(ns, acc_total)

Get the model normalization resulting in ns, given the total acceptance.

Parameters:

ns (float) – number of signal events
acc_total (float) – total signal acceptance

Returns:

model normalization

Return type:

float

Todo

Check units of acc_total

This is just a semantically meaningful way of expressing ns / acc_total.

to_ns(flux, acc_total, E0=1, unit=1, use_E2dNdE=True)

Get number of signal events corresponding to a particular flux.

Parameters:

flux (float) – dN/dE or E0^2 * dN/dE(see argument E2dNdE)
acc_total (float) – total signal acceptance
E0 (float) – reference energy in units unit * GeV
unit (float) – energy unit for E0 as well as in denominator of return value
E2dNdE (bool) – if true(default), flux gives E0^2 * dN/dE. otherwise, flux gives dN/dE.

Returns:

number of signal events

Return type:

float

Todo

Check units of acc_total

class csky.hyp.SeyfertCoreCoronaFlux(psp_table, log_xray_lumin, src_dist, norm, lumin_scale=1.0, crit_log_energy_flux=-50, crit_log_nu_energy_lower=2.0, crit_log_nu_energy_upper=7.0)[source]

Bases: SplineFluxModel

Implements the Core-Corona Seyfert Galaxy neutrino flux model of A. Kheirandish et al., Astrophys.J. 922 (2021) 45 by means of B-spline interpolation.

Methods:

`__call__`(E)	The flux value dN/dE at energy E.
`__init__`(psp_table, log_xray_lumin, ...[, ...])
`to_E2dNdE`(ns, acc_total[, E0, unit])	Get E0^2 * dN/dE for particular reference energy and units.
`to_dNdE`(ns, acc_total[, E0, unit])	Get dN/dE for particular reference energy and units.
`to_model_norm`(ns, acc_total)	Get the model normalization resulting in ns, given the total acceptance.
`to_ns`(flux, acc_total[, E0, unit, use_E2dNdE])	Get number of signal events corresponding to a particular flux.

Attributes:

`as_params`	Convert to str->value fit parameter kwargs.
`crit_log_energy_flux`	defines when the flux is considered to be 0
`crit_log_nu_energy_lower`	The lower bound of the support of the spline interpolator
`crit_log_nu_energy_upper`	The upper bound of the support of the spline interpolator
`energy_range`	Get the energy range containing nonzero flux, (Emin,Emax), in GeV.
`log_xray_lumin`	The log10 of the intrinsic source luminosity in 2-10keV x-ray band
`lumin_scale`	correct normalization of model flux by relative factor
`math_function_str`	We have used hash of this representation in skyllh.core.source_hypothesis.get_fluxmodel_to_source_map() to map fluxes with KDE PDFs.
`norm`	The relative flux normalization.
`psp_table`	The photospline.SplineTable object that describes the neutrino flux as function of neutrino energy via B-spline interpolation.
`src_dist`	The distance to the source in units of Mpc

__call__(E)[source]

The flux value dN/dE at energy E. :param E: Evaluation energy [GeV] :type E: float | 1d numpy.ndarray of float

Returns:: flux – Flux at energy E in units of GeV^-1 cm^-2 s^-1.
Return type:: float | 1d ndarray of float

__init__(psp_table, log_xray_lumin, src_dist, norm, lumin_scale=1.0, crit_log_energy_flux=-50, crit_log_nu_energy_lower=2.0, crit_log_nu_energy_upper=7.0)[source]

property as_params: Convert to str->value fit parameter kwargs.

property crit_log_energy_flux: defines when the flux is considered to be 0

property crit_log_nu_energy_lower: The lower bound of the support of the spline interpolator

property crit_log_nu_energy_upper: The upper bound of the support of the spline interpolator

property energy_range: Get the energy range containing nonzero flux, (Emin,Emax), in GeV.

property log_xray_lumin: The log10 of the intrinsic source luminosity in 2-10keV x-ray band

property lumin_scale: correct normalization of model flux by relative factor

property math_function_str: We have used hash of this representation in skyllh.core.source_hypothesis.get_fluxmodel_to_source_map() to map fluxes with KDE PDFs. So far KDEs only depend on log_xray_lumin value and not on src_dist.

property norm: The relative flux normalization. norm=1 corresponds to the nominal model flux.

property psp_table: The photospline.SplineTable object that describes the neutrino flux as function of neutrino energy via B-spline interpolation.

property src_dist: The distance to the source in units of Mpc

to_E2dNdE(ns, acc_total, E0=1, unit=1)

Get E0^2 * dN/dE for particular reference energy and units.

Just like Flux.to_dNdE(), except multiplied by E0^2.

to_dNdE(ns, acc_total, E0=1, unit=1)

Get dN/dE for particular reference energy and units.

Parameters:

ns (float) – number of signal events
acc_total (float) – total signal acceptance
E0 (float) – reference energy in units unit * GeV
unit (float) – energy unit for E0 as well as in denominator of return value

Returns:

dN/dE in units 1/(unit * GeV)/cm2/s

Return type:

float

For example, to get dN/dE at 100 TeV in units 1/TeV/cm2/s, use to_dNdE(ns, acc_total, E0=100, unit=1000).

Todo

Check units of acc_total

to_model_norm(ns, acc_total)

Get the model normalization resulting in ns, given the total acceptance.

Parameters:

ns (float) – number of signal events
acc_total (float) – total signal acceptance

Returns:

model normalization

Return type:

float

Todo

Check units of acc_total

This is just a semantically meaningful way of expressing ns / acc_total.

to_ns(flux, acc_total, E0=1, unit=1, use_E2dNdE=True)

Get number of signal events corresponding to a particular flux.

Parameters:

flux (float) – dN/dE or E0^2 * dN/dE(see argument E2dNdE)
acc_total (float) – total signal acceptance
E0 (float) – reference energy in units unit * GeV
unit (float) – energy unit for E0 as well as in denominator of return value
E2dNdE (bool) – if true(default), flux gives E0^2 * dN/dE. otherwise, flux gives dN/dE.

Returns:

number of signal events

Return type:

float

Todo

Check units of acc_total

class csky.hyp.SplineFluxModel(norm, psp_table, crit_log_nu_energy_lower, crit_log_nu_energy_upper)[source]

Bases: Flux

Spline flux spectrum.

Methods:

`__call__`(E)	Get dN/dE given E.
`__init__`(norm, psp_table, ...)
`to_E2dNdE`(ns, acc_total[, E0, unit])	Get E0^2 * dN/dE for particular reference energy and units.
`to_dNdE`(ns, acc_total[, E0, unit])	Get dN/dE for particular reference energy and units.
`to_model_norm`(ns, acc_total)	Get the model normalization resulting in ns, given the total acceptance.
`to_ns`(flux, acc_total[, E0, unit, use_E2dNdE])	Get number of signal events corresponding to a particular flux.

Attributes:

`as_params`	Convert to str->value fit parameter kwargs.
`crit_log_nu_energy_lower`	The lower bound of the support of the spline interpolator
`crit_log_nu_energy_upper`	The upper bound of the support of the spline interpolator
`energy_range`	Get the energy range containing nonzero flux, (Emin,Emax), in GeV.
`norm`	The relative flux normalization.
`psp_table`	The photospline.SplineTable object that describes the neutrino flux as function of neutrino energy via B-spline interpolation.

abstract __call__(E)

Get dN/dE given E.

Parameters:: E (float) – energy in GeV
Returns:: dN/dE in units 1/GeV/cm2/s
Return type:: float

__init__(norm, psp_table, crit_log_nu_energy_lower, crit_log_nu_energy_upper)[source]

abstract property as_params: Convert to str->value fit parameter kwargs.

property crit_log_nu_energy_lower: The lower bound of the support of the spline interpolator

property crit_log_nu_energy_upper: The upper bound of the support of the spline interpolator

abstract property energy_range: Get the energy range containing nonzero flux, (Emin,Emax), in GeV.

property norm: The relative flux normalization. norm=1 corresponds to the nominal model flux.

property psp_table: The photospline.SplineTable object that describes the neutrino flux as function of neutrino energy via B-spline interpolation.

to_E2dNdE(ns, acc_total, E0=1, unit=1)

Get E0^2 * dN/dE for particular reference energy and units.

Just like Flux.to_dNdE(), except multiplied by E0^2.

to_dNdE(ns, acc_total, E0=1, unit=1)

Get dN/dE for particular reference energy and units.

Parameters:

ns (float) – number of signal events
acc_total (float) – total signal acceptance
E0 (float) – reference energy in units unit * GeV
unit (float) – energy unit for E0 as well as in denominator of return value

Returns:

dN/dE in units 1/(unit * GeV)/cm2/s

Return type:

float

For example, to get dN/dE at 100 TeV in units 1/TeV/cm2/s, use to_dNdE(ns, acc_total, E0=100, unit=1000).

Todo

Check units of acc_total

to_model_norm(ns, acc_total)

Get the model normalization resulting in ns, given the total acceptance.

Parameters:

ns (float) – number of signal events
acc_total (float) – total signal acceptance

Returns:

model normalization

Return type:

float

Todo

Check units of acc_total

This is just a semantically meaningful way of expressing ns / acc_total.

to_ns(flux, acc_total, E0=1, unit=1, use_E2dNdE=True)

Get number of signal events corresponding to a particular flux.

Parameters:

flux (float) – dN/dE or E0^2 * dN/dE(see argument E2dNdE)
acc_total (float) – total signal acceptance
E0 (float) – reference energy in units unit * GeV
unit (float) – energy unit for E0 as well as in denominator of return value
E2dNdE (bool) – if true(default), flux gives E0^2 * dN/dE. otherwise, flux gives dN/dE.

Returns:

number of signal events

Return type:

float

Todo

Check units of acc_total

inj

Event injectors and randomization methods.

Classes:

`DataInjector`(ana, data, keep[, randomizers, ...])	Inject real data events.
`DataOnOffInjector`(ana, time_model, keep[, ...])	Inject off-time events into on-time windows.
`DecBandSelector`(src[, cut_n_sigma])	Select events based on proximity in declination to source(s).
`DecRandomizer`([dec_rand_method, ...])	Randomize events in declination.
`EnergyRandomizer`(sigma_log10)	Randomize event energies.
`Injector`()	Base class for event injectors.
`MCBackgroundInjector`(ana, bg_weight_names, keep)	Inject background events using weighted MC events to model the background instead of scrambling data in RA.
`MJDGRLRandomizer`(grl[, presample, rates_by])	Randomize event times according to GRL.
`MJDShuffleRandomizer`()	Shuffle event times.
`PointSourceBinnedTimeInjector`(ana, src, ...)	Inject events distributed according to a binned lightcurve.
`PointSourceInjector`(ana, src, flux, keep[, ...])	Inject one or more point-like sources.
`PointSourceRedNoiseInjector`(ana, src, flux, ...)	Inject red noise as signal.
`PointSourceTimeDepInjector`(ana, src, flux, keep)	Base class for for time-dependent injection from point-like sources.
`PointSourceTimeWindowInjector`(ana, src, ...)	Inject Gaussian or time box flares.
`PoleRandomizer`(pole_radius[, shuffle, ...])	Randomize declinations of events near the North and South poles.
`RARandomizer`()	Randomize events uniformly in right ascension.
`Randomizer`()	Base class for event randomizers.
`Selector`()	Base class for event selectors.
`SigInjector`()	Base class for signal injectors.
`SinDecRandomizer`(sindec_jitter)	Randomize events in sin declination.
`SpatialPriorInjector`(ana, src, flux, keep, ...)	Injection of one or more point-like sources, whose position is not exactly known, but given in the form of a spatial prior (healpix map).
`TemplateBGInjector`(ana, template_model, ...)	Inject typical types of background along with a highly-extended source according to a spatial template.
`TemplateInjector`(ana, template_model, keep)	Inject a highly-extended source according to a spatial template.
`TimeWindowSelector`(mjd_min, mjd_max)	Select events within a given time window.
`TransientInjector`(ana, time_model, flux, keep)	Inject events from transient sources.

class csky.inj.DataInjector(ana, data, keep, randomizers=None, bg_replace=False)[source]

Bases: Injector

Inject real data events.

This class injects real events from the data, possibly(in fact, almost always) with one or more randomizations applied.

Methods:

`__init__`(ana, data, keep[, randomizers, ...])	Construct a DataInjector.
`inject`([seed])	Get some number of events.

__init__(ana, data, keep, randomizers=None, bg_replace=False)[source]

Construct a DataInjector.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
data (utils.Events) – the data events to inject
keep (list of str) – names of event features to keep
randomizers (list of Randomizer) – the randomizers to apply; if not given, RARandomizer is used.
bg_replace (bool) – If true, sampling with replacement will be used.

inject(seed=None)[source]

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

class csky.inj.DataOnOffInjector(ana, time_model, keep, randomizers=None, cut_n_sigma=5, full_sky=False, bg_replace=True)[source]

Bases: Injector

Inject off-time events into on-time windows.

This class performs the sort of background injection typically used in transient analyses, where the background is taken from off-time data while the analysis is performed in specific on-time windows.

Methods:

`__init__`(ana, time_model, keep[, ...])	Construct a DataOnOffInjector.
`inject`([seed])	Get some number of events.

__init__(ana, time_model, keep, randomizers=None, cut_n_sigma=5, full_sky=False, bg_replace=True)[source]

Construct a DataOnOffInjector.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
time_model (pdf.TimePDFRatioModel) – the time PDF ratio model
keep (list of str) – names of event features to keep
randomizers (list of Randomizer) – the randomizers to apply. you probably don’t want to apply any(TODO: reconsider, why exactly is this argument here?)
cut_n_sigma (float) – number of sigmas away from a source an event must be within to be included
bg_replace (bool) – If true, sampling with replacement will be used.

inject(seed=None)[source]

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

class csky.inj.DecBandSelector(src, cut_n_sigma=5)[source]

Bases: Selector

Select events based on proximity in declination to source(s).

This class allows pre-emptive application of the so-called “band” cut(see PointSourceSpacePDFRatioModel for more).

Methods:

`__call__`(arr)	Get a downselected copy of `arr`.
`__init__`(src[, cut_n_sigma])	Construct a DecBandSelector.
`mask`(arr)	Get a boolean mask for selected elements of `arr`.
`where`(arr)	Get the indices of selected elements of `arr`.

__call__(arr)[source]

Get a downselected copy of arr.

Parameters:: arr (utils.Arrays) – the Arrays to mask
Return type:: same as arr
Returns:: downselected copy(or potentially view) of arr

__init__(src, cut_n_sigma=5)[source]

Construct a DecBandSelector.

Parameters:

src (utis.Sources) – the source list. required keys: ra and dec. optional keys: weight, extension
cut_n_sigma (float) – number of sigmas away from a source an event must be within to be included

mask(arr)[source]

Get a boolean mask for selected elements of arr.

Parameters:: arr (utils.Arrays) – the Arrays to mask
Return type:: ndarray of bool
Returns:: the mask

where(arr)[source]

Get the indices of selected elements of arr.

Parameters:: arr (utils.Arrays) – the Arrays to mask
Return type:: ndarray of bool
Returns:: the indices

class csky.inj.DecRandomizer(dec_rand_method='gaussian_sigma', dec_rand_kwargs=None, dec_rand_pole_exlusion=None)[source]

Bases: Randomizer

Randomize events in declination.

This class applies randomization to the event declinations. Multiple randomization methods are available, which need to be specified via the dec_rand_method keyword. Additionally, events near the poles can be excluded from the randomization by setting the pole exclusion radius via the keyword dec_rand_pole_exlusion.

Methods:

`__call__`(ev[, seed])	Randomize events
`__init__`([dec_rand_method, dec_rand_kwargs, ...])	Construct a DecRandomizer

__call__(ev, seed=None)[source]: Randomize events

__init__(dec_rand_method='gaussian_sigma', dec_rand_kwargs=None, dec_rand_pole_exlusion=None)[source]

Construct a DecRandomizer

Parameters:

dec_rand_method (str) –
The randomization method to use. ‘gaussian_sigma’:

Applies Gaussian-distributed noise to the event declinations. The width of the Gaussians is based on the per-event angular uncertainties sigma. Additional kwargs (to be passed via dec_rand_kwargs):

’sigma_scaling’:
Scaling factor that will be applied to the per-event angular uncertainties prior to randomization. This only affects the randomization - event angular uncertainties will not be modified.

’uniform_sigma’:
Applies uniform jitter to the event declinations. The width of the jitter is based on the per-event angular uncertainties sigma. Additional kwargs (to be passed via dec_rand_kwargs):

’sigma_scaling’:
Scaling factor that will be applied to the per-event angular uncertainties prior to randomization. This only affects the randomization - event angular uncertainties will not be modified.

’gaussian_fixed’:
Applies Gaussian-distributed noise to the event declinations. The width of the Gaussians is set to randomization_width. Additional kwargs (to be passed via dec_rand_kwargs):

’randomization_width’:
The sigma of the Gaussian.

’uniform_fixed’:
Applies uniform jitter to the event declinations. The width of the Gaussians is set to randomization_width. Additional kwargs (to be passed via dec_rand_kwargs):

’randomization_width’:
The width for the jitter. Events will be randomized to dec +- width.
dec_rand_kwargs (dict) – Additional keyword arguments provided as a dictionary. See dec_rand_method for details on which kwargs are necessary for the specified randomization method.
dec_rand_pole_exlusion (float) – Radius of the polar regions(in radians), which will be ignored. Events with declinations in this pole region will not be randomized. If `None`(default), then all events are randomized.

class csky.inj.EnergyRandomizer(sigma_log10)[source]

Bases: Randomizer

Randomize event energies.

This class randomizes log10(E/GeV). This may be useful if some events are causing intense spikes in the background distributions at particular declinations. Note that this feature has not been thoroughly tested, if at all.

Methods:

`__call__`(ev[, seed])	Randomize events in place.
`__init__`(sigma_log10)	Construct an EnergyRandomizer.

__call__(ev, seed=None)[source]

Randomize events in place.

Parameters:

ev (utils.Events) – the events
seed (int or numpy.random.RandomState) – seed spec

__init__(sigma_log10)[source]

Construct an EnergyRandomizer.

Parameters:: sigma_log10 (float) – width of the Gaussian noise to apply to log10(E/GeV)

class csky.inj.Injector[source]

Bases: object

Base class for event injectors.

Injectors produce event ensembles along with counts of events that could have been injected but were instead pre-emptively excluded using a Selector.

TODO: revisit hierarchy. Background injectors don’t seem to want an n_inj argument at all, in the Injector.inject() method.

Methods:

`__init__`()
`inject`(n_inj[, seed])	Get some number of events.

__init__()

abstract inject(n_inj, seed=None)[source]

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

class csky.inj.MCBackgroundInjector(ana, bg_weight_names, keep, randomizers=None)[source]

Bases: Injector

Inject background events using weighted MC events to model the background instead of scrambling data in RA.

This class injects simulated events from the monte carlo, possibly(in fact, almost always) with one or

more randomizations applied.

Methods:

`__init__`(ana, bg_weight_names, keep[, ...])	Construct a PointSourceInjector.
`inject`([replace, seed])	Get some number of events.

__init__(ana, bg_weight_names, keep, randomizers=None)[source]

Construct a PointSourceInjector.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
mc (utils.Events) – the monte carlo events to inject
bg_weight_names (list of str) – names of the monte carlo weights saved in the signal events. Weights should be in Hz and already have flux incorperated.
keep (list of str) – names of event features to keep
randomizers (list of Randomizer) – the randomizers to apply to the monte carlo events

inject(replace=False, seed=None)[source]

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

class csky.inj.MJDGRLRandomizer(grl, presample=0, rates_by='events')[source]

Bases: Randomizer

Randomize event times according to GRL.

Methods:

`__call__`(ev[, seed])	Randomize events in place.
`__init__`(grl[, presample, rates_by])

__call__(ev, seed=None)[source]

Randomize events in place.

Parameters:

ev (utils.Events) – the events
seed (int or numpy.random.RandomState) – seed spec

__init__(grl, presample=0, rates_by='events')[source]

class csky.inj.MJDShuffleRandomizer[source]

Bases: Randomizer

Shuffle event times.

This class shuffles the arrival times(in MJD) of the events.

Methods:

`__call__`(ev[, seed])	Randomize events in place.
`__init__`()

__call__(ev, seed=None)[source]

Randomize events in place.

Parameters:

ev (utils.Events) – the events
seed (int or numpy.random.RandomState) – seed spec

__init__()

class csky.inj.PointSourceBinnedTimeInjector(ana, src, flux, keep, lcs, threshs=None, lags=None, sindec_bandwidth=0.03490658503988659)[source]

Bases: PointSourceTimeDepInjector

Inject events distributed according to a binned lightcurve.

Methods:

`__init__`(ana, src, flux, keep, lcs[, ...])	Construct a PointSourceBinnedTimeInjector.
`get_frac_during_livetime`()	Get the fraction of the total lightcurve covered by a given dataset.
`get_time_model`()	Get the time PDF ratio model.
`inject`(n_inj[, seed])	Get some number of events.
`sample_mjd`(src_n_injs[, seed])	Get per-source `mjd` arrays.

Attributes:

acc_total

Total signal acceptance.

__init__(ana, src, flux, keep, lcs, threshs=None, lags=None, sindec_bandwidth=0.03490658503988659)[source]

Construct a PointSourceBinnedTimeInjector.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
src (utis.Sources) – the source list. required keys: ra and dec. optional keys: weight, extension
flux (hyp.Flux) – the signal spectrum
keep (list of str) – names of event features to keep
lcs (list of histlite.Hist) – the per-source lightcurves(see pdf.BinnedTimePDFRatioModel)
threshs (ndarray of float) – per-source thresholds
sindec_bandwidth (float) – width(not half-width) of the sin(dec) bands from which to draw simulated signal events

property acc_total

Total signal acceptance.

This is similar to LLHEvaluatorBase.get_acc_total(), but evaluated in light of all specified signal injection details and typically in terms of the actual simulated signal events that will be used.

get_frac_during_livetime()[source]

Get the fraction of the total lightcurve covered by a given dataset.

See pdf.TimePDFRatioModel.get_frac_during_livetime().

get_time_model()[source]

Get the time PDF ratio model.

Return type:: pdf.TimePDFRatioModel

inject(n_inj, seed=None)

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

sample_mjd(src_n_injs, seed=None)[source]

Get per-source mjd arrays.

Parameters:

src_n_injs (ndarray of int) – per-source number of events to inject
seed (int or numpy.random.RandomState) – seed spec

class csky.inj.PointSourceInjector(ana, src, flux, keep, sindec_bandwidth=0.03490658503988659)[source]

Bases: SigInjector

Inject one or more point-like sources.

This class implements the signal injection for classic time-integrated analyses. It transparently supports stacking of multiple sources as well as per-source Gaussian spatial extension.

Methods:

`__init__`(ana, src, flux, keep[, ...])	Construct a PointSourceInjector.
`inject`(n_inj[, seed])	Get some number of events.

Attributes:

acc_total

Total signal acceptance.

__init__(ana, src, flux, keep, sindec_bandwidth=0.03490658503988659)[source]

Construct a PointSourceInjector.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
src (utis.Sources) – the source list. required keys: ra and dec. optional keys: weight, extension
flux (hyp.Flux) – the signal spectrum
keep (list of str) – names of event features to keep
sindec_bandwidth (float) – width(not half-width) of the sin(dec) bands from which to draw simulated signal events

property acc_total

Total signal acceptance.

This is similar to LLHEvaluatorBase.get_acc_total(), but evaluated in light of all specified signal injection details and typically in terms of the actual simulated signal events that will be used.

inject(n_inj, seed=None)[source]

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

class csky.inj.PointSourceRedNoiseInjector(ana, src, flux, keep, dt, beta=2, sindec_bandwidth=0.03490658503988659)[source]

Bases: PointSourceTimeDepInjector

Inject red noise as signal.

This class implements white/pink/red/black noise injector as a signal light-curve for one source.

Methods:

`__init__`(ana, src, flux, keep, dt[, beta, ...])	Construct a PointSourceRedNoiseInjector.
`generate_rednoise`(N, dt, beta, random[, ...])	Generate a power-law lightcurve This uses the method from Timmer & Koenig [1]
`get_frac_during_livetime`()	Get the fraction of the total lightcurve covered by a given dataset.
`get_time_model`()	Get the time PDF ratio model.
`inject`(n_inj[, seed])	Get some number of events.
`sample_mjd`(src_n_injs[, seed])	Get per-source `mjd` arrays.

Attributes:

acc_total

Total signal acceptance.

__init__(ana, src, flux, keep, dt, beta=2, sindec_bandwidth=0.03490658503988659)[source]

Construct a PointSourceRedNoiseInjector.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
src (utils.Sources) – the source list. required keys: ra and dec. optional keys: weight, extension
flux (hyp.Flux) – the signal spectrum
keep (list of str) – names of event features to keep
dt (float (single number)) – spacing between time-steps/noise-curve in units of days (recommended < 50 days)
beta (integer (single number)) – power-law index of noise (default=2 for red noise)
sindec_bandwidth (float) – width(not half-width) of the sin(dec) bands from which to draw simulated signal events

property acc_total

Total signal acceptance.

This is similar to LLHEvaluatorBase.get_acc_total(), but evaluated in light of all specified signal injection details and typically in terms of the actual simulated signal events that will be used.

generate_rednoise(N, dt, beta, random, generate_complex=False)[source]

Generate a power-law lightcurve This uses the method from Timmer & Koenig [1]

:type N : integer :param N : number of equal-spaced time steps to generate :type dt : float :param dt : spacing between time-steps :type beta : float :param beta : power-law index, the spectrum will be (1 / f)^beta

white noise is beta=0 pink noise is beta=1 red noise is beta=2 black noise is beta>2

:type generate_complex : boolean (optional) :param generate_complex : if True, generate a complex time series rather

than a real time series

References

[1] Timmer, J. & Koenig, M. On Generating Power Law Noise. A&A 300:707

get_frac_during_livetime()[source]

Get the fraction of the total lightcurve covered by a given dataset.

See pdf.TimePDFRatioModel.get_frac_during_livetime().

get_time_model()[source]

Get the time PDF ratio model.

Return type:: pdf.TimePDFRatioModel

inject(n_inj, seed=None)

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

sample_mjd(src_n_injs, seed=None)[source]

Get per-source mjd arrays.

Parameters:

src_n_injs (ndarray of int) – per-source number of events to inject
seed (int or numpy.random.RandomState) – seed spec

class csky.inj.PointSourceTimeDepInjector(ana, src, flux, keep, sindec_bandwidth=0.03490658503988659, spatial_prior=False)[source]

Bases: SigInjector

Base class for for time-dependent injection from point-like sources.

This is an abstract base class for time-dependent injection from point-like sources. Most of the injection work is done in terms of per-source instances of PointSourceInjector; time-dependence is handled by derived classes’ implementations of the abstract methods PointSourceTimeDepInjector.get_time_model(), PointSourceTimeDepInjector.get_frac_during_livetime(), and PointSourceTimeDepInjector.sample_mjd().

Methods:

`__init__`(ana, src, flux, keep[, ...])	Construct a PointSourceTimeDepInjector.
`get_frac_during_livetime`()	Get the fraction of the total lightcurve covered by a given dataset.
`get_time_model`()	Get the time PDF ratio model.
`inject`(n_inj[, seed])	Get some number of events.
`sample_mjd`(src_n_injs[, seed])	Get per-source `mjd` arrays.

Attributes:

acc_total

Total signal acceptance.

__init__(ana, src, flux, keep, sindec_bandwidth=0.03490658503988659, spatial_prior=False)[source]

Construct a PointSourceTimeDepInjector.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
src (utis.Sources) – the source list. required keys: ra and dec. optional keys: weight, extension
flux (hyp.Flux) – the signal spectrum
keep (list of str) – names of event features to keep
sindec_bandwidth (float) – width(not half-width) of the sin(dec) bands from which to draw simulated signal events
spatial_prior (bool) – option to use a spatial prior healpix map for source localisations

property acc_total

Total signal acceptance.

This is similar to LLHEvaluatorBase.get_acc_total(), but evaluated in light of all specified signal injection details and typically in terms of the actual simulated signal events that will be used.

abstract get_frac_during_livetime()[source]

Get the fraction of the total lightcurve covered by a given dataset.

See pdf.TimePDFRatioModel.get_frac_during_livetime().

abstract get_time_model()[source]

Get the time PDF ratio model.

Return type:: pdf.TimePDFRatioModel

inject(n_inj, seed=None)[source]

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

abstract sample_mjd(src_n_injs, seed=None)[source]

Get per-source mjd arrays.

Parameters:

src_n_injs (ndarray of int) – per-source number of events to inject
seed (int or numpy.random.RandomState) – seed spec

class csky.inj.PointSourceTimeWindowInjector(ana, src, flux, keep, t0, dt, box=False, box_mode='center', endpoints=False, sindec_bandwidth=0.03490658503988659)[source]

Bases: PointSourceTimeDepInjector

Inject Gaussian or time box flares.

This class implements Gaussian or top-hat time window flare injection for one or more(possibly spatially-extended) sources.

Methods:

`__init__`(ana, src, flux, keep, t0, dt[, ...])	Construct a PointSourceTimeInjector.
`get_frac_during_livetime`()	Get the fraction of the total lightcurve covered by a given dataset.
`get_time_model`()	Get the time PDF ratio model.
`inject`(n_inj[, seed])	Get some number of events.
`sample_mjd`(src_n_injs[, seed])	Get per-source `mjd` arrays.

Attributes:

acc_total

Total signal acceptance.

__init__(ana, src, flux, keep, t0, dt, box=False, box_mode='center', endpoints=False, sindec_bandwidth=0.03490658503988659)[source]

Construct a PointSourceTimeInjector.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
src (utis.Sources) – the source list. required keys: ra and dec. optional keys: weight, extension
flux (hyp.Flux) – the signal spectrum
keep (list of str) – names of event features to keep
t0 (ndarray of float) – per-source anchor time-of-flare
dt (ndarray of float) – per-source width of flare(in days)
box (bool) – whether to use box rather than Gaussian flare
box_mode (str) – one of ‘center’, ‘pre’, or ‘post’
endpoints (bool) – injects and flare start and stop times, rather than throughout flare
sindec_bandwidth (float) – width(not half-width) of the sin(dec) bands from which to draw simulated signal events

property acc_total

Total signal acceptance.

This is similar to LLHEvaluatorBase.get_acc_total(), but evaluated in light of all specified signal injection details and typically in terms of the actual simulated signal events that will be used.

get_frac_during_livetime()[source]

Get the fraction of the total lightcurve covered by a given dataset.

See pdf.TimePDFRatioModel.get_frac_during_livetime().

get_time_model()[source]

Get the time PDF ratio model.

Return type:: pdf.TimePDFRatioModel

inject(n_inj, seed=None)

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

sample_mjd(src_n_injs, seed=None)[source]

Get per-source mjd arrays.

Parameters:

src_n_injs (ndarray of int) – per-source number of events to inject
seed (int or numpy.random.RandomState) – seed spec

class csky.inj.PoleRandomizer(pole_radius, shuffle=True, sample_dec=False)[source]

Bases: Randomizer

Randomize declinations of events near the North and South poles.

This class randomizes the declinations of events near the poles. This is mainly useful for cascade analysis, where it is non-obvious how much of the polar regions should simply be neglected. It may also be useful in track analyses if the user is specifically interested in sources near one or both of the poles.

Methods:

`__call__`(ev[, seed])	Randomize events in place.
`__init__`(pole_radius[, shuffle, sample_dec])	Construct a PoleRandomizer.

__call__(ev, seed=None)[source]

Randomize events in place.

Parameters:

ev (utils.Events) – the events
seed (int or numpy.random.RandomState) – seed spec

__init__(pole_radius, shuffle=True, sample_dec=False)[source]

Construct a PoleRandomizer.

Parameters:

pole_radius (float) – radius of the polar regions(in radians)
shuffle (bool) – if false, scramble polar events uniformly throughout the poles in which they each actually originated; if true(default), then instead randomize polar events’ declinations by applying a random permutation of these. if sample_dec is set to true, then the declination values are additionally smeared by a Gaussian (see description for sample_dec).
sample_dec (bool) – if true, the event declinations at the poles are additionally randomized by sampling from a Gaussian with the width set to the event sigma. Note that this setting may push events from the poles away towards the equator and therefore modify the joint PDFs! If false(default), only the permutation of event declination across events at the pole is performed. This setting is only relevant if shuffle is set to True.

class csky.inj.RARandomizer[source]

Bases: Randomizer

Randomize events uniformly in right ascension.

Methods:

`__call__`(ev[, seed])	Randomize events in place.
`__init__`()

__call__(ev, seed=None)[source]

Randomize events in place.

Parameters:

ev (utils.Events) – the events
seed (int or numpy.random.RandomState) – seed spec

__init__()

class csky.inj.Randomizer[source]

Bases: object

Base class for event randomizers.

Classes which randomize(or scramble) data events should derive from this class.

Methods:

`__call__`(ev[, seed])	Randomize events in place.
`__init__`()

abstract __call__(ev, seed=None)[source]

Randomize events in place.

Parameters:

ev (utils.Events) – the events
seed (int or numpy.random.RandomState) – seed spec

__init__()[source]

class csky.inj.Selector[source]

Bases: object

Base class for event selectors.

This is an interface for pre-emptively selecting events prior to [any randomizations as well as] likelihood evaluation. For per-trial masking, see PDFRatioEvaluator.

TODO: actually go and add documentation for masking functionality under PDFRatioEvaluator

Methods:

`__call__`(arr)	Get a downselected copy of `arr`.
`__init__`()
`mask`(arr)	Get a boolean mask for selected elements of `arr`.
`where`(arr)	Get the indices of selected elements of `arr`.

abstract __call__(arr)[source]

Get a downselected copy of arr.

Parameters:: arr (utils.Arrays) – the Arrays to mask
Return type:: same as arr
Returns:: downselected copy(or potentially view) of arr

__init__()

abstract mask(arr)[source]

Get a boolean mask for selected elements of arr.

Parameters:: arr (utils.Arrays) – the Arrays to mask
Return type:: ndarray of bool
Returns:: the mask

abstract where(arr)[source]

Get the indices of selected elements of arr.

Parameters:: arr (utils.Arrays) – the Arrays to mask
Return type:: ndarray of bool
Returns:: the indices

class csky.inj.SigInjector[source]

Bases: Injector

Base class for signal injectors.

Signal injectors, unlike background injectors, must also be able to provide their total acceptance. This class enforces that requirement using SigInjector.acc_total().

Methods:

`__init__`()
`inject`(n_inj[, seed])	Get some number of events.

Attributes:

acc_total

Total signal acceptance.

__init__()

abstract property acc_total

Total signal acceptance.

This is similar to LLHEvaluatorBase.get_acc_total(), but evaluated in light of all specified signal injection details and typically in terms of the actual simulated signal events that will be used.

abstract inject(n_inj, seed=None)

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

class csky.inj.SinDecRandomizer(sindec_jitter)[source]

Bases: Randomizer

Randomize events in sin declination.

This class applies a jitter of specified strength to the sin(dec) values. The jitter strength should not be chosen much larger than the sin(dec) bin-widths in the PDFs. Randomization should not change the joint PDFs of the observables, since the background PDF is only created once at initialization (and optionally updated with the injected signal via the update_bg parameter). If binning is performed uniformly in sin(dec), a jitter strength of the order of the bin-width is probably fine as long as there are no big jumps in the PDF between consecutive bins. Otherwise, a smaller jitter value should be chosen.

Because the SinDecRandomizer modifies the sin(dec) values instead of the declination directly, it applies a stronger randomization in terms of modified declination at the poles.

Methods:

`__call__`(ev[, seed])	Randomize events in place.
`__init__`(sindec_jitter)	Construct a SinDecRandomizer

__call__(ev, seed=None)[source]

Randomize events in place.

Parameters:

ev (utils.Events) – the events
seed (int or numpy.random.RandomState) – seed spec

__init__(sindec_jitter)[source]

Construct a SinDecRandomizer

Parameters:

sindec_jitter (float) –

jitter strength in sin(dec). The sin(dec) values will be modified via:

sin(dec) += uniform(-sindec_jitter, sindec_jitter)

class csky.inj.SpatialPriorInjector(ana, src, flux, keep, hp_idx, sindec_bandwidth=0.03490658503988659)[source]

Bases: SigInjector

Injection of one or more point-like sources, whose position is not exactly known, but given in the form of a spatial prior (healpix map).

Methods:

`__init__`(ana, src, flux, keep, hp_idx[, ...])	Construct a SpatialPriorInjector.
`inject`(n_inj[, seed, ra, dec])	Get some number of events.

Attributes:

acc_total

Total signal acceptance.

__init__(ana, src, flux, keep, hp_idx, sindec_bandwidth=0.03490658503988659)[source]

Construct a SpatialPriorInjector.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
src (utis.Sources) – the source list. required keys: ra, dec, prior. optional keys: weight, extension
flux (hyp.Flux) – the signal spectrum
keep (list of str) – names of event features to keep

property acc_total

Total signal acceptance.

This is similar to LLHEvaluatorBase.get_acc_total(), but evaluated in light of all specified signal injection details and typically in terms of the actual simulated signal events that will be used.

inject(n_inj, seed=None, ra=None, dec=None)[source]

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

class csky.inj.TemplateBGInjector(ana, template_model, keep, bg_weight_names, flux=None, sig_replace=True, randomizers=None, mcbg=False)[source]

Bases: Injector

Inject typical types of background along with a highly-extended source according to a spatial template.

This class implements addition injection to the typical background (scrambled or MC) according to an all-sky source template(such as for the Galactic “plane”). If per-pixel energy-binned spectra are available then they are used.

Methods:

`__init__`(ana, template_model, keep, ...[, ...])	Construct a TemplateInjector.
`inject`(gp_norm[, n_inj, seed, replace, ...])	Get some number of events.

__init__(ana, template_model, keep, bg_weight_names, flux=None, sig_replace=True, randomizers=None, mcbg=False)[source]

Construct a TemplateInjector.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
template_model (pdf.TemplateSpacePDFRatioModel) – the template space PDF ratio model
keep (list of str) – names of event features to keep
flux (hyp.Flux) – if given, the signal spectrum. otherwise, inject according to per-pixel spectra if available, or use template_model.flux
sig_replace (bool) – If true, sampling with replacement will be used, means an element can be reused multiple times.

inject(gp_norm, n_inj=None, seed=None, replace=False, truth=False, debug=False)[source]

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

class csky.inj.TemplateInjector(ana, template_model, keep, flux=None, sig_replace=True)[source]

Bases: SigInjector

Inject a highly-extended source according to a spatial template.

This class implements injection according to an all-sky source template(such as for the Galactic “plane”). If per-pixel energy-binned spectra are available then they are used, in line with the principle that “you can do whatever you want in your likelihood but you’d better use everything you’ve got for signal injection”, despite the fact that this is relatively slow computationally.

Methods:

`__init__`(ana, template_model, keep[, flux, ...])	Construct a TemplateInjector.
`inject`(n_inj[, seed])	Get some number of events.

Attributes:

acc_total

Total signal acceptance.

__init__(ana, template_model, keep, flux=None, sig_replace=True)[source]

Construct a TemplateInjector.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
template_model (pdf.TemplateSpacePDFRatioModel) – the template space PDF ratio model
keep (list of str) – names of event features to keep
flux (hyp.Flux) – if given, the signal spectrum. otherwise, inject according to per-pixel spectra if available, or use template_model.flux
sig_replace (bool) – If true, sampling with replacement will be used, means an element can be reused multiple times.

property acc_total

Total signal acceptance.

This is similar to LLHEvaluatorBase.get_acc_total(), but evaluated in light of all specified signal injection details and typically in terms of the actual simulated signal events that will be used.

inject(n_inj, seed=None)[source]

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

class csky.inj.TimeWindowSelector(mjd_min, mjd_max)[source]

Bases: Selector

Select events within a given time window.

You probably don’t want to use this class; you might be looking for DataOnOffInjector.

Methods:

`__call__`(arr)	Get a downselected copy of `arr`.
`__init__`(mjd_min, mjd_max)
`mask`(arr)	Get a boolean mask for selected elements of `arr`.
`where`(arr)	Get the indices of selected elements of `arr`.

__call__(arr)[source]

Get a downselected copy of arr.

Parameters:: arr (utils.Arrays) – the Arrays to mask
Return type:: same as arr
Returns:: downselected copy(or potentially view) of arr

__init__(mjd_min, mjd_max)[source]

mask(arr)[source]

Get a boolean mask for selected elements of arr.

Parameters:: arr (utils.Arrays) – the Arrays to mask
Return type:: ndarray of bool
Returns:: the mask

where(arr)[source]

Get the indices of selected elements of arr.

Parameters:: arr (utils.Arrays) – the Arrays to mask
Return type:: ndarray of bool
Returns:: the indices

class csky.inj.TransientInjector(ana, time_model, flux, keep, sindec_bandwidth=0.03490658503988659, spatial_prior=False)[source]

Bases: PointSourceTimeDepInjector

Inject events from transient sources.

Methods:

`__init__`(ana, time_model, flux, keep[, ...])	Construct a TransientInjector.
`get_frac_during_livetime`()	Get the fraction of the total lightcurve covered by a given dataset.
`get_time_model`()	Get the time PDF ratio model.
`inject`(n_inj[, seed])	Get some number of events.
`sample_mjd`(src_n_injs[, seed])	Get per-source `mjd` arrays.

Attributes:

acc_total

Total signal acceptance.

__init__(ana, time_model, flux, keep, sindec_bandwidth=0.03490658503988659, spatial_prior=False)[source]

Construct a TransientInjector.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
time_model (pdf.TimePDFRatioModel) – the time PDF ratio model
flux (hyp.Flux) – the signal spectrum
keep (list of str) – names of event features to keep
sindec_bandwidth (float) – width(not half-width) of the sin(dec) bands from which to draw simulated signal events
spatial_prior (bool) – option to use a spatial prior healpix map for source localisations

property acc_total

Total signal acceptance.

This is similar to LLHEvaluatorBase.get_acc_total(), but evaluated in light of all specified signal injection details and typically in terms of the actual simulated signal events that will be used.

get_frac_during_livetime()[source]

Get the fraction of the total lightcurve covered by a given dataset.

See pdf.TimePDFRatioModel.get_frac_during_livetime().

get_time_model()[source]

Get the time PDF ratio model.

Return type:: pdf.TimePDFRatioModel

inject(n_inj, seed=None)

Get some number of events.

Parameters:

n_inj (int) – number of events to inject
seed (int or numpy.random.RandomState) – seed spec

Return type:

(list of utils.Events, int)

Returns:

(event ensembles, number of pre-emptively excluded events)

sample_mjd(src_n_injs, seed=None)[source]

Get per-source mjd arrays.

Parameters:

src_n_injs (ndarray of int) – per-source number of events to inject
seed (int or numpy.random.RandomState) – seed spec

inspect

Find csky objects inside other csky objects.

Functions:

`get_energy_eval`(L[, key, i])	Get the energy `csky.pdf.PDFRatioEvaluator`.
`get_energy_model`(L[, key])	Get the energy `csky.pdf.PDFRatioModel`.
`get_llh_eval`(L[, key])	Get a specific `csky.llh.LLHEvaluator` from a `csky.llh.MultiLLHEvaluator`.
`get_llh_model`(L[, key])	Get a `csky.llh.LLHModel` from a `csky.llh.MultiLLHEvaluator` or `csky.trial.TrialRunner`.
`get_pdf_ratio_evals`(L[, key])	Get the overall `csky.pdf.PDFRatioEvaluator`.
`get_pdf_ratio_model`(L[, key])	Get the overall `csky.pdf.PDFRatioModel`.
`get_space_eval`(L[, key, i])	Get the spatial `csky.pdf.PDFRatioEvaluator`.
`get_space_model`(L[, key])	Get the spatial `csky.pdf.PDFRatioModel`.
`get_time_eval`(L[, key, i])	Get the temporal `csky.pdf.PDFRatioEvaluator`.
`get_time_model`(L[, key])	Get the temporal `csky.pdf.PDFRatioModel`.

csky.inspect.get_energy_eval(L, key=0, i=0)[source]: Get the energy csky.pdf.PDFRatioEvaluator.

csky.inspect.get_energy_model(L, key=0)[source]: Get the energy csky.pdf.PDFRatioModel.

csky.inspect.get_llh_eval(L, key=0)[source]: Get a specific csky.llh.LLHEvaluator from a csky.llh.MultiLLHEvaluator.

csky.inspect.get_llh_model(L, key=0)[source]: Get a csky.llh.LLHModel from a csky.llh.MultiLLHEvaluator or csky.trial.TrialRunner.

csky.inspect.get_pdf_ratio_evals(L, key=0)[source]: Get the overall csky.pdf.PDFRatioEvaluator.

csky.inspect.get_pdf_ratio_model(L, key=0)[source]: Get the overall csky.pdf.PDFRatioModel.

csky.inspect.get_space_eval(L, key=0, i=0)[source]: Get the spatial csky.pdf.PDFRatioEvaluator.

csky.inspect.get_space_model(L, key=0)[source]: Get the spatial csky.pdf.PDFRatioModel.

csky.inspect.get_time_eval(L, key=0, i=0)[source]: Get the temporal csky.pdf.PDFRatioEvaluator.

csky.inspect.get_time_model(L, key=0)[source]: Get the temporal csky.pdf.PDFRatioModel.

llh

Log likelihood ratio modeling, evaluation, and fitting.

This module provides the log likelihood modeling and evaluation framework for csky. The structure is very similar to that of pdf. Specifically, log likelihood ratios are modeled using LLHModel, and they are evaluated and used for maximum-likelihood parameter fitting using LLHEvaluator or MultiLLHEvaluator, each of which derive from LLHEvaluatorBase which does most of the heavy lifting for parameter fitting.

See LLHModel for more on the likelihood formalism.
See LLHEvaluatorBase.fit for more on the parameter fitting implementation.

Classes:

`LLHEvaluator`(llh_model, evs, n_excluded)	Single-dataset LLH evaluator.
`LLHEvaluatorBase`()	LLHEvaluator base class.
`LLHModel`(ana, pdf_ratio_model[, N_bg, ...])	Model that defines a likelihood.
`MultiLLHEvaluator`(llhs)	Multi-dataset LLH evaluator.

class csky.llh.LLHEvaluator(llh_model, evs, n_excluded)[source]

Bases: LLHEvaluatorBase

Single-dataset LLH evaluator.

This class implements the LLHEvaluatorBase interface for a single dataset.

Methods:

`__init__`(llh_model, evs, n_excluded)	Construct an LLHEvaluator.
`fit`([prior, seeder])	Optimize the test statistic.
`get_acc_total`(**params)	Get the total acceptance for the hypothesis specified by `**params`.
`get_counts_Xs_nsmax`([_masks])	Get total event counts, \(X\) arrays, and maximum \(n_s\) values.
`get_minus_llh_ratio`([ns, prior])	Evaluate the negative of the log likelihood ratio.
`get_raw_calc`([prior, taylor_tol, _masks])	Get a C++ log likelihood ratio calculator.
`get_ts`([ns, prior])	Evaluate the test statistic(for fixed ns and other parameters).
`scan_ts`(ns[, get_ns_from_arg, prior])	Scan -2ln[L(ns=0)/L(ns,...)] over a grid.

__init__(llh_model, evs, n_excluded)[source]

Construct an LLHEvaluator.

Parameters:

llh_model (LLHModel) – the LLH model
evs (list of utils.Events) – the event ensembles(usually data, signal[, signal_2,,, …])
n_excluded (int) – number of events masked out of evs before even looking at the specific per-event pdf ratios. For example, in a point source analysis it is an extremely effective optimization to apply a declination band cut prior to generating, let alone analyzing, randomized realizations of the data; this argument tracks how many events have already been dropped by such cuts pre-emtive cuts.

fit(prior=None, seeder=None, **params)

Optimize the test statistic.

Parameters:

prior (callable) – if given, the prior as a function, of the form f(**params)->float
seeder (callable) – if given, function which generates fit seeds, of the form f(llh_evaluator,bounds_dict)->list(params dicts)
params –
fit parameter specifications, and fit options. fit options are prefixed with _. fit parameters can take one of two forms:
- name=float, e.g. gamma=2.19 (fixed value)
- name=[min_value, seed0, [seed1, ...], max_value], e.g. ``gamma=np.arange(1, 4.01, .25)
  (fit \(\gamma\in[1,4]\), using seeds [1.25, 1.5, …, 4])
fit options include:
- _ns_min (default: 0)
- _ts_min (default: 0)
- _mllh_max (default: -_ts_min/2)
- _ns_tol (default: 1e-4)
- _taylor_tol (default: 1e-3)
- _fit_null (default: False) – whether to try to optimize parameters other than
  \(n_s\) when the best fit is \(n_s=0\)
- _full_output (default: False) – see below

Returns:

at least, (ts[float], fitted_params[dict(str->float)]). if: _full_output=True, additional elements are (fixed_params[dict(str->float)], options[dict(str->value)])

Return type:

tuple

This method is arguably the most important one in csky. It provides flexible and (we think) efficient likelihood optimization. The optimization is performed in two major steps. First, a grid of seeds is computed by using cycler to iterate over all combinations of seed values included in **params. In this step, \(n_s\) optimizations are performed(in C++) for each seed in the cycle, holding other parameters fixed each time. The set of parameters yielding the best log likelihood ratio in that procedure is then taken as the seed for a multi-dimensional optimization using scipy.optimize.fmin_l_bfgs_b(). The results of that fit are finally collated and returned as described above.

get_acc_total(**params)[source]: Get the total acceptance for the hypothesis specified by **params.

get_counts_Xs_nsmax(_masks=None, **params)[source]

Get total event counts, \(X\) arrays, and maximum \(n_s\) values.

Returns:

counts: list of event counts, including masked-out events unless extended=True for the LLHModel.
Xs: list of \(X\) arrays, where \(X\) is the coefficient to \(n_s\) in the log likelihood ratio. For example, when sigsub=False and extended=False, \(X = (S/B - 1) / N\)
nsmax: maximum allowed \(n_s\), which is equal to the sum of the lengths of the \(X\) arrays.

Return type:

tuple of (list of int, list of ndarray of float, int)

get_minus_llh_ratio(ns=0, prior=None, **params): Evaluate the negative of the log likelihood ratio.

get_raw_calc(prior=None, taylor_tol=0.001, _masks=None, **params)

Get a C++ log likelihood ratio calculator.

Parameters:

prior (callable) – the prior function
taylor_tol (float) – distance from the singularity at which the log likelihood ratio should be Taylor-expanded.

Returns:

the log likelihood ratio calculator

Return type:

cext.LogLCalc

get_ts(ns=0, prior=None, **params): Evaluate the test statistic(for fixed ns and other parameters).

scan_ts(ns, get_ns_from_arg=None, prior=None, **params)

Scan -2ln[L(ns=0)/L(ns,…)] over a grid.

Todo

clarify this in detail!

class csky.llh.LLHEvaluatorBase[source]

Bases: object

LLHEvaluator base class.

This abstract base class implements the core logic for likelihood evaluation. The primary functionality is accessed through LLHEvaluatorBase.fit(). Internally, the likelihood is evaluated using the C++ extension, which enables extremely efficient numerical optimization of \(n_s\). When other fit parameters are defined (typically, at least \(\gamma\)), these are optimized numerically using scipy.optimize.fmin_l_bfgs_b(), with the \(n_s\) optimization nested inside each function call performed by fmin_l_bfgs_b.

This class provides all of its functionality via two abstract methods: LLHEvaluatorBase.get_counts_Xs_nsmax() and LLHEvaluatorBase.get_acc_total(). The former provides event counts, per-event \(X\) arrays (see docstring for details), and maximum allowed values for \(n_s\). The latter specifies the total signal acceptance, given the fit parameters (again, typically at least \(\gamma\)). These abstract methods are implemented by both LLHEvaluator (for single-dataset analysis) and MultiLLHEvaluator (for multi-dataset analysis).

All user-facing methods in this class accept kwargs listed as **params. Unless otherwise specified (see e.g. LLHEvaluatorBase.fit()) these kwargs are the fit parameters, e.g. gamma=2.19.

Methods:

`__init__`()
`fit`([prior, seeder])	Optimize the test statistic.
`get_acc_total`(**params)	Get the total signal acceptance.
`get_counts_Xs_nsmax`(**params)	Get total event counts, \(X\) arrays, and maximum \(n_s\) values.
`get_minus_llh_ratio`([ns, prior])	Evaluate the negative of the log likelihood ratio.
`get_raw_calc`([prior, taylor_tol, _masks])	Get a C++ log likelihood ratio calculator.
`get_ts`([ns, prior])	Evaluate the test statistic(for fixed ns and other parameters).
`scan_ts`(ns[, get_ns_from_arg, prior])	Scan -2ln[L(ns=0)/L(ns,...)] over a grid.

__init__()

fit(prior=None, seeder=None, **params)[source]

Optimize the test statistic.

Parameters:

prior (callable) – if given, the prior as a function, of the form f(**params)->float
seeder (callable) – if given, function which generates fit seeds, of the form f(llh_evaluator,bounds_dict)->list(params dicts)
params –
fit parameter specifications, and fit options. fit options are prefixed with _. fit parameters can take one of two forms:
- name=float, e.g. gamma=2.19 (fixed value)
- name=[min_value, seed0, [seed1, ...], max_value], e.g. ``gamma=np.arange(1, 4.01, .25)
  (fit \(\gamma\in[1,4]\), using seeds [1.25, 1.5, …, 4])
fit options include:
- _ns_min (default: 0)
- _ts_min (default: 0)
- _mllh_max (default: -_ts_min/2)
- _ns_tol (default: 1e-4)
- _taylor_tol (default: 1e-3)
- _fit_null (default: False) – whether to try to optimize parameters other than
  \(n_s\) when the best fit is \(n_s=0\)
- _full_output (default: False) – see below

Returns:

at least, (ts[float], fitted_params[dict(str->float)]). if: _full_output=True, additional elements are (fixed_params[dict(str->float)], options[dict(str->value)])

Return type:

tuple

This method is arguably the most important one in csky. It provides flexible and (we think) efficient likelihood optimization. The optimization is performed in two major steps. First, a grid of seeds is computed by using cycler to iterate over all combinations of seed values included in **params. In this step, \(n_s\) optimizations are performed(in C++) for each seed in the cycle, holding other parameters fixed each time. The set of parameters yielding the best log likelihood ratio in that procedure is then taken as the seed for a multi-dimensional optimization using scipy.optimize.fmin_l_bfgs_b(). The results of that fit are finally collated and returned as described above.

abstract get_acc_total(**params)[source]

Get the total signal acceptance.

Returns:: the total signal acceptance
Return type:: float

This method returns the total, as opposed to potentially per-source, signal acceptance, given any required fit parameters. Ultimately, the calculation is performed by whatever AccModel is activated in the PDF ratio model that is being used.

abstract get_counts_Xs_nsmax(**params)[source]

Get total event counts, \(X\) arrays, and maximum \(n_s\) values.

Returns:

counts: list of event counts, including masked-out events unless extended=True for the LLHModel.
Xs: list of \(X\) arrays, where \(X\) is the coefficient to \(n_s\) in the log likelihood ratio. For example, when sigsub=False and extended=False, \(X = (S/B - 1) / N\)
nsmax: maximum allowed \(n_s\), which is equal to the sum of the lengths of the \(X\) arrays.

Return type:

tuple of (list of int, list of ndarray of float, int)

get_minus_llh_ratio(ns=0, prior=None, **params)[source]: Evaluate the negative of the log likelihood ratio.

get_raw_calc(prior=None, taylor_tol=0.001, _masks=None, **params)[source]

Get a C++ log likelihood ratio calculator.

Parameters:

prior (callable) – the prior function
taylor_tol (float) – distance from the singularity at which the log likelihood ratio should be Taylor-expanded.

Returns:

the log likelihood ratio calculator

Return type:

cext.LogLCalc

get_ts(ns=0, prior=None, **params)[source]: Evaluate the test statistic(for fixed ns and other parameters).

scan_ts(ns, get_ns_from_arg=None, prior=None, **params)[source]

Scan -2ln[L(ns=0)/L(ns,…)] over a grid.

Todo

clarify this in detail!

class csky.llh.LLHModel(ana, pdf_ratio_model, N_bg=None, sigsub=False, extended=False, update_bg=False, cut_n_sigma=5, keep_pdfs=False)[source]

Bases: object

Model that defines a likelihood.

This class is used to define a likelihood of roughly the form

\[L(\vec\mu) = \prod_\text{events} \left[ \frac{n_s}{N} \cdot S(\text{event}|\vec\mu) + \left( 1 - \frac{n_s}{N} \right) \cdot B(\text{event}|\vec\mu) \right].\]

The fit parameters \(\vec\mu\) always start with \(n_s\) and typically also include a spectral index \(\gamma\). Three significant modifications are possible:

signal subtracting likelihood: If the baseline background PDF is just the data PDF \(D\), then the signal subtracting likelihood redefines it such that \(D = (n_s/N)\cdot S + (1 - n_s/N)\cdot B\). This modification is enabled using sigsub=True.
extended likelihood: The so-called extended likelihood replaces \(N\to\langle n_b\rangle\) in the above, and then multiplies the above by a Poisson factor \(P(N|\langle n_b\rangle)\). This modification is enabled using extended=True.
update background PDFs: If the baseline background PDF is just the data PDF \(D\), then in principle it should be updated during signal injection trials so that, in those trials, the background PDF is constructed the way it would have been if those events had really occurred. This can make trials more computationally expensive for little gain in accuracy in most cases; however, you may need it in cases where the number of simulated signal events is large. This modification is enabled using update_bg=True.

Methods:

`__call__`(evs, n_excluded)	Obtain a `LLHEvaluator`.
`__init__`(ana, pdf_ratio_model[, N_bg, ...])	Construct a LLHModel.

__call__(evs, n_excluded)[source]

Obtain a LLHEvaluator.

Parameters:

evs (list of utils.Events) – the event ensembles
n_excluded (int) – the number of events excluded from the evs

__init__(ana, pdf_ratio_model, N_bg=None, sigsub=False, extended=False, update_bg=False, cut_n_sigma=5, keep_pdfs=False)[source]

Construct a LLHModel.

Parameters:

pdf_ratio_model (pdf.AccWeightedPDFRatioModel) – the PDF ratio model
N_bg (int) – the number of background events
sigsub (bool) – whether to use the signal subtracting likelihood
extended (bool) – whether to use the extended likelihood
update_bg (bool) – whether to update the background PDFs in signal-injected trials
cut_n_sigma (float) – number of sigmas away from a source an event must be within to be included

class csky.llh.MultiLLHEvaluator(llhs)[source]

Bases: LLHEvaluatorBase

Multi-dataset LLH evaluator.

This class implements the LLHEvaluatorBase interface for multiple datasets. Most of the “heavy lifting” is still done by LLHEvaluator, but this class enables combined analysis by renormalizing the \(X\) arrays according to the relative signal acceptance for each dataset, as well as by ensuring that nsmax is valid for all datasets given those relative signal acceptances.

Methods:

`__init__`(llhs)	Construct a MultiLLHEvaluator.
`fit`([prior, seeder])	Optimize the test statistic.
`get_acc_total`(**params)	Get the total signal acceptance.
`get_counts_Xs_nsmax`([_masks])	Get total event counts, \(X\) arrays, and maximum \(n_s\) values.
`get_minus_llh_ratio`([ns, prior])	Evaluate the negative of the log likelihood ratio.
`get_raw_calc`([prior, taylor_tol, _masks])	Get a C++ log likelihood ratio calculator.
`get_ts`([ns, prior])	Evaluate the test statistic(for fixed ns and other parameters).
`scan_ts`(ns[, get_ns_from_arg, prior])	Scan -2ln[L(ns=0)/L(ns,...)] over a grid.

__init__(llhs)[source]

Construct a MultiLLHEvaluator.

Parameters:: llhs (list of LLHEvaluator) – the single-dataset LLH evaluators

fit(prior=None, seeder=None, **params)

Optimize the test statistic.

Parameters:

prior (callable) – if given, the prior as a function, of the form f(**params)->float
seeder (callable) – if given, function which generates fit seeds, of the form f(llh_evaluator,bounds_dict)->list(params dicts)
params –
fit parameter specifications, and fit options. fit options are prefixed with _. fit parameters can take one of two forms:
- name=float, e.g. gamma=2.19 (fixed value)
- name=[min_value, seed0, [seed1, ...], max_value], e.g. ``gamma=np.arange(1, 4.01, .25)
  (fit \(\gamma\in[1,4]\), using seeds [1.25, 1.5, …, 4])
fit options include:
- _ns_min (default: 0)
- _ts_min (default: 0)
- _mllh_max (default: -_ts_min/2)
- _ns_tol (default: 1e-4)
- _taylor_tol (default: 1e-3)
- _fit_null (default: False) – whether to try to optimize parameters other than
  \(n_s\) when the best fit is \(n_s=0\)
- _full_output (default: False) – see below

Returns:

at least, (ts[float], fitted_params[dict(str->float)]). if: _full_output=True, additional elements are (fixed_params[dict(str->float)], options[dict(str->value)])

Return type:

tuple

This method is arguably the most important one in csky. It provides flexible and (we think) efficient likelihood optimization. The optimization is performed in two major steps. First, a grid of seeds is computed by using cycler to iterate over all combinations of seed values included in **params. In this step, \(n_s\) optimizations are performed(in C++) for each seed in the cycle, holding other parameters fixed each time. The set of parameters yielding the best log likelihood ratio in that procedure is then taken as the seed for a multi-dimensional optimization using scipy.optimize.fmin_l_bfgs_b(). The results of that fit are finally collated and returned as described above.

get_acc_total(**params)[source]

Get the total signal acceptance.

Returns:: the total signal acceptance
Return type:: float

This method returns the total, as opposed to potentially per-source, signal acceptance, given any required fit parameters. Ultimately, the calculation is performed by whatever AccModel is activated in the PDF ratio model that is being used.

get_counts_Xs_nsmax(_masks=None, **params)[source]

Get total event counts, \(X\) arrays, and maximum \(n_s\) values.

Returns:

counts: list of event counts, including masked-out events unless extended=True for the LLHModel.
Xs: list of \(X\) arrays, where \(X\) is the coefficient to \(n_s\) in the log likelihood ratio. For example, when sigsub=False and extended=False, \(X = (S/B - 1) / N\)
nsmax: maximum allowed \(n_s\), which is equal to the sum of the lengths of the \(X\) arrays.

Return type:

tuple of (list of int, list of ndarray of float, int)

get_minus_llh_ratio(ns=0, prior=None, **params): Evaluate the negative of the log likelihood ratio.

get_raw_calc(prior=None, taylor_tol=0.001, _masks=None, **params)

Get a C++ log likelihood ratio calculator.

Parameters:

prior (callable) – the prior function
taylor_tol (float) – distance from the singularity at which the log likelihood ratio should be Taylor-expanded.

Returns:

the log likelihood ratio calculator

Return type:

cext.LogLCalc

get_ts(ns=0, prior=None, **params): Evaluate the test statistic(for fixed ns and other parameters).

scan_ts(ns, get_ns_from_arg=None, prior=None, **params)

Scan -2ln[L(ns=0)/L(ns,…)] over a grid.

Todo

clarify this in detail!

pdf

PDF and PDF ratio implementation.

This module provides the PDF modeling and evaluation framework for csky. Here is an overview of the key concepts:

PDF ratio models

In csky, it is computationally convenient and efficient to work in terms of log likelihood ratios rather than directly in terms of the likelihood itself. This enables us to represent all PDF calculations in terms of PDF ratios. In the standard case, it looks like:

\[\Lambda(\vec\mu) = -2 \ln \left[ \frac{L(n_s=0)}{L(\vec\mu)} \right] = 2\sum_\text{events} \ln \left[ \left(\frac{S(\text{event}|\vec\mu)}{B(\text{event}|\vec\mu)} - 1 \right) \cdot \frac{n_s}{N} + 1 \right],\]

where \(S\) and \(B\) are the signal and background PDFs, and the fit parameters \(\vec\mu\) always start with \(n_s\) and typically also include a spectral index \(\gamma\). Because we have only the ratio \(S/B\), we can define PDF ratios \(W(\text{event}|\vec\mu)=S(\text{event}|\vec\mu)/B(\text{event}|\vec\mu)\).

Because the \(\text{events}\) vary with each trial, and \(\vec\mu\) changes with each log likelihood ratio evaluation, we want a piecemeal prescription for evaluating \(W\). The first step is the PDF ratio model, which you can think of as \(W(?|?)\). In csky, classes providing such a definition are derived from PDFRatioModel. These are the classes that the user is most likely to interact with directly, as they must be configured prior to running any trials.

PDF ratio evaluators

We have discussed trials so far without providing a definition. In csky, a trial is simply a list of one or more event ensembles (structured as utils.Events) along with a count of events that are present in principle but masked out of the detailed calculations in practice.

Once we have a trial, we want to be able to evaluate \(\Lambda(\vec\mu)\). In general, many evaluations will be required to find the parameters \(\hat{\vec\mu}\) which maximize \(\Lambda\). Usually, significant optimizations can be achieved by caching certain values in between evaluations. Therefore, we do not directly ask PDF ratio models to evaluate PDF ratios. Instead we ask them for PDF ratio evaluators, which you can think of as \(W(\text{events}|?)\). In csky, classes providing this functionality are derived from PDFRatioEvaluator.

Acceptance Models

One of the crucial observations during the design of this framework is that we need to care about variable signal acceptance at many steps in the analysis process:

For source stacking, \(S(\text{event}|\vec\mu)\) becomes a weighted average over the values it would take for each individual source. The weighting in this average includes a factor \(w_\text{acc}\), the relative signal acceptance for each source. This typically depends on declination, but in the case of time-dependent analysis it also depends on the relative flux expectations for ecah source at the event arrival times.
For multi-dataset analysis, the overall likelihood is simply a product over the per-dataset likelihoods. However, we do not fit \((n_s)_d\) independently for each dataset \(d\). Instead, these values are constrained such that the ratios \((n_s)_1/(n_s)_2\), etc. are equal to the ratios of the absolute signal acceptances for each dataset, \((w_\text{acc})_1/(w_\text{acc})_2\) and so on.
For counts <–> flux/fluence conversion, we need a factor \(n_s/w_\text{acc}\). If we have a signal injector for the specific spectrum under consideration, this \(w_\text{acc}\) can be obtained directly from the signal MC. However, in general \(w_\text{acc}=w_\text{acc}(\vec\mu)\), and we want to be able to calculate its value for arbitrary \(\hat{\vec\mu}\).

Because variable signal acceptance is relevant in so many places, in csky we treat it as roughly on par with PDF ratio models and evaluators. Specifically, we only evaluate likelihoods in terms of acceptance weighted PDF ratio model/evaluator pairs: AccWeightedPDFRatioModel and AccWeightedPDFRatioEvaluator. Each class implementing AccWeightedPDFRatioModel must provide an acceptance parameterization, which by convention is an instance of a nested class (e.g. PointSourceSpacePDFRatioModel.AccModel) that implements the AccModel interface.

Parameterizations

Often, we can pre-configure some of the functions needed by PDF ratio evaluators, so that they can be re-used for many different analyses. The most common examples are the background space PDF (e.g. BgSinDecParameterization) and the declination- and spectrum-dependent signal acceptance (e.g. SinDecAccParameterization). These parameterizations do not constitute PDF ratio models or evaluators themselves, but they are defined in this module because they are tools used by those classes. Note that for most applications we can also pre-compute an EnergyPDFRatioModel, which is indeed a proper PDFRatioModel.

TODO:

discuss signal subtraction
document TransientTimePDFRatioModel and TransientTimePDFRatioEvaluator

Classes:

`AccModel`()	Acceptance model base class.
`AccWeightedGenericPDFRatioModel`(func, ...[, ...])
`AccWeightedPDFRatioEvaluator`()
`AccWeightedPDFRatioModel`()	Base class for acceptance-weighted PDF ratio models.
`BgAzimuthSinDecParameterization`(ev[, hkw2, ...])	Background space PDF that accounts for azimuthal uncertainty.
`BgSinDecKDEParameterization`(ev[, hkw, skw, ...])	KDE-like parameterization of the background space PDF.
`BgSinDecParameterization`(ev[, dOmega_corr, ...])	Traditional parameterization of the background space PDF.
`BinnedTimePDFRatioEvaluator`(ev, model)
`BinnedTimePDFRatioModel`(ana, src, lcs[, ...])	Binned lightcurve time PDF ratio model.
`CustomFluxEnergyPDFRatioEvaluator`(ev, model)
`CustomFluxEnergyPDFRatioModel`(bg_ev, sig_ev, ...)	Model of energy PDF ratios for a custom flux.
`EnergyPDFRatioEvaluator`(ev, model)
`EnergyPDFRatioModel`(bg_ev, sig_ev[, ...])	Model of energy PDF ratios.
`FitPointSourceSpacePDFRatioEvaluator`(ev, model)
`FitPointSourceSpacePDFRatioModel`(ana, src[, ...])	Space PDF ratio model that fits the sources.
`Gauss2DAngResParameterization`(sig[, flux, ...])	Parameterization of the dec and ra angular resolution in reco sin(dec) and log(energy) bins.
`Gauss2DPointSourceSpacePDFRatioEvaluator`(ev, ...)
`GenericPDFRatioEvaluator`(ev, model[, i])
`GenericPDFRatioModel`(func, features, fits[, ...])	Generic PDF ratio model.
`MCPointSourceSpacePDFRatioEvaluator`(ev, model)
`MultiPDFRatioEvaluator`(ev, model)
`MultiPDFRatioModel`(acc_weighted_model, ...)	Model of a product of separable PDF ratio models.
`NullBgSpaceParameterization`()	Uniform background space PDF.
`NullEnergyPDFRatioEvaluator`(ev, model)
`NullEnergyPDFRatioModel`()	Uniform energy PDF ratio model.
`PDFRatioEvaluator`()	PDF ratio evaluator base class.
`PDFRatioModel`()	PDF ratio model base class.
`PointSourceSpacePDFRatioEvaluator`(ev, model)
`PointSourceSpacePDFRatioModel`(ana, src[, ...])	Space PDF ratio model for one or more point or point-like sources.
`PriorSpacePDFRatioEvaluator`(ev, model[, i])
`PriorSpacePDFRatioModel`(ana, src[, ...])	Space PDF ratio model for one or more point or point-like sources whose
`SinDecAccParameterization`(sig_ev[, hkw, ...])	Signal acceptance parameterization vs sin(dec) over a range of power law spectra.
`SinDecCustomFluxAccParameterization`(sig_ev, flux)	Signal acceptance parameterization vs sin(dec) for a specific spectrum.
`SpaceTimePDFRatioEvaluator`(ev, model)
`SpaceTimePDFRatioModel`(space_model, time_model)	Space * Time acceptance weighted PDF ratio model.
`TemplateSpacePDFRatioEvaluator`(ev, model)
`TemplateSpacePDFRatioModel`(ana, template[, ...])	Template space PDF ratio model.
`TimePDFRatioModel`()	Base class for time PDF ratio models.
`TransientTimePDFRatioEvaluator`(ev, model)
`TransientTimePDFRatioModel`(ana, src[, ...])
`UntriggeredTimePDFRatioEvaluator`(ev, model)
`UntriggeredTimePDFRatioModel`(ana, src[, ...])	Untriggered flare time PDF ratio model.

class csky.pdf.AccModel[source]

Bases: object

Acceptance model base class.

This is a base class for acceptance models. Any acceptance-weighted PDFRatioModel must implement an acceptance model that derives from this class.

Methods:

`__init__`()
`get_acc_per_source`(**params)	Get the absolute acceptance on a per-source basis.
`get_acc_total`(**params)	Get the absolute acceptance, summing over sources.

__init__()

abstract get_acc_per_source(**params)[source]

Get the absolute acceptance on a per-source basis.

Parameters:: params – fit arguments

get_acc_total(**params)[source]

Get the absolute acceptance, summing over sources.

Parameters:: params – fit arguments

class csky.pdf.AccWeightedGenericPDFRatioModel(func, features, fits, ana=None, src=None, acc_param=None, bg_param=None, kent_min=0.12217304763960307, cut_n_sigma=5, sigsub=False)[source]

Bases: AccWeightedPDFRatioModel

Classes:

AccModel(src, acc_param, livetime)

Acceptance model for PointSourceSpacePDFRatioModel.

Methods:

`__call__`(ev[, i])	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(func, features, fits[, ana, src, ...])	Construct an AccGenericPDFRatioModel.
`get_updated`(evs)	Update the GenericPDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

Attributes:

acc_model

The acceptance model (just a formality in this case).

class AccModel(src, acc_param, livetime)[source]

Bases: AccModel

Acceptance model for PointSourceSpacePDFRatioModel. This acceptance model evaluates the time-independent acceptance for each source, weights these by the per-source intrinsic weights, and finally multiplies by the livetime.

Methods:

`__init__`(src, acc_param, livetime)
`get_acc_per_source`(**params)	Get the absolute acceptance on a per-source basis.
`get_acc_total`(**params)	Get the absolute acceptance, summing over sources.

__init__(src, acc_param, livetime)[source]

get_acc_per_source(**params)[source]

Get the absolute acceptance on a per-source basis.

Parameters:: params – fit arguments

get_acc_total(**params)

Get the absolute acceptance, summing over sources.

Parameters:: params – fit arguments

__call__(ev, i=(None, None))[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(func, features, fits, ana=None, src=None, acc_param=None, bg_param=None, kent_min=0.12217304763960307, cut_n_sigma=5, sigsub=False)[source]

Construct an AccGenericPDFRatioModel. Intended for use with funcs that already take acceptance weighting into account.

Parameters:

func (callable) – python object that is callable and returns S/B. If update_bg setting is set to True, then the object must also implement a get_updated method.
features (str->str mapping) – mapping of func argname => Events column name
fits (str->number-or-sequence mapping) – mapping with str keys and number (for fixed values) or sequence (for fittable values) pairs. sequence form should be like (min_value, seed1, seed2, …, max_value)

abstract property acc_model: The acceptance model (just a formality in this case).

get_updated(evs)[source]

Update the GenericPDFRatioModel in light of a set of events.

Parameters:: evs (utils.Events) – the event ensembles

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.AccWeightedPDFRatioEvaluator[source]

Bases: PDFRatioEvaluator

Methods:

`__call__`([mask])	Call self as a function.
`__init__`()

abstract __call__(mask=None, **params): Call self as a function.

__init__()

class csky.pdf.AccWeightedPDFRatioModel[source]

Bases: PDFRatioModel

Base class for acceptance-weighted PDF ratio models.

This is the base class for PDF ratio models which track acceptance weighting. Examples include:

PointSourceSpacePDFRatioModel: Tracks declination-dependent signal acceptance.
SpaceTimePDFRatioModel: Integrates declination-dependent signal acceptance over per-source time PDFs.
MultiPDFRatioModel: Combines exactly one AccWeightedPDFRatioModel with one or more additional non-acceptance-weighted PDFRatioModels`(such as :class:`EnergyPDFRatioModel).

Methods:

`__call__`(ev, i)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`()
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

Attributes:

acc_model

Acceptance model.

abstract __call__(ev, i)

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__()

abstract property acc_model

Acceptance model.

Returns:: the acceptance model
Return type:: AccModel

get_updated(evs)

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.BgAzimuthSinDecParameterization(ev, hkw2={}, skw2={}, hkw={}, skw={}, smooth=None, _other=None)[source]

Bases: object

Background space PDF that accounts for azimuthal uncertainty.

This class implements a zenith-dependentn azimuthal asymmetry that modulates the expectation measured by a BgSinDecParameterization.

Methods:

`__call__`(ev)	Call self as a function.
`__init__`(ev[, hkw2, skw2, hkw, skw, smooth, ...])

__call__(ev)[source]: Call self as a function.

__init__(ev, hkw2={}, skw2={}, hkw={}, skw={}, smooth=None, _other=None)[source]

class csky.pdf.BgSinDecKDEParameterization(ev, hkw={}, skw={}, _other=None)[source]

Bases: BgSinDecParameterization

KDE-like parameterization of the background space PDF.

This class implements a background space PDF using a pseudo-KDE, where the kernel widths are determined by the event angular uncertainties.

Methods:

`__call__`([ev, sindec])	Evaluate for given events or sindecs.
`__init__`(ev[, hkw, skw, _other])	Construct a BgSinDecParameterization.

__call__(ev=[], sindec=None)

Evaluate for given events or sindecs.

Parameters:

ev (utils.Events) – the events
sindec (array-like) – the sin(dec) array)

Returns:

ndarray

__init__(ev, hkw={}, skw={}, _other=None)[source]

Construct a BgSinDecParameterization.

Parameters:

ev (utils.Events) – the background-like events
hkw (mapping) –
skw (mapping) – histlite.Hist.spline_fit() kwargs
_other (BgSinDecParameterization) – another similarly-constructed BgSinDecKDEParameterization to be added to the first one.

class csky.pdf.BgSinDecParameterization(ev, dOmega_corr=[], hkw={}, skw={}, bg_mc_weight='', _other=None)[source]

Bases: object

Traditional parameterization of the background space PDF.

This class implements the traditional background space PDF as a function of the sin of the declination. The PDF satisfies the condition \(\int d\Omega f = \int d\alpha \int d(\sin(\delta)) f = 1\).

Methods:

`__call__`([ev, sindec])	Evaluate for given events or sindecs.
`__init__`(ev[, dOmega_corr, hkw, skw, ...])	Construct a BgSinDecParameterization.

__call__(ev=[], sindec=None)[source]

Evaluate for given events or sindecs.

Parameters:

ev (utils.Events) – the events
sindec (array-like) – the sin(dec) array)

Returns:

ndarray

__init__(ev, dOmega_corr=[], hkw={}, skw={}, bg_mc_weight='', _other=None)[source]

Construct a BgSinDecParameterization.

Parameters:

ev (utils.Events) – the background-like events
hkw (mapping) –
skw (mapping) – histlite.Hist.spline_fit() kwargs
_other (BgSinDecParameterization) – another similarly-constructed BgSinDecKDEParameterization to be added to the first one.

class csky.pdf.BinnedTimePDFRatioEvaluator(ev, model)[source]

Bases: PDFRatioEvaluator

Methods:

`__call__`(**params)	Call self as a function.
`__init__`(ev, model)

__call__(**params)[source]: Call self as a function.

__init__(ev, model)[source]

class csky.pdf.BinnedTimePDFRatioModel(ana, src, lcs, n_seeds_thresh=10, n_seeds_lag=10, range_lag=(-1, 1), thresh_seeding='quantiles', lag_seeding='uniform', use_grl=True, cut_n_sigma=None, thresh_seeds=None, lag_seeds=None)[source]

Bases: TimePDFRatioModel

Binned lightcurve time PDF ratio model.

This class implements a binned lightcurve signal time PDF. The PDFs are specified as histlite histograms. By default the flux threshold is a fit parameter. In principle, stacking of multiple sources should be supported, though this has not been well tested, if at all. Note that in the stacking case, this implementation fits the threshold for each source independently; some refactoring would be required to constrain the thresholds to match.

Todo

probably want to (optionally?) correct
BinnedTimePDFRatioModel.get_frac_during_livetime() to integrate only over true livetime.

Methods:

`__call__`(ev)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(ana, src, lcs[, n_seeds_thresh, ...])	Construct a BinnedTimePDFRatioModel
`get_frac_during_livetime`(**params)	Get the fraction of the total lightcurve covered by a given dataset.
`get_lags`(**params)	Extract per-source lagolds.
`get_lc_pdfs`([threshs_array, lags_array])	Generate normalized per-source time PDFs from lightcurves, given the per-source thresholds.
`get_lc_pdfs_cropped`(**params)	Get cropped signal time PDFs spanning only the livetime of the originally specified sub analysis.
`get_threshs`(**params)	Extract per-source thresholds.
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

__call__(ev)[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(ana, src, lcs, n_seeds_thresh=10, n_seeds_lag=10, range_lag=(-1, 1), thresh_seeding='quantiles', lag_seeding='uniform', use_grl=True, cut_n_sigma=None, thresh_seeds=None, lag_seeds=None)[source]

Construct a BinnedTimePDFRatioModel

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
src (utils.Sources) – the source list(this class just needs to know it’s length)
lcs (array of histlite.Hist) – the lightcurves. values are flux; bin edges are in units of MJD
n_seeds_thresh (int) – number of flux threshold seeds to try
n_seeds_lag (int) – number of lag seeds to try
range_lag (tuple of (float, float)) – range of lag times to fit
thresh_seeding (str) – one of ‘quantiles’ (recommended), ‘log’, ‘linear’,’custom’; determines how threshold seeds are selected
lag_seeding (str) – one of ‘uniform’,’custom’; determines how lag seeds are selected
cut_n_sigma (float) – ignored (but present for consistency with other time models)
thresh_seeds (list) – list of threshold seeds if thresh_seeding is ‘custom’
lag_seeds (list) – list of lag seeds if lag_seeding is ‘custom’

get_frac_during_livetime(**params)[source]

Get the fraction of the total lightcurve covered by a given dataset.

Returns:

fraction of the total integrated signal lightcurve that is covered by: active livetime for a given analysis.SubAnalysis

Return type:

float

get_lags(**params)[source]

Extract per-source lagolds.

Returns:: the flux lags for each source
Return type:: ndarray of float

get_lc_pdfs(threshs_array=None, lags_array=None, **params)[source]

Generate normalized per-source time PDFs from lightcurves, given the per-source thresholds.

Parameters:

threshs_array (ndarray of float) – the thresholds
lags_array (ndarray of float) – the lags in MJD
params – if threshs_array and/or lags_array are not given, they are extracted from the fit parameters which are given as kwargs

Returns:

the signal time PDFs

Return type:

ndarray of histlite.Hist

get_lc_pdfs_cropped(**params)[source]

Get cropped signal time PDFs spanning only the livetime of the originally specified sub analysis.

Crop to sub analysis, or edge of lc pdf

Returns:: the cropped signal time PDFs
Return type:: ndarray of histlite.Hist

The PDFs are not renormalized after being cropped.

get_threshs(**params)[source]

Extract per-source thresholds.

Returns:: the flux thresholds for each source
Return type:: ndarray of float

get_updated(evs)

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.CustomFluxEnergyPDFRatioEvaluator(ev, model)[source]

Bases: PDFRatioEvaluator

Methods:

`__call__`(**params)	Call self as a function.
`__init__`(ev, model)

__call__(**params)[source]: Call self as a function.

__init__(ev, model)[source]

class csky.pdf.CustomFluxEnergyPDFRatioModel(bg_ev, sig_ev, flux, features=['sindec', 'log10energy'], normalize_axes=[1], hkw={}, skw={}, bg_mc_weight='', keep_pdfs=False)[source]

Bases: PDFRatioModel

Model of energy PDF ratios for a custom flux.

This class implements the energy PDF ratio similar to EnergyPDFRatioModel, except for a specific flux which need not be a simple unbroken power law.

Methods:

`__call__`(ev)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(bg_ev, sig_ev, flux[, features, ...])	Construct a CustomFluxEnergyPDFRatioModel.
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

__call__(ev)[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(bg_ev, sig_ev, flux, features=['sindec', 'log10energy'], normalize_axes=[1], hkw={}, skw={}, bg_mc_weight='', keep_pdfs=False)[source]

Construct a CustomFluxEnergyPDFRatioModel.

Parameters:

flux (hyp.Flux) – the spectral hypothesis.
parameters (Other) – see EnergyPDFRatioModel.__init__().

get_updated(evs)[source]

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.EnergyPDFRatioEvaluator(ev, model)[source]

Bases: PDFRatioEvaluator

Methods:

`__call__`([_mask])	Call self as a function.
`__init__`(ev, model)

__call__(_mask=None, **params)[source]: Call self as a function.

__init__(ev, model)[source]

class csky.pdf.EnergyPDFRatioModel(bg_ev, sig_ev, features=['sindec', 'log10energy'], normalize_axes=[1], hkw={}, skw={}, gammas=array([1., 1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 1.875, 2., 2.125, 2.25, 2.375, 2.5, 2.625, 2.75, 2.875, 3., 3.125, 3.25, 3.375, 3.5, 3.625, 3.75, 3.875, 4.]), bg_mc_weight='')[source]

Bases: PDFRatioModel

Model of energy PDF ratios.

This class implements the classic energy PDF ratio. Specifically, it computes 1D or 2D histograms for an ensemble of background-like and signal-like events, and then fits a spline (by default, in terms of log-value) tto the resulting histogram ratio. This is done on a grid of spectral indices \(\gamma`(see also :class:`CustomFluxEnergyPDFRatioModel\) for a fixed, non-power-law spectrum). The dimensions in practice are typically \(\sin(\delta,\log_{10)(E)})\). The axes along which to normalize can be specified; typically this is done along the \(\log_{10}(E)\) axis.

Formally, we can say that \(S=S(\log_{10}(E)|\vec{\mu})\), and similarly for \(B\).

In principle, this class can be used for arbitrary 1D or 2D spline-of-histogram PDF ratios, and it could be readily generalized to a signal weighting parameter other than \(\gamma\). If these cases arise in practice, this class should be refactored accordingly.

The 2D limitation is driven by difficulties fitting higher dimensional splines with scipy. A motivated developer could in principle migrate these fits to photospline, in which case arbitrary dimensionality should be possible.

Methods:

`__call__`(ev)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(bg_ev, sig_ev[, features, ...])	Construct an EnergyPDFRatioModel.
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

__call__(ev)[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(bg_ev, sig_ev, features=['sindec', 'log10energy'], normalize_axes=[1], hkw={}, skw={}, gammas=array([1., 1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 1.875, 2., 2.125, 2.25, 2.375, 2.5, 2.625, 2.75, 2.875, 3., 3.125, 3.25, 3.375, 3.5, 3.625, 3.75, 3.875, 4.]), bg_mc_weight='')[source]

Construct an EnergyPDFRatioModel.

Parameters:

bg_ev (utils.Events) – background-like events
sig_ev (utils.Events) – signal-like events
features (list of str) – the event features to histogram
normalize_axes (list of int) – the axes along which to normalize
hkw (mapping) – histlite.hist() kwargs
skw (mapping) – histlite.Hist.spline_fit() kwargs
gammas (array of float) – spectral indices to consider
bg_mc_weights (str) – if given, background histogram is weighted by bg_ev[bg_mc_weight]. otherwise, the default is equal weights per event.

get_updated(evs)[source]

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.FitPointSourceSpacePDFRatioEvaluator(ev, model, i=(None, None))[source]

Bases: AccWeightedPDFRatioEvaluator

Methods:

`__call__`(**params)	Call self as a function.
`__init__`(ev, model[, i])

__call__(**params)[source]: Call self as a function.

__init__(ev, model, i=(None, None))[source]

class csky.pdf.FitPointSourceSpacePDFRatioModel(ana, src, bg_param=None, acc_param=None, cut_n_sigma=5, cut_extension=0.05235987755982989, sigma=None, sigsub=False)[source]

Bases: AccWeightedPDFRatioModel

Space PDF ratio model that fits the sources.

This model, which is implemented internally in terms of PointSourceSpacePDFRatioModel, adds per source cooordinates \((lpha,\delta)\) to the fitter parameters. This is mainly useful for finding the true hottest spot in a small, already-identified region, in which case the user can ask for such refinement without mentioning this class explicitly by name.

Do not expect good performance for multiple sources. This class may be broken if a time PDF ratio model is used.

Classes:

AccModel(model)

Acceptance model for FitPointSourceSpacePDFRatioModel.

Methods:

`__call__`(ev[, i])	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(ana, src[, bg_param, acc_param, ...])	Construct a FitPointSourceSpacePDFRatioModel
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`params_to_src`(src, **params)	Extract source coordinate fit parameters into a source list.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.
`src_to_params`(src)	Convert a source list into a fit parameters mapping.

Attributes:

acc_model

Acceptance model.

class AccModel(model)[source]

Bases: AccModel

Acceptance model for FitPointSourceSpacePDFRatioModel.

This acceptance model functions similarly to PointSourceSpacePDFRatioModel.AccModel, but it requires the per-source coordinates to be specified in order to perform the evaluation.

Methods:

`__init__`(model)
`get_acc_per_source`(**params)	Get the absolute acceptance on a per-source basis.
`get_acc_total`(**params)	Get the absolute acceptance, summing over sources.

__init__(model)[source]

get_acc_per_source(**params)[source]

Get the absolute acceptance on a per-source basis.

Parameters:: params – fit arguments

get_acc_total(**params)

Get the absolute acceptance, summing over sources.

Parameters:: params – fit arguments

__call__(ev, i=(None, None))[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(ana, src, bg_param=None, acc_param=None, cut_n_sigma=5, cut_extension=0.05235987755982989, sigma=None, sigsub=False)[source]

Construct a FitPointSourceSpacePDFRatioModel

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
src (utils.Sources) – the source list. required keys: ra and dec. optional keys: weight, extension
bg_param – the background space PDF parameterization
acc_param – the signal acceptance parameterization
cut_n_sigma (float) – number of sigmas away from a source an event must be within to be included
cut_extension (float) – additional Gaussian smearing width(in radians) to apply before pre-computing band and/or box cuts
sigma (float) – offset(in radians) of additional seed grid points in (ra,dec)
sigsub (bool) – whether signal subtraction will be enabled

property acc_model

Acceptance model.

Returns:: the acceptance model
Return type:: AccModel

get_updated(evs)

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

static params_to_src(src, **params)[source]

Extract source coordinate fit parameters into a source list.

Parameters:: src (utils.Sources) – base source list(may include weights and extensions)
Returns:: the source list
Return type:: utils.Sources

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

static src_to_params(src)[source]

Convert a source list into a fit parameters mapping.

Parameters:: src (utils.Sources) – source list
Returns:: the fit parameters mapping
Return type:: dict

class csky.pdf.Gauss2DAngResParameterization(sig, flux=PowerLawFlux(gamma=2.5), hkw={}, smooth=None, smooth_bins=None, fit=True)[source]

Bases: object

Parameterization of the dec and ra angular resolution in reco sin(dec) and log(energy) bins.

This is mainly here for reproducing the MESC-2yr analysis. You probably don’t want or need to use it.

Methods:

`__call__`(ev)	Call self as a function.
`__init__`(sig[, flux, hkw, smooth, ...])

__call__(ev)[source]: Call self as a function.

__init__(sig, flux=PowerLawFlux(gamma=2.5), hkw={}, smooth=None, smooth_bins=None, fit=True)[source]

class csky.pdf.Gauss2DPointSourceSpacePDFRatioEvaluator(ev, model, i=(None, None))[source]

Bases: PointSourceSpacePDFRatioEvaluator

Methods:

`__call__`([_mask])	Call self as a function.
`__init__`(ev, model[, i])

__call__(_mask=None, **params): Call self as a function.

__init__(ev, model, i=(None, None))[source]

class csky.pdf.GenericPDFRatioEvaluator(ev, model, i=(None, None))[source]

Bases: PDFRatioEvaluator

Methods:

`__call__`(**params)	Call self as a function.
`__init__`(ev, model[, i])

__call__(**params)[source]: Call self as a function.

__init__(ev, model, i=(None, None))[source]

class csky.pdf.GenericPDFRatioModel(func, features, fits, ana=None, src=None)[source]

Bases: PDFRatioModel

Generic PDF ratio model.

This class allows an arbitrary callable to be used as a PDFRatioModel.

Methods:

`__call__`(ev)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(func, features, fits[, ana, src])	Construct a GenericPDFRatioModel.
`get_updated`(evs)	Update the GenericPDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

__call__(ev)[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(func, features, fits, ana=None, src=None)[source]

Construct a GenericPDFRatioModel.

Parameters:

func (callable) – python object that is callable and returns S/B. If update_bg setting is set to True, then the object must also implement a get_updated method.
features (str->str mapping) – mapping of func argname => Events column name
fits (str->number-or-sequence mapping) – mapping with str keys and number (for fixed values) or sequence (for fittable values) pairs. sequence form should be like (min_value, seed1, seed2, …, max_value)

get_updated(evs)[source]

Update the GenericPDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.MCPointSourceSpacePDFRatioEvaluator(ev, model, i=(None, None))[source]

Bases: AccWeightedPDFRatioEvaluator

Methods:

`__call__`([_mask])	Call self as a function.
`__init__`(ev, model[, i])

__call__(_mask=None, **params)[source]: Call self as a function.

__init__(ev, model, i=(None, None))[source]

class csky.pdf.MultiPDFRatioEvaluator(ev, model)[source]

Bases: AccWeightedPDFRatioEvaluator

Methods:

`__call__`([_mask])	Call self as a function.
`__init__`(ev, model)

__call__(_mask=None, **params)[source]: Call self as a function.

__init__(ev, model)[source]

class csky.pdf.MultiPDFRatioModel(acc_weighted_model, *other_models, **kwargs)[source]

Bases: AccWeightedPDFRatioModel

Model of a product of separable PDF ratio models.

Methods:

`__call__`(ev)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(acc_weighted_model, *other_models, ...)	Construct a MultiPDFRatioModel.
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

Attributes:

acc_model

Acceptance model.

__call__(ev)[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(acc_weighted_model, *other_models, **kwargs)[source]

Construct a MultiPDFRatioModel.

Parameters:

acc_weighted_model (AccWeightedPDFRatioModel) – the acceptance-weighted model
*other_models (PDFRatioModel) – one or more additional models

property acc_model

Acceptance model.

Returns:: the acceptance model
Return type:: AccModel

get_updated(evs)[source]

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)[source]

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.NullBgSpaceParameterization[source]

Bases: object

Uniform background space PDF.

This class is intended to implement a uniform background space PDF. It has not been tested thoroughly… if at all.

Methods:

`__call__`(ev, **kw)	Call self as a function.
`__init__`()

__call__(ev, **kw)[source]: Call self as a function.

__init__()

class csky.pdf.NullEnergyPDFRatioEvaluator(ev, model)[source]

Bases: PDFRatioEvaluator

Methods:

`__call__`(**params)	Call self as a function.
`__init__`(ev, model)

__call__(**params)[source]: Call self as a function.

__init__(ev, model)[source]

class csky.pdf.NullEnergyPDFRatioModel[source]

Bases: PDFRatioModel

Uniform energy PDF ratio model.

This class is intended to implement a uniform energy PDF ratio, i.e. simply \(S/B=1\). It has not been tested thoroughly… if at all. It also is likely useless; if you don’t want an energy PDF to be used, you should build your likelihood without one.

Methods:

`__call__`(ev)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`()
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

__call__(ev)[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__()

get_updated(evs)

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.PDFRatioEvaluator[source]

Bases: object

PDF ratio evaluator base class.

This is the base class for PDF ratio evaluators. Evaluators based on e.g. space, energy, and/or time information must derive from this class. Subclasses apply corresponding PDFRatioModels to a given set of events such that, once further given a set of fit parameter kwargs, they can provide the S/B PDF ratios, as well as the S-tilde/B scrambled signal PDF ratios used in the signal subtracting likelihood.

Users should not need to make direct, explicit use of any PDFRatioEvaluators; they are created automatically by the trial running machinery in trial.

Methods:

`__call__`([mask])	Call self as a function.
`__init__`()

abstract __call__(mask=None, **params)[source]: Call self as a function.

__init__()

class csky.pdf.PDFRatioModel[source]

Bases: object

PDF ratio model base class.

This is the base class for PDF ratio models. Models based on e.g. space, energy, and/or time information must derive from this class. These models essentially define the signal/background discrimination that will be used in the likelihood, and they should be aware of everything that is needed to compute the PDF ratios, once given a set of events. PDFRatioEvaluator instances, created automatically by the trial running machinery in trial, are responsible for the actual evaluation based on a given set of events.

Methods:

`__call__`(ev, i)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`()
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

abstract __call__(ev, i)[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__()

get_updated(evs)[source]

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)[source]

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.PointSourceSpacePDFRatioEvaluator(ev, model, i=(None, None))[source]

Bases: AccWeightedPDFRatioEvaluator

Methods:

`__call__`([_mask])	Call self as a function.
`__init__`(ev, model[, i])

__call__(_mask=None, **params)[source]: Call self as a function.

__init__(ev, model, i=(None, None))[source]

class csky.pdf.PointSourceSpacePDFRatioModel(ana, src, bg_param=None, acc_param=None, cut_n_sigma=5, sigsub=False, mc_pdf=False, kent_min=0.12217304763960307)[source]

Bases: AccWeightedPDFRatioModel

Space PDF ratio model for one or more point or point-like sources.

This class implements the space PDF ratio; it is one of the most important workhorses of csky. The source list should give the coordinates \((\alpha,\delta)\) for one or more sources; multiple sources will result in a stacking analysis. Per-source extensions and intrinsic weights may be given as well.

A declination band cut is applied at a given \(n\sigma\) from each event, using the per-event angular error estimates. If signal subtraction is not enabled, an additional right ascension cut is applied(this is the so-called “box” cut).

Classes:

AccModel(src, acc_param, livetime)

Acceptance model for PointSourceSpacePDFRatioModel.

Methods:

`__call__`(ev[, i])	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(ana, src[, bg_param, acc_param, ...])	Construct a PointSourceSpacePDFRatioModel.
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

Attributes:

acc_model

Acceptance model.

class AccModel(src, acc_param, livetime)[source]

Bases: AccModel

Acceptance model for PointSourceSpacePDFRatioModel.

This acceptance model evaluates the time-independent acceptance for each source, weights these by the per-source intrinsic weights, and finally multiplies by the livetime.

Methods:

`__init__`(src, acc_param, livetime)
`get_acc_per_source`(**params)	Get the absolute acceptance on a per-source basis.
`get_acc_total`(**params)	Get the absolute acceptance, summing over sources.

__init__(src, acc_param, livetime)[source]

get_acc_per_source(**params)[source]

Get the absolute acceptance on a per-source basis.

Parameters:: params – fit arguments

get_acc_total(**params)

Get the absolute acceptance, summing over sources.

Parameters:: params – fit arguments

__call__(ev, i=(None, None))[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(ana, src, bg_param=None, acc_param=None, cut_n_sigma=5, sigsub=False, mc_pdf=False, kent_min=0.12217304763960307)[source]

Construct a PointSourceSpacePDFRatioModel.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
src (utils.Sources) – the source list. required keys: ra and dec. optional keys: weight, extension
bg_param – the background space PDF parameterization
acc_param – the signal acceptance parameterization
cut_n_sigma (float) – number of sigmas away from a source an event must be within to be included
sigsub (bool) – whether signal subtraction will be enabled

property acc_model

Acceptance model.

Returns:: the acceptance model
Return type:: AccModel

get_updated(evs)[source]

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)[source]

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.PriorSpacePDFRatioEvaluator(ev, model, i=(None, None))[source]

Bases: PointSourceSpacePDFRatioEvaluator

Methods:

`__call__`([_mask])	Call self as a function.
`__init__`(ev, model[, i])

__call__(_mask=None, **params)[source]: Call self as a function.

__init__(ev, model, i=(None, None))[source]

class csky.pdf.PriorSpacePDFRatioModel(ana, src, bg_param=None, acc_param=None, cut_n_sigma=5, sigsub=False, mc_pdf=False, kent_min=0.12217304763960307)[source]

Bases: PointSourceSpacePDFRatioModel

Space PDF ratio model for one or more point or point-like sources whose: position is described by a spatial prior

The source list should contain the healpix prior for each of the sources.

Classes:

AccModel(src, acc_param, livetime)

Acceptance model for PriorSpacePDFRatioModel.

Methods:

`__call__`(ev[, i])	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(ana, src[, bg_param, acc_param, ...])	Construct a PointSourceSpacePDFRatioModel.
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

Attributes:

acc_model

Acceptance model.

class AccModel(src, acc_param, livetime)[source]

Bases: AccModel

Acceptance model for PriorSpacePDFRatioModel.

This acceptance model evaluates the time-independent acceptance for each source, weights these by the per-source weights given the fitted positions, and finally multiplies by the livetime.

Methods:

`__init__`(src, acc_param, livetime)
`get_acc_per_source`(**params)	Get the absolute acceptance on a per-source basis.
`get_acc_total`(**params)	Get the absolute acceptance, summing over sources.

__init__(src, acc_param, livetime)[source]

get_acc_per_source(**params)[source]

Get the absolute acceptance on a per-source basis.

Parameters:: params – fit arguments

get_acc_total(**params)

Get the absolute acceptance, summing over sources.

Parameters:: params – fit arguments

__call__(ev, i=(None, None))[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(ana, src, bg_param=None, acc_param=None, cut_n_sigma=5, sigsub=False, mc_pdf=False, kent_min=0.12217304763960307)[source]

Construct a PointSourceSpacePDFRatioModel.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
src (utils.Sources) – the source list. required keys: ra and dec. optional keys: weight, extension
bg_param – the background space PDF parameterization
acc_param – the signal acceptance parameterization
cut_n_sigma (float) – number of sigmas away from a source an event must be within to be included
sigsub (bool) – whether signal subtraction will be enabled

property acc_model

Acceptance model.

Returns:: the acceptance model
Return type:: AccModel

get_updated(evs)

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.SinDecAccParameterization(sig_ev, hkw={}, skw={}, gammas=array([1., 1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 1.875, 2., 2.125, 2.25, 2.375, 2.5, 2.625, 2.75, 2.875, 3., 3.125, 3.25, 3.375, 3.5, 3.625, 3.75, 3.875, 4.]))[source]

Bases: object

Signal acceptance parameterization vs sin(dec) over a range of power law spectra.

This class implements a spline lookup of the signal acceptance \(A(\sin(\delta),\gamma)\). This is the time-independent acceptance – it must be integrated over time to obtain, e.g., the relative acceptance ratios across multiple data taking configurations.

Methods:

`__call__`(a, **params)	Compute acceptance.
`__init__`(sig_ev[, hkw, skw, gammas])	Construct a SinDecAccParameterization.

__call__(a, **params)[source]

Compute acceptance.

Parameters:

a (utils.Arrays) – object with a .dec property
params – fit arguments(should include gamma)

Returns:

the acceptance values

Return type:

ndarray

__init__(sig_ev, hkw={}, skw={}, gammas=array([1., 1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 1.875, 2., 2.125, 2.25, 2.375, 2.5, 2.625, 2.75, 2.875, 3., 3.125, 3.25, 3.375, 3.5, 3.625, 3.75, 3.875, 4.]))[source]

Construct a SinDecAccParameterization.

Parameters:

sig_ev (utils.Events) – the signal-like events
hkw (mapping) – histlite.hist() kwargs
skw (mapping) – histlite.Hist.spline_fit() kwargs
gammas (array-like) – set of gammas to evaluate

class csky.pdf.SinDecCustomFluxAccParameterization(sig_ev, flux, hkw={}, skw={})[source]

Bases: object

Signal acceptance parameterization vs sin(dec) for a specific spectrum.

Parameters:

sig_ev (utils.Events) – the signal-like events
flux (hyp.Flux) – the spectrum of interest
hkw (mapping) – histlite.hist() kwargs
skw (mapping) – histlite.Hist.spline_fit() kwargs

Methods:

`__call__`(a, **params)	Call self as a function.
`__init__`(sig_ev, flux[, hkw, skw])

__call__(a, **params)[source]: Call self as a function.

__init__(sig_ev, flux, hkw={}, skw={})[source]

class csky.pdf.SpaceTimePDFRatioEvaluator(ev, model)[source]

Bases: PDFRatioEvaluator

Methods:

`__call__`([_mask])	Call self as a function.
`__init__`(ev, model)

__call__(_mask=None, **params)[source]: Call self as a function.

__init__(ev, model)[source]

class csky.pdf.SpaceTimePDFRatioModel(space_model, time_model, transient=False)[source]

Bases: AccWeightedPDFRatioModel

Space * Time acceptance weighted PDF ratio model.

This class implements a space * time PDF ratio model, where the acceptance weight is computed by integrating the declination-dependent signal acceptance over the per-source time PDFs.

Classes:

AccModel(spacetime_model)

Acceptance model for SpaceTimePDFRatioModel.

Methods:

`__call__`(ev)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(space_model, time_model[, transient])	Consstruct a SpaceTimePDFRatioModel.
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

Attributes:

acc_model

Acceptance model.

class AccModel(spacetime_model)[source]

Bases: AccModel

Acceptance model for SpaceTimePDFRatioModel.

This acceptance model corrects the acceptance given by the space model’s acc_model, accounting for the fraction of the total livetime covered by the signal time PDF.

Methods:

`__init__`(spacetime_model)
`get_acc_per_source`(**params)	Get the absolute acceptance on a per-source basis.
`get_acc_total`(**params)	Get the absolute acceptance, summing over sources.

__init__(spacetime_model)[source]

get_acc_per_source(**params)[source]

Get the absolute acceptance on a per-source basis.

Parameters:: params – fit arguments

get_acc_total(**params)

Get the absolute acceptance, summing over sources.

Parameters:: params – fit arguments

__call__(ev)[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(space_model, time_model, transient=False)[source]

Consstruct a SpaceTimePDFRatioModel.

Parameters:

space_model (AccWeightedPDFRatioModel) – the space model(should be a point-like source model, not a template model)
time_model (TimePDFRatioModel) – the time model
transient (bool) – whether the time model should be understood as a transient source model. TODO: why is this here? probably a never-finished idea for some optimization?

Todo

why is transient here? probably a never-finished idea for some optimization?

property acc_model

Acceptance model.

Returns:: the acceptance model
Return type:: AccModel

get_updated(evs)[source]

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)[source]

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.TemplateSpacePDFRatioEvaluator(ev, model)[source]

Bases: PDFRatioEvaluator

Methods:

`__call__`([_mask])	Call self as a function.
`__init__`(ev, model)

__call__(_mask=None, **params)[source]: Call self as a function.

__init__(ev, model)[source]

class csky.pdf.TemplateSpacePDFRatioModel(ana, template, bg_param=None, acc_param=None, flux=None, bins_energy=None, sigmas=None, cut_n_sigma=5, sigsub=False, fast_weight=False, sindec_bandwidth=0.017453292519943295, dir='', quiet=False, floor=1e-12, keep_pdfs=False)[source]

Bases: AccWeightedPDFRatioModel

Template space PDF ratio model.

This class implements a template space PDF ratio using healpy [1]. The template should be specified either as a 1D array of shape (npix,) or as a 2D array of shape (npix,nEbins). In the former case, a spectrum must be given, wheres in the latter case the energy bin edges must also be given. Because the Gaussian smoothing calculations can be computationally intense(despite the relatively efficient implementation in healpy), it is possible to specify a directory in which to keep a cache.

If per-pixel binned spectra are given, the multiplication of the map by the signal acceptance is computed in declination bands using the signal MC directly unless requested otherwise(in which case the pixel-averaged spectrum is used).

[1] https://healpy.readthedocs.io/en/latest/

Classes:

AccModel(src, acc_param, livetime, sindec_range)

Methods:

`__call__`(ev)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(ana, template[, bg_param, ...])	Construct a TemplateSpacePDFRatioModel.
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

Attributes:

acc_model

Acceptance model.

class AccModel(src, acc_param, livetime, sindec_range)[source]

Bases: AccModel

Methods:

`__init__`(src, acc_param, livetime, sindec_range)
`get_acc_per_source`(**params)	Get the absolute acceptance on a per-source basis.
`get_acc_total`(**params)	Get the absolute acceptance, summing over sources.

__init__(src, acc_param, livetime, sindec_range)[source]

get_acc_per_source(**params)[source]

Get the absolute acceptance on a per-source basis.

Parameters:: params – fit arguments

get_acc_total(**params)

Get the absolute acceptance, summing over sources.

Parameters:: params – fit arguments

__call__(ev)[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(ana, template, bg_param=None, acc_param=None, flux=None, bins_energy=None, sigmas=None, cut_n_sigma=5, sigsub=False, fast_weight=False, sindec_bandwidth=0.017453292519943295, dir='', quiet=False, floor=1e-12, keep_pdfs=False)[source]

Construct a TemplateSpacePDFRatioModel.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
template (ndarray) – (npix,) array of per-pixel weights, or (npix,nEbins) array of per-pixel and per-energy-bin dN/dE values
bg_param – the background space PDF parameterization
acc_param – the signal acceptance parameterization
flux (hyp.Flux) – the signal spectrum
bins_energy (kind) – the energy bin left edges
sigmas (ndarray) – if given, the Gaussian smoothing grid to consider. if not given, reasonable defaults are chosen depending on whether ana includes tracks or cascades
cut_n_sigma (float) – number of sigmas away from a source an event must be within to be included
sigsub (bool) – whether signal subtraction will be enabled
fast_weight (bool) – whether to compute the acceptance-weighted maps very fast using the acc_param; otherwise, by default small declination bands of signal MC are used for each row of constant declination pixels
sindec_bandwidth (float) – width of declination bands for computing the acceptance-weighted map.
dir (str) – path to a directory in which to cache this template’s state(filename will include ana.key)
quiet (bool) – if true, suppress status printouts during construction
floor (float) – minimum allowed value for the normalized template and smoothed signal space PDFs
keep_pdfs (bool) – when have more pdfs have option of keep bkg pdfs for custom fluxes

property acc_model

Acceptance model.

Returns:: the acceptance model
Return type:: AccModel

get_updated(evs)[source]

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.TimePDFRatioModel[source]

Bases: PDFRatioModel

Base class for time PDF ratio models.

Time PDF ratio models must provide an additional method, TimePDFRatioModel.get_frac_during_livetime(), which gives the fraction of the total integrated signal lightcurve that is covered by active livetime for a given analysis.SubAnalysis.

Methods:

`__call__`(ev, i)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`()
`get_frac_during_livetime`(**params)	Get the fraction of the total lightcurve covered by a given dataset.
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

abstract __call__(ev, i)

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__()

abstract get_frac_during_livetime(**params)[source]

Get the fraction of the total lightcurve covered by a given dataset.

Returns:

fraction of the total integrated signal lightcurve that is covered by: active livetime for a given analysis.SubAnalysis

Return type:

float

get_updated(evs)

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.TransientTimePDFRatioEvaluator(ev, model)[source]

Bases: PDFRatioEvaluator

Methods:

`__call__`(**params)	Call self as a function.
`__init__`(ev, model)

__call__(**params)[source]: Call self as a function.

__init__(ev, model)[source]

class csky.pdf.TransientTimePDFRatioModel(ana, src, use_grl=True, cut_n_sigma=4)[source]

Bases: TimePDFRatioModel

Methods:

`__call__`(ev)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(ana, src[, use_grl, cut_n_sigma])
`get_frac_during_livetime`(**params)	Get the fraction of the total lightcurve covered by a given dataset.
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

__call__(ev)[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(ana, src, use_grl=True, cut_n_sigma=4)[source]

get_frac_during_livetime(**params)[source]

Get the fraction of the total lightcurve covered by a given dataset.

Returns:

fraction of the total integrated signal lightcurve that is covered by: active livetime for a given analysis.SubAnalysis

Return type:

float

get_updated(evs)

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

class csky.pdf.UntriggeredTimePDFRatioEvaluator(ev, model)[source]

Bases: PDFRatioEvaluator

Methods:

`__call__`(**params)	Call self as a function.
`__init__`(ev, model)

__call__(**params)[source]: Call self as a function.

__init__(ev, model)[source]

class csky.pdf.UntriggeredTimePDFRatioModel(ana, src, t0_min=None, t0_max=None, dt_min=1.1574074074074074e-11, dt_max=200, livetime=None, n_seeds_t0=6, n_seeds_dt=5, box=False, box_mode='center', use_grl=True, cut_n_sigma=5)[source]

Bases: TimePDFRatioModel

Untriggered flare time PDF ratio model.

This class implements a Gaussian or box-time-window flare signal hypothesis for which the peak time and flare duration are free to float. Gaussian is the default, but box=True can easily be selected. The anchor time is the center of the flare by default, but box_mode='pre' and box_mode='post' allow precursor and afterglow fits, respectively. The user can specify the allowed range for each of these values.

In principle it is possible to fit for flares from multiple source candidates simultaneously in a stacking analysis, though this is not well tested if at all.

Methods:

`__call__`(ev)	Apply the model to a set of events to obtain a PDFRatioEvaluator.
`__init__`(ana, src[, t0_min, t0_max, dt_min, ...])	Construct a UntriggeredTimePDFRatioModel.
`get_frac_during_livetime`([t0_array, dt_array])	Get the fraction of the total lightcurve covered by a given dataset.
`get_t0_dt`([self, n_src, box_mode])	Extract timing parameters as aligned arrays `(t0,dt)`
`get_updated`(evs)	Update the PDFRatioModel in light of a set of events.
`set_ra`(ra)	Set the right ascension of the(presumably one-and-only) source.

__call__(ev)[source]

Apply the model to a set of events to obtain a PDFRatioEvaluator.

Parameters:

ev (utils.Events) – the events
i (tuple of ([array of int], [array of int])) – i_ev, i_src: aligned arrays indicating pairs of events and sources to which the evaluator should apply

Returns:

PDFRatioEvaluator

__init__(ana, src, t0_min=None, t0_max=None, dt_min=1.1574074074074074e-11, dt_max=200, livetime=None, n_seeds_t0=6, n_seeds_dt=5, box=False, box_mode='center', use_grl=True, cut_n_sigma=5)[source]

Construct a UntriggeredTimePDFRatioModel.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
src (utils.Sources) – the source list(this class just needs to know it’s length)
t0_min (float) – earliest allowed t0 (default: time of first event in the sub subanalysis
t0_max (float) – latest allowed t0 (default: time of first event in the sub subanalysis
dt_min (float) – smallest allowed dt
dt_max (float) – largest allowed dt
n_seed_t0 (int) – number of t0 seeds to try
n_seeds_dt (int) – number of dt seeds to try
cut_n_sigma (float) – number of sigmas away from a source an event must be within to be included

get_frac_during_livetime(t0_array=None, dt_array=None, **params)[source]

Get the fraction of the total lightcurve covered by a given dataset.

Returns:

fraction of the total integrated signal lightcurve that is covered by: active livetime for a given analysis.SubAnalysis

Return type:

float

static get_t0_dt(self=None, n_src=None, box_mode='center', **params)[source]

Extract timing parameters as aligned arrays (t0,dt)

Returns:: the t0 and dt arrays
Return type:: tuple of (ndarray,ndarray)

get_updated(evs)

Update the PDFRatioModel in light of a set of events.

Parameters:: evs (sequence of utils.Events) – the event ensembles

Users do not need to call this method explicitly.

This method updates the PDFRatioModel, given a sequence of event ensembles in the form of utils.Events. This allows an optional mode in which the background PDFs should be recomputed to account not only for data events(evs[0]) but also for injected signal events(evs[1:]).

set_ra(ra)

Set the right ascension of the(presumably one-and-only) source.

Parameters:: ra (float) – the right ascension

Users do not need to call this method explicitly.

This method sets the right ascension of a single source, enabling a powerful optimization when computing skymaps. Specifically, space PDFs can be recomputed at a new coordinate without repeating “band” and/or “box” masking calculations. The default implementation here does nothing, but any space PDFRatioModel’s should update when this method is called.

plotting

Plotting of skymaps, test statistic distributions, PDFs, etc.

TODO: This module is a complete mess. Starting point for things to clean up:

abandon colormaps.py; everyone should be using modern matplotlib by now.
get rid of Plot and SkyPlot
make SkyPlotter more flexible and generally tidy its interface
probably get rid of plot_energy_pdf and plot_gauss_2d_angres_param
add docstrings to everything that remains

Classes:

`Plot`([fig])	Base class for plots.
`SkyPlot`([fig, m, rot, coord])	Skymap plotter.
`SkyPlotter`([coord, projection, pc_kw, cb_kw])	Skymap plotter using matplotlib directly for projections.

class csky.plotting.Plot(fig=None)[source]

Base class for plots.

Methods:

__init__([fig])

__init__(fig=None)[source]

class csky.plotting.SkyPlot(fig=None, m=None, rot=None, coord='C', *a, **kw)[source]

Skymap plotter.

Methods:

`__init__`([fig, m, rot, coord])
`colorbar`([unit])	Draw a colorbar.

__init__(fig=None, m=None, rot=None, coord='C', *a, **kw)[source]

colorbar(unit='')[source]: Draw a colorbar.

class csky.plotting.SkyPlotter(coord='C', projection='aitoff', pc_kw={}, cb_kw={})[source]

Skymap plotter using matplotlib directly for projections.

Methods:

__init__([coord, projection, pc_kw, cb_kw])

__init__(coord='C', projection='aitoff', pc_kw={}, cb_kw={})[source]

seeding

Seeding algorithms for tricky likelihood optimizations.

Time-dependent analyses with fitted flare profiles are notoriously challenging minimizer problems. This module provides UTFSeeder, which basically works and is recommended in the tutorial notebooks. It also provides some other classes that you should ignore.

Classes:

`BoxUTFSeeder`([n_test, threshold, gammas, ...])	Seeder for box-profile untriggered flare search.
`GaussianUTFSeeder`([n_flare, ...])	Seeder for Gaussian-profile untriggered flare search.
`UTFSeeder`([threshold, gammas, test_gammas, ...])	Seeder for box-profile untriggered flare search.

class csky.seeding.BoxUTFSeeder(n_test=100, threshold=1, gammas=array([1., 1.5, 2., 2.5, 3., 3.5, 4.]), weight_boundary_only=True, re_fit=False, debug=False)[source]

Seeder for box-profile untriggered flare search.

NOTE: you probably want UTFSeeder.

Methods:

`__call__`(ML, bounds, fixed_params)	Call self as a function.
`__init__`([n_test, threshold, gammas, ...])

__call__(ML, bounds, fixed_params)[source]: Call self as a function.

__init__(n_test=100, threshold=1, gammas=array([1., 1.5, 2., 2.5, 3., 3.5, 4.]), weight_boundary_only=True, re_fit=False, debug=False)[source]

class csky.seeding.GaussianUTFSeeder(n_flare=array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 25, 50]), remove_duplicates=True, threshold=0.5)[source]

Seeder for Gaussian-profile untriggered flare search.

NOTE: you probably want UTFSeeder.

Methods:

`__call__`(ML, bounds, fixed_params)	Call self as a function.
`__init__`([n_flare, remove_duplicates, threshold])

__call__(ML, bounds, fixed_params)[source]: Call self as a function.

__init__(n_flare=array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 25, 50]), remove_duplicates=True, threshold=0.5)[source]

class csky.seeding.UTFSeeder(threshold=1000, gammas=array([0]), test_gammas=array([1, 2, 3, 4]), n_test=array([1, 2, 3, 4, 5, 10, 15, 20, 25, 30, 35, 40]), perc_best=1, reseed=True, refit=None, debug=False)[source]

Seeder for box-profile untriggered flare search.

This class provides some heuristics for choosing reasonable seeds for Gaussian or box-profile untriggered flare searches. It’s absolutely disgusting, but it does appear to basically work - even for multi-dataset analysis. “Inspired” by some code in psLab.

Methods:

`__call__`(ML, bounds, fixed_params[, _masks])	Call self as a function.
`__init__`([threshold, gammas, test_gammas, ...])

__call__(ML, bounds, fixed_params, _masks=None)[source]: Call self as a function.

__init__(threshold=1000, gammas=array([0]), test_gammas=array([1, 2, 3, 4]), n_test=array([1, 2, 3, 4, 5, 10, 15, 20, 25, 30, 35, 40]), perc_best=1, reseed=True, refit=None, debug=False)[source]

selections

Dataset wrangling module.

Classes:

`Dataset`(key, livetime, sig[, data, bg, ...])	Background and signal data.
`Repository`([local_root, remote_root, username])	Data file repository.

Exceptions:

StopModifications(*a, **kw)

class csky.selections.Dataset(key, livetime, sig, data=None, bg=None, cascades=False, grl=None, range_sindec=(-1, 1))[source]

Bases: object

Background and signal data.

These normally should not be constructed directly; see csky.utils.get_analysis().

Methods:

__init__(key, livetime, sig[, data, bg, ...])

Construct a Dataset.

__init__(key, livetime, sig, data=None, bg=None, cascades=False, grl=None, range_sindec=(-1, 1))[source]

Construct a Dataset.

Parameters:

key (str) – a unique identifier for these datasets
livetime (float) – The livetime of this analysis.
sig (mapping or object) – an object containing at least the following members or keys: dec, ra, true_dec, true_ra, true_energy, oneweight (OneWeight/ngen). Corresponding values should be per-event numpy arrays.
bg (mapping or object, optional) – like data, but must include weight.
data (mapping or object, optional) – an object containing at least the following members or keys: dec, ra, weight. Corresponding values should be per-event numpy arrays, except for weight which may be None.
cascades (bool) – whether this is cascade data as opposed to tracks
grl (csky.utils.Arrays) – the good run list

One of bg or data must be provided. sig, bg, and data will generally need additional keys such as energy (for use with energy PDFs) and/or sigma (for per-event Gaussian signal space PDFs).

class csky.selections.Repository(local_root=None, remote_root=None, username=None)[source]

Bases: object

Data file repository.

This class provides a mechanism for loading data from the UW-Madison stores. If working elsewhere, it is possible to mirror the required data on request.

Methods:

`__init__`([local_root, remote_root, username])	Construct a Repository.
`ensure`(path)	Ensure that the specified file is locally accessible.
`load`(path[, backend, target_class, quiet, ...])	Load a file.

__init__(local_root=None, remote_root=None, username=None)[source]

Construct a Repository.

Parameters:: local_root (str) – local directory in which to read/write files

ensure(path)[source]: Ensure that the specified file is locally accessible.

load(path, backend=None, target_class=<class 'csky.utils.Arrays'>, quiet=False, kw={'convert': True}, mmap=False)[source]: Load a file.

exception csky.selections.StopModifications(*a, **kw)[source]

Bases: Exception

__init__(*a, **kw)[source]

add_note(): Exception.add_note(note) – add a note to the exception

with_traceback(): Exception.with_traceback(tb) – set self.__traceback__ to tb and return self.

timing

Simple task timing utility.

Classes:

Timer()

Simple timer.

class csky.timing.Timer[source]

Simple timer.

Example usage:

# set up timer
import csky as cy
timer = cy.time.Timer()
time = timer.time

# time some things
with time('some lengthy calculation'):
    [do some work here]
with time('another lengthy calculation'):
    [do some more work here]

# get summary
print(timer)

Methods:

__init__()

__init__()[source]

trial

User-level interface for csky trial operations.

This module gives the user-level interface for generating actual or randomized datasets and evaluating and/or maximizing likelihoods. It also includes tools for performing batches of trials; estimating threshold quantities such as sensitivities, discovery potentials, and upper limits; and distributing trials over multiple local cores.

Trial Runner

Programmatically, a trial runner is anything which implements the interface declared in the abstract base class, TrialRunnerBase. Broadly speaking, trial runners are objects which:

Know how to generate actual or randomized datasets, possibly including simulated signal injections.

Know how to apply a csky.llh.LLHModel to the data in order to obtain a csky.llh.LLHEvaluator.

Know how to map operations 1-2 onto multiple datasets to perform multi-dataset analysis.

Know how to associate fit results (e.g. [3.7, 12.36, 2.5]) with column names (e.g. [‘TS, ns, gamma’])

Know how to report per-dataset signal acceptance for use in counts <–> flux calculations.

Thru this interface, TrialRunnerBase provides trial batching, sens/disc/UL calculation, and embarassingly parallel trial multiprocessing.

This module provides the following trial runners:

TrialRunner – for performing trials that involve a single maximum likelihood fit. This includes single source analysis, multi-source stacking, and template analysis. Other trial runners are generally implemented in terms of this one.
SkyScanTrialRunner – for performing trials that search for the hottest spot over much or all of the sky, on a grid of healpix pixels.
SpatialPriorTrialRunner – for performing trials that search for the hottest spot associated with each of one or more spatial priors, again expressed in terms of healpix maps.
MultiTrialRunner – the most generic trial runner, for performing trials that search for the maximum test statistic that can be obtained using the likelihood fits from one or more other arbitrary trial runners.

Classes:

`MultiTrialRunner`(ana, tr_inj, trs[, bgs, ...])
`MultiflareTrialRunner`(ana, srclist, tr_all)
`SkyScanTrialRunner`(ana, get_tr, get_selector)
`SpatialPriorTrial`(evss, n_exs, ra, dec)
`SpatialPriorTrialRunner`(ana, get_tr, ...[, ...])
`Trial`(evss, n_exs)
`TrialRunner`(ana, get_llh, get_injs[, ...])	Main useful trial runner.
`TrialRunnerBase`(ana)	Base class for trial runners.

class csky.trial.MultiTrialRunner(ana, tr_inj, trs, bgs=None, mp_cpus=1, mp_nice=0)[source]

Bases: TrialRunnerBase

Methods:

`__init__`(ana, tr_inj, trs[, bgs, mp_cpus, ...])	Construct a TrialRunnerBase.
`find_n_sig`(ts, beta[, tol, n_batches, ...])	Find the signal strength that exceeds `ts` in `beta` fraction of trials.
`get_acc_total`(**params)	Sum over `TrialRunnerBase.get_acc_total_accs()`.
`get_accs`(**params)	Per-SubAnalysis signal acceptances given fit parameters.
`get_many_fits`(n_trials[, n_sig, poisson, ...])	Get fits for a batch of new trials.
`get_one_fit`([n_sig, poisson, seed, TRUTH, flat])	Get a fit for one new trial.
`get_one_fit_from_trial`(evss_ns_excluded_tuple)	Get a fit, given a trial.
`get_one_trial`([n_sig, poisson, seed, TRUTH])	Get a trial.
`sig_inj_acc_total`()	Sum over `TrialRunnerBase.sig_inj_accs`.
`to_E2dNdE`(ns[, flux])	Get E0^2 * dN/dE.
`to_dNdE`(ns[, flux])	Get dN/dE.
`to_dict`(x)	Convert an ordered per-SubAnalysis quantity into a key->value dict.
`to_model_norm`(ns[, flux])	Get the model normalization.
`to_ns`(E2dNdE[, flux])	Get number of signal events corresponding to a particular flux.

Attributes:

`keys`	The Analysis keys.
`n_ana`	The number of SubAnalysis datasets.
`sig_inj_accs`	Per-SubAnalysis SigInjector total acceptances.

__init__(ana, tr_inj, trs, bgs=None, mp_cpus=1, mp_nice=0)[source]

Construct a TrialRunnerBase.

Parameters:: ana (csky.analysis.Analysis) – the(multi-dataset) Analysis

find_n_sig(ts, beta, tol=0.03, n_batches=6, batch_size=50, coverage=2, n_sig_step=3, n_sig_start=None, first_batch_size=50, min_batch_size=0, max_batch_size=inf, tss=None, trials=None, n_bootstrap=100, mp_cpus=None, seed=None, logging=True, **kw)

Find the signal strength that exceeds ts in beta fraction of trials.

This method estimates sensitivity or discovery potential in terms of a Poisson mean expected number of events n_sig.

Under typical interactive usage, this method estimates sensitivity or discovery potential in two steps. First a preliminary scan finds a very coarse estimate used to define a grid of relevant signal strengths. Then, batches of trials are performed for each n_inj in that grid, and the TS-passing-fraction curve (as a function of n_inj) is interpolated to find a more precise n_sig estimate.

The trial batches are grown iteratively in chunks of size batch_size (or min_batch_size for the first iteration, if that is larger) until one of two stopping conditions is reached: either the batches have reached max_batch_size, or the estimated relative error (n_sig_error/n_sig ratio) is less than tol.

Pre-computed trials can be used if either tss or trials is passed. Then, if tol=np.inf is passed, no new trials will be performed. If n_bootstrap=1 is passed, the error estimate will also be effectively skipped, saving time if these estimates are not needed.

Parameters:

ts (float) – test statistic threshold
beta (float) – fraction of signal trials that should exceed ts
tol (float) – stop main trial loop if relative error (n_sig_error/n_sig) reaches this value or less
n_batches (int) – number of n_inj grid points (columns in the printed log)
batch_size (int) – number of trials per grid point per iteration
coverage (float) – factor to multiply by the initial coarse scan n_sig estimate to obtain the maximum n_inj grid point
n_sig_step (float) – n_sig step size for initial coarse scan
n_sig_start (float) – n_sig starting point for initial coarse scan
first_batch_size (int) – number of trials per iteration for initial coarse scan
min_batch_size (int) – number of trials per n_inj for first iteration of main trial batches
max_batch_size (int) – stop main trial loop if batches reach this size
tss (dict(float,numpy.ndarray)) – TS values from pre-computed trials, with items like n_inj: ts_array
trials (dict(float, csky.utils.Arrays )) – pre-computed trials, with items like n_inj: trials
n_bootstrap (int) – number of times to resample from trial arrays to estimate n_sig_error
mp_cpus (int) – number of cores to use for trial batches
seed – random number generator seed
logging (bool) – whether to print a log as trials are performed
**kw – additional keyword arguments are passed to TrialRunnerBase.get_many_fits().

get_acc_total(**params): Sum over TrialRunnerBase.get_acc_total_accs().

get_accs(**params)[source]: Per-SubAnalysis signal acceptances given fit parameters.

get_many_fits(n_trials, n_sig=0, poisson=False, seed=None, mp_cpus=None, logging=True, **fitter_args)

Get fits for a batch of new trials.

Parameters:

n_trials (int) – number of trials to perform
n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
mp_cpus (int) – number of cores to use
logging (bool) – whether to print a log as trials are performed
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

Returns:

fit results – ts, ns, and any other fit parameters

Return type:

csky.utils.Arrays

get_one_fit(n_sig=0, poisson=False, seed=None, TRUTH=False, flat=True, **fitter_args)

Get a fit for one new trial.

Parameters:

n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
TRUTH (bool) – whether to use background-like(“scrambled”) or real data
flat (bool) – whether to flatten the fit result
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

get_one_fit_from_trial(evss_ns_excluded_tuple, flat=True, **fitter_args)[source]

Get a fit, given a trial.

Parameters:

evss_ns_excluded_tuple (list of lists of csky.utils.Events, list of int) – same as return from TrialRunnerBase.get_one_trial()
flat (bool) – whether to flatten the fit result
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

get_one_trial(n_sig=0, poisson=False, seed=None, TRUTH=False)[source]

Get a trial.

Parameters:

n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
TRUTH (bool) – whether to use background-like(“scrambled”) or real data

Returns:

events(for each SubAnalysis(for each injection origin))
ns_excluded(number of events pre-emptively masked out for each SubAnalysis)

Return type:

list of lists of csky.utils.Events, list of int

property keys: The Analysis keys.

property n_ana: The number of SubAnalysis datasets.

sig_inj_acc_total(): Sum over TrialRunnerBase.sig_inj_accs.

property sig_inj_accs: Per-SubAnalysis SigInjector total acceptances.

to_E2dNdE(ns, flux=None, *a, **kw)

Get E0^2 * dN/dE.

See also hyp.Flux.to_dNdE().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_dNdE(ns, flux=None, *a, **kw)

Get dN/dE.

See also hyp.Flux.to_dNdE().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_dict(x): Convert an ordered per-SubAnalysis quantity into a key->value dict.

to_model_norm(ns, flux=None, **params)

Get the model normalization.

See also hyp.Flux.to_model_norm().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_ns(E2dNdE, flux=None, *a, **kw)

Get number of signal events corresponding to a particular flux.

See also hyp.Flux.to_ns().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

class csky.trial.MultiflareTrialRunner(ana, srclist, tr_all, conf_doc={}, filter_srcs=True, muonflag=False, threshold=1, dt_max=inf, only_max=False, mp_cpus=1, fitter_args={})[source]

Bases: TrialRunnerBase

Methods:

`__init__`(ana, srclist, tr_all[, conf_doc, ...])	Construct a MultiflareTrialRunner.
`find_n_sig`(ts, beta[, tol, n_batches, ...])	Find the signal strength that exceeds `ts` in `beta` fraction of trials.
`get_acc_total`(**params)	Sum over `TrialRunnerBase.get_acc_total_accs()`.
`get_accs`(**params)	Per-SubAnalysis signal acceptances given fit parameters.
`get_many_fits`(n_trials[, n_sig, poisson, ...])	Get fits for a batch of new trials.
`get_one_fit`([Am, Io, alpha, deltaT, ...])	Get a fit for one new trial.
`get_one_fit_from_trial`(evss_ns_excluded_tuple)	Get a fit, given a trial.
`get_one_trial`([Am, Io, alpha, injflares, ...])	Get a trial.
`sig_inj_acc_total`()	Sum over `TrialRunnerBase.sig_inj_accs`.
`to_E2dNdE`(ns[, flux])	Get E0^2 * dN/dE.
`to_dNdE`(ns[, flux])	Get dN/dE.
`to_dict`(x)	Convert an ordered per-SubAnalysis quantity into a key->value dict.
`to_model_norm`(ns[, flux])	Get the model normalization.
`to_ns`(E2dNdE[, flux])	Get number of signal events corresponding to a particular flux.

Attributes:

`keys`	The Analysis keys.
`n_ana`	The number of SubAnalysis datasets.
`sig_inj_accs`	Per-SubAnalysis SigInjector total acceptances.

__init__(ana, srclist, tr_all, conf_doc={}, filter_srcs=True, muonflag=False, threshold=1, dt_max=inf, only_max=False, mp_cpus=1, fitter_args={})[source]: Construct a MultiflareTrialRunner.

find_n_sig(ts, beta, tol=0.03, n_batches=6, batch_size=50, coverage=2, n_sig_step=3, n_sig_start=None, first_batch_size=50, min_batch_size=0, max_batch_size=inf, tss=None, trials=None, n_bootstrap=100, mp_cpus=None, seed=None, logging=True, **kw)

Find the signal strength that exceeds ts in beta fraction of trials.

This method estimates sensitivity or discovery potential in terms of a Poisson mean expected number of events n_sig.

Under typical interactive usage, this method estimates sensitivity or discovery potential in two steps. First a preliminary scan finds a very coarse estimate used to define a grid of relevant signal strengths. Then, batches of trials are performed for each n_inj in that grid, and the TS-passing-fraction curve (as a function of n_inj) is interpolated to find a more precise n_sig estimate.

The trial batches are grown iteratively in chunks of size batch_size (or min_batch_size for the first iteration, if that is larger) until one of two stopping conditions is reached: either the batches have reached max_batch_size, or the estimated relative error (n_sig_error/n_sig ratio) is less than tol.

Pre-computed trials can be used if either tss or trials is passed. Then, if tol=np.inf is passed, no new trials will be performed. If n_bootstrap=1 is passed, the error estimate will also be effectively skipped, saving time if these estimates are not needed.

Parameters:

ts (float) – test statistic threshold
beta (float) – fraction of signal trials that should exceed ts
tol (float) – stop main trial loop if relative error (n_sig_error/n_sig) reaches this value or less
n_batches (int) – number of n_inj grid points (columns in the printed log)
batch_size (int) – number of trials per grid point per iteration
coverage (float) – factor to multiply by the initial coarse scan n_sig estimate to obtain the maximum n_inj grid point
n_sig_step (float) – n_sig step size for initial coarse scan
n_sig_start (float) – n_sig starting point for initial coarse scan
first_batch_size (int) – number of trials per iteration for initial coarse scan
min_batch_size (int) – number of trials per n_inj for first iteration of main trial batches
max_batch_size (int) – stop main trial loop if batches reach this size
tss (dict(float,numpy.ndarray)) – TS values from pre-computed trials, with items like n_inj: ts_array
trials (dict(float, csky.utils.Arrays )) – pre-computed trials, with items like n_inj: trials
n_bootstrap (int) – number of times to resample from trial arrays to estimate n_sig_error
mp_cpus (int) – number of cores to use for trial batches
seed – random number generator seed
logging (bool) – whether to print a log as trials are performed
**kw – additional keyword arguments are passed to TrialRunnerBase.get_many_fits().

get_acc_total(**params): Sum over TrialRunnerBase.get_acc_total_accs().

get_accs(**params)[source]: Per-SubAnalysis signal acceptances given fit parameters.

get_many_fits(n_trials, n_sig=0, poisson=False, seed=None, mp_cpus=None, logging=True, **fitter_args)

Get fits for a batch of new trials.

Parameters:

n_trials (int) – number of trials to perform
n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
mp_cpus (int) – number of cores to use
logging (bool) – whether to print a log as trials are performed
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

Returns:

fit results – ts, ns, and any other fit parameters

Return type:

csky.utils.Arrays

get_one_fit(Am=0, Io=2.0, alpha=3.0, deltaT=100.0, poisson=False, seed=None, TRUTH=False, flat=True, replace_evts=True, **fitter_args)[source]

Get a fit for one new trial.

Parameters:

n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
TRUTH (bool) – whether to use background-like(“scrambled”) or real data
flat (bool) – whether to flatten the fit result
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

get_one_fit_from_trial(evss_ns_excluded_tuple, flat=True, decorrelate=True, logging=False, **fitter_args)[source]

Get a fit, given a trial.

Parameters:

evss_ns_excluded_tuple (list of lists of csky.utils.Events, list of int) – same as return from TrialRunnerBase.get_one_trial()
flat (bool) – whether to flatten the fit result
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

get_one_trial(Am=0, Io=2.0, alpha=3.0, injflares=None, deltaT=None, t0=None, n_sig=0, poisson=False, seed=None, TRUTH=False, replace_evts=False, save_evts=False)[source]

Get a trial.

Parameters:

n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
TRUTH (bool) – whether to use background-like(“scrambled”) or real data

Returns:

events(for each SubAnalysis(for each injection origin))
ns_excluded(number of events pre-emptively masked out for each SubAnalysis)

Return type:

list of lists of csky.utils.Events, list of int

property keys: The Analysis keys.

property n_ana: The number of SubAnalysis datasets.

sig_inj_acc_total(): Sum over TrialRunnerBase.sig_inj_accs.

property sig_inj_accs: Per-SubAnalysis SigInjector total acceptances.

to_E2dNdE(ns, flux=None, *a, **kw)

Get E0^2 * dN/dE.

See also hyp.Flux.to_dNdE().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_dNdE(ns, flux=None, *a, **kw)

Get dN/dE.

See also hyp.Flux.to_dNdE().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_dict(x): Convert an ordered per-SubAnalysis quantity into a key->value dict.

to_model_norm(ns, flux=None, **params)

Get the model normalization.

See also hyp.Flux.to_model_norm().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_ns(E2dNdE, flux=None, *a, **kw)

Get number of signal events corresponding to a particular flux.

See also hyp.Flux.to_ns().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

class csky.trial.SkyScanTrialRunner(ana, get_tr, get_selector, fitter_args={}, param_names=None, nside=128, get_pix=None, pixmask=None, scan_ra=None, scan_dec=None, min_dec=-1.5707963267948966, max_dec=1.5707963267948966, refine_max=False, ts_to_p=None, llh_kw={}, inj_kw={}, src_kw={}, src=None, sig_src=None, mp_cpus=1, mp_scan_cpus=1, mp_nice=0, inj=False)[source]

Bases: TrialRunnerBase

Methods:

`__init__`(ana, get_tr, get_selector[, ...])	Construct a TrialRunnerBase.
`find_n_sig`(ts, beta[, tol, n_batches, ...])	Find the signal strength that exceeds `ts` in `beta` fraction of trials.
`get_acc_total`(**params)	Sum over `TrialRunnerBase.get_acc_total_accs()`.
`get_accs`(**params)	Per-SubAnalysis signal acceptances given fit parameters.
`get_many_fits`(n_trials[, n_sig, poisson, ...])	Get fits for a batch of new trials.
`get_one_fit`([n_sig, poisson, seed, TRUTH, flat])	Get a fit for one new trial.
`get_one_fit_from_trial`(evss_ns_excluded_tuple)	Get a fit, given a trial.
`get_one_trial`([n_sig, poisson, seed, TRUTH, ...])	Get a trial.
`sig_inj_acc_total`()	Sum over `TrialRunnerBase.sig_inj_accs`.
`sig_inj_accs`()	Per-SubAnalysis SigInjector total acceptances.
`to_E2dNdE`(ns[, flux])	Get E0^2 * dN/dE.
`to_dNdE`(ns[, flux])	Get dN/dE.
`to_dict`(x)	Convert an ordered per-SubAnalysis quantity into a key->value dict.
`to_model_norm`(ns[, flux])	Get the model normalization.
`to_ns`(E2dNdE[, flux])	Get number of signal events corresponding to a particular flux.

Attributes:

`keys`	The Analysis keys.
`n_ana`	The number of SubAnalysis datasets.

__init__(ana, get_tr, get_selector, fitter_args={}, param_names=None, nside=128, get_pix=None, pixmask=None, scan_ra=None, scan_dec=None, min_dec=-1.5707963267948966, max_dec=1.5707963267948966, refine_max=False, ts_to_p=None, llh_kw={}, inj_kw={}, src_kw={}, src=None, sig_src=None, mp_cpus=1, mp_scan_cpus=1, mp_nice=0, inj=False)[source]

Construct a TrialRunnerBase.

Parameters:: ana (csky.analysis.Analysis) – the(multi-dataset) Analysis

find_n_sig(ts, beta, tol=0.03, n_batches=6, batch_size=50, coverage=2, n_sig_step=3, n_sig_start=None, first_batch_size=50, min_batch_size=0, max_batch_size=inf, tss=None, trials=None, n_bootstrap=100, mp_cpus=None, seed=None, logging=True, **kw)

Find the signal strength that exceeds ts in beta fraction of trials.

This method estimates sensitivity or discovery potential in terms of a Poisson mean expected number of events n_sig.

Under typical interactive usage, this method estimates sensitivity or discovery potential in two steps. First a preliminary scan finds a very coarse estimate used to define a grid of relevant signal strengths. Then, batches of trials are performed for each n_inj in that grid, and the TS-passing-fraction curve (as a function of n_inj) is interpolated to find a more precise n_sig estimate.

The trial batches are grown iteratively in chunks of size batch_size (or min_batch_size for the first iteration, if that is larger) until one of two stopping conditions is reached: either the batches have reached max_batch_size, or the estimated relative error (n_sig_error/n_sig ratio) is less than tol.

Pre-computed trials can be used if either tss or trials is passed. Then, if tol=np.inf is passed, no new trials will be performed. If n_bootstrap=1 is passed, the error estimate will also be effectively skipped, saving time if these estimates are not needed.

Parameters:

ts (float) – test statistic threshold
beta (float) – fraction of signal trials that should exceed ts
tol (float) – stop main trial loop if relative error (n_sig_error/n_sig) reaches this value or less
n_batches (int) – number of n_inj grid points (columns in the printed log)
batch_size (int) – number of trials per grid point per iteration
coverage (float) – factor to multiply by the initial coarse scan n_sig estimate to obtain the maximum n_inj grid point
n_sig_step (float) – n_sig step size for initial coarse scan
n_sig_start (float) – n_sig starting point for initial coarse scan
first_batch_size (int) – number of trials per iteration for initial coarse scan
min_batch_size (int) – number of trials per n_inj for first iteration of main trial batches
max_batch_size (int) – stop main trial loop if batches reach this size
tss (dict(float,numpy.ndarray)) – TS values from pre-computed trials, with items like n_inj: ts_array
trials (dict(float, csky.utils.Arrays )) – pre-computed trials, with items like n_inj: trials
n_bootstrap (int) – number of times to resample from trial arrays to estimate n_sig_error
mp_cpus (int) – number of cores to use for trial batches
seed – random number generator seed
logging (bool) – whether to print a log as trials are performed
**kw – additional keyword arguments are passed to TrialRunnerBase.get_many_fits().

get_acc_total(**params): Sum over TrialRunnerBase.get_acc_total_accs().

get_accs(**params)[source]: Per-SubAnalysis signal acceptances given fit parameters.

get_many_fits(n_trials, n_sig=0, poisson=False, seed=None, mp_cpus=None, logging=True, **fitter_args)

Get fits for a batch of new trials.

Parameters:

n_trials (int) – number of trials to perform
n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
mp_cpus (int) – number of cores to use
logging (bool) – whether to print a log as trials are performed
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

Returns:

fit results – ts, ns, and any other fit parameters

Return type:

csky.utils.Arrays

get_one_fit(n_sig=0, poisson=False, seed=None, TRUTH=False, flat=True, **fitter_args)

Get a fit for one new trial.

Parameters:

n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
TRUTH (bool) – whether to use background-like(“scrambled”) or real data
flat (bool) – whether to flatten the fit result
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

get_one_fit_from_trial(evss_ns_excluded_tuple, flat=True, logging=False, mp_cpus=None, **fitter_args)[source]

Get a fit, given a trial.

Parameters:

evss_ns_excluded_tuple (list of lists of csky.utils.Events, list of int) – same as return from TrialRunnerBase.get_one_trial()
flat (bool) – whether to flatten the fit result
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

get_one_trial(n_sig=0, poisson=False, seed=None, TRUTH=False, mp_cpus=None)[source]

Get a trial.

Parameters:

n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
TRUTH (bool) – whether to use background-like(“scrambled”) or real data

Returns:

events(for each SubAnalysis(for each injection origin))
ns_excluded(number of events pre-emptively masked out for each SubAnalysis)

Return type:

list of lists of csky.utils.Events, list of int

property keys: The Analysis keys.

property n_ana: The number of SubAnalysis datasets.

sig_inj_acc_total(): Sum over TrialRunnerBase.sig_inj_accs.

sig_inj_accs()[source]: Per-SubAnalysis SigInjector total acceptances.

to_E2dNdE(ns, flux=None, *a, **kw)

Get E0^2 * dN/dE.

See also hyp.Flux.to_dNdE().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_dNdE(ns, flux=None, *a, **kw)

Get dN/dE.

See also hyp.Flux.to_dNdE().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_dict(x): Convert an ordered per-SubAnalysis quantity into a key->value dict.

to_model_norm(ns, flux=None, **params)

Get the model normalization.

See also hyp.Flux.to_model_norm().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_ns(E2dNdE, flux=None, *a, **kw)

Get number of signal events corresponding to a particular flux.

See also hyp.Flux.to_ns().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

class csky.trial.SpatialPriorTrial(evss, n_exs, ra, dec)

Bases: tuple

Methods:

`__init__`()
`count`(value, /)	Return number of occurrences of value.
`index`(value[, start, stop])	Return first index of value.

Attributes:

`dec`	Alias for field number 3
`evss`	Alias for field number 0
`n_exs`	Alias for field number 1
`ra`	Alias for field number 2

__init__()

count(value, /): Return number of occurrences of value.

dec: Alias for field number 3

evss: Alias for field number 0

index(value, start=0, stop=9223372036854775807, /)

Return first index of value.

Raises ValueError if the value is not present.

n_exs: Alias for field number 1

ra: Alias for field number 2

class csky.trial.SpatialPriorTrialRunner(ana, get_tr, get_selector, llh_priors, inj_priors=None, acc_param='acc_param', prior_scan_each=False, prior_mode='stack', fitter_args={}, param_names=None, get_pixmask=None, min_dec=-1.5707963267948966, max_dec=1.5707963267948966, refine_max=False, mp_cpus=1, mp_scan_cpus=1, mp_nice=0, llh_kw={}, inj_kw={}, src_kw={}, inj_src_kw={}, sig_src=None)[source]

Bases: TrialRunnerBase

Methods:

`__init__`(ana, get_tr, get_selector, llh_priors)	Construct a TrialRunnerBase.
`find_n_sig`(ts, beta[, tol, n_batches, ...])	Find the signal strength that exceeds `ts` in `beta` fraction of trials.
`get_acc_total`(**params)	Sum over `TrialRunnerBase.get_acc_total_accs()`.
`get_accs`(**params)	Per-SubAnalysis signal acceptances given fit parameters.
`get_many_fits`(n_trials[, n_sig, poisson, ...])	Get fits for a batch of new trials.
`get_one_fit`([n_sig, poisson, seed, TRUTH, flat])	Get a fit for one new trial.
`get_one_fit_from_trial`(trial[, flat, ...])	Get a fit, given a trial.
`get_one_trial`([n_sig, poisson, seed, TRUTH, ...])	Get a trial.
`sig_inj_acc_total`()	Sum over `TrialRunnerBase.sig_inj_accs`.
`sig_inj_accs`()	Per-SubAnalysis SigInjector total acceptances.
`to_E2dNdE`(ns[, flux])	Get E0^2 * dN/dE.
`to_dNdE`(ns[, flux])	Get dN/dE.
`to_dict`(x)	Convert an ordered per-SubAnalysis quantity into a key->value dict.
`to_model_norm`(ns[, flux])	Get the model normalization.
`to_ns`(E2dNdE[, flux])	Get number of signal events corresponding to a particular flux.

Attributes:

`keys`	The Analysis keys.
`n_ana`	The number of SubAnalysis datasets.

__init__(ana, get_tr, get_selector, llh_priors, inj_priors=None, acc_param='acc_param', prior_scan_each=False, prior_mode='stack', fitter_args={}, param_names=None, get_pixmask=None, min_dec=-1.5707963267948966, max_dec=1.5707963267948966, refine_max=False, mp_cpus=1, mp_scan_cpus=1, mp_nice=0, llh_kw={}, inj_kw={}, src_kw={}, inj_src_kw={}, sig_src=None)[source]

Construct a TrialRunnerBase.

Parameters:: ana (csky.analysis.Analysis) – the(multi-dataset) Analysis

find_n_sig(ts, beta, tol=0.03, n_batches=6, batch_size=50, coverage=2, n_sig_step=3, n_sig_start=None, first_batch_size=50, min_batch_size=0, max_batch_size=inf, tss=None, trials=None, n_bootstrap=100, mp_cpus=None, seed=None, logging=True, **kw)

Find the signal strength that exceeds ts in beta fraction of trials.

This method estimates sensitivity or discovery potential in terms of a Poisson mean expected number of events n_sig.

Under typical interactive usage, this method estimates sensitivity or discovery potential in two steps. First a preliminary scan finds a very coarse estimate used to define a grid of relevant signal strengths. Then, batches of trials are performed for each n_inj in that grid, and the TS-passing-fraction curve (as a function of n_inj) is interpolated to find a more precise n_sig estimate.

The trial batches are grown iteratively in chunks of size batch_size (or min_batch_size for the first iteration, if that is larger) until one of two stopping conditions is reached: either the batches have reached max_batch_size, or the estimated relative error (n_sig_error/n_sig ratio) is less than tol.

Pre-computed trials can be used if either tss or trials is passed. Then, if tol=np.inf is passed, no new trials will be performed. If n_bootstrap=1 is passed, the error estimate will also be effectively skipped, saving time if these estimates are not needed.

Parameters:

ts (float) – test statistic threshold
beta (float) – fraction of signal trials that should exceed ts
tol (float) – stop main trial loop if relative error (n_sig_error/n_sig) reaches this value or less
n_batches (int) – number of n_inj grid points (columns in the printed log)
batch_size (int) – number of trials per grid point per iteration
coverage (float) – factor to multiply by the initial coarse scan n_sig estimate to obtain the maximum n_inj grid point
n_sig_step (float) – n_sig step size for initial coarse scan
n_sig_start (float) – n_sig starting point for initial coarse scan
first_batch_size (int) – number of trials per iteration for initial coarse scan
min_batch_size (int) – number of trials per n_inj for first iteration of main trial batches
max_batch_size (int) – stop main trial loop if batches reach this size
tss (dict(float,numpy.ndarray)) – TS values from pre-computed trials, with items like n_inj: ts_array
trials (dict(float, csky.utils.Arrays )) – pre-computed trials, with items like n_inj: trials
n_bootstrap (int) – number of times to resample from trial arrays to estimate n_sig_error
mp_cpus (int) – number of cores to use for trial batches
seed – random number generator seed
logging (bool) – whether to print a log as trials are performed
**kw – additional keyword arguments are passed to TrialRunnerBase.get_many_fits().

get_acc_total(**params): Sum over TrialRunnerBase.get_acc_total_accs().

get_accs(**params)[source]: Per-SubAnalysis signal acceptances given fit parameters.

get_many_fits(n_trials, n_sig=0, poisson=False, seed=None, mp_cpus=None, logging=True, **fitter_args)

Get fits for a batch of new trials.

Parameters:

n_trials (int) – number of trials to perform
n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
mp_cpus (int) – number of cores to use
logging (bool) – whether to print a log as trials are performed
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

Returns:

fit results – ts, ns, and any other fit parameters

Return type:

csky.utils.Arrays

get_one_fit(n_sig=0, poisson=False, seed=None, TRUTH=False, flat=True, **fitter_args)

Get a fit for one new trial.

Parameters:

n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
TRUTH (bool) – whether to use background-like(“scrambled”) or real data
flat (bool) – whether to flatten the fit result
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

get_one_fit_from_trial(trial, flat=True, logging=False, mp_cpus=None, **fitter_args)[source]

Get a fit, given a trial.

Parameters:

evss_ns_excluded_tuple (list of lists of csky.utils.Events, list of int) – same as return from TrialRunnerBase.get_one_trial()
flat (bool) – whether to flatten the fit result
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

get_one_trial(n_sig=0, poisson=False, seed=None, TRUTH=False, mp_cpus=None)[source]

Get a trial.

Parameters:

n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
TRUTH (bool) – whether to use background-like(“scrambled”) or real data

Returns:

events(for each SubAnalysis(for each injection origin))
ns_excluded(number of events pre-emptively masked out for each SubAnalysis)

Return type:

list of lists of csky.utils.Events, list of int

property keys: The Analysis keys.

property n_ana: The number of SubAnalysis datasets.

sig_inj_acc_total(): Sum over TrialRunnerBase.sig_inj_accs.

sig_inj_accs()[source]: Per-SubAnalysis SigInjector total acceptances.

to_E2dNdE(ns, flux=None, *a, **kw)

Get E0^2 * dN/dE.

See also hyp.Flux.to_dNdE().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_dNdE(ns, flux=None, *a, **kw)

Get dN/dE.

See also hyp.Flux.to_dNdE().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_dict(x): Convert an ordered per-SubAnalysis quantity into a key->value dict.

to_model_norm(ns, flux=None, **params)

Get the model normalization.

See also hyp.Flux.to_model_norm().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_ns(E2dNdE, flux=None, *a, **kw)

Get number of signal events corresponding to a particular flux.

See also hyp.Flux.to_ns().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

class csky.trial.Trial(evss, n_exs)

Bases: tuple

Methods:

`__init__`()
`count`(value, /)	Return number of occurrences of value.
`index`(value[, start, stop])	Return first index of value.

Attributes:

`evss`	Alias for field number 0
`n_exs`	Alias for field number 1

__init__()

count(value, /): Return number of occurrences of value.

evss: Alias for field number 0

index(value, start=0, stop=9223372036854775807, /)

Return first index of value.

Raises ValueError if the value is not present.

n_exs: Alias for field number 1

class csky.trial.TrialRunner(ana, get_llh, get_injs, fitter_args={}, param_names=None, mp_cpus=1, mp_nice=0, llh_kw={}, inj_kw={}, concat_evs=False)[source]

Bases: TrialRunnerBase

Main useful trial runner.

This class performs trials that involve a single maximum likelihood fit. This includes single source analysis, multi-source stacking, and template analysis. Other trial runners are generally implemented in terms of this one.

Methods:

`__init__`(ana, get_llh, get_injs[, ...])	Construct a TrialRunner.
`find_n_sig`(ts, beta[, tol, n_batches, ...])	Find the signal strength that exceeds `ts` in `beta` fraction of trials.
`get_acc_total`(**params)	Sum over `TrialRunnerBase.get_acc_total_accs()`.
`get_accs`(**params)	Per-SubAnalysis signal acceptances given fit parameters.
`get_many_fits`(n_trials[, n_sig, poisson, ...])	Get fits for a batch of new trials.
`get_one_fit`([n_sig, poisson, seed, TRUTH, flat])	Get a fit for one new trial.
`get_one_fit_from_trial`(evss_ns_excluded_tuple)	Get a fit, given a trial.
`get_one_trial`([n_sig, poisson, seed, TRUTH, ...])	Get a trial.
`modify_ra`(ra)	Modify source right ascension in place.
`sig_n_injs`([n_sig, poisson, seed])	Multinomial sample from `TrialRunner.sig_inj_probs` .
`to_E2dNdE`(ns[, flux])	Get E0^2 * dN/dE.
`to_dNdE`(ns[, flux])	Get dN/dE.
`to_dict`(x)	Convert an ordered per-SubAnalysis quantity into a key->value dict.
`to_model_norm`(ns[, flux])	Get the model normalization.
`to_ns`(E2dNdE[, flux])	Get number of signal events corresponding to a particular flux.

Attributes:

`keys`	The Analysis keys.
`n_ana`	The number of SubAnalysis datasets.
`sig_inj_acc_total`	Sum over `TrialRunnerBase.sig_inj_accs`.
`sig_inj_accs`	Per-SubAnalysis SigInjector total acceptances.
`sig_inj_probs`	Relative weight of each signal injector (by signal acceptance).
`sig_injs`	The signal injectors.

__init__(ana, get_llh, get_injs, fitter_args={}, param_names=None, mp_cpus=1, mp_nice=0, llh_kw={}, inj_kw={}, concat_evs=False)[source]

Construct a TrialRunner.

Parameters:

ana (analysis.SubAnalysis) – the sub analysis
get_llh (callable) – get_llh(ana, **kwargs) -> LLHEvaluator
get_injs (callable) – get_injs(ana, **kwargs) -> (truth, bg, sig), each of which must be an instance of inj.Injector, or should be None. Specifically, as an optimization, in cases where we know we don’t need the signal injector (e.g. when a TrialRunner is merely an intermediary for a SkyScanTrialRunner), the call will look like get_injs(..., inj=False)`, in which case ``sig can be None.
fitter_args (dict) – arguments to pass to LLHEvaluatorBase.fit()
param_names (list of str) – ordered list of parameter names. This ensures that fit results can be represented sequentially in a reliable way.
mp_cpus (int) – default number of cores to use for trial batches (can be overridden at the specific call sites, e.g. tr.get_many_fits(..., mp_cpus=2))
mp_nice (int) – OS-level niceness (TODO: actually this is currently ignored, but it should be easy to respect if we want to)
llh_kw (dict) – additional kwargs to pass to get_llh
inj_kw (dict) – additional kwargs to pass to get_injs

Todo

mp_nice is currently ignored, but it should
be easy to respect if we want to

find_n_sig(ts, beta, tol=0.03, n_batches=6, batch_size=50, coverage=2, n_sig_step=3, n_sig_start=None, first_batch_size=50, min_batch_size=0, max_batch_size=inf, tss=None, trials=None, n_bootstrap=100, mp_cpus=None, seed=None, logging=True, **kw)

Find the signal strength that exceeds ts in beta fraction of trials.

This method estimates sensitivity or discovery potential in terms of a Poisson mean expected number of events n_sig.

Under typical interactive usage, this method estimates sensitivity or discovery potential in two steps. First a preliminary scan finds a very coarse estimate used to define a grid of relevant signal strengths. Then, batches of trials are performed for each n_inj in that grid, and the TS-passing-fraction curve (as a function of n_inj) is interpolated to find a more precise n_sig estimate.

The trial batches are grown iteratively in chunks of size batch_size (or min_batch_size for the first iteration, if that is larger) until one of two stopping conditions is reached: either the batches have reached max_batch_size, or the estimated relative error (n_sig_error/n_sig ratio) is less than tol.

Pre-computed trials can be used if either tss or trials is passed. Then, if tol=np.inf is passed, no new trials will be performed. If n_bootstrap=1 is passed, the error estimate will also be effectively skipped, saving time if these estimates are not needed.

Parameters:

ts (float) – test statistic threshold
beta (float) – fraction of signal trials that should exceed ts
tol (float) – stop main trial loop if relative error (n_sig_error/n_sig) reaches this value or less
n_batches (int) – number of n_inj grid points (columns in the printed log)
batch_size (int) – number of trials per grid point per iteration
coverage (float) – factor to multiply by the initial coarse scan n_sig estimate to obtain the maximum n_inj grid point
n_sig_step (float) – n_sig step size for initial coarse scan
n_sig_start (float) – n_sig starting point for initial coarse scan
first_batch_size (int) – number of trials per iteration for initial coarse scan
min_batch_size (int) – number of trials per n_inj for first iteration of main trial batches
max_batch_size (int) – stop main trial loop if batches reach this size
tss (dict(float,numpy.ndarray)) – TS values from pre-computed trials, with items like n_inj: ts_array
trials (dict(float, csky.utils.Arrays )) – pre-computed trials, with items like n_inj: trials
n_bootstrap (int) – number of times to resample from trial arrays to estimate n_sig_error
mp_cpus (int) – number of cores to use for trial batches
seed – random number generator seed
logging (bool) – whether to print a log as trials are performed
**kw – additional keyword arguments are passed to TrialRunnerBase.get_many_fits().

get_acc_total(**params): Sum over TrialRunnerBase.get_acc_total_accs().

get_accs(**params)[source]: Per-SubAnalysis signal acceptances given fit parameters.

get_many_fits(n_trials, n_sig=0, poisson=False, seed=None, mp_cpus=None, logging=True, **fitter_args)

Get fits for a batch of new trials.

Parameters:

n_trials (int) – number of trials to perform
n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
mp_cpus (int) – number of cores to use
logging (bool) – whether to print a log as trials are performed
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

Returns:

fit results – ts, ns, and any other fit parameters

Return type:

csky.utils.Arrays

get_one_fit(n_sig=0, poisson=False, seed=None, TRUTH=False, flat=True, **fitter_args)

Get a fit for one new trial.

Parameters:

n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
TRUTH (bool) – whether to use background-like(“scrambled”) or real data
flat (bool) – whether to flatten the fit result
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

get_one_fit_from_trial(evss_ns_excluded_tuple, flat=True, **fitter_args)[source]

Get a fit, given a trial.

Parameters:

evss_ns_excluded_tuple (list of lists of csky.utils.Events, list of int) – same as return from TrialRunnerBase.get_one_trial()
flat (bool) – whether to flatten the fit result
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

get_one_trial(n_sig=0, poisson=False, seed=None, TRUTH=False, debug=False)[source]

Get a trial.

Parameters:

n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
TRUTH (bool) – whether to use background-like(“scrambled”) or real data

Returns:

events(for each SubAnalysis(for each injection origin))
ns_excluded(number of events pre-emptively masked out for each SubAnalysis)

Return type:

list of lists of csky.utils.Events, list of int

property keys: The Analysis keys.

modify_ra(ra)[source]

Modify source right ascension in place.

This method is used internally for efficiently computing sky scans. Users should not need to call this method.

property n_ana: The number of SubAnalysis datasets.

property sig_inj_acc_total: Sum over TrialRunnerBase.sig_inj_accs.

property sig_inj_accs: Per-SubAnalysis SigInjector total acceptances.

property sig_inj_probs: Relative weight of each signal injector (by signal acceptance).

property sig_injs: The signal injectors.

sig_n_injs(n_sig=0, poisson=False, seed=None)[source]: Multinomial sample from TrialRunner.sig_inj_probs .

to_E2dNdE(ns, flux=None, *a, **kw)

Get E0^2 * dN/dE.

See also hyp.Flux.to_dNdE().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_dNdE(ns, flux=None, *a, **kw)

Get dN/dE.

See also hyp.Flux.to_dNdE().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_dict(x): Convert an ordered per-SubAnalysis quantity into a key->value dict.

to_model_norm(ns, flux=None, **params)

Get the model normalization.

See also hyp.Flux.to_model_norm().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_ns(E2dNdE, flux=None, *a, **kw)

Get number of signal events corresponding to a particular flux.

See also hyp.Flux.to_ns().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

class csky.trial.TrialRunnerBase(ana)[source]

Bases: object

Base class for trial runners.

This is the base class for csky’s somewhat idiosyncratic concept of trial runners. In short, a trial runner is an object which knows how to:

generate ensembles of events(optionally including signal injection)
find the parameters that optimize the likelihood for given event ensembles
repeat 1-2 en masse

In most work with csky, classes derived from this one do account for the vast majority of user-library interaciton. This base class provides a number of methods for performing(1-2) or (3) very easily; converting counts <–> flux-or-fluence; conveniently estimating threshold quantities such as sensitivity, discovery potential, or upper limits; and otherwise inspecting the source data. Derived classes are responsible only for implementing(1-2); providing the the source data. Derived classes are responsible only for implementing(1-2); providing the named array mapping for (3); and of course tracking the signal acceptance of each dataset for use in those counts <–> flux-or-fluence calculations.

Methods:

`__init__`(ana)	Construct a TrialRunnerBase.
`find_n_sig`(ts, beta[, tol, n_batches, ...])	Find the signal strength that exceeds `ts` in `beta` fraction of trials.
`get_acc_total`(**params)	Sum over `TrialRunnerBase.get_acc_total_accs()`.
`get_accs`(**params)	Per-SubAnalysis signal acceptances given fit parameters.
`get_many_fits`(n_trials[, n_sig, poisson, ...])	Get fits for a batch of new trials.
`get_one_fit`([n_sig, poisson, seed, TRUTH, flat])	Get a fit for one new trial.
`get_one_fit_from_trial`(evss_ns_excluded_tuple)	Get a fit, given a trial.
`get_one_trial`([n_sig, poisson, seed, TRUTH])	Get a trial.
`sig_inj_acc_total`()	Sum over `TrialRunnerBase.sig_inj_accs`.
`to_E2dNdE`(ns[, flux])	Get E0^2 * dN/dE.
`to_dNdE`(ns[, flux])	Get dN/dE.
`to_dict`(x)	Convert an ordered per-SubAnalysis quantity into a key->value dict.
`to_model_norm`(ns[, flux])	Get the model normalization.
`to_ns`(E2dNdE[, flux])	Get number of signal events corresponding to a particular flux.

Attributes:

`keys`	The Analysis keys.
`n_ana`	The number of SubAnalysis datasets.
`sig_inj_accs`	Per-SubAnalysis SigInjector total acceptances.

__init__(ana)[source]

Construct a TrialRunnerBase.

Parameters:: ana (csky.analysis.Analysis) – the(multi-dataset) Analysis

find_n_sig(ts, beta, tol=0.03, n_batches=6, batch_size=50, coverage=2, n_sig_step=3, n_sig_start=None, first_batch_size=50, min_batch_size=0, max_batch_size=inf, tss=None, trials=None, n_bootstrap=100, mp_cpus=None, seed=None, logging=True, **kw)[source]

Find the signal strength that exceeds ts in beta fraction of trials.

This method estimates sensitivity or discovery potential in terms of a Poisson mean expected number of events n_sig.

Under typical interactive usage, this method estimates sensitivity or discovery potential in two steps. First a preliminary scan finds a very coarse estimate used to define a grid of relevant signal strengths. Then, batches of trials are performed for each n_inj in that grid, and the TS-passing-fraction curve (as a function of n_inj) is interpolated to find a more precise n_sig estimate.

The trial batches are grown iteratively in chunks of size batch_size (or min_batch_size for the first iteration, if that is larger) until one of two stopping conditions is reached: either the batches have reached max_batch_size, or the estimated relative error (n_sig_error/n_sig ratio) is less than tol.

Pre-computed trials can be used if either tss or trials is passed. Then, if tol=np.inf is passed, no new trials will be performed. If n_bootstrap=1 is passed, the error estimate will also be effectively skipped, saving time if these estimates are not needed.

Parameters:

ts (float) – test statistic threshold
beta (float) – fraction of signal trials that should exceed ts
tol (float) – stop main trial loop if relative error (n_sig_error/n_sig) reaches this value or less
n_batches (int) – number of n_inj grid points (columns in the printed log)
batch_size (int) – number of trials per grid point per iteration
coverage (float) – factor to multiply by the initial coarse scan n_sig estimate to obtain the maximum n_inj grid point
n_sig_step (float) – n_sig step size for initial coarse scan
n_sig_start (float) – n_sig starting point for initial coarse scan
first_batch_size (int) – number of trials per iteration for initial coarse scan
min_batch_size (int) – number of trials per n_inj for first iteration of main trial batches
max_batch_size (int) – stop main trial loop if batches reach this size
tss (dict(float,numpy.ndarray)) – TS values from pre-computed trials, with items like n_inj: ts_array
trials (dict(float, csky.utils.Arrays )) – pre-computed trials, with items like n_inj: trials
n_bootstrap (int) – number of times to resample from trial arrays to estimate n_sig_error
mp_cpus (int) – number of cores to use for trial batches
seed – random number generator seed
logging (bool) – whether to print a log as trials are performed
**kw – additional keyword arguments are passed to TrialRunnerBase.get_many_fits().

get_acc_total(**params)[source]: Sum over TrialRunnerBase.get_acc_total_accs().

abstract get_accs(**params)[source]: Per-SubAnalysis signal acceptances given fit parameters.

get_many_fits(n_trials, n_sig=0, poisson=False, seed=None, mp_cpus=None, logging=True, **fitter_args)[source]

Get fits for a batch of new trials.

Parameters:

n_trials (int) – number of trials to perform
n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
mp_cpus (int) – number of cores to use
logging (bool) – whether to print a log as trials are performed
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

Returns:

fit results – ts, ns, and any other fit parameters

Return type:

csky.utils.Arrays

get_one_fit(n_sig=0, poisson=False, seed=None, TRUTH=False, flat=True, **fitter_args)[source]

Get a fit for one new trial.

Parameters:

n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
TRUTH (bool) – whether to use background-like(“scrambled”) or real data
flat (bool) – whether to flatten the fit result
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

abstract get_one_fit_from_trial(evss_ns_excluded_tuple, flat=True, **fitter_args)[source]

Get a fit, given a trial.

Parameters:

evss_ns_excluded_tuple (list of lists of csky.utils.Events, list of int) – same as return from TrialRunnerBase.get_one_trial()
flat (bool) – whether to flatten the fit result
fitter_args – same as params in csky.llh.LLHEvaluatorBase.fit()

abstract get_one_trial(n_sig=0, poisson=False, seed=None, TRUTH=False)[source]

Get a trial.

Parameters:

n_sig (int) – number of signal events to inject
poisson (bool) – whether n_sig is specific or just a Poissonian mean
seed (int) – seed spec
TRUTH (bool) – whether to use background-like(“scrambled”) or real data

Returns:

events(for each SubAnalysis(for each injection origin))
ns_excluded(number of events pre-emptively masked out for each SubAnalysis)

Return type:

list of lists of csky.utils.Events, list of int

property keys: The Analysis keys.

property n_ana: The number of SubAnalysis datasets.

sig_inj_acc_total()[source]: Sum over TrialRunnerBase.sig_inj_accs.

abstract property sig_inj_accs: Per-SubAnalysis SigInjector total acceptances.

to_E2dNdE(ns, flux=None, *a, **kw)[source]

Get E0^2 * dN/dE.

See also hyp.Flux.to_dNdE().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_dNdE(ns, flux=None, *a, **kw)[source]

Get dN/dE.

See also hyp.Flux.to_dNdE().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_dict(x)[source]: Convert an ordered per-SubAnalysis quantity into a key->value dict.

to_model_norm(ns, flux=None, **params)[source]

Get the model normalization.

See also hyp.Flux.to_model_norm().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

to_ns(E2dNdE, flux=None, *a, **kw)[source]

Get number of signal events corresponding to a particular flux.

See also hyp.Flux.to_ns().

Todo

clarify sig_inj_acc_total vs get_acc_total situation in both code and docs

utils

Utility classes and functions.

Classes:

`Arrays`([init, names, convert, _require_array])	Optimized DataFrame-like data structure.
`Events`(a, *kw)	Subclass of `Arrays` automatically providing certain event properties.
`Sources`(a, *kw)	Subclass of `Arrays` automatically providing certain source properties.

Functions:

`ensure_dir`(dirname)	Make sure `dirname` exists and is a directory.
`get_pixmask`(nside, events[, radius])	Query disc of given radius around selected events.
`get_pixval`(map, ra, dec)	Get the value of a healpix map at specific equatorial coordinates.
`get_random`([seed])	Get numpy RandomState.
`intersection_time_intervals`(intervals_1, ...)	Calculate intersections between two sorted lists containing continuous, non-overlapping intervals
`sources`(ra, dec, **kw)	Convenience wrapper for constructing `Sources` instances.
`union_continuous_intervals`(intervals)	Calculate the union of a given list of continuous intervals

class csky.utils.Arrays(init={}, names=None, convert=False, _require_array=True, **kw)[source]

Bases: object

Optimized DataFrame-like data structure.

This is the baseline data structure for all event data in csky. It acts a lot like a pandas DataFrame, and can be converted to one with the .as_dataframe property. Conversion to astropy Table and to numpy structured array are supported by the .as_table and .as_array properties.

Each stored array is treated as a column and can be accessed like, e.g., data['ra'] or data.ra. Assignment of new columns is only available by the [] syntax, e.g. data['something_else'] = some_new_array.

Additional dictionary-like access is possible through the .keys(), .values(), and .items() methods.

When the [] operand is not a string, slicing/masking is performed, e.g. src_0_1_2 = srcs[:3] or north_data = data[data.dec > 0]. Iteration over an Arrays yields per-row instances, that is, new Arrays instances with column lengths of 1.

This class is easier than DataFrames to use for computationally-intense work because it does very little to impose and maintain column alignment. This means it can expose direct references to the underlying numpy arrays rather than routing through some kind of Series objects.

Methods:

__init__([init, names, convert, _require_array])

Construct an Arrays collection.

__init__(init={}, names=None, convert=False, _require_array=True, **kw)[source]

Construct an Arrays collection.

Parameters:

init (numpy.ndarray or object with keys() or ._dict.keys() method) – the data
names (list of str) – the column names
convert (bool) – if True, add references converting from skylab to csky naming convention

class csky.utils.Events(*a, **kw)[source]

Bases: Arrays

Subclass of Arrays automatically providing certain event properties.

Bonus properties include sindec, log10energy, ra_deg, dec_deg, sigma_deg, and (for MC where true values are known) dpsi_deg.

Methods:

__init__(*a, **kw)

Construct an Arrays collection.

__init__(*a, **kw)[source]

Construct an Arrays collection.

Parameters:

init (numpy.ndarray or object with keys() or ._dict.keys() method) – the data
names (list of str) – the column names
convert (bool) – if True, add references converting from skylab to csky naming convention

class csky.utils.Sources(*a, **kw)[source]

Bases: Arrays

Subclass of Arrays automatically providing certain source properties.

Bonus properties include weight (normalization: sum to unity) and extension (zero by default). If deg=True on construction, then (ra,dec,extension) are converted from degrees to radians.

Methods:

__init__(*a, **kw)

Construct an Arrays collection.

__init__(*a, **kw)[source]

Construct an Arrays collection.

Parameters:

init (numpy.ndarray or object with keys() or ._dict.keys() method) – the data
names (list of str) – the column names
convert (bool) – if True, add references converting from skylab to csky naming convention

csky.utils.ensure_dir(dirname)[source]: Make sure dirname exists and is a directory.

csky.utils.get_pixmask(nside, events, radius=3)[source]

Query disc of given radius around selected events.

Parameters:

events (np.ndarray) – Array of data events
nside (int) – Resolution of the healpix map to use
radius (float) – Radius around each event which is used to select pixels

Returns:

list of pixels within ‘radius’ of the data events

Return type:

src_pix (np.ndarray)

csky.utils.get_pixval(map, ra, dec)[source]

Get the value of a healpix map at specific equatorial coordinates.

Parameters:

map (ndarray) – the healpix map
ra (float or ndarray) – the right ascension(s)
dec (float or ndarray) – the declination(s)

Returns:

the value(s)

Return type:

float or ndarray

csky.utils.get_random(seed=None)[source]: Get numpy RandomState.

csky.utils.intersection_time_intervals(intervals_1, intervals_2)[source]: Calculate intersections between two sorted lists containing continuous, non-overlapping intervals

csky.utils.sources(ra, dec, **kw)[source]: Convenience wrapper for constructing Sources instances.

csky.utils.union_continuous_intervals(intervals)[source]: Calculate the union of a given list of continuous intervals