{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Tutorial: fitting a BL Lac broad-band SED using angpy and sherpa\n", "\n", "In order to perform a fit of the broad-band SED of a jetted AGN, `agnpy` includes a `sherpa` wrapper.\n", "The wrapper defines a custom [model.RegriddableModel1D](https://sherpa.readthedocs.io/en/latest/model_classes/api/sherpa.models.model.RegriddableModel1D.html), representing the emission due to a combination of radiative processes. The model can be used to fit flux points.\n", "\n", "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:\n", "\n", " * `SycnhrotronSelfComptonModel`, representing the sum of synchrotron and synchrotron self Compton (SSC) radiation. This scenario is commonly considered to model BL Lac sources;\n", " \n", " * `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).\n", "\n", "In this tutorial, we will show how to use the `SynchrotronSelfComptonModel` to fit the broad-band SED of Mrk 421, measured by a MWL campaign in 2009 [(Abdo et al. 2011)](https://ui.adsabs.harvard.edu/abs/2011ApJ...736..131A/abstract).\n", "\n", "[sherpa](https://sherpa.readthedocs.io/en/latest/index.html) is required to run this notebook." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# import numpy, astropy and matplotlib for basic functionalities\n", "import numpy as np\n", "import astropy.units as u\n", "import matplotlib.pyplot as plt\n", "import pkg_resources\n", "\n", "# import agnpy classes\n", "from agnpy.spectra import BrokenPowerLaw\n", "from agnpy.fit import SynchrotronSelfComptonModel, load_sherpa_flux_points\n", "from agnpy.utils.plot import load_mpl_rc, sed_y_label\n", "\n", "load_mpl_rc()\n", "\n", "# import sherpa classes\n", "from sherpa.fit import Fit\n", "from sherpa.stats import Chi2\n", "from sherpa.optmethods import LevMar" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### `sherpa` wrapper of `agnpy`\n", "\n", "\n", "The `SynchrotronSelfComptonModel` wraps the `agnpy` functions to compute synchrotron and SSC radiation and returns a `model.RegriddableModel1D`. To initialise the model, only the electron distribution has to be specified, the remaining parameters (the ones of the emission region) will be initialised automatically and can be modified at a later stage.\n", "\n", "The `SynchrotronSelfComptonModel` class provides both the `sherpa` and `gammapy` wrappers. You should specify, through the `backend` argument, which package you want to use." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# electron energy distribution\n", "n_e = BrokenPowerLaw(\n", " k=1e-8 * u.Unit(\"cm-3\"),\n", " p1=2.02,\n", " p2=3.43,\n", " gamma_b=1e5,\n", " gamma_min=500,\n", " gamma_max=1e6,\n", ")\n", "\n", "# initialise the sherpa model\n", "ssc_model = SynchrotronSelfComptonModel(n_e, backend=\"sherpa\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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\n", "\n", "\\begin{equation}\n", "R_{\\rm b} = \\frac{c \\delta_{\\rm D} t_{\\rm var}}{1 + z},\n", "\\end{equation}\n", "\n", "where $c$ is the speed of light, $\\delta_{\\rm D}$ the Doppler factor, and $z$ the redshift." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "ssc_model.z = 0.0308\n", "ssc_model.delta_D = 18\n", "ssc_model.t_var = (1 * u.d).to_value(\"s\")\n", "ssc_model.t_var.freeze()\n", "ssc_model.log10_B = -1.3" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
| Component | Parameter | Thawed | Value | Min | Max | Units |
|---|---|---|---|---|---|---|
| ssc | log10_k | -8.0 | -10.0 | 10.0 | ||
| p1 | 2.02 | 1.0 | 5.0 | |||
| p2 | 3.43 | 1.0 | 5.0 | |||
| log10_gamma_b | 5.0 | 2.0 | 6.0 | |||
| log10_gamma_min | 2.6989700043360187 | 0.0 | 4.0 | |||
| log10_gamma_max | 6.0 | 4.0 | 8.0 | |||
| z | 0.0308 | 0.001 | 10.0 | |||
| delta_D | 18.0 | 1.0 | 100.0 | |||
| log10_B | -1.3 | -4.0 | 2.0 | |||
| t_var | 86400.0 | 10.0 | 31415926.535897933 | s |