Retrieving Minos errors in Light-Curve analysis

[piemmegi]

Feature Request

When studying gamma-ray sources in low-counts regimes (e.g., faint sources at very high energies) it would be very useful to have access to the Minos errors on the flux points in a Light-Curve. These are the de-facto standard used in the FGL catalogs when reporting the flux points in Light-Curves corresponding to a very low TS (usually TS < 1), since this is a regime where the standard errors estimated with the quadratic approximation would significantly underestimate the real error on each point.

At the time of writing, fermipy already provides Minos errors on the points in a SED (see the docs), but not in a Light-Curve. It could be useful thus to extend the gta.lightcurve method, in order to allow a user to retrieve and save the Minos errors on the points in a Light-Curve.

A possible implementation of this feature follows what is already done for the SED, specifically in the sed.py code (L430-520):

ref_flux = self.like[name].flux(emin, emax)
(...)
lnlp = self.profile_norm(name, logemin=logemin, logemax=logemax,
                                     savestate=True, reoptimize=True,
                                     npts=npts, optimizer=config['optimizer'])
(...)
ul_data = utils.get_parameter_limits(lnlp['flux'],
                                                 lnlp['dloglike'],
                                                 cl_limit=ul_confidence)

o['norm_err_hi'][i] = ul_data['err_hi'] / ref_flux
o['norm_err_lo'][i] = ul_data['err_lo'] / ref_flux
(...)
for t in ['flux', 'eflux', 'dnde', 'e2dnde']:
    o['%s_err' % t] = o['norm_err'] * o['ref_%s' % t]
    o['%s_err_hi' % t] = o['norm_err_hi'] * o['ref_%s' % t]
    o['%s_err_lo' % t] = o['norm_err_lo'] * o['ref_%s' % t]

The idea is thus to use the the gta.lightcurve output to get the equivalent of ref_flux, lnlp['flux'] and lnlp['dloglike'], and then reproduce the same code from there. Note that the SED ref_flux is the flux obtained with the global model but integrating over the restricted energy range of a single energy bin. In the Light-Curve case, ref_flux would be simply the global flux of the source, because from time bin to time bin the basic assumption is that the same global model holds.

A snippet of code for a hypothetical implementation after the of the call to gta.lightcurve, from a users perspective, would be:

lightOut = gta.lightcurve(...)
flux_err = lightOut['flux_err']
flux_scan = lightOut['flux_scan']
dloglike_scan = lightOut['dloglike_scan']
ref_flux = gta.roi[srcName].data['flux']
ref_eflux = gta.roi[srcName].data['eflux']
ref_dnde = gta.roi[srcName].data['dnde']

flux_err_hi = np.zeros_like(flux_err)
flux_err_lo = np.zeros_like(flux_err)
eflux_err_hi = np.zeros_like(flux_err)
eflux_err_lo = np.zeros_like(flux_err)
dnde_err_hi = np.zeros_like(flux_err)
dnde_err_lo = np.zeros_like(flux_err)

for j in range(nbins):
    # This loop is probably unavoidable, because utils.get_parameter_limits() seem to accept only 1D inputs
    ul_data  = utils.get_parameter_limits(flux_scan[j,:],dloglike_scan[j,:])       
    norm_err_hi = ul_data['err_hi'] / ref_flux
    norm_err_lo = ul_data['err_lo'] / ref_flux
    flux_err_hi[j] = norm_err_hi * ref_flux
    flux_err_lo[j] = norm_err_lo * ref_flux
    eflux_err_hi[j] = norm_err_hi * ref_eflux
    eflux_err_lo[j] = norm_err_lo * ref_eflux
    dnde_err_hi[j] = norm_err_hi * ref_dnde
    dnde_err_lo[j] = norm_err_lo * ref_dnde

This is obviuously a patch solution from a user's perspective, not a direct suggestion on how to actually implement it in the fermipy source code. I'll work on how to do that and, if I put together something that can work eventually I'll open a dedicated new branch and do the PR.