# Licensed under a 3-clause BSD style license - see LICENSE.rst
from collections import OrderedDict
import numpy as np
import astropy.units as u
from astropy.table import Table
from astropy.time import Time
from ..spectrum.utils import SpectrumEvaluator
from ..stats.poisson import excess_error, excess_ul_helene
from ..utils.scripts import make_path
from ..utils.table import table_from_row_data
from ..stats.poisson import significance_on_off
__all__ = ["LightCurve", "LightCurveEstimator"]
[docs]class LightCurve:
"""Lightcurve container.
The lightcurve data is stored in ``table``.
For now we only support times stored in MJD format!
TODO: specification of format is work in progress
See https://github.com/open-gamma-ray-astro/gamma-astro-data-formats/pull/61
Usage: :ref:`time-lc`
Parameters
----------
table : `~astropy.table.Table`
Table with lightcurve data
"""
def __init__(self, table):
self.table = table
def __repr__(self):
return "{}(len={})".format(self.__class__.__name__, len(self.table))
@property
def time_scale(self):
"""Time scale (str).
Taken from table "TIMESYS" header.
Common values: "TT" or "UTC".
Assumed default is "UTC".
"""
return self.table.meta.get("TIMESYS", "utc")
@property
def time_format(self):
"""Time format (str)."""
return "mjd"
# @property
# def time_ref(self):
# """Time reference (`~astropy.time.Time`)."""
# return time_ref_from_dict(self.table.meta)
def _make_time(self, colname):
val = self.table[colname].data
scale = self.time_scale
format = self.time_format
return Time(val, scale=scale, format=format)
@property
def time(self):
"""Time (`~astropy.time.Time`)."""
return self._make_time("time")
@property
def time_min(self):
"""Time bin start (`~astropy.time.Time`)."""
return self._make_time("time_min")
@property
def time_max(self):
"""Time bin end (`~astropy.time.Time`)."""
return self._make_time("time_max")
@property
def time_mid(self):
"""Time bin center (`~astropy.time.Time`)."""
return self.time_min + 0.5 * self.time_delta
@property
def time_delta(self):
"""Time bin width (`~astropy.time.TimeDelta`)."""
return self.time_max - self.time_min
[docs] @classmethod
def read(cls, filename, **kwargs):
"""Read from file.
Parameters
----------
filename : str
Filename
kwargs : dict
Keyword arguments passed to `astropy.table.Table.read`.
"""
filename = make_path(filename)
table = Table.read(str(filename), **kwargs)
return cls(table=table)
[docs] def write(self, filename, **kwargs):
"""Write to file.
Parameters
----------
filename : str
Filename
kwargs : dict
Keyword arguments passed to `astropy.table.Table.write`.
"""
filename = make_path(filename)
self.table.write(str(filename), **kwargs)
[docs] def compute_fvar(self):
r"""Calculate the fractional excess variance.
This method accesses the the ``FLUX`` and ``FLUX_ERR`` columns
from the lightcurve data.
The fractional excess variance :math:`F_{var}`, an intrinsic
variability estimator, is given by
.. math::
F_{var} = \sqrt{\frac{S^{2} - \bar{\sigma^{2}}}{\bar{x}^{2}}}.
It is the excess variance after accounting for the measurement errors
on the light curve :math:`\sigma`. :math:`S` is the variance.
Returns
-------
fvar, fvar_err : `~numpy.ndarray`
Fractional excess variance.
References
----------
.. [Vaughan2003] "On characterizing the variability properties of X-ray light
curves from active galaxies", Vaughan et al. (2003)
http://adsabs.harvard.edu/abs/2003MNRAS.345.1271V
"""
flux = self.table["flux"].data.astype("float64")
flux_err = self.table["flux_err"].data.astype("float64")
flux_mean = np.mean(flux)
n_points = len(flux)
s_square = np.sum((flux - flux_mean) ** 2) / (n_points - 1)
sig_square = np.nansum(flux_err ** 2) / n_points
fvar = np.sqrt(np.abs(s_square - sig_square)) / flux_mean
sigxserr_a = np.sqrt(2 / n_points) * (sig_square / flux_mean) ** 2
sigxserr_b = np.sqrt(sig_square / n_points) * (2 * fvar / flux_mean)
sigxserr = np.sqrt(sigxserr_a ** 2 + sigxserr_b ** 2)
fvar_err = sigxserr / (2 * fvar)
return fvar, fvar_err
[docs] def compute_chisq(self):
"""Calculate the chi-square test for `LightCurve`.
Chisquare test is a variability estimator. It computes
deviations from the expected value here mean value
Returns
-------
ChiSq, P-value : tuple of float or `~numpy.ndarray`
Tuple of Chi-square and P-value
"""
import scipy.stats as stats
flux = self.table["flux"]
yexp = np.mean(flux)
yobs = flux.data
chi2, pval = stats.chisquare(yobs, yexp)
return chi2, pval
[docs] def plot(self, ax=None, time_format="mjd", flux_unit="cm-2 s-1", **kwargs):
"""Plot flux points.
Parameters
----------
ax : `~matplotlib.axes.Axes`, optional.
The `~matplotlib.axes.Axes` object to be drawn on.
If None, uses the current `~matplotlib.axes.Axes`.
time_format : {'mjd', 'iso'}, optional
If 'iso', the x axis will contain Matplotlib dates.
For formatting these dates see: https://matplotlib.org/gallery/ticks_and_spines/date_demo_rrule.html
flux_unit : str, `~astropy.units.Unit`, optional
Unit of the flux axis
kwargs : dict
Keyword arguments passed to :func:`matplotlib.pyplot.errorbar`
Returns
-------
ax : `~matplotlib.axes.Axes`
Axis object
"""
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter
if ax is None:
ax = plt.gca()
y, yerr = self._get_fluxes_and_errors(flux_unit)
is_ul, yul = self._get_flux_uls(flux_unit)
x, xerr = self._get_times_and_errors(time_format)
# length of the ul arrow
ul_arr = (
np.nanmax(np.concatenate((y[~is_ul], yul[is_ul])))
- np.nanmin(np.concatenate((y[~is_ul], yul[is_ul])))
) * 0.1
# join fluxes and upper limits for the plot
y[is_ul] = yul[is_ul]
yerr[0][is_ul] = ul_arr
# set plotting defaults and plot
kwargs.setdefault("marker", "+")
kwargs.setdefault("ls", "None")
ax.errorbar(x=x, y=y, xerr=xerr, yerr=yerr, uplims=is_ul, **kwargs)
ax.set_xlabel("Time ({})".format(time_format.upper()))
ax.set_ylabel("Flux ({:FITS})".format(u.Unit(flux_unit)))
if time_format == "iso":
ax.xaxis.set_major_formatter(DateFormatter("%Y-%m-%d %H:%M:%S"))
plt.setp(
ax.xaxis.get_majorticklabels(),
rotation=30,
ha="right",
rotation_mode="anchor",
)
return ax
def _get_fluxes_and_errors(self, unit="cm-2 s-1"):
"""Extract fluxes and corresponding errors
Helper function for the plot method.
Parameters
----------
unit : str, `~astropy.units.Unit`, optional
Unit of the returned flux and errors values
Returns
-------
y : `numpy.ndarray`
Flux values
(yn, yp) : tuple of `numpy.ndarray`
Flux error values
"""
y = self.table["flux"].quantity.to(unit)
if all(k in self.table.colnames for k in ["flux_errp", "flux_errn"]):
yp = self.table["flux_errp"].quantity.to(unit)
yn = self.table["flux_errn"].quantity.to(unit)
elif "flux_err" in self.table.colnames:
yp = self.table["flux_err"].quantity.to(unit)
yn = self.table["flux_err"].quantity.to(unit)
else:
yp, yn = np.zeros_like(y), np.zeros_like(y)
return y.value, (yn.value, yp.value)
def _get_flux_uls(self, unit="cm-2 s-1"):
"""Extract flux upper limits
Helper function for the plot method.
Parameters
----------
unit : str, `~astropy.units.Unit`, optional
Unit of the returned flux upper limit values
Returns
-------
is_ul : `numpy.ndarray`
Is flux point is an upper limit? (boolean array)
yul : `numpy.ndarray`
Flux upper limit values
"""
try:
is_ul = self.table["is_ul"].data.astype("bool")
except KeyError:
is_ul = np.zeros_like(self.table["flux"]).data.astype("bool")
if is_ul.any():
yul = self.table["flux_ul"].quantity.to(unit)
else:
yul = np.zeros_like(self.table["flux"]).quantity
yul[:] = np.nan
return is_ul, yul.value
def _get_times_and_errors(self, time_format="mjd"):
"""Extract times and corresponding errors
Helper function for the plot method.
Parameters
----------
time_format : {'mjd', 'iso'}, optional
Time format of the times. If 'iso', times and errors will be
returned as `~datetime.datetime` and `~datetime.timedelta` objects
Returns
-------
x : `~numpy.ndarray` or of `~datetime.datetime`
Time values or `~datetime.datetime` instances if 'iso' is chosen
as time format
(xn, xp) : tuple of `numpy.ndarray` of `~datetime.timedelta`
Tuple of time error values or `~datetime.timedelta` instances if
'iso' is chosen as time format
"""
from matplotlib.dates import num2timedelta
try:
x = self.time
except KeyError:
x = self.time_mid
try:
xn, xp = x - self.time_min, self.time_max - x
except KeyError:
xn, xp = x - x, x - x
# convert to given time format
if time_format == "iso":
x = x.datetime
# TODO: In astropy version >3.1 the TimeDelta.to_datetime() method
# should start working, for now i will use num2timedelta.
xn = np.asarray(num2timedelta(xn.to("d").value))
xp = np.asarray(num2timedelta(xp.to("d").value))
elif time_format == "mjd":
x = x.mjd
xn, xp = xn.to("d").value, xp.to("d").value
else:
raise ValueError("Invalid time_format: {}".format(time_format))
return x, (xn, xp)
[docs]class LightCurveEstimator:
"""Light curve estimator.
For a usage example see :gp-notebook:`light_curve`.
Parameters
----------
spec_extract : `~gammapy.spectrum.SpectrumExtraction`
Contains statistics, IRF and event lists
"""
def __init__(self, spec_extract):
self.observations = spec_extract.observations
self.obs_spec = spec_extract.spectrum_observations
self.off_evt_list = self._get_off_evt_list(spec_extract)
self.on_evt_list = self._get_on_evt_list(spec_extract)
@staticmethod
def _get_off_evt_list(spec_extract):
"""
Returns list of OFF events for each observations
"""
off_evt_list = []
for bg in spec_extract.bkg_estimate:
off_evt_list.append(bg.off_events)
return off_evt_list
@staticmethod
def _get_on_evt_list(spec_extract):
"""
Returns list of ON events for each observations
"""
on_evt_list = []
for obs in spec_extract.bkg_estimate:
on_evt_list.append(obs.on_events)
return on_evt_list
[docs] @staticmethod
def make_time_intervals_fixes(time_step, spectrum_extraction):
"""Create time intervals of fixed size.
Parameters
----------
time_step : float
Size of the light curve bins in seconds
spectrum_extraction : `~gammapy.spectrum.SpectrumExtraction`
Contains statistics, IRF and event lists
Returns
-------
table : `~astropy.table.Table`
Table of time intervals
Examples
--------
To extract intervals for light curve::
intervals = list(zip(table['t_start'], table['t_stop']))
"""
rows = []
time_start = Time(100000, format="mjd")
time_end = Time(0, format="mjd")
time_step = time_step / (24 * 3600)
for obs in spectrum_extraction.observations:
time_events = obs.events.time
if time_start > time_events.min():
time_start = time_events.min()
if time_end < time_events.max():
time_end = time_events.max()
time = time_start.value
while time < time_end.value:
time += time_step
rows.append(
dict(
t_start=Time(time - time_step, format="mjd", scale="tt"),
t_stop=Time(time, format="mjd", scale="tt"),
)
)
return Table(rows=rows)
def _create_and_filter_onofflists(
self, t_index, energy_range=None, interval=None, extra=False
):
"""Extract on and off events list from an observation and apply energy and time filters
Helper function for compute_flux_point and make_time_intervals_min_significance
Parameters
----------
t_index : int
index in self of the observation to use
energy_range : `~astropy.units.Quantity`
True energy range to filter the events
interval : `~astropy.time.Time`
Time interval (2-element array)
extra : boolean
Define if we want spec and e_reco in output
Returns
-------
on : `~gammapy.data.EventList`
List of on events
off : `~gammapy.data.EventList`
List of on events
"""
spec = self.obs_spec[t_index]
# get ON and OFF evt list
off = self.off_evt_list[t_index]
on = self.on_evt_list[t_index]
# introduce the e_reco binning here
e_reco = spec.counts.energy.edges
if energy_range is not None:
emin = e_reco[e_reco.searchsorted(max(spec.lo_threshold, energy_range[0]))]
emax = e_reco[
e_reco.searchsorted(min(spec.hi_threshold, energy_range[1])) - 1
]
# filter the event list with energy
on = on.select_energy([emin, emax])
on = on.select_energy(energy_range)
off = off.select_energy([emin, emax])
off = off.select_energy(energy_range)
if interval is not None:
# filter the event list with time
tmin = interval[0]
tmax = interval[1]
on = on.select_time([tmin, tmax])
off = off.select_time([tmin, tmax])
if extra:
return on, off, spec, e_reco
return on, off
@staticmethod
def _alpha(time_holder, obs_properties, n, istart, i):
"""Helper function for make_time_intervals_min_significance.
Parameters
----------
time_holder : list of float and flag
Contains a list of a time and a flag in 2-element arrays
obs_properties : `~astropy.table.Table`
Contains the dead time fraction and ratio of the on/off region
n : int
First observation to use
istart : int
index of the first event of the interval
i : int
index of the last event of the interval
"""
def in_list(item, L):
o, j = np.where(L == item)
for index in o:
yield index
alpha = 0
tmp = 0
time = 0
xm1 = istart
# loop over observations
for x in in_list("end", time_holder[istart : i + 1]):
if tmp == 0:
alpha += (
(1 - obs_properties["deadtime"][n])
* (
float(time_holder[x][0])
- (float(time_holder[xm1][0]) + float(time_holder[xm1 - 1][0]))
/ 2
)
* obs_properties["A_off"][n + tmp]
)
time += (1 - obs_properties["deadtime"][n]) * (
float(time_holder[x][0])
- (float(time_holder[xm1][0]) + float(time_holder[xm1 - 1][0])) / 2
)
xm1 = x + 1
tmp += 1
else:
alpha += (
(1 - obs_properties["deadtime"][n + tmp])
* (float(time_holder[x][0]) - float(time_holder[xm1][0]))
* obs_properties["A_off"][n + tmp]
)
time += (1 - obs_properties["deadtime"][n + tmp]) * (
float(time_holder[x][0]) - float(time_holder[xm1][0])
)
xm1 = x + 1
tmp += 1
alpha += (
(1 - obs_properties["deadtime"][n + tmp])
* (
(float(time_holder[i][0]) + float(time_holder[i][0])) / 2
- float(time_holder[xm1][0])
)
* obs_properties["A_off"][n + tmp]
)
time += (1 - obs_properties["deadtime"][n + tmp]) * (
(float(time_holder[i][0]) + float(time_holder[i][0])) / 2
- float(time_holder[xm1][0])
)
alpha = time / alpha
return alpha
[docs] def make_time_intervals_min_significance(
self,
significance,
significance_method,
energy_range,
spectrum_extraction,
separators=None,
):
"""
Create time intervals such that each bin of a light curve reach a given significance
The function work event by event to create an interval containing enough statistic and then starting a new one
Parameters
----------
significance : float
Target significance for each light curve point
significance_method : {'lima', 'simple'}
Significance method (see `~gammapy.stats.significance_on_off`)
energy_range : `~astropy.units.Quantity`
True energy range to evaluate integrated flux (true energy)
spectrum_extraction : `~gammapy.spectrum.SpectrumExtraction`
Contains statistics, IRF and event lists
separators : list of `~astropy.time.Time`
Contains a list of time to stop the current point creation (not saved) and start a new one
Mostly useful between observations separated by a large time gap
Returns
-------
table : `~astropy.table.Table`
Table of time intervals and information about their content : on/off events, alpha, significance
Examples
--------
extract intervals for light curve :
intervals = list(zip(table['t_start'], table['t_stop']))
"""
# The function create a list of time associated with identifiers called time_holder.
# The identifiers can be 'on' for the on events, 'off' for the off events, 'start' for the start of an
# observation, 'end for the end of an observation and 'break' for a separator.
# The function then loop other all the elements of this list
time_holder = []
obs_properties = []
n_obs = 0
# extract the separators
if separators is not None:
for time in separators:
time_holder.append([time.tt.mjd, "break"])
# recovers the starting and ending time of each observations and useful properties
for obs in spectrum_extraction.observations:
time_holder.append([obs.events.observation_time_start.tt.mjd, "start"])
time_holder.append([obs.events.observation_time_end.tt.mjd, "end"])
obs_properties.append(
dict(
deadtime=obs.observation_dead_time_fraction,
A_off=spectrum_extraction.bkg_estimate[n_obs].a_off,
)
)
n_obs += 1
obs_properties = Table(rows=obs_properties)
# prepare the on and off photon list as in the flux point computation -> should be updated accordingly
for t_index, obs in enumerate(self.observations):
on, off = self._create_and_filter_onofflists(
t_index=t_index, energy_range=energy_range
)
for time in on.time.tt.mjd:
time_holder.append([time, "on"])
for time in off.time.tt.mjd:
time_holder.append([time, "off"])
# sort all elements in the table by time
time_holder = sorted(time_holder, key=lambda item: item[0])
time_holder = np.asarray(time_holder)
rows = []
istart = 1
i = 1
n = 0
while time_holder[i][0] < time_holder[-1][0]:
i += 1
if time_holder[i][1] == "break":
while time_holder[i + 1][1] != "on" and time_holder[i + 1][1] != "off":
i += 1
n += np.sum(time_holder[istart:i] == "end")
istart = i
continue
if time_holder[i][1] != "on" and time_holder[i][1] != "off":
continue
# compute alpha
alpha = self._alpha(time_holder, obs_properties, n, istart, i)
non = np.sum(time_holder[istart : i + 1] == "on")
noff = np.sum(time_holder[istart : i + 1] == "off")
# check the significance
signif = significance_on_off(non, noff, alpha, method=significance_method)
if signif > significance:
rows.append(
dict(
t_start=(
float(time_holder[istart - 1][0])
+ float(time_holder[istart][0])
)
/ 2,
t_stop=(float(time_holder[i][0]) + float(time_holder[i + 1][0]))
/ 2,
n_on=non,
n_off=noff,
alpha=alpha,
significance=signif,
)
)
# start the next interval
while (
time_holder[i + 1][0] < time_holder[-1][0]
and time_holder[i + 1][1] != "on"
and time_holder[i + 1][1] != "off"
):
i += 1
n += np.sum(time_holder[istart : i + 1] == "end")
istart = i + 1
i = istart
table = Table(rows=rows)
table["t_start"] = Time(table["t_start"], format="mjd", scale="tt")
table["t_stop"] = Time(table["t_stop"], format="mjd", scale="tt")
return table
[docs] def light_curve(
self, time_intervals, spectral_model, energy_range, ul_significance=3
):
"""Compute light curve.
Implementation follows what is done in:
http://adsabs.harvard.edu/abs/2010A%26A...520A..83H.
To be discussed: assumption that threshold energy in the
same in reco and true energy.
Parameters
----------
time_intervals : list of `~astropy.time.Time`
List of time intervals
spectral_model : `~gammapy.spectrum.models.SpectralModel`
Spectral model
energy_range : `~astropy.units.Quantity`
True energy range to evaluate integrated flux (true energy)
ul_significance : float
Upper limit confidence level significance
Returns
-------
lc : `~gammapy.time.LightCurve`
Light curve
"""
rows = []
for time_interval in time_intervals:
useinterval, row = self.compute_flux_point(
time_interval, spectral_model, energy_range, ul_significance
)
if useinterval:
rows.append(row)
table = table_from_row_data(rows)
# TODO: check why stripping of the units here is needed
table["time_min"] = table["time_min"].data
table["time_max"] = table["time_max"].data
return LightCurve(table)
[docs] def compute_flux_point(
self, time_interval, spectral_model, energy_range, ul_significance=3
):
"""Compute one flux point for one time interval.
Parameters
----------
time_interval : `~astropy.time.Time`
Time interval (2-element array, or a tuple of Time objects)
spectral_model : `~gammapy.spectrum.models.SpectralModel`
Spectral model
energy_range : `~astropy.units.Quantity`
True energy range to evaluate integrated flux (true energy)
ul_significance : float
Upper limit confidence level significance
Returns
-------
useinterval : bool
Is True if the time_interval produce a valid flux point
measurements : dict
Dictionary with flux point measurement in the time interval
"""
tmin, tmax = time_interval[0], time_interval[1]
livetime = 0
alpha_mean = 0.0
alpha_mean_backup = 0.0
measured_excess = 0
predicted_excess = 0
n_on = 0
n_off = 0
useinterval = False
# Loop on observations
for t_index, obs in enumerate(self.observations):
# discard observations not matching the time interval
obs_start = obs.events.observation_time_start
obs_stop = obs.events.observation_time_end
if (tmin < obs_start and tmax < obs_start) or (tmin > obs_stop):
continue
useinterval = True
# get ON and OFF evt list
on, off, spec, e_reco = self._create_and_filter_onofflists(
t_index=t_index,
energy_range=energy_range,
interval=[tmin, tmax],
extra=True,
)
n_on_obs = len(on.table)
n_off_obs = len(off.table)
# compute effective livetime (for the interval)
if tmin >= obs_start and tmax <= obs_stop:
# interval included in obs
livetime_to_add = (tmax - tmin).to("s")
elif tmin >= obs_start and tmax >= obs_stop:
# interval min above tstart from obs
livetime_to_add = (obs_stop - tmin).to("s")
elif tmin <= obs_start and tmax <= obs_stop:
# interval min below tstart from obs
livetime_to_add = (tmax - obs_start).to("s")
elif tmin <= obs_start and tmax >= obs_stop:
# obs included in interval
livetime_to_add = (obs_stop - obs_start).to("s")
else:
livetime_to_add = u.Quantity(0, "s")
# Take into account dead time
livetime_to_add *= 1.0 - obs.observation_dead_time_fraction
# Compute excess
obs_measured_excess = n_on_obs - spec.alpha * n_off_obs
# Compute the expected excess in the range given by the user
# but must respect the energy threshold of the observation
# (to match the energy range of the measured excess)
# We use the effective livetime and the right energy threshold
e_idx = np.where(
np.logical_and.reduce(
(
e_reco >= spec.lo_threshold, # threshold
e_reco <= spec.hi_threshold, # threshold
e_reco >= energy_range[0], # user
e_reco <= energy_range[-1],
) # user
)
)[0]
counts_predictor = SpectrumEvaluator(
livetime=livetime_to_add,
aeff=spec.aeff,
edisp=spec.edisp,
model=spectral_model,
)
counts_predicted_excess = counts_predictor.compute_npred().data.data[
e_idx[:-1]
]
obs_predicted_excess = np.sum(counts_predicted_excess)
# compute effective normalisation between ON/OFF (for the interval)
livetime += livetime_to_add
alpha_mean += spec.alpha * n_off_obs
alpha_mean_backup += spec.alpha * livetime_to_add
measured_excess += obs_measured_excess
predicted_excess += obs_predicted_excess
n_on += n_on_obs
n_off += n_off_obs
# Fill time interval information
if useinterval:
int_flux = spectral_model.integral(energy_range[0], energy_range[1])
if n_off > 0.0:
alpha_mean /= n_off
if livetime > 0.0:
alpha_mean_backup /= livetime
if alpha_mean == 0.0: # use backup if necessary
alpha_mean = alpha_mean_backup
flux = measured_excess / predicted_excess.value
flux *= int_flux
# Gaussian errors, TODO: should be improved
flux_err = int_flux / predicted_excess.value
delta_excess = excess_error(n_on=n_on, n_off=n_off, alpha=alpha_mean)
flux_err *= delta_excess
sigma = significance_on_off(
n_on=n_on, n_off=n_off, alpha=alpha_mean, method="lima"
)
is_ul = sigma <= 3
flux_ul = excess_ul_helene(
n_on - alpha_mean * n_off, delta_excess, ul_significance
)
flux_ul *= int_flux / predicted_excess.value
else:
flux = u.Quantity(0, "cm-2 s-1")
flux_err = u.Quantity(0, "cm-2 s-1")
flux_ul = u.Quantity(-1, "cm-2 s-1")
# Store measurements in a dict and return that
return (
useinterval,
OrderedDict(
[
# TODO: check why time is stored without units
("time_min", Time(tmin, format="mjd").value),
("time_max", Time(tmax, format="mjd").value),
("flux", flux),
("flux_err", flux_err),
("flux_ul", flux_ul),
("is_ul", is_ul),
("livetime", livetime),
("alpha", alpha_mean),
("n_on", n_on),
("n_off", n_off),
("measured_excess", measured_excess),
("expected_excess", predicted_excess),
]
),
)