Tutorial: fitting a FSRQ broad-band SED using angpy and Gammapy

In order to perform a fit of the broad-band SED of a jetted AGN, agnpy includes a gammapy wrapper. The wrapper defines a custom SpectralModel, representing the emission due to a combination of radiative processes. The SpectralModel can be used either to fit flux points or to perform a forward-folding likelihood fit (if the instrument response is available in a format compatible with gammapy).

Several combination of radiative processes can be considered to model the broad-band emission of jetted AGN. For simplicity, we provide wrappers for the two scenarios most-commonly considered:

  • SycnhrotronSelfComptonModel, representing the sum of synchrotron and synchrotron self Compton (SSC) radiation. This scenario is commonly considered to model BL Lac sources;

  • ExternalComptonModel, representing the sum of synchrotron and synchrotron self Compton radiation along with an external Compton (EC) component. EC scattering can be computed considering a list of target photon fields. This scenario is commonly considered to model flat spectrum radio quasars (FSRQs).

In this tutorial we will show how to use ExternalComptonSpectralModel to perform the fit of the broad-band SED of PKS1510-089, measured during its gamma-ray flaring activity in 2015 (Ahnen et al. 2017). We select the MWL SED corresponding to the period identified in the paper as “Period B” (MJD 57164-57166).

Gammapy is required to run this notebook.

[1]:
# import numpy, astropy and matplotlib for basic functionalities
import numpy as np
import astropy.units as u
from astropy.constants import c, G, M_sun
import matplotlib.pyplot as plt
import pkg_resources

# import agnpy classes
from agnpy.spectra import LogParabola
from agnpy.fit import ExternalComptonModel, load_gammapy_flux_points
from agnpy.utils.plot import load_mpl_rc, sed_y_label

load_mpl_rc()

# import gammapy classes
from gammapy.modeling.models import SkyModel
from gammapy.modeling import Fit

The ExternalComptonModel wraps the agnpy functions to compute synchrotron, SSC, and EC radiation. It returns a gammapy.modeling.SpectralModel. To initialise the model, the electron distribution and a list of targets for EC has to be specified, the remaining parameters (the ones of the emission region and of the line / thermal emitters) will be initialised automatically and can be modified at a later stage.

The ExternalComptonModel class provides both the sherpa and gammapy wrappers. You should specify, through the backend argument, which package you want to use.

In this case we choose to consider only EC on the DT radiation.

[2]:
# electron energy distribution
n_e = LogParabola(
    k=1 * u.Unit("cm-3"), p=2.0, q=0.2, gamma_0=1e2, gamma_min=1, gamma_max=3e4
)

# initialise the sherpa model, consider only EC on the DT radiation
ec_model = ExternalComptonModel(n_e, ["dt"], ssa=True, backend="gammapy")

Let us set appropriate parameters for the emission region. The size of the blob, \(R_{\rm b}\), is set by the variability timescale, \(t_{\rm var}\), via

\begin{equation} R_{\rm b} = \frac{c \delta_{\rm D} t_{\rm var}}{1 + z}, \end{equation}

where \(c\) is the speed of light, \(\delta_{\rm D}\) the Doppler factor, and \(z\) the redshift.

The parameters of the accretion disk and the DT are taken from (Ahnen et al. 2017).

[3]:
z = 0.361

# - blob parameters
Gamma = 20
delta_D = 25
Beta = np.sqrt(1 - 1 / np.power(Gamma, 2))  # jet relativistic speed
mu_s = (1 - 1 / (Gamma * delta_D)) / Beta  # viewing angle
B = 0.35 * u.G
# size and location of the emission region
t_var = 0.5 * u.d
r = 6e17 * u.cm

# - disk parameters
L_disk = 6.7e45 * u.Unit("erg s-1")  # disk luminosity
M_BH = 5.71 * 1e7 * M_sun
eta = 1 / 12
m_dot = (L_disk / (eta * c ** 2)).to("g s-1")
R_g = ((G * M_BH) / c ** 2).to("cm")
R_in = 6 * R_g
R_out = 3e4 * R_g

