\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.15?urlpath=lab/tree/analysis_mwl.ipynb)\n",
"- You can contribute with your own notebooks in this\n",
"[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n",
"- **Source files:**\n",
"[analysis_mwl.ipynb](../_static/notebooks/analysis_mwl.ipynb) |\n",
"[analysis_mwl.py](../_static/notebooks/analysis_mwl.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Joint modeling, fitting, and serialization\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This tutorial illustrates how to perfom a joint modeling and fitting of the Crab Nebula spectrum using different datasets.\n",
"We look at the gamma-ray emission from the Crab nebula between 10 GeV and 100 TeV.\n",
"The spectral parameters are optimized by combining a 3D analysis of Fermi-LAT data, a ON/OFF spectral analysis of HESS data, and flux points from HAWC.\n",
"\n",
"In this tutorial we are going to use pre-made datasets. We prepared maps of the Crab region as seen by Fermi-LAT using the same event selection than the [3FHL catalog](https://arxiv.org/abs/1702.00664) (7 years of data with energy from 10 GeV to 2 TeV). For the HESS ON/OFF analysis we used two observations from the [first public data release](https://arxiv.org/abs/1810.04516) with a significant signal from energy of about 600 GeV to 10 TeV. These observations have an offset of 0.5° and a zenith angle of 45-48°. The HAWC flux points data are taken from a [recent analysis](https://arxiv.org/pdf/1905.12518.pdf) based on 2.5 years of data with energy between 300 Gev and 300 TeV. \n",
"\n",
"More details on how to prepare datasets with the high and low level interfaces are available in these tutorials: \n",
"\n",
"- https://docs.gammapy.org/0.14/notebooks/fermi_lat.html\n",
"- https://docs.gammapy.org/dev/notebooks/hess.html\n",
"- https://docs.gammapy.org/dev/notebooks/spectrum_analysis.html\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from astropy import units as u\n",
"import matplotlib.pyplot as plt\n",
"from gammapy.modeling import Fit, Datasets\n",
"from gammapy.spectrum import (\n",
" FluxPoints,\n",
" FluxPointsEstimator,\n",
" FluxPointsDataset,\n",
" SpectrumDatasetOnOff,\n",
")\n",
"from gammapy.maps import MapAxis\n",
"from pathlib import Path"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Data and models files\n",
"\n",
"\n",
"The datasets serialization produce YAML files listing the datasets and models. In the following cells we show an example containning only the Fermi-LAT dataset and the Crab model. \n",
"\n",
"Fermi-LAT-3FHL_datasets.yaml:\n",
"\n",
"```yaml\n",
"datasets:\n",
"- name: Fermi-LAT\n",
" type: MapDataset\n",
" likelihood: cash\n",
" models:\n",
"- Crab Nebula\n",
" background: background\n",
" filename: $GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_data_Fermi-LAT.fits\n",
"```\n",
"\n",
"We used as model a point source with a log-parabola spectrum. The initial parameters were taken from the latest Fermi-LAT catalog [4FGL](https://arxiv.org/abs/1902.10045), then we have re-optimized the spectral parameters for our dataset in the 10 GeV - 2 TeV energy range (fixing the source position).\n",
"\n",
"Fermi-LAT-3FHL_models.yaml:\n",
"\n",
"```yaml\n",
"components:\n",
"- name: Crab Nebula\n",
" type: SkyModel\n",
" spatial:\n",
" type: PointSpatialModel\n",
" frame: icrs\n",
" parameters:\n",
" - name: lon_0\n",
" value: 83.63310241699219\n",
" unit: deg\n",
" min: .nan\n",
" max: .nan\n",
" frozen: true\n",
" - name: lat_0\n",
" value: 22.019899368286133\n",
" unit: deg\n",
" min: -90.0\n",
" max: 90.0\n",
" frozen: true\n",
" spectral:\n",
" type: LogParabolaSpectralModel\n",
" parameters:\n",
" - name: amplitude\n",
" value: 0.3415498620816483\n",
" unit: cm-2 s-1 TeV-1\n",
" min: .nan\n",
" max: .nan\n",
" frozen: false\n",
" - name: reference\n",
" value: 5.054833602905273e-05\n",
" unit: TeV\n",
" min: .nan\n",
" max: .nan\n",
" frozen: true\n",
" - name: alpha\n",
" value: 2.510798031388936\n",
" unit: ''\n",
" min: .nan\n",
" max: .nan\n",
" frozen: false\n",
" - name: beta\n",
" value: -0.022476498188855533\n",
" unit: ''\n",
" min: .nan\n",
" max: .nan\n",
" frozen: false\n",
"- name: background\n",
" type: BackgroundModel\n",
" parameters:\n",
" - name: norm\n",
" value: 0.9544383244743555\n",
" unit: ''\n",
" min: 0.0\n",
" max: .nan\n",
" frozen: false\n",
" - name: tilt\n",
" value: 0.0\n",
" unit: ''\n",
" min: .nan\n",
" max: .nan\n",
" frozen: true\n",
" - name: reference\n",
" value: 1.0\n",
" unit: TeV\n",
" min: .nan\n",
" max: .nan\n",
" frozen: true\n",
"\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reading different datasets\n",
"\n",
"\n",
"### Fermi-LAT 3FHL: map dataset for 3D analysis\n",
"For now we let's use the datasets serialization only to read the 3D `MapDataset` associated to Fermi-LAT 3FHL data and models."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"path = \"$GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL\"\n",
"filedata = Path(path + \"_datasets.yaml\")\n",
"filemodel = Path(path + \"_models.yaml\")\n",
"datasets = Datasets.read(filedata=filedata, filemodel=filemodel)\n",
"dataset_fermi = datasets[0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We get the Crab model in order to share it with the other datasets"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"LogParabolaSpectralModel\n",
"\n",
" name value error unit min max frozen\n",
"--------- ---------- ----- -------------- --- --- ------\n",
"amplitude 3.415e-01 nan cm-2 s-1 TeV-1 nan nan False\n",
"reference 5.055e-05 nan TeV nan nan True\n",
" alpha 2.511e+00 nan nan nan False\n",
" beta -2.248e-02 nan nan nan False\n"
]
}
],
"source": [
"crab_model = dataset_fermi.models[\"Crab Nebula\"]\n",
"crab_spec = crab_model.spectral_model\n",
"print(crab_spec)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### HESS-DL3: 1D ON/OFF dataset for spectral fitting\n",
"\n",
"The ON/OFF datasets can be read from PHA files following the [OGIP standards](https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node5.html).\n",
"We read the PHA files from each observation, and compute a stacked dataset for simplicity.\n",
"Then the Crab spectral model previously defined is added to the dataset."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"datasets = []\n",
"\n",
"for obs_id in [23523, 23526]:\n",
" dataset = SpectrumDatasetOnOff.from_ogip_files(\n",
" f\"$GAMMAPY_DATA/joint-crab/spectra/hess/pha_obs{obs_id}.fits\"\n",
" )\n",
" datasets.append(dataset)\n",
"\n",
"dataset_hess = Datasets(datasets).stack_reduce()\n",
"dataset_hess.name = \"HESS\"\n",
"dataset_hess.models = crab_model"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### HAWC: 1D dataset for flux point fitting\n",
"\n",
"The HAWC flux point are taken from https://arxiv.org/pdf/1905.12518.pdf. Then these flux points are read from a pre-made FITS file and passed to a `FluxPointsDataset` together with the source spectral model.\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# read flux points from https://arxiv.org/pdf/1905.12518.pdf\n",
"filename = \"$GAMMAPY_DATA/hawc_crab/HAWC19_flux_points.fits\"\n",
"flux_points_hawc = FluxPoints.read(filename)\n",
"dataset_hawc = FluxPointsDataset(crab_model, flux_points_hawc, name=\"HAWC\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Datasets serialization\n",
"\n",
"The `datasets` object contains each dataset previously defined. \n",
"It can be saved on disk as datasets.yaml, models.yaml, and several data files specific to each dataset. Then the `datasets` can be rebuild later from these files."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"datasets = Datasets([dataset_fermi, dataset_hess, dataset_hawc])\n",
"path = Path(\"crab-3datasets\")\n",
"path.mkdir(exist_ok=True)\n",
"\n",
"datasets.write(path=path, prefix=\"crab_10GeV_100TeV\", overwrite=True)\n",
"filedata = path / \"crab_10GeV_100TeV_datasets.yaml\"\n",
"filemodel = path / \"crab_10GeV_100TeV_models.yaml\"\n",
"datasets = Datasets.read(filedata=filedata, filemodel=filemodel)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Joint analysis\n",
"\n",
"We run the fit on the `Datasets` object that include a dataset for each instrument\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OptimizeResult\n",
"\n",
"\tbackend : minuit\n",
"\tmethod : minuit\n",
"\tsuccess : True\n",
"\tmessage : Optimization terminated successfully.\n",
"\tnfev : 866\n",
"\ttotal stat : -14820.65\n",
"\n",
" name value error unit min max frozen\n",
"--------- --------- --------- -------------- ---------- --------- ------\n",
" lon_0 8.363e+01 0.000e+00 deg nan nan True\n",
" lat_0 2.202e+01 0.000e+00 deg -9.000e+01 9.000e+01 True\n",
"amplitude 3.721e-03 6.347e-04 cm-2 s-1 TeV-1 nan nan False\n",
"reference 5.055e-05 0.000e+00 TeV nan nan True\n",
" alpha 1.246e+00 3.935e-02 nan nan False\n",
" beta 6.247e-02 2.178e-03 nan nan False\n",
" norm 9.826e-01 3.029e-01 0.000e+00 nan False\n",
" tilt 0.000e+00 0.000e+00 nan nan True\n",
"reference 1.000e+00 0.000e+00 TeV nan nan True\n",
"CPU times: user 10.6 s, sys: 118 ms, total: 10.7 s\n",
"Wall time: 10.8 s\n"
]
}
],
"source": [
"%%time\n",
"fit_joint = Fit(datasets)\n",
"results_joint = fit_joint.run()\n",
"print(results_joint)\n",
"print(results_joint.parameters.to_table())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's display only the parameters of the Crab spectral model"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"LogParabolaSpectralModel\n",
"\n",
" name value error unit min max frozen\n",
"--------- --------- --------- -------------- --- --- ------\n",
"amplitude 3.721e-03 6.347e-04 cm-2 s-1 TeV-1 nan nan False\n",
"reference 5.055e-05 0.000e+00 TeV nan nan True\n",
" alpha 1.246e+00 3.935e-02 nan nan False\n",
" beta 6.247e-02 2.178e-03 nan nan False\n"
]
}
],
"source": [
"crab_spec = datasets[0].models[\"Crab Nebula\"].spectral_model\n",
"crab_spec.parameters.covariance = results_joint.parameters.get_subcovariance(\n",
" crab_spec.parameters\n",
")\n",
"print(crab_spec)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can compute flux points for Fermi-LAT and HESS datasets in order plot them together with the HAWC flux point."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# compute Fermi-LAT and HESS flux points\n",
"e_edges = MapAxis.from_bounds(\n",
" 0.01, 2.0, nbin=6, interp=\"log\", unit=\"TeV\"\n",
").edges\n",
"\n",
"flux_points_fermi = FluxPointsEstimator(\n",
" datasets=[dataset_fermi], e_edges=e_edges, source=\"Crab Nebula\"\n",
").run()\n",
"\n",
"\n",
"e_edges = MapAxis.from_bounds(1, 15, nbin=6, interp=\"log\", unit=\"TeV\").edges\n",
"flux_points_hess = FluxPointsEstimator(\n",
" datasets=[dataset_hess], e_edges=e_edges\n",
").run()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, Let's plot the Crab spectrum fitted and the flux points of each instrument.\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAF3CAYAAABE0Ck1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3zV1f3H8dfJzc1eQEII2RCEsMGwUXDhQJSq1VoVxK2AdbXVX+usFmtbrQWqYlWcaB1FXCgOVIZAUHbYM4wQEghkkXV+fwQoIJAQcvNN7n0/H4/7gPu9936/7yDyued8zzDWWkRERMQ3+DkdQERERBqOCr+IiIgPUeEXERHxISr8IiIiPkSFX0RExIeo8IuIiPgQf6cDNITo6GibkpLidAwREZEGs3Dhwl3W2pijj/tE4U9JSSEzM9PpGCIiIg3GGLPpWMfV1S8iIuJDVPhFRER8iAq/iIiID/GJe/wiIuK88vJysrOzKS0tdTqKVwkKCiIhIQG3212r96vwi4hIg8jOziY8PJyUlBSMMU7H8QrWWvLy8sjOziY1NbVWn1FXv4iINIjS0lJatGihol+PjDG0aNHipHpRVPhFRKTBqOjXv5P9M1XhFxERn2GM4brrrjv0vKKigpiYGC6++OKTOk9KSgq7du065fc4odEXfmNMG2PMS8aY9w47FmqMedUY86Ix5hon84mISNMRGhrKsmXLKCkpAWDGjBnEx8c7nKphebTwG2NeNsbsNMYsO+r4BcaYVcaYtcaY+090DmvtemvtjUcdvgx4z1p7M3BJPccWEREvduGFF/LJJ58AMGXKFK6++upDr+Xn5zN8+HC6du1K3759WbJkCQB5eXkMGTKEHj16cOutt2KtPfSZN954g969e9O9e3duvfVWKisrG/YHOkmeHtU/GZgAvHbwgDHGBUwEzgOygQXGmGmACxh31OdvsNbuPMZ5E4ClB37fuP+ERUTkZx79aDkrtu2t13N2bB3Bw8M61fi+X/3qVzz22GNcfPHFLFmyhBtuuIHvv/8egIcffpgePXowdepUvv76a0aMGMGiRYt49NFHGThwIA899BCffPIJkyZNAiArK4t33nmH2bNn43a7ueOOO3jzzTcZMWJEvf5s9cmjhd9a+50xJuWow72Btdba9QDGmLeBS62144Da3mTJprr4L6IJ3K44ad+Mg7MecDqFiIhX6tq1Kxs3bmTKlClcdNFFR7w2a9Ys3n//fQDOPvts8vLyKCgo4LvvvuODDz4AYOjQoTRr1gyAr776ioULF9KrVy8ASkpKaNmyZQP+NCfPiXn88cCWw55nA32O92ZjTAvgCaCHMeaBA18QPgAmGGOGAh8d53O3ALcAJCUl1VP0BvLtkyr8IuLVatMy96RLLrmE++67j5kzZ5KXl3fo+OFd+AcdHDV/rNHz1lpGjhzJuHFHd1g3Xk4U/mPNO/j5n/TBF6zNA2476lgRMOpEF7HWTgImAWRkZBz3/HX2ytB6PyUAO5Z47vyjPqn/c4qINEE33HADkZGRdOnShZkzZx46fuaZZ/Lmm2/y4IMPMnPmTKKjo4mIiDh0/I9//COfffYZu3fvBuCcc87h0ksv5e6776Zly5bk5+ezb98+kpOTHfrJauZE4c8GEg97ngBscyBH47JnExQc1hGyaVb1r5GJENV4/wKJiDRFCQkJ/OY3v/nZ8UceeYRRo0bRtWtXQkJCePXVV4Hqe/9XX301PXv2ZNCgQYd6kjt27Mjjjz/OkCFDqKqqwu12M3HixEZd+M2xujXq9QLV9/g/ttZ2PvDcH1gNnANsBRYAv7bWLvdUhoyMDJuZmemp0x9XZWUlFRUVVFZWHvGoqqo69OvBR1l5JaUVlZSWV5HxXi++GjqHiipLeaWlvLKKKlvdpVQF2AOdJi5jcPkZ/PwMLmMIdLsI8HcR6HYR6O8iJNBFSICbkEB//P39cblchx7+/v74+Xnf8AgRabyysrJIT093OoZXOtafrTFmobU24+j3erTFb4yZAgwGoo0x2cDD1tqXjDFjgM+pHsn/sieLvidUVFRQVFREeXn5oUdFRQXl5eUUlpSRW1hGXlEZBfur2Hvgsa+siqIyy76yKgrLqigut5RUWErKqyipsFRU/e/8G4Pgxvc31mvmIJch2F39CPH3I8RtCA/0IzzQRWSwP1HBAbQIdRMdHkRMeBAtI4JpGRVCcGAAAQEB+pIgIuIlPD2q/+rjHP8U+NST1/aUkrJKvlq6mWXrt7KruJJdxZXklVSRX1LJ7pIqiiuO3YMS6DKEBRhCA/wIcxtahroI9j9QjP39CPI3BLoMAf6GObt+zb3xzXH7G9wuPwJcfrj8DC4/P/z8DMYYDNUDIw72BFRUQUVVFWUVlooqS2lFFfvLKymtqKK0ooqSsiqKyqsoLquipKL6S8iuknIKy8ooLKvuUTiaHxAZ5EfzYD+aB/vTKjyAVpGBtI4KJrF5CCnR4bRuHkZQUJC+GIiINBFevTufMWYYMCwtLa3ezrm7uIwx72YB1YWxWbAf0SEukiL96RbrolmwH82CXEQF+REZ6KJ5qJvo8EBCgwLx96/ucne73Ud0vR/8vZ+fHy6XC2POoH+9JT6StfbQrYaDj4qKCvYU7Sd3Xym79pWyq7CMXYX7yS082HtRQU5RBSty91NUfuS82wAXxIb6ExfuT2JUECnRIbSJCaNdq0iSYyIJDAz00E8iIiJ14dWF31r7EfBRRkbGzfV1ztiIIF76dSfKC3YSE+omJDiIwMBAAgICjni43e5a743ckIwxh75wHC4qClJO8LmDtzLy9xWzJa+QLfnFbMkvJntPCdsKyti2r5xFO/ZQVrnn0GeC/Q0JEf6kNAukbUwo7VtF0DmxGSktoxrln42IiC/w6sLvCS4/w+BOiUDiz4qnNzvYWxEfHEx8yxb0Per1qqoqSkpL2bJrH2t2FLAut5D1u4rZkFfKgq3FzFhXCOQAEBZgSI0KoF3LEDrFRdAtqTkdE5oTEhzU4D+XiIivUeGvA18q+LXl5+dHaEgIHZJC6JAUe8Rr5eXl5Ozex7Ls3WRtK2D1zkLW7irl46zdfLBsN7CJABe0aRZAx1ahdEuIJCM1mvbxLfD3119REZ92cE0TrUNSb/Svqnic2+0moWVzElo254Ke/zteVFxCVnY+i7fks2zbXlbmFPPxioNfBjYSEWDo0DKI7vER9EptQe+0WCLDQhz7OUSk6XO5XHTp0uXQ86lTp5KSklKv15g2bRorVqzg/vt/vgddSkoKmZmZREdH/+y1Sy+9lJ07dzJ37lwAnnjiCd59910Ali5deij3DTfcwJ133lnnfB6fx98YODWPX05eYXEJyzbnkbkxj8XZBazIKWbrvup9mPwNpLUIoHt8OH3btGBA+zhiosLqfjG1JEQaVJ3m8b8ytHqBs7uX1fzeWggLC6OwsPCkP1dZWVkvvb3HK/x79uyhS5cuhIWF8emnn5KamnrE6zXlbjTz+EVOVlhIMH07JNC3QwJQPQth6649/LBmJ5mb8vlpayHvLsnj7cV5+LGatBYBZCSG0z8thjM6tCYyLPjkLrhnkwd+ChGpV4evauoBlZWV3H///cycOZP9+/czevRobr31VmbOnMmjjz5KXFwcixYt4tNPP+WCCy5g4MCB/PDDD3Tr1o1Ro0bx8MMPs3PnTt5880169+7N5MmTyczMZMKECbXO8P777zNs2DBiY2N5++23eeABz+3X4tWF3xPT+aRhGWNIiGnGFTHNuOLAHMfde4v4Ye0O5qzdReaWvbyzOI+3FuXh77eSji2D6JcSxVkdYjm9bSvc7hr+inv4HxQROUlH71NyvP1L6thTV1JSQvfu3QFITU3lv//9Ly+99BKRkZEsWLCA/fv3M2DAAIYMGQLA/PnzWbZsGampqWzcuJG1a9fy7rvvMmnSJHr16sVbb73FrFmzmDZtGn/+85+ZOnVqnXJNmTKFhx9+mNjYWK644goV/rryxHQ+ccBR/8M3Ay488KA57ItwMa+oJXOKE5i1K5kXdpTywg87aO43hzOCNzIodBODwrJp4a448rye3BAJdAtB5FR4aP+S4OBgFi1adMSxL774giVLlvDee+8BUFBQwJo1awgICKB3795HdLunpqYeutfeqVMnzjnnHIwxdOnShY0bN9YpU05ODmvXrmXgwIEYY/D392fZsmV07ty5bj9kDby68ItvCPev5NzI7ZwbuR3LfHLKgvh2XzzfFyfxbUkqHxZ1wm9nFd0CtnF26AauNjOILl77vxNoQySRxuPoL8yvDK3+f/SRAo9d0lrL+PHjOf/88484PnPmTEJDQ484dviiZH5+foee+/n5UVFxZOOisrKS008/HajeBvixxx475vXfeecddu/efegLxt69e3n77bd5/PHHT+0HOw4Vfmn8TqLlbIBWwFUHHvvLypm7ehtfrdjB7A2B/H13An/nDFqFuTjfzOPR8qcp/2M+bn9N0RTxVeeffz7PPfccZ599Nm63m9WrVxMfH3/K53W5XD/rXTiWKVOmMH36dPr16wfAhg0bOO+881T4ReoiMMDN4M7JDO5c3ZJfs3UXny/dyszVeby1rTuPBkHGn75gYGo4F3SO47yuSQQFaFVBkUYlMrHm95yCm266iY0bN9KzZ0+stcTExNT5Xn1tdO3a9dD+Jr1792bz5s307fu/ZdFSU1OJiIhg3rx59OnTp96vr+l8curqebpNQyl48RK25eYzrvnjzM8uprTCEhZgOCMlgqFd4zhXXwJE6lWdp/OBxszUQNP5pOE1wdHxkf6VRMZF8tqowRQWlzJ98WY+W7aDmev38tnqAsKnrWJQm0iG94xncMdE/HU7QKThqeDXO68u/JrOdxRPjV5vqqPjDztvWEgQV/Q7jSv6nXboS8DHS7YzY00BH6/cQ4vgLM49rRm/7JVMRlqcZ/KIiDQAry78ms7nYR6abuO0w78E7CksYdrCDXy8NIf3luTxzuI8UqLcDO0UzZV92pDcMsrpuCIiJ0X3+OXUNcB0m8Zge95e3p2/no+X5bI6rww/oGd8MFf0jOeSjDaEBGo8gMiJ1Okev9TKydzj92uwVOK9fOQeXFyLCO68sDtf/PY8pt2Wwa97RLMhfz/3f7SW3k/M4J435rJw3Q6nY4qInJBXd/VLAxr0812ovFnXlFi6psTycEUl0xdt5N3MbKatyOeDZfl0iA7k8h5xXNWvLREhQU5HFWnSRk0fBcArF7zicBLvoRa/1I+zPLeudGPm9ncxLKMtr902iO/vO5OxA1tTVFbJEzM20nfc19z1+hwWbchxOqaIHBAWduSOnpMnT2bMmDEAPPLII8THx9O9e/dDjz179lBcXMw111xDly5d6Ny5MwMHDjy0U94TTzxBp06d6Nq1K927d2fevHkN/jOdLLX4RepJXPNw7r24B3dfVMW3y7N544eNfJy1m6nLM+kcG8TVGfFc3ieNoAD9bydyMrYWbm2wa919993cd999RxwbN24csbGxLF26FIBVq1bhdruZO3cuH3/8MT/++COBgYHs2rWLsrKyBstaV/oXSKSe+fn5cVaXJM7qksT2/H28PmsNHyzZyR8+WcffvtrAZV1jGHXmaSRERzgdVaRJ2F603dnrb99OcvL/Ziq1b9/+0PHo6OhD6/VHR0c7ku9kaVS/SAOorKzi05828NrczSzYWoy/gUFtwrn5zDT6tm/tdDyRBlGbUf0H7+kftDJ/JYXlhWTEHjk4va73/F0u16Hd9QDy8/O55JJLmDBhAo888ggvvvgiMTExADRr1oxvvvmGRYsWMWTIENq2bcs555zDyJEjadeuHYWFhQwcOJDi4mLOPfdcrrrqKgYNGlSnXKdKK/cdoAV8pLFwufwYltGWYRltWb45l5e+W8tnK3fz1Ss/0anlCq7vl8wverfF36VhNyJQ3b1/eEs/M6e68RYXGkd8WN030Dl6W97JkydzeMPwWF393bt3Z/369XzxxRd8+eWX9OrVi7lz55Kens7ChQv5/vvv+eabb7jqqqt48sknuf766+ucryF4deHXAj7SGHVKiuHpa2P4474SXv1uFVN+3MFvP1zN375cx68zWnPD4PaEBwfWfCIRL3R0S37U9FFk5mSydORShxJVCwsL47LLLuOyyy7Dz8+PTz/9lPT0dFwuF4MHD2bw4MF06dKFV199tdEXfjUvRBzSPDyYu4d2Z87/DeEvl7SjeYg/z3y7hX7jvuLh9xawLX+f0xFFBJg9eza7d+8GoKysjBUrVpCcnMyqVatYs2bNofctWrToiLEAjZVXt/hFmgJ/lx9X9T+Nq/qfxvcrsnn+27W8lrmTt37cyQUdmnHnuem0a93M6ZgijokLbbj9MZ555hneeOONQ8+nTp3KunXruP3227HWUlVVxdChQ7n88sv58ccfGTt2LHv27MHf35+0tDQmTZrUYFnrSoP7RBqhFVt2MfGrVXy+eg9VVXBWWgRjz2lP99SW2qZUmqy6LNmrBXxqR4P7RJq4jonRTLw+ms25e5n4ZRYfLs/jqxcW0C8pjN/ZFvQIyXM6okiDUMGvf7rHL9KIJcVE8Jer+/D9bwdzY+9WLN1RxC+2XMV1my5g1sqGW9RERLyHCr9IExATGcKDl53OrN+fw73NZrGsrBXXTl7EL/75Nd9n6QuAiNSeCr9IExIVGsjYlkuY1fZ1xp6RwPq8/Vz36iIuH/8Nc1Ztczreib0y9H/jE0TEMSr8Ik1QqF8l9w7txuwHzmHMwHjW7Crl16/8xJUTZpK5fqfT8USkEfPqwm+MGWaMmVRQUOB0FJH6tWcTAGFBAdx3cXdm/f5s7hgQz4qdxVwxaQEjXviOFdn5DocUOXWbrhvBputGOB3Dq3h14bfWfmStvSUyMtLpKCL1q2DLEU8jQgL53bDqLwA39G7F/C2FXDxxLre9MpuNuXsdCinS+JxoW96DunXrxtVXX33o+eLFi+nevfuh51OmTCEkJITy8nIAli5dSteuXQEoLy/n/vvvp127dnTu3JnevXvz2WefeerHqRNN5xPxFE/dz96x5LjnjwIeAkYnBzIxtxtvrOrOV6u+4ZqIxYyNWUQL//01n1/rA0gjU7a14QawZmVlUVVVxXfffUdRURGhoaF06dKFTZs2sW/fPsLDw5kzZw4dOnTgp59+onfv3syZM4cBAwYA8OCDD7J9+3aWLVtGYGAgOTk5fPvttw2Wvza8usUv4lX2bIJNs2D/gRb8plnVjwPd/odr4b+fh+Lm83XqG1wStpLX9vZg0PpreXZnN4qrXA0cXOTUVGxruIGrb731Ftdddx1Dhgxh2rRpQPVW27169WLevHkALFy4kNGjRzNnzhwA5syZQ//+/SkuLubFF19k/Pjxh7bqjY2N5corr2yw/LWhFr+Ip3iq5fzK0OqC/0jNY1cSgL8DozbvYtwny3lm0wBeLT2TsYNTGHFGe1x+xjMZRRqpkpKSI7rtD27Le9A777zDjBkzWLVqFRMmTDjU5d+/f3/mzJlDv3798PPzY/DgwTzwwAPcddddzJkzh4cffpi1a9eSlJREREREg/9cJ0MtfhEf0DkpmjdvH8RrI7vTMszNo5+t45y/fsnnS7bU/GERB+SOn0BWh3SKFywAIKtDOlkd0skdP+GUzntwW96Dj8cee+zQawsWLCAmJobk5GTOOeccfvzxx0Ob8wwYMIA5c+Ywf/58evXqRdu2bVm7di25ubkUFhbSpk2bU8rVkFT4RZqiyMQ6fezM9Hg+u+dsnhrentLyKm59awlXTpypGQDS6MSMHUP6yixCevUCIH1lFukrs4gZO6aGT9bdlClTWLlyJSkpKbRt25a9e/fy/vvvA9C3b18WLFjArFmz6NevHwAJCQm8/fbb9O/fH4C0tDQ2b97Mvn2Ne2dNFX6Rpiiq7lt/GmO4sm8aM39/DncNSmL5jmIunjiXe96cx659pfUYUqTpqKqq4t1332XJkiVs3LiRjRs38uGHHzJlyhQAwsPDSUxMZPLkyYcKf79+/fjHP/5xqPCHhIRw4403cuedd1JWVgbA9u3bj9jtrzFQ4RfxUUFuf+66sAvf/PYshneOZuqyXQx66msmfLGc8soqp+OJHOLfurXHr/Hdd98RHx9PfHz8oWNnnnkmK1asYPv27UB1d//+/ftJTKzucevXrx/r168/VPgBHn/8cWJiYujYsSOdO3dm+PDhxMTEeDz/ydC2vCJNjYe25V2Rnc9DU5eQmV1EUlQADw1N59wuCfV3AW0n7PPqsi3vwcV7kl9/zRORvMbJbMurFr+IANAxoTnvjh7EP3/ZibKKKm56czHXvvA9G3Mb9/1K8W7Jr7+mol/PVPhF5BBjDJecnsK395/LmDMSydyyjyH/+I5x0xZRWl7pdDwRqQcq/CLyM4H+Lu4b2pUv7xnEwNRIXpizlbOe+pJPF212OpqInCKvLvzapEe80qhPGuw+eULzUF6+aSAvX9cdt5/hjreXcu0L37Mlr7BuJzzGKoPiW3xhXFlDO9k/U68u/NqkR6R+nN0pni9/ey53DkpiweZ9nPv0dzwzfdnJj/4v0IJBviwoKIi8vDwV/3pkrSUvL4+goKBaf0aj+kXkSDVsLrSlPIxHtvfnq5I00ty7+HOrmfQO2VnzeXcsqd5nIHlgPQU9imYLNHrl5eVkZ2dTWqr1IupTUFAQCQkJuN3uI44fb1S/1uoXkZOS6C7kpaQv+GJvFg/vHMSVW67gqvClPBA7jyhX9aIlubPyiRnYvPoDezYd2dLfNKv618jEU1qISJoet9tNamqq0zF8ngq/iBypli3nIcCA/RU8MW0xb/8IX5R3549D23P6Px+mZMF6Yl6ae+QHTmJzIRHxHK++xy8inhUa6M+ff3k6H9zWl5gwN/e+v4LfN+tf8wdFxDFq8YvIKeue3ILXA1aQP/Vfh45ldaheRSx69Oj/baxSx82FRKT+qMUvIvUi9s6xpK/MInhW9TaqFw7/G7+//V8UXDXyf2/SPX0Rx6nwi0i9SokOA+DxS9qzPq+Ui579nvEzVlBpjUeuN2r6KEZNH+WRc4t4IxV+Eal30aNHc23/NL66dzB9kiP4+1cb+MXG4azeH+V0NBGfp8IvIvXu4D392MhgXrt5AH+7rCObKyIZuvFKnvl8ORXa9lfEMSr8IuJRxhiu6J3KjNS3OTtkA89+s5Hu417nqg/GOh1NxCep8ItIg4jxL+X5hBn87fJOlJaEM3/Bufzrq5VUVXn/6qEijYkKv4g0GGPgil4p9O49g6ioHTw1Yx2XT/yOLflFp3TerYVb6ymhiPdT4ReRBhcYWEq3rj/w50vTWbmziCHPfMubc9fXefOW7UXb6zmhiPfSAj4i4ohtRVv59YVtOKN9LHe+uYA/fJjF58u28/TVGUSHBR56X01T9Vbmr6zV+wBeueCVUwst4gXU4hcRRxxspSc2D+X90YP47bltmbNhD+f9/RtmLN9W4+e3Fm4lMyeTwvJCADJzMsnMyVS3v0gNvHpbXmPMMGBYWlrazWvWrHE6johPOG7Le8fS6l9bdWFl/koKywvJiD1yx9DCwgiWr+hFcXEz4uJW0y5tOa8N/XeN18vMyWTpyKX1EV/EaxxvW16vbvFbaz+y1t4SGRnpdBQRAbZSccJWeljYXjJO/4b41ivZvv00FmSezbLsPQ2aUSsBirfTPX4RqVfHvY/+ytDqX6//pFat9Jkrd3DvfxYz/F+zuffcNG4dfBp+fsde9jcuNO5UY4v4DK9u8YtI0zW4Qytm3HsW/VOj+MsXa7lm0mxy9+0/9Hru+AmHfh8fFu9ERJEmSYVfRBxRm1Z689AAXr2pPw9d1J6FW/Yy5Olv+DprBwC7Jk70dEQRr6SufhFxRG1b6cYYbjgzjf7tWnLH6wu44dWF/GrvCkbW/FEROQYVfhFpGKM+OaWPd4iL4M3A5eyZ+tyhY1kd0vkdMPv8eLjgFPOJ+Ah19YtIkxH3mztJX5nF6re+AOCXVz7Nff93GbMvTHQ4mUjToRa/iDQ5l/ZMJAtIbBbE8hX96R4aT1lFFQH+9dOW0SJA4s1U+EWkwdXH0rnRo0fz4e2DeGzaEt6Yv5WfNu/mhZF9SGwecsrn1tr/4s28euW+gzIyMmxmZqbTMUTEQz5dspXfvrcEA/z1l924sEvrY76vNgvzHG9VwWPR2v/SmPnkyn0i4hsu6hrPZ78ZREJUELe/+ROPfbiEisqqkzqH1v4XX6EWv4h4jdLySh7672L+8+N2eiSE88KI3rSMCDqpc2jtf/EWavGLiNcLcrt46sqe/O3yzqzYUcj5z8xk7tpcp2OJNCoq/CLida7olcyHowcSHujPNS/N5/lvVnMyvZta+1+8mQq/iHilDnERfHLXIM4+rQVPfr6GW1+bT9H+ilp9Vmv/izdT4RcRrxUe5ObF6/tw33lpfLlyF0P/8S3rcgudjiXiKBV+EfFqxhjGnNOe127oze6Scob983s+X6Z5+uK7VPhFxCcMbBfDZ3cNIrl5ELe+8SN//WwFVVXeP6tJ5Ggq/CLiM1pHBfPfMWdySddYJn67gVEvz2VvabnTsUQalAq/iPiUILeLZ68+nQeHtmfWut1c/Kxn7/uPmj6qVisGijQUFX4R8TnGGG48I403b+pDQUkFl4z/nq+zcoDqZXgb81K8+iIhp0qFX0R8Vt+20XzymzNpHRnIja9m8q+v/zffP3f8BIfTiXiGVxd+Y8wwY8ykgoICp6OISCOV0CyED8eeyXkdonnqizWMeTOT0vJKdk2c6HQ0EY/w6sJvrf3IWntLZGSk01FEpBELCfDnhZG9uevstnyybCfDf/eG05FEPMbf6QAiIo2BMYZrVn3J+VP/19LP6pAOQPTo0cSMHeNUNJF65dUtfhGRkxEzdgzpK7Oo+uoHAIZf/nfWvfOlir54FRV+EZGjdIqvvj3YPjaUMVMW8Y8vVp7UJj8ijZkKv4jIMUSPHs27dwzk4s4t+cfX67jzrYXsr6h0OpbIKVPhFxE5hpixYwj0dzH+mgzuPKsNHy3N4arnZpNfVOZ0NJFTUqvCb4zpbowZa4wZZ4x5yBhzmTFGQ+VFxOsZY7jn/Pgj9r4AACAASURBVHSevaorK3YUcvGz37J258mt9Le1cKuH0omcvBMWfmPMtcaYhcCjQDNgE7AXOBeYaYx5yRiT4PmYIiLOurRHIm/f0o+S8kqGT/ieOWtza/3Z7UXaDVAaj5qm87UAzrTWFh3rRWNMBpAOZNd3MBGRxqZncjOmjT2T6/49lxEvzyet3Xzi4jaf8DMr81cC1GqZ3ca8VLB4jxMWfmvtszW8nlm/cUREGrfE5tUr/d30yg8sWNWb4pJQ2qRmYcyR79tauPWIln5mTvU/l3GhccSHxTdkZJEj1LiAjzEmALgIOANoDZQAy4BPrbUrPRtPRKTxiQx28+YtA/j9uz/x38XQLfIs/n7V6QT4//zu6ajpo8jMyWTpyKX1dn2NGZBTUdM9/j8C84CzgMXAq8A0qr8wPGOMmW6M6ezxlCIijUyAvx9P/6onvzmrDR8t3ck1k2azt7S8Qa6tMQNyKmpq8S+11j5+nNeeMsbEAYn1nElEpEkwxnD3+enENwvhgf8uY/j473jj5v60jgo+4n1xoXG1Ol9txgFozICcqhO2+K21Hx59zFQLPfD6dmvtfE+FExFpCq7snczkUb3YsbeMS8Z/x4ptR+4IWh/39LcWbiUzJ5PC8uqphJk5mWTmZKrbX05arTbpMca8BowBKoBMINoY86S19mlPhhMRaSrOOK0l79/RnxH//oErnpvDiyMyGNAu5qTOUZsWuifGDIhvqe3KfV2stXuB4cAXQAJwvadCiYg0RelxkXw49kxiwwO4/pX5fPjTFqcjifxMbQt/gDHGH7gUmGqtLQOqPBdLRKRpah0VzNQxZ9I5LpzfvLOEzZvT6v0atR0zIHIstS38/wY2U71637fGmCTg5NasFBHxEZEhbqbcNoBz27dg3frurFnbmaqq+tvdT+sAyKmoVeG31j5jrW1trR1iq/emzAbO9mw0EZGmK8jt4oWRfYiLW0N2dgd+MyWT8kp1lIrz6rQ7n7W26kB3v4iIHIfLz9D+tMWkJC/jo6U7GfXSXIrLKpyOJT5O2/KKiHiQMZCaupLHLklnzoY9XPncLPYUq90kzlHhFxFpACP6t2HC1d1ZlVPMLyZ8z46CUqcjiY+qsfAbY0KNMcnHON7JM5FERLzTRV3jeWVUBjn79nPp+O9Yn6sx0tLwalqr/3JgLfCJMWapMabnYS+/7tFkIiJeaGC7lky5pR+lFVVcNnEWy5/4q9ORxMfU1OJ/EMiw1nYGbgWmGGMuOfCaOf7HREQEqlfjO3pFvm6JzXj/jgEE+vvh9/rL/LAu16F04otqKvx+1tqtANbaOVRP4XvUGDMaqL9JqSIiPiatZTgT1rwPwIiXF/DVCu24Jw2jpsJfZIxJPfjkwJeAwcAvgY4ezCUi4rVyx08gq0M6YfO+A+DD9++l9WVn8/3/PeFwMvEFNRX+0Ry1kY+1tgAYQnXXv4iInKSYsWNIX5lF+sosAOIXLeH+O55jRFV3Xp+9zuF04u1q2pb3R2vtmoPPjTEhxpgIIAj4r6fDiYj4goggN2/fdgZ9UyJ58KOVPPf1KqcjiRer1Tx+Y8xNxpjtwGpgGbD8wK8iInIKokePBiA4wMWrN/XnnNOa85cv1vK3z5ZTvUK6SP2q7QI+vwe6WWsTrLVJ1tpEa22SJ4OJiPiCmLFjDv0+wN+PF0b24eJOMUz4diN/mrZExV/qnX/NbwFgPbDXk0FERAT8XX7885peBL+7kJfnZlO0v4JxV/TEz696BvXRUwNFTlZtC//9wGxjzA/A/oMHrbX3eCSViIgP8/MzPHXl6YQELObVeVspLV/A01f3wuWn5VPk1NW28D8PzAaWAtpXUkTEw4wxPDK8G8FuF8/P2kxZxTzGX9sbf5e2WJFTU9vCX2WtvdOjSURE5AjGGO6/uAuB/n48O3Mjt0z+gReu74tbxV9OQW3/9nxljLnBGBNjjIk4+PBosgOMMW2MMS8ZY9470TEREW919wWd+N15bfl6zW5ueGkO+ysqnY4kTVhtC/9I4FHgR6qn8tVqOp8x5mVjzE5jzLKjjl9gjFlljFlrjLn/ROew1q631t5Y0zEREW92xzkd+OMFaXy/voDrX5xDabmKv9RNrbr6rbWJdTz/ZGAC8NrBA8YYFzAROA/IBhYYY6YBLmDcUZ+/wVq7s47XFhHxKjcNbo/b349HPl7NiEmzeO2WgQS5XU7Hkiamtgv43GaMiTrseTNjzC01fc5a+x2Qf9Th3sDaA632MuBt4FJr7VJr7cVHPepc9I0xtxhjMo0xmbm52vlKRLzDyIHt+NOw9izYUsg1L8yipEwtfzk5te3qv81au+fgE2vtbuD2Ol4zHthy2PPsA8eOyRjTwhjzPNDDGPPA8Y4dzVo7yVqbYa3NiImJqWNUEZHG59oBafz50g78tLWQq5//nuKyCqcjSRNS21H9R/QlGWP8AHcdr3msiajHXZrKWpsH3FbTMRERX3J1v7b4GcMDH2bx6+e/563bziAkoLb/pIsvq22Lf4YxZooxZpAx5kzgTeDLOl4zGzh8zEACsK2O5xIR8VlX9W3DX4ans2R7MVc/p5a/1E5tC/9vqV7A527gXmAWcF8dr7kAaGeMSTXGBAC/AqbV8VwiIj7tl32qi//SHdXFv2h/udORpJGrVeG31lZaaydYa4dbay+11k601tb41dIYMwWYC7Q3xmQbY2488LkxwOdAFvAfa+3yU/khRER82S/7tOHJA8X/18/PUvGXEzIn2vnJGDMVeAGYcXShN8YkUz2/P9ta+7JHU9aRMWYYMCwtLe3mNWvWOB1HRMSj/jNvPfdPzaJrXAhv3XoGIYG65+/LjDELrbUZRx+vqcU/mur59quNMXONMdOMMV8YY9YCrwDLG2vRB7DWfmStvSUyMtLpKCIiHndlnzaMu7QDS7YXc+2kWZTonr8cwwm/DlprtwL3APcYY9KAOKAEWGWt3dcA+URE5CRc1bctFZVV/PGj1Vw3aRZv3DKQII32l8PUeqcHa+1aa+331tpMFX0RkcbrmgHteHRoOxZmFzHyxdmUlleQO36C07GkkdAWTyIiXmjEGafx4IVtmbelkJH3/JtdEyc6HUkaCfX/iIh4qRsGdaCiooo/z6h+XlFZhb+29PV5J/wbYIz5lzFmmDEmpKECiYhI/cgdP4EzRl/OZ1Orl11Z06kTWR3S1e3v42qazjcQuAA4Gyikeu799KYy717T+UREYNN1IyhesIALh/+NC9pH8a+R/fHzO9bq6eJN6jSdz1o7y1r7R2ttf+BaYCfwB2PMj8aYScaYyzyUt15oOp+ICCS/Xr0z+o194pi+ag93vfEDJ2r0iXc7mVH9O621r1trfw2cDrwEdPFYMhERqTfRo0fzx+E9uLZnS6atyOf+dxao+PuoOg3us9V/W+YdeIiISCMXM3YMAH/6ZQbF5fN4Z1EuIe6FPHz5z3qCxctpeKeIiA8xxvC3q/twYfsoXlmQw18/+snpSNLAahrV72qoICIi0jD8/AzjR/TjrLbhTJy9jQmfL3U6kjSgmlr824wxzxljzmyQNCIi0iD8XX68MGoA/ZJC+ds3m5k8M8vpSNJAair8XYBlwBPGmM3GmL8ZY05vgFwiIuJhAf4uXr5xAN3jgnls+nrenatpz76gpul8O621E621ZwADgO3A88aY1caYRxsk4Sk4sPjQpIKCAqejiIg0SsGBbl67eQCnRQfywEer+ezHDU5HEg87mel8W4DngGeAIqq37G3UNI9fRKRmESGBvHFLfxIi3Nz1fhYzl29xOpJ4UI2F3xgTYIz5hTHmP8AGYCjwENVb9IqIiBeIjgjhzZv70iLExR1vLyVz7XanI4mH1DSq/zUgGxgJfACkWmuvOdCSLm+IgCIi0jDiW0Tw2g29CPb346Y3FrFiS67TkcQDamrxfwucZq0dbq1921pb3BChRETEGWlxzXl5ZA8qq+D6yQvZmLPH6UhSz2oa3PeStXaPMSbGGPOCMeYTAGNMR2PM9Q2SUEREGlS3lFiev7oze/dXMeLleeTsLnQ6ktSj2g7um0x16z/hwPM1wL2eCCQiIs4bkJ7IM5e1Z9u+Ckb8ey4FRSVOR5J6UtvC39Ja+xZQBXDg/n6lx1KJiIjjLuzZlicuasuavDJGvTSH4tIypyNJPaht4S8yxjQHLIAxphewz2Op6onm8YuInJqrBnbgt2cl8OO2Uu54dQ7l5RVOR5JTVNvCfx/wEdDGGPMtMAUY67FU9UTz+EVETt3t53fjpt4xzNxQxH1TfqCyUh2+TdkJt+U1xvS11v5grc00xpwFpAMGWGGtVZ+PiIiP+L/hGeQXzeWD5Xto/kEmD13RG2OM07GkDmpq8f/r4G+stWXW2sXW2kUq+iIivsXPz4+//KoPZ7UJ5ZWFu/jX9EVOR5I6qvWSvSIi4tvcbn8mjuhHz7gg/v7tNt7+foXTkaQOTtjVT/U9/WnHe9Fae0k95xERkUYsJCiQF6/vw69emMNDn22geWgAQ3qmAbDpuhGUbd1Ku6+/cjilnEhNhT8X+HtDBBERkaahRWQYL43M4KoX53PPf1fzSmggvdonAlCxbZvD6aQmNRX+Qmvttw2SREREmoyk2Ob8+9ruXDP5J257exlTbggk0OlQUis13ePXxswiInJMnVJaMeGKdIYtmk7leYMoXrAAgKwO6WR1SCd3/ASHE8qx1NTif/pELxpjIoAka+2y+oskIiJNxcDOKeT+7jYunjaETkXb+MsXT5O+MsvpWHICNRX+y40xTwHTgYVU3/MPAtKAs4BktGa/iIhP+0W/dHbuLWXcN9XPS0v3ExSkjv/GylhrT/wGY5oBVwADgDigBMgCPrHWzvJ4wlNgjBkGDEtLS7t5zZo1TscREfFaVVVV/OmWP1GWm0fJ5b/kr9f0x+VyOR3LpxljFlprM352vKbC7w0yMjJsZmam0zFERLxaeXk5d78xl49X7eP2PtH8brhW93PS8Qq/FvAREZF64Xa7eepXvemTEMTz83bx2teLnY4kx6DCLyIi9SYkOIiJ1/aiXQs3T3y1lc8WrHQ6khxFhV9EROpVdFQEz/26O82DXfzuo/XMz9rodCQ5TK0KvzEmxBjzoDHmxQPP2xljLvZsNBER8bRR00cxavqoej9v2/iW/POKDhjgzvdWsi47p96vIXVT2xb/K8B+oN+B59nA4x5JJCIiXqF3egrjhqaSX1LJHVMWk5u/x+lIQu0Lf1tr7VNAOYC1tgTQUE0RETmhi3p34P7BcazOK+c3U36kqLjY6Ug+r7aFv8wYEwxYAGNMW6p7AERERI7LGMPIc7oxqmcz5mwp4aH3MikvL3c6lk+rbeF/mOrV+xKNMW8CXwG/81gqERHxGi6Xi99d0pML0kJ5f8U+nv1kIZWVlU7H8lm1KvzW2hnAZcD1wBQgw1o703OxRETEmwQFBfHklafTs1Ug/5qXxzvfLcUXFpBrjE5Y+I0xPQ8+qF6XfzuwDUg6cExERKRWoiLCefaqbiRF+POnr7bxzU+rnY7kk2rapOfvB34NAjKAxVQP6usKzAMGei6aiIh4m8S4GP5xWXtGTcnidx+vZ3JEMJ3TkpyO5VNO2OK31p5lrT0L2AT0tNZmWGtPB3oAaxsi4KkwxgwzxkwqKChwOoqIiBzQvX0K4y5MYt9+y13vr2RrTq7TkXxKbQf3dbDWLj34xFq7DOjumUj1x1r7kbX2lsjISKejiIjIYYb0Suf3Z7Zk3e5y7n13CXv37XM6ks+obeHPMsb82xgz2Bgz6MAKflmeDCYiIt7Lz8+P687uxqgeUfyQXcpj//2J0tJSp2P5hNoW/lHAcuA3wF3AigPHRERE6sTtdnPvxd0Z0iaE91bs44XPF1FRUeF0LK9X2+l8pdbaZ6y1vzjweMZaq69mIiJySkJDQ3n8sq50aRnA+B/ymDp7mab5eVhtN+kZYIyZYYxZbYxZf/Dh6XAiIuL9Wka34K/D29My1MWjX21j9mJN8/Ok2nb1vwQ8TfX0vV6HPURERE5ZhzZJPDk0BWvh959sYNWGzU5H8lq1LfwF1trPrLU7rbV5Bx8eTSYiIj5lYLfTePDsVuworOT3U1exc1cem64bwabrRjgdzavUtvB/Y4z5qzGm31Gr+YmIiNQLPz8/fjGgM7f1asainDIe/XAJlZWVlG3d6nQ0r1LTyn0H9Tnwa8Zhxyxwdv3GERERXxYQEMAd53dj4+75fLKmmDaBSVy07UenY3mVWhX+A6v3iYiIeFxoaCgPDuvM1jcX8S/bn4uYSmVlJS6Xy+loXuGEhd8Yc8+JXrfWPl2/cURERMD1zjuMe2HioeerO3UGIHr0aGLGjnEqlleo6R5/+IFHBnA7EH/gcRvQ0bPRRETEV8WMHUOHrBVs69EPgJtH/APzyWcq+vWgpk16HrXWPgpEU71Jz73W2nuB04GEhggoIiK+yRhDW1c5ANv2VfCHj9eyI2enw6mavtqO6k8Cyg57Xgak1HsaERGRwxhjcMXFcWOPSBZu389Tn61gnzb0OSW1HdX/OjDfGPNfqkfz/wJ41WOpREREDghMSOD2czuycc8iPlhZRMrXS7n1gtMJDAx0OlqTVNtR/U8YYz4DzjhwaJS19ifPxaofxphhwLC0tDSno4iISB0kv/7aod//4YJ2ZP8niwnzdpPUbDnDBnTTSP86qG1XP9baH621zx54NPqiD2Ct/chae0tkZKTTUURE5BSlJCXyp/OTiAzy40/f5JC5bLU29KmDWhd+ERERp/Xo2I4Hz4plX1kVD3+xhY2btzgdqclR4RcRkSbD5XIxpFdH7uzTjJV55TwxfS25ublOx2pSVPhFRKRJCQoK4tozO3JZh1C+3FDCi9+spLCw0OlYTYYKv4iINDlRUVHcdU4ap8cF8tJPe5k2dwVlZWU1f1BU+EVEpGlKSkzgwXMTiA1z8eSsPOYuXkVVVZXTsRo9FX4REWmyOrdP48HBsZRXWh77ehur121wOlKjp8IvIiJNlr+/P4N7duCuvs1Yt7uCJ7/cyPbt252O1aip8IuISJMWHBzML/u358qOoczcVMqL366hoKDA6ViNlgq/iIg0ec2bN+eOQan0ah3I5MX7mPbDSkpLS52O1Sip8IuIiFdITkrigbPiaRXm4q9z8vlhySoqKyudjtXoqPCLiIhXMMbQpUMafzgzhrJKy+Pf5LBqzTqnYzU6KvwiIuI13G43g3q0587eUazdXc7T32azbds2p2M1Kir8IiLiVcLCwri8T9tDK/u9Nnsde/bscTpWo6HCLyIiXic2NpbbBybSLTaAf/+0l88zV2mw3wEq/CIi4pXapKZw/6BYmge5eGp2PguXrdZgP1T4RUTES7lcLrqnt+P3A5uzb38V477NYe269U7HcpwKv4iIeK3g4GDO6p7GzadHsCy3jOdmb/X5wX4q/CIi4tWaN2/OVb2SGdImmKmrivjvgg0+PdhPhV9ERLxeYmIio/u3Iq2ZmwkLCpi9ZI3PDvZT4RcREa9njCH9tDR+P7A5bj8YNyuf5avW+uQ2vir8IiI+bmvhVqcjNIiAgAB6d27HXX2j2Lq3gmdm5bBhg+9t4+vVhd8YM8wYM0m7NImIHN/2It/ZxjYyMpLzuyXzq85hzNpSypTMreTk5Dgdq0EZa63TGTwuIyPDZmZmOh1DRKRBjZo+qsb3rMxfSWF5IRmxGTW+95ULXqmPWI6z1rJq9Wr+MH0Li3fs5/GzWzC0byfCw8OdjlavjDELrbU/+w/r1S1+ERE5tq2FW8nMyaSwvBCAzJxMMnMyfaLb3xhD2zZtuKd/C5oHu/jrnN38tGIN5eXlTkdrEP5OBxAREc+oTQt91PRRZOZksnTk0gZI1Hi43W66tG/Lbwv2839f5/H3OXn8OXwtHdM7YIxxOp5HqcUvIiI+KSIiggEdk7ixRwSLc8p4LXMn2dnZTsfyOBV+EREfFxca53QEx7Ru3ZpfdG3J4OQg/rOikBlLt5Cfn+90LI9S4RcR8XHxYfFOR3BUmzZtuKNPCxIi/PnHD3tYmLWekpISp2N5jAq/iIj4NLfbTXq7ttzXL4qyKvjbnHxWr13ntTv5qfCLiIjPCw8Pp1f7RG4/PYJVeeW8tCCXTZs2OR3LIzSqX0REhOr7/Rd2KiRrVznTVhfTMWY7vwgNJTY21ulo9UotfhERkQNSU1O56fRmtInyZ/z8An5avYnCwkKnY9UrFX4REZED3G43p6W14d5+UVgLf5+7h1Vr1lFRUeF0tHqjwi8iInKYiIgIuqclcEevSNbkl/Pyj/msX7/e6Vj1Rvf4RUREjtK6dWvO61DIitwyPllTTMfoXMLCwmjdurXT0U6ZCr+IiMhRjDGkpqYyqrCY1XnlTMwsICVqC2FhYURERDgd75Soq19EROQYAgICaNc2lXv6RQLwzA8FrF67nrKyMoeTnRoVfhERkeOIioqiS2prxvSKZO3ucib/tJv169fTlLe0V1e/iIjICSQkJHBWYSHLdpbx8ZpiOsXkERa2lYSEBKej1YkKv4iIyHFsum4EAKkvTuL6ohJW7ipjwoICUqK2ERYWRlRUlMMJT566+kVERGoQFBRE29Rk7ukXRZWFZ37Yw9r1G5rk/X4VfhERkVpo0aIFnZJacntGBKvzy3ljcUGTvN+vrn4REZFaSkpK4uyiIpbklDF1VRFdWuYTGppNYmKi09FqTYVfRESkllwuF6mpqdxYvJ9VeeX8c34ByVH+hIeHN5n7/erqFxEROQmhoaGkJsVzb78oSiqqeHZeAes3NJ37/Sr8IiIiJ6lVq1Z0jG/GjT0iWLqzjPeW720y9/vV1S8iInICZVu3HvN4amoq5xcXszSnjLeXF9IpJqBJzO9X4RcRETmBim3bjnnc7XaTkpLCraXlrM4v5x/z9pAYUX2/PzIysoFT1p66+kVEROooKiqK5Nax3NM3kvySKp5bWMCGRn6/X4VfRETkKLnjJ5DVIZ3iBQsAyOqQTlaHdHLHT/jZexMSEugaH8mvO4cxN3s/09fsY8OGDY32fr+6+kVERI4SM3YMMWPHsOm6ERQvWED6yqzjvtfPz4/U1FQuLSlhyc4yXl60lw7RAYSHb6d169YNmLp21OIXERE5RcHBwSQmJDC2dyRB/n48PXcPG7O3sW/fPqej/YwKv4iIyAn417LVHhsbS3LLZtzZO5LNeyt4fUl1l39FRYWHE54cdfWLiIicQEB8fK3fm5KSQnFxMRe328/Ha4rpHhtASMhG0tLSPJjw5KjFLyIiUk/cbjfJyclc2yWclEh/Ji4oYOOOfHJyco77mU3XjTi0/W9DUOEXERGpR1FRUcTFxnB33yhKKy3/nF/AluxsiouLnY4GqPCLiIjUu8TERNrGhDKqewRLdpYxbVUR69evp6qqyuloTaPwG2PaGGNeMsa8d9ix4caYF40xHxpjhjiZT0RE5HAHp/gNaRNCn/hA3lq6jxU7Ctm8ebPT0Txf+I0xLxtjdhpjlh11/AJjzCpjzFpjzP0nOoe1dr219sajjk211t4MXA9cVe/BRURETkFoaCitW7fm9oxIIoL8+Me8PWzNySU/P9/RXA3R4p8MXHD4AWOMC5gIXAh0BK42xnQ0xnQxxnx81KNlDef/44FziYiINCqtWrWiVbNw7uwdyfZ9lUxevI/Nmzezf/9+xzJ5fDqftfY7Y0zKUYd7A2uttesBjDFvA5daa8cBF9fmvMYYAzwJfGat/bH+EouIiNQPYwypqamUlJRwafsypq4qokerIoKCNtC+fXuqS1nDcmoefzyw5bDn2UCf473ZGNMCeALoYYx54MAXhLHAuUCkMSbNWvv8UZ+5BbgFICkpqZ7ji4iIL0h+/bVTPkdgYCCJiYn8qqKSJTn7eS6zgHbN3YSHbyP+JNYIqC9OFf5jfcU57m4G1to84Lajjv0T+OcJPjMJmASQkZHROHdKEBERnxAdHU1BQQF39a3ktzPy+Of8Ah4K8iMiIqLBszg1qj8bSDzseQJw7A2PRUREvEBycjIpzYMZ1T2cpTvL+GRNsSO7+DlV+BcA7YwxqcaYAOBXwDSHsoiIiHicv78/ycnJnJsaTK/WgbyxdB9rcosp3V/aoDkaYjrfFGAu0N4Yk22MudFaWwGMAT4HsoD/WGuXezqLiIiIkyIjI4mJieH2jEjCAvx4dl4BxZWW8gbcyMfjhd9ae7W1Ns5a67bWJlhrXzpw/FNr7WnW2rbW2ic8nUNERKQxSExMpGVEMGN6Ve/i92p0BtY23Ip+TWLlvroyxgwzxkwqKChwOoqIiAhQvapfSkoKPeOCuDAthA+bd2FhUFzDXb/BruQAa+1H1tpbIiMjnY4iIiJySFhYGLGxsVzXNZzE/btZl1fK7qKyBrm2U9P5REREfFrr1q0pKCjgvu0zOW3NbLJ3l9AsNMDj11XhFxERccDBVf2q9lev3d8loWF6p726q19ERKSxyh0/gY09euLKygIgq0M6WR3SyR0/waPXVYtfRETEATFjxxAzdgybrhtB8YIFpK/MapDrenWLX6P6RUREjuTVhV+j+kVEpCnwb926wa7l1YVfRESkKQhowF36VPhFRER8iAq/iIiID1HhFxER8SEq/CIiIj5EhV9ERMSHeHXh1zx+ERGRI3l14dc8fhERkSN5deEXERGRI6nwi4iI+BAVfhERER+iwi8iIuJDVPhFRER8iAq/iIiID/n/9u4+VNKyjOP492cmyhob+BKhpUlaWULmZpRFL5YohOZLWNmLshn9ofRHEUZChIJG/SWZZehuSaxvabb2YpFiiGK6apu6CGKmy2KtRoGiQdvVH/OcnD3NnDPn7MyZZ+f5fmDYOfdz3/dcOxf3uZ57Zs4zFn5Jkjpkpgu/F/CRJGlne047gEmqqo3AxjVr1pw77VgkqY3Wnbhu2iFohc30jl+SJO3Mwi9JUodY+CVJv6nWGgAABq1JREFU6pCZfo9fkqS2O+SaH6/o47njlySpQyz8kiR1iIVfkqQOsfBLktQhM134vXKfJEk7m+nCX1Ubq+oLq1evnnYokiS1wkwXfkmStDMLvyRJHWLhlySpQyz8kiR1iIVfkqQOsfBLktQhFn5JkjrEwi9JUoekqqYdw8Ql2Q78ZV7zamDQJf0GtQ9q2x94diwBLs2wuCc9xyhjFuuz0PFRn/dB7btzLpYzz6j9l5sP18b4x7g2VmYe18bODqmqA/6vtao6eQOuHLV9SNv9bYp70nOMMmaxPgsdH/V5H9S+O+diOfOM2n+5+XBtjH+Ma2Nl5nFtjHbr8kv9G5fQPqzvNIwjluXMMcqYxfosdHwpz3tb8jGuOJY6z6j9l5sP18b4x7g2VmYe18YIOvFS/yQkub+q1kw7DpmLtjEf7WEu2qUt+ejyjn9XXTntAPQ/5qJdzEd7mIt2aUU+3PFLktQh7vglSeoQC78kSR1i4ZckqUMs/GOW5GNJfpjkliQnTDuerktyWJKrktw47Vi6KMmqJD9q1sRZ046n61wP7THNWmHh75Pk6iR/S/LwvPYTkzyW5PEkFyw0R1X9rKrOBc4GzpxguDNvTPl4oqrWTjbSblliXk4DbmzWxMkrHmwHLCUfrofJWmIuplYrLPw7Ww+c2N+Q5BXA5cBJwJHAJ5McmeSoJLfOux3YN/TCZpyWbz3jy4fGZz0j5gU4GHi66bZjBWPskvWMng9N1nqWnosVrxV7ruSDtV1V/T7JofOajwUer6onAJJcC5xSVZcAH50/R5IAlwK/qqoHJhvxbBtHPjR+S8kLsJVe8X8INxoTscR8PLqy0XXLUnKRZAtTqhUuxMUdxMs7Fuj9Ijtogf7nAx8GzkjyxUkG1lFLykeS/ZJ8Hzg6ydcmHVyHDcvLTcDpSa6gxZcwnUED8+F6mIpha2NqtcId/+IyoG3oVY+q6jLgssmF03lLzcdzgCdgkzcwL1X1AnDOSgejoflwPay8YbmYWq1wx7+4rcDr+n4+GNg2pVhkPtrKvLSL+WiP1uXCwr+4+4DDk7whyV7AJ4CfTzmmLjMf7WRe2sV8tEfrcmHh75NkA3AP8KYkW5Osrap/A+cBtwFbgOur6pFpxtkV5qOdzEu7mI/22F1y4Zf0SJLUIe74JUnqEAu/JEkdYuGXJKlDLPySJHWIhV+SpA6x8EuS1CEWfmnGJdmR5KG+24JfZbySktzYfEf8vU1sTyXZ3hfroUPGXZzkonlta5Jsbu7/Lsnqyf8PpN2Pf8cvzbgkz1fVvmOec8/mwiS7MsdbgYur6tS+trOBNVV13ghjb66qI/ravgM8V1WXJFkL7F9V39qVGKVZ5I5f6qgkTyb5ZpIHkvwpyZub9lVJrk5yX5IHk5zStJ+d5IYkG4HfJNkjyfeSPJLk1iS/THJGkuOT3Nz3OB9JctOAEM4CbhkhzpOS3NPEeV2SVc2Vz15KckzTJ8DHgWubYbcAn9qV50eaVRZ+afbtM++l/jP7jj1bVe8ArgC+0rR9Hbi9qt4JfBD4dpJVzbF3A5+rqg8BpwGHAkcBn2+OAdwOvCXJAc3P5wDrBsR1HLBpocCTHAhcABzfxLkZ+FJzeAO9657PzbWtqv4MUFXPAq9K8uqF5pe6yK/llWbfi1X19iHH5nbim+gVcoATgJOTzJ0I7A28vrn/26r6e3P/vcANVfUf4Jkkd0Dv+0aTXAN8Osk6eicEnx3w2K8Fti8S+3uAI4G7e5t69gLuao5tAO5M8lV6JwAb5o3d3jzGPxZ5DKlTLPxSt/2r+XcHL/8+CHB6VT3W3zHJu4AX+psWmHcdsBF4id7JwaDPA7xI76RiIQF+XVWfmX+gqp5Msg14H3AqcMy8Lns3jyGpjy/1S5rvNuD85n1zkhw9pN9dwOnNe/2vAT4wd6CqttH7zvELgfVDxm8B3rhILHcD709yWBPLqiSH9x3fAFwGbKmqZ+Yak+wB7A88vcj8UudY+KXZN/89/ksX6X8R8Epgc5KHm58H+SmwFXgY+AFwL/DPvuM/AZ6uqkeHjP8FfScLg1TVX4G1wHVJ/kjvROCIvi7XA2/j5Q/1zTkWuKuqdiw0v9RF/jmfpGVLsm9VPZ9kP+APwHFzO+8k3wUerKqrhozdB7ijGTPWAp3kcnrfe37nOOeVZoHv8UvaFbc2n5zfC7ior+hvovd5gC8PG1hVLyb5BnAQ8NSY43rQoi8N5o5fkqQO8T1+SZI6xMIvSVKHWPglSeoQC78kSR1i4ZckqUMs/JIkdch/Abh9IXvIJY/sAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# display spectrum and flux points\n",
"energy_range = [0.01, 120] * u.TeV\n",
"plt.figure(figsize=(8, 6))\n",
"ax = crab_spec.plot(energy_range=energy_range, energy_power=2, label=\"Model\")\n",
"crab_spec.plot_error(ax=ax, energy_range=energy_range, energy_power=2)\n",
"flux_points_fermi.plot(ax=ax, energy_power=2, label=\"Fermi-LAT\")\n",
"flux_points_hess.plot(ax=ax, energy_power=2, label=\"HESS\")\n",
"flux_points_hawc.plot(ax=ax, energy_power=2, label=\"HAWC\")\n",
"plt.legend();"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.0"
},
"nbsphinx": {
"orphan": true
}
},
"nbformat": 4,
"nbformat_minor": 4
}