# - DT parameters
xi_dt = 0.6
T_dt = 2e3 * u.K
# Ghisellini and Tavecchio's Formula
# relating the DT radius to the disk luminosity (check the reference)
R_dt = 2.5 * 1e18 * np.sqrt(L_disk.to_value("erg s-1") / 1e45) * u.cm
[4]:
# emission regions parameters
ec_model.parameters["z"].value = z
ec_model.parameters["delta_D"].value = delta_D
ec_model.parameters["log10_B"].value = np.log10(0.35)
ec_model.parameters["mu_s"].value = mu_s
ec_model.parameters["t_var"].value = t_var.to_value("s")
ec_model.parameters["t_var"].frozen = True
ec_model.parameters["log10_r"].value = np.log10(r.to_value("cm"))

# we freeze the reference energy of the log parabola
ec_model.parameters["log10_gamma_0"].frozen = True

# target parameters
# - disk
ec_model.parameters["log10_L_disk"].value = np.log10(L_disk.to_value("erg s-1"))
ec_model.parameters["M_BH"].value = M_BH.to_value("g")
ec_model.parameters["m_dot"].value = m_dot.to_value("g s-1")
ec_model.parameters["R_in"].value = R_in.to_value("cm")
ec_model.parameters["R_out"].value = R_out.to_value("cm")
# - DT
ec_model.parameters["xi_dt"].value = xi_dt
ec_model.parameters["T_dt"].value = T_dt.to_value("K")
ec_model.parameters["R_dt"].value = R_dt.to_value("cm")

We can display all the parameters of the model at once

[5]:
ec_model.parameters.to_table()
[5]:
Table length=21
typenamevalueuniterrorminmaxfrozenis_normlink
str8str15float64str5int64float64float64boolboolstr1
spectrallog10_k0.0000e+000.000e+00-1.000e+011.000e+01FalseFalse
spectralp2.0000e+000.000e+001.000e+005.000e+00FalseFalse
spectralq2.0000e-010.000e+001.000e-031.000e+00FalseFalse
spectrallog10_gamma_02.0000e+000.000e+002.000e+006.000e+00TrueFalse
spectrallog10_gamma_min0.0000e+000.000e+000.000e+004.000e+00TrueFalse
spectrallog10_gamma_max4.4771e+000.000e+004.000e+008.000e+00TrueFalse
spectralz3.6100e-010.000e+001.000e-031.000e+01TrueFalse
spectraldelta_D2.5000e+010.000e+001.000e+001.000e+02FalseFalse
spectrallog10_B-4.5593e-010.000e+00-4.000e+002.000e+00FalseFalse
spectralt_var4.3200e+04s0.000e+001.000e+013.142e+07TrueFalse
spectralmu_s9.9925e-010.000e+000.000e+001.000e+00TrueFalse
spectrallog10_r1.7778e+010.000e+001.600e+012.200e+01TrueFalse
spectrallog10_L_disk4.5826e+010.000e+004.200e+014.800e+01TrueFalse
spectralM_BH1.1354e+41g0.000e+001.000e+321.000e+45TrueFalse
spectralm_dot8.9457e+25g s-10.000e+001.000e+241.000e+30TrueFalse
spectralR_in5.0589e+13cm0.000e+001.000e+121.000e+16TrueFalse
spectralR_out2.5295e+17cm0.000e+001.000e+121.000e+19TrueFalse
spectralxi_dt6.0000e-010.000e+000.000e+001.000e+00TrueFalse
spectralT_dt2.0000e+03K0.000e+001.000e+021.000e+04TrueFalse
spectralR_dt6.4711e+18cm0.000e+001.000e+171.000e+20TrueFalse
spectralnorm1.0000e+000.000e+001.000e-011.000e+01TrueTrue

or display separately the parameters related to the electrons energy distribution…

[6]:
ec_model.spectral_parameters.to_table()
[6]:
Table length=6
typenamevalueuniterrorminmaxfrozenis_normlink
str8str15float64str1int64float64float64boolboolstr1
spectrallog10_k0.0000e+000.000e+00-1.000e+011.000e+01FalseFalse
spectralp2.0000e+000.000e+001.000e+005.000e+00FalseFalse
spectralq2.0000e-010.000e+001.000e-031.000e+00FalseFalse
spectrallog10_gamma_02.0000e+000.000e+002.000e+006.000e+00TrueFalse
spectrallog10_gamma_min0.0000e+000.000e+000.000e+004.000e+00TrueFalse
spectrallog10_gamma_max4.4771e+000.000e+004.000e+008.000e+00TrueFalse

and the parameters related to the emission region, the blob in this case

[7]:
ec_model.emission_region_parameters.to_table()
[7]:
Table length=6
typenamevalueuniterrorminmaxfrozenis_normlink
str8str7float64str1int64float64float64boolboolstr1
spectralz3.6100e-010.000e+001.000e-031.000e+01TrueFalse
spectraldelta_D2.5000e+010.000e+001.000e+001.000e+02FalseFalse
spectrallog10_B-4.5593e-010.000e+00-4.000e+002.000e+00FalseFalse
spectralt_var4.3200e+04s0.000e+001.000e+013.142e+07TrueFalse
spectralmu_s9.9925e-010.000e+000.000e+001.000e+00TrueFalse
spectrallog10_r1.7778e+010.000e+001.600e+012.200e+01TrueFalse

in case of the ExternalComptonSpectralModel, we can display also the parameters of the targets (i.e. of the line and thermal emitters producing the photon targets for EC scattering.

[8]:
ec_model.targets_parameters.to_table()
[8]:
Table length=8
typenamevalueuniterrorminmaxfrozenis_normlink
str8str12float64str5int64float64float64boolboolstr1
spectrallog10_L_disk4.5826e+010.000e+004.200e+014.800e+01TrueFalse
spectralM_BH1.1354e+41g0.000e+001.000e+321.000e+45TrueFalse
spectralm_dot8.9457e+25g s-10.000e+001.000e+241.000e+30TrueFalse
spectralR_in5.0589e+13cm0.000e+001.000e+121.000e+16TrueFalse
spectralR_out2.5295e+17cm0.000e+001.000e+121.000e+19TrueFalse
spectralxi_dt6.0000e-010.000e+000.000e+001.000e+00TrueFalse
spectralT_dt2.0000e+03K0.000e+001.000e+021.000e+04TrueFalse
spectralR_dt6.4711e+18cm0.000e+001.000e+171.000e+20TrueFalse

Fit with gammapy

Here we start the procedure to fit with Gammapy.

1) load the MWL flux points, add systematics

A function is provided in agnpy.fit to directly load flux points in a gammapy.estimators.FluxPoints object. It reads the data from a file, included in the package, containing a MWL SED following these specifications.

The same function allows to add a systematic error on the flux points, this is done with a dictionary specifying the instrument name and the systematic error, expressed as a relative error on the flux. The systematic error is summed in quadrature to the statistical error.

In this example, we use a very rough and conservative estimate of the systematic errors (\(30\%\) of the flux for VHE instruments, \(10\%\) for HE and X-ray instruments, \(5\%\) for all the other instruments).

Specifying the systematic errors through the dictionary is optional.

We can also set the minimum and maximum energy to be used in the fit. We exclude points below \(10^{11}\,{\rm Hz}\), as they are measured in the radio band with large integration regions. They hence include the extended emission of the jet, while in our model we are considering the emission from a finite region of the jet, the blob.

[9]:
systematics_dict = {
    "Fermi-LAT": 0.10,
    "KVA1": 0.05,
    "KVA2": 0.05,
    "MAGIC": 0.30,
    "Metsahovi": 0.05,
    "NICS1": 0.05,
    "NICS2": 0.05,
    "SMART1": 0.05,
    "Swift-XRT": 0.10,
    "TCS1": 0.05,
    "TCS2": 0.05,
    "TCS3": 0.05,
    "TCS4": 0.05,
    "TCS5": 0.05,
    "TCS6": 0.05,
    "UVOT": 0.05,
}

# min and max energy to be considered in the fit
E_min = (1e11 * u.Hz).to("eV", equivalencies=u.spectral())
E_max = 100 * u.TeV

datasets = load_gammapy_flux_points(
    "../../agnpy/data/mwl_seds/PKS1510-089_2015b.ecsv", E_min, E_max, systematics_dict
)
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.

For the gammapy wrapper, we have to define a spectral model and associate it to the list of datasets we obtained from the flux points.

[10]:
sky_model = SkyModel(spectral_model=ec_model, name="PKS1510-089")
datasets.models = [sky_model]
/Users/cosimo/work/agnpy/agnpy/synchrotron/synchrotron.py:66: RuntimeWarning: divide by zero encountered in true_divide
  u = 1 / 2 + np.exp(-tau) / tau - (1 - np.exp(-tau)) / np.power(tau, 2)
/Users/cosimo/work/agnpy/agnpy/synchrotron/synchrotron.py:66: RuntimeWarning: invalid value encountered in true_divide
  u = 1 / 2 + np.exp(-tau) / tau - (1 - np.exp(-tau)) / np.power(tau, 2)
/Users/cosimo/software/miniconda3/envs/gammapy-0.20.1/lib/python3.8/site-packages/astropy/units/quantity.py:613: RuntimeWarning: divide by zero encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/cosimo/software/miniconda3/envs/gammapy-0.20.1/lib/python3.8/site-packages/astropy/units/quantity.py:613: RuntimeWarning: overflow encountered in expm1
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/cosimo/software/miniconda3/envs/gammapy-0.20.1/lib/python3.8/site-packages/astropy/units/quantity.py:613: RuntimeWarning: overflow encountered in multiply
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)

Let us plot all the flux points and the initial model

[11]:
fig, ax = plt.subplots(figsize=(8, 6))

for dataset in datasets:
    dataset.data.plot(ax=ax, label=dataset.name)

ec_model.plot(
    ax=ax,
    energy_bounds=[1e-6, 1e14] * u.eV,
    energy_power=2,
    label="EC model",
    color="k",
    lw=1.6,
)

ax.set_ylabel(sed_y_label)
ax.set_xlabel(r"$E\,/\,{\rm eV}$")
ax.set_xlim([1e-5, 1e12])
ax.set_ylim([1e-13, 1e-8])
ax.legend(ncol=4, fontsize=9)

plt.show()
../_images/tutorials_ec_dt_gammapy_fit_20_0.png

3) run the fit

[12]:
%%time
# define the fitter
fitter = Fit()
results = fitter.run(datasets)
/Users/cosimo/work/agnpy/agnpy/synchrotron/synchrotron.py:67: RuntimeWarning: overflow encountered in true_divide
  return np.where(tau < 1e-3, 1, 3 * u / tau)
/Users/cosimo/work/agnpy/agnpy/synchrotron/synchrotron.py:66: RuntimeWarning: divide by zero encountered in true_divide
  u = 1 / 2 + np.exp(-tau) / tau - (1 - np.exp(-tau)) / np.power(tau, 2)
/Users/cosimo/work/agnpy/agnpy/synchrotron/synchrotron.py:66: RuntimeWarning: invalid value encountered in true_divide
  u = 1 / 2 + np.exp(-tau) / tau - (1 - np.exp(-tau)) / np.power(tau, 2)
/Users/cosimo/software/miniconda3/envs/gammapy-0.20.1/lib/python3.8/site-packages/astropy/units/quantity.py:613: RuntimeWarning: divide by zero encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/cosimo/software/miniconda3/envs/gammapy-0.20.1/lib/python3.8/site-packages/astropy/units/quantity.py:613: RuntimeWarning: overflow encountered in expm1
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/cosimo/software/miniconda3/envs/gammapy-0.20.1/lib/python3.8/site-packages/astropy/units/quantity.py:613: RuntimeWarning: overflow encountered in multiply
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
CPU times: user 4min 48s, sys: 14.9 s, total: 5min 2s
Wall time: 5min 4s
[13]:
print(results)
print(ec_model.parameters.to_table())
OptimizeResult

        backend    : minuit
        method     : migrad
        success    : True
        message    : Optimization terminated successfully.
        nfev       : 317
        total stat : 220.71

CovarianceResult

        backend    : minuit
        method     : hesse
        success    : True
        message    : Hesse terminated successfully.

  type         name         value     unit   error      min        max    frozen is_norm link
-------- --------------- ----------- ----- --------- ---------- --------- ------ ------- ----
spectral         log10_k  2.4521e-01       1.577e-01 -1.000e+01 1.000e+01  False   False
spectral               p  2.1681e+00       1.243e-01  1.000e+00 5.000e+00  False   False
spectral               q  3.3617e-01       3.368e-02  1.000e-03 1.000e+00  False   False
spectral   log10_gamma_0  2.0000e+00       0.000e+00  2.000e+00 6.000e+00   True   False
spectral log10_gamma_min  0.0000e+00       0.000e+00  0.000e+00 4.000e+00   True   False
spectral log10_gamma_max  4.4771e+00       0.000e+00  4.000e+00 8.000e+00   True   False
spectral               z  3.6100e-01       0.000e+00  1.000e-03 1.000e+01   True   False
spectral         delta_D  2.5823e+01       1.788e+00  1.000e+00 1.000e+02  False   False
spectral         log10_B -4.0655e-01       3.696e-02 -4.000e+00 2.000e+00  False   False
spectral           t_var  4.3200e+04     s 0.000e+00  1.000e+01 3.142e+07   True   False
spectral            mu_s  9.9925e-01       0.000e+00  0.000e+00 1.000e+00   True   False
spectral         log10_r  1.7778e+01       0.000e+00  1.600e+01 2.200e+01   True   False
spectral    log10_L_disk  4.5826e+01       0.000e+00  4.200e+01 4.800e+01   True   False
spectral            M_BH  1.1354e+41     g 0.000e+00  1.000e+32 1.000e+45   True   False
spectral           m_dot  8.9457e+25 g s-1 0.000e+00  1.000e+24 1.000e+30   True   False
spectral            R_in  5.0589e+13    cm 0.000e+00  1.000e+12 1.000e+16   True   False
spectral           R_out  2.5295e+17    cm 0.000e+00  1.000e+12 1.000e+19   True   False
spectral           xi_dt  6.0000e-01       0.000e+00  0.000e+00 1.000e+00   True   False
spectral            T_dt  2.0000e+03     K 0.000e+00  1.000e+02 1.000e+04   True   False
spectral            R_dt  6.4711e+18    cm 0.000e+00  1.000e+17 1.000e+20   True   False
spectral            norm  1.0000e+00       0.000e+00  1.000e-01 1.000e+01   True    True

Now let us plot the final model and the flux points

[14]:
fig, ax = plt.subplots(figsize=(8, 6))

for dataset in datasets:
    dataset.data.plot(ax=ax, label=dataset.name)

ec_model.plot(
    ax=ax,
    energy_bounds=[1e-6, 1e14] * u.eV,
    energy_power=2,
    label="model",
    color="k",
    lw=1.6,
)

ax.set_ylabel(sed_y_label)
ax.set_xlabel(r"$E\,/\,{\rm eV}$")
ax.set_xlim([1e-5, 1e12])
ax.set_ylim([1e-13, 1e-8])
ax.legend(ncol=4, fontsize=9)

plt.show()
/Users/cosimo/work/agnpy/agnpy/synchrotron/synchrotron.py:66: RuntimeWarning: divide by zero encountered in true_divide
  u = 1 / 2 + np.exp(-tau) / tau - (1 - np.exp(-tau)) / np.power(tau, 2)
/Users/cosimo/work/agnpy/agnpy/synchrotron/synchrotron.py:66: RuntimeWarning: invalid value encountered in true_divide
  u = 1 / 2 + np.exp(-tau) / tau - (1 - np.exp(-tau)) / np.power(tau, 2)
/Users/cosimo/software/miniconda3/envs/gammapy-0.20.1/lib/python3.8/site-packages/astropy/units/quantity.py:613: RuntimeWarning: divide by zero encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/cosimo/software/miniconda3/envs/gammapy-0.20.1/lib/python3.8/site-packages/astropy/units/quantity.py:613: RuntimeWarning: overflow encountered in expm1
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/Users/cosimo/software/miniconda3/envs/gammapy-0.20.1/lib/python3.8/site-packages/astropy/units/quantity.py:613: RuntimeWarning: overflow encountered in multiply
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
../_images/tutorials_ec_dt_gammapy_fit_25_1.png