\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.17?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": [
"## Prerequisites\n",
"\n",
"- Handling of Fermi-LAT data with gammapy [see the corresponding tutorial](fermi_lat.ipynb)\n",
"- Knowledge of spectral analysis to produce 1D On-Off datasets, [see the following tutorial](spectrum_analysis.ipynb)\n",
"- Using flux points to directly fit a model (without forward-folding) [see the SED fitting tutorial](sed_fitting.ipynb)\n",
"\n",
"## Context\n",
"\n",
"Some science studies require to combine heterogeneous data from various instruments to extract physical informations. In particular, it is often useful to add flux measurements of a source at different energies to an analysis to better constrain the wide-band spectral parameters. This can be done using a joint fit of heterogeneous datasets.\n",
" \n",
"**Objectives: Constrain the spectral parameters of the gamma-ray emission from the Crab nebula between 10 GeV and 100 TeV, using a 3D Fermi dataset, a H.E.S.S. reduced spectrum and HAWC flux points.**\n",
"\n",
"## Proposed approach\n",
"\n",
"This tutorial illustrates how to perfom a joint modeling and fitting of the Crab Nebula spectrum using different datasets.\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",
"## The setup\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\n",
"from gammapy.datasets import Datasets, FluxPointsDataset, SpectrumDatasetOnOff\n",
"from gammapy.estimators import (\n",
" FluxPoints,\n",
" FluxPointsEstimator,\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": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Models\n",
"\n",
"Component 0: SkyModel\n",
"\n",
" Name : Crab Nebula\n",
" Datasets names : None\n",
" Spectral model type : LogParabolaSpectralModel\n",
" Spatial model type : PointSpatialModel\n",
" Temporal model type : None\n",
" Parameters:\n",
" amplitude : 3.42e-01 1 / (cm2 s TeV)\n",
" reference (frozen) : 0.000 TeV \n",
" alpha : 2.511 \n",
" beta : -0.022 \n",
" lon_0 (frozen) : 83.633 deg \n",
" lat_0 (frozen) : 22.020 deg \n",
"\n",
"Component 1: BackgroundModel\n",
"\n",
" Name : background\n",
" Datasets names : ['Fermi-LAT']\n",
" Parameters:\n",
" norm : 0.954 \n",
" tilt (frozen) : 0.000 \n",
" reference (frozen) : 1.000 TeV \n",
"\n",
"\n"
]
}
],
"source": [
"print(datasets[0].models)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We get the Crab model in order to share it with the other datasets"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"LogParabolaSpectralModel\n",
"\n",
" name value unit min max frozen error \n",
"--------- ---------- -------------- --- --- ------ ---------\n",
"amplitude 3.415e-01 cm-2 s-1 TeV-1 nan nan False 0.000e+00\n",
"reference 5.055e-05 TeV nan nan True 0.000e+00\n",
" alpha 2.511e+00 nan nan False 0.000e+00\n",
" beta -2.248e-02 nan nan False 0.000e+00\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": 5,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/terrier/Code/gammapy-dev/gammapy-docs/build/v0.17/gammapy/gammapy/utils/interpolation.py:163: RuntimeWarning: overflow encountered in log\n",
" return np.log(values)\n"
]
}
],
"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(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": 6,
"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": 7,
"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": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Models\n",
"\n",
"Component 0: SkyModel\n",
"\n",
" Name : Crab Nebula\n",
" Datasets names : None\n",
" Spectral model type : LogParabolaSpectralModel\n",
" Spatial model type : PointSpatialModel\n",
" Temporal model type : None\n",
" Parameters:\n",
" amplitude : 3.42e-01 1 / (cm2 s TeV)\n",
" reference (frozen) : 0.000 TeV \n",
" alpha : 2.511 \n",
" beta : -0.022 \n",
" lon_0 (frozen) : 83.633 deg \n",
" lat_0 (frozen) : 22.020 deg \n",
"\n",
"Component 1: BackgroundModel\n",
"\n",
" Name : background\n",
" Datasets names : ['Fermi-LAT']\n",
" Parameters:\n",
" norm : 0.954 \n",
" tilt (frozen) : 0.000 \n",
" reference (frozen) : 1.000 TeV \n",
"\n",
"\n"
]
}
],
"source": [
"print(datasets[\"HESS\"].models)"
]
},
{
"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": 9,
"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 : 829\n",
"\ttotal stat : -14819.11\n",
"\n",
" name value unit min max frozen error \n",
"--------- --------- -------------- ---------- --------- ------ ---------\n",
"amplitude 4.015e-03 cm-2 s-1 TeV-1 nan nan False 1.072e-04\n",
"reference 5.055e-05 TeV nan nan True 0.000e+00\n",
" alpha 1.260e+00 nan nan False 7.101e-04\n",
" beta 6.185e-02 nan nan False 1.962e-04\n",
" lon_0 8.363e+01 deg nan nan True 0.000e+00\n",
" lat_0 2.202e+01 deg -9.000e+01 9.000e+01 True 0.000e+00\n",
" norm 9.839e-01 0.000e+00 nan False 3.239e-01\n",
" tilt 0.000e+00 nan nan True 0.000e+00\n",
"reference 1.000e+00 TeV nan nan True 0.000e+00\n",
"CPU times: user 12.6 s, sys: 105 ms, total: 12.7 s\n",
"Wall time: 12.7 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": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"LogParabolaSpectralModel\n",
"\n",
" name value unit min max frozen error \n",
"--------- --------- -------------- --- --- ------ ---------\n",
"amplitude 4.015e-03 cm-2 s-1 TeV-1 nan nan False 1.072e-04\n",
"reference 5.055e-05 TeV nan nan True 0.000e+00\n",
" alpha 1.260e+00 nan nan False 7.101e-04\n",
" beta 6.185e-02 nan nan False 1.962e-04\n"
]
}
],
"source": [
"crab_spec = datasets[0].models[\"Crab Nebula\"].spectral_model\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": 11,
"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",
" e_edges=e_edges, source=\"Crab Nebula\"\n",
").run([dataset_fermi])\n",
"\n",
"\n",
"e_edges = MapAxis.from_bounds(1, 15, nbin=6, interp=\"log\", unit=\"TeV\").edges\n",
"flux_points_hess = FluxPointsEstimator(\n",
" e_edges=e_edges, source=\"Crab Nebula\"\n",
").run([dataset_hess])"
]
},
{
"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": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAF3CAYAAABE0Ck1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd3xUVf7/8dcnnQQSSkKHEIoQmpQEaYIVRUWxIrKIIKIu4C7uul/dtaA/Xdd11VXEggWsYC/YG0gVCIiABKQjvYQWSkKS8/sjmAUEEsLM3GTm/Xw85oFz5s6dd7gOn5x77znHnHOIiIhIaAjzOoCIiIgEjgq/iIhICFHhFxERCSEq/CIiIiFEhV9ERCSEqPCLiIiEkAivAwRCYmKia9CggdcxREREAmbu3LnbnHNJR7eHROFv0KABGRkZXscQEREJGDNbc6x2neoXEREJISr8IiIiIUSFX0REJISExDV+ERHx3sGDB1m3bh0HDhzwOkpQiYmJoW7dukRGRpZoexV+EREJiHXr1lGpUiUaNGiAmXkdJyg459i+fTvr1q0jJSWlRO/RqX4REQmIAwcOUK1aNRV9HzIzqlWrdlJnUcp84Tezhmb2kpm9e1hbnJm9YmYvmFk/L/OJiEjJqej73sn+nfq18JvZy2a2xcwWHdV+oZktNbPlZnbnifbhnFvpnLvxqOYrgHedczcBl/o4toiIBCkzo3///kXP8/LySEpK4pJLLjmp/TRo0IBt27ad8jZe8HePfxxw4eENZhYOjAZ6As2BvmbW3MxamdknRz2qH2e/dYFfD/13vp+yi4hIkImLi2PRokXs378fgK+//po6dep4nCqw/Fr4nXNTgKyjmjsAyw/15HOBCcBlzrmFzrlLjnpsOc6u11FY/KEcXK4QEZGyo2fPnnz66acAjB8/nr59+xa9lpWVRe/evWndujUdO3ZkwYIFAGzfvp0ePXrQtm1bbr75ZpxzRe95/fXX6dChA23atOHmm28mP79s90e9uKu/Dv/rrUNhET/jeBubWTXgIaCtmd3lnHsYeB942swuBiYe531DgCEA9evX91F0ERHxhfsn/sziDbt9us/mteO5r1eLYre79tpreeCBB7jkkktYsGABgwYNYurUqQDcd999tG3blg8//JDvvvuO66+/nvnz53P//ffTtWtX7r33Xj799FPGjBkDQGZmJm+99RbTp08nMjKSP/7xj7zxxhtcf/31Pv3ZfMmLwn+suxDcMdoKX3BuO3DLUW17gYEn+hDn3BhgDEBaWtpx918mTXoYzr7L6xQiIkGpdevWrF69mvHjx3PRRRcd8dq0adN47733ADjnnHPYvn07u3btYsqUKbz//vsAXHzxxVSpUgWAb7/9lrlz55Keng7A/v37qV79eFepywYvCv86oN5hz+sCGzzIUWY45ygoKMA5h3OOyO//xYFOI4peO9pvd3CaWdEjLCys6E8RkbKuJD1zf7r00kv561//yuTJk9m+fXtRe3H/5h7NOceAAQN4+OGH/RfWx7wo/HOAJmaWAqwHrgWu8yDHqRl78RFPHQ5XUFi49+aHkZUfw468GHbkx5CVH8Pugmj25EexuyCa7IJo9rpI9hVEss9Fsd9FkuMiyHER5OUXMDUaOj0wkXxnFBBGPoY77ERJ4TNHOI5IyyecAiKtgCjLI8byiCKPCmF5xFousWEHiQ07SCXLxZpdRHyFSOJjIqkcG0W1StEkVowhKb4CCXExREZGaqiNiISEQYMGkZCQQKtWrZg8eXJRe7du3XjjjTe45557mDx5MomJicTHxxe133333Xz++efs2LEDgHPPPZfLLruMESNGUL16dbKystizZw/Jycke/WTF82vhN7PxwFlAopmtA+5zzr1kZsOAL4Fw4GXn3M/+zOFLBw7mM+OXjWzd1oDN+RXZlBfHlryKbCuIY3t+HNsLYtnnoo/7/gjyiQ87QKzlUsEKi3IFO8jAsE/pzeSiI/Jj1GAAvrOOTInoRNihqyG//QJQ4Iw8wsg79OdBF06uiyDHhXPARXDARbKzIJ59B6PY66LYWxBFzpytx80VGQYJMWFUrRBB1dgIEuMiqREfTa2ECtSuHEvdanEkJ1YiPq6CfjkQkXKvbt26/OlPf/pd+8iRIxk4cCCtW7cmNjaWV155BSi89t+3b1/atWtH9+7di+4da968OQ8++CA9evSgoKCAyMhIRo8eXaYLvx3rtEawSUtLcxkZGT7Z1/qd++nyr++KnsdGGFUrhFGlQjiVY8KpGhtBtbgoqsRGUjUuiipxUVSNi6ZyXBSVY6OJjY4kIiKC8PBwwsLCiv48wsgEGLnLJ3kLCgooKCggPz+f/TkH2bkvh137ctmxN4ft2Tls35vDjr0H2b43l6y9uWzfl8eOfXlk7c9nd+7v/9+oHBNGjbgIasZHUa9KBVIS42iYVInGNROoU62SLjWIyHFlZmaSmprqdYygdKy/WzOb65xLO3rboJ6r38x6Ab0aN27ss33WqBTNy/3bEGcHqV0llsoVKxAZGUlkZCTh4eE++xxfCQsLIywsjIiICKKjo6kcX7FE78vPzyd73wHW79jLuqy9rN+xj/U79rN+53427s5lyZb9TF2dTYH731mE2EijbnwkyVVjaJgYR7Oa8bSoV5VGNSoTHq5fCEREyoKgLvzOuYnAxLS0tJt8tc+I8DDOaeHnyR66n3Ayw4AIDw8noVIcCZXiaH6c0ZD7DuSwestuVmzZzcqt2azcupc1Ow4wb102Xy/bDWwEIDrcSK4cSeOkWJrXiqdl3cqcnpxElYoxgfuBREQECPLCX26Vk6F8sTHRNK+fRPP6Sb97LWv3PhavzyJzw06WbtrD8m37mbZqN58t2QmsBaBWpQiaVY+lZe142iZXJa1RDeIrRAX4pxARCS0q/OIXVeNj6RofS9fUukVtzjnWbt3NT2u3s2j9ThZv3MOijXuZtGI3TF1HGFC/ciStalekbf0qdGxcg9Q6VXQzoYiID6nwS8CYGcnVE0iunsClh91usmlHNnNWbuHHNTtYsH433y3bycTFO4CVVIoOo1XNWNrXr0znJtVJa1SDSN0vICJSair84rmaVSrSq31FerUvfF5Q4FiyPosfVmxhzuosflqfzYw16xg1dR0xEUarmrGkN6jMmU1rkt6wOhH6RUAkeP02Z8rAT73NEURU+KXMCQszmterRvN61Rh0qG3Tzr1MXbKRmSu2MffX3TwzbT3PTFtPbKTRpnYcXRoncm6LOjStlaBLAyJyXOHh4bRq1aro+YcffkiDBg18+hkff/wxixcv5s47f3+jdoMGDcjIyCAxMfF3r1122WVs2bKFmTNnAvDQQw/xzjvvALBw4cKi3IMGDeK2224rdb6gHsd/2HC+m5YtW+Z1HPGhzbv28X3mBqb+spU5a3ezKTsPgOpxEXRskMDZzWpwTss6JOhmQZEyo1Tj+MdeDDvXwIhFPslQsWJFsrOzT/p9+fn5PhmyfbzCv3PnTlq1akXFihX57LPPSElJOeL14nKfzDj+oD5H6pyb6JwbkpCQ4HUU8bEaCbFc07Exo67vxA93X8B3t5/J33s0pGn1Cny9NIsR7y2m3QNf0/upSTz15SKWbfLNhEgi4oFdvxa/zSnIz8/njjvuID09ndatW/P8888DMHnyZM4++2yuu+46WrVqxerVq2nWrBmDBw+mZcuW9OvXj2+++YYuXbrQpEkTZs+eDcC4ceMYNmzYSWV477336NWrF9deey0TJkzw+c94OJ3ql6DQsHo8Q86JZ8g5qRzML+CHZZv5+ucNTF2exeOT1vD4pDXUiY/krCZVuej0enRsXJ3wMF0SEClzjloHhU0Ljt1eymv++/fvp02bNgCkpKTwwQcf8NJLL5GQkMCcOXPIycmhS5cu9OjRA4DZs2ezaNEiUlJSWL16NcuXL+edd95hzJgxpKen8+abbzJt2jQ+/vhj/vnPf/Lhhx+WKtf48eO57777qFGjBldddRV33eW/Yd0q/BJ0IsPDOLNZLc5sVguAtduy+fyntXy7ZCtv/biZN+ZuJiEmjDMbVuGSNnU4p3kdoiKC+uSXSPmzc82RPf010wr/TKgHlUs/D36FChWYP3/+EW1fffUVCxYs4N133wVg165dLFu2jKioKDp06HDEafeUlJSia+0tWrTg3HPPxcyKzgiUxubNm1m+fDldu3bFzIiIiGDRokW0bNmydD9kMVT4JejVT6zIzec25+ZzYff+XL5auI4vFm7ku2VZfLJ4O3GRCzmzUWV6tanLeS3rEB1R9qZeFgkZR/fkx15cWPR9tH7JsTjnGDVqFBdccMER7ZMnTyYuLu6Itujo/y3CFhYWVvQ8LCyMvLy8I7bNz8+nffvC4UqXXnopDzzwwDE//6233mLHjh1Fv2Ds3r2bCRMm8OCDD57aD3YcKvwSUuIrRHFVh4Zc1aEhOXn5fPfzBj75aR3fL9/JF0t2EBu5iO6Nq3B5+/qcnVpLcwaIhIALLriAZ599lnPOOYfIyEh++eUX6tQ59anZw8PDf3d24VjGjx/PF198QadOnQBYtWoV559/vgq/iK9FR4TT8/R69Dy9HgfzC5icuZEP5/3K5GU7+Dwzi0rRCzi/aTWu6ZBCh4aJhOmeABFvJNTz6+4HDx7M6tWradeuHc45kpKSSn2tviRat25dtJJphw4dWLt2LR07dix6PSUlhfj4eGbNmsUZZ5zh888P6uF8v/HlsrxyDD4ebuO1nLx8vl20gffnrWXqil3k5DtqVozgopY1uK5TQxrXiPc6oki5VOrhfKAJfIqhZXkP8ceyvHIcfh5uE0jREeFc1KYeF7WpR3ZOHhPnreGDH9cxbtZ6Xv5hPc1rVODKdnW5Kj2FhNhIr+OKBDcVfJ9Tjz+UHD0cxlc2LYCc3ZDc1T/7LyNf/M279vPWrJW8/+NGVu/IITIMzm5ShT90bkTXJtV1KUCkGKXq8UuJqMcvgeGn4TZlVY2ECtzWowXDz2/OT7/u4I0ZK/l88Va+WppBrUqRXNG2Fv27NKFmQozXUUVEjkuFP5T4q+ccgOE2ZYmZ0aZ+VdrUr8r/O5jPxB/XMmH2WkZPWcuzU9dyZsPKXN+lIWc3q6mzACJS5qjwy6kb+CmMDM1pkWMiw7m6QwpXd0hh9bZsxk1dzgc/beL7V+dRs1IkfdLqMKBrE6rGac0AESkbNEhZfKP771ehCjUNEisy8vI2zL77fJ64uiW14qN5ctJqznjoG4a++gM/rsnyOqJIuTPwi4EM/GKg1zGCigq/+MbZ/ptXuryJjgjn8vbJfDC8O1/++Ux6t67Od79kcfmzM+n5xCTez1jDwfwCr2OKhKSKFSse8fzwBXVGjhxJnTp1aNOmTdFj586d7Nu3j379+tGqVStatmxJ165di1bKe+ihh2jRogWtW7emTZs2zJo1K+A/08nSqX4RP2paM55Hr03j3gMHmfDDKl6ZuYbb313EQ58t4boOdRl4pi4DiBRnffb6gH3WiBEj+Otf/3pE28MPP0yNGjVYuHAhAEuXLiUyMpKZM2fyySefMG/ePKKjo9m2bRu5ubkBy1paQV34NY5fyopKMZHcdNZp3NitCd9mbuKF75cxavJqxkxdQ69W1bn1nKY0ql7J65giZdLGvRu9/fyNG0lO/t9IpaZNmxa1JyYmFs3Xn5iY6Em+k6Vx/CIeWbppN89+u5RPft5CfgF0bViZoec25YyG1TDTaAAJPiUZx3/09fwlWUvIPphNWo0jh6OPvXBsqTKEh4cXra4HkJWVxaWXXsrTTz/NyJEjeeGFF0hKSgKgSpUqTJo0ifnz59OjRw8aNWrEueeey4ABA2jSpAnZ2dl07dqVffv2cd5559GnTx+6d+9eqlyn6mTG8esav4hHmtaM57/90pl513kM6Vqfn9bv4doXZtHrye/5fOEGCgqC/5dykeNZn72ejM0ZZB8svJaesTmDjM0Zp3za/7dleX97HL1i3ogRI4pemzRpEgBt2rRh5cqV3HHHHWRlZZGenk5mZiYVK1Zk7ty5jBkzhqSkJPr06cO4ceNOKV8gBPWpfpHyIKlSNHdd0oo/9UjljZmreHHqKm5940eSqyzm1rMac2Vafa0SKCHj6J78wC8GkrE5g4UDFnqUqFDFihW54ooruOKKKwgLC+Ozzz4jNTWV8PBwzjrrLM466yxatWrFK6+8wg033OBp1uLoXxORMiI2KoKbujdh2l3n8dhVrYgIM+784GfO/Nc3jJ26ggMH872OKBKSpk+fzo4dOwDIzc1l8eLFJCcns3TpUpYtW1a03fz584+4F6CsUo9fpIyJDA/jyrT6XNG+Ht8s3sSTXy/l/k+X8PSk5QzumsINXRtRISrc65giAVMrrlbAPuuJJ57g9ddfL3r+4YcfsmLFCm699VaccxQUFHDxxRdz5ZVXMm/ePIYPH87OnTuJiIigcePGjBkzJmBZS0s394mUcc45pi/fxhNfZTL31z1UqRDO4K4pDDyzEbFR+t1dyo/SLNLz281+pb2ZL1RokR6RIGJmdG2SRNcmSfywYhuPf5nJo18v54Wpq7jpzBQGndlYZwAkaKng+56u8YuUIx0bJfL2H8/k3Vs60axGHI9+vZzOD3/DmMnLdA+AiJSICr9IOZTWoCoTbj2Tt2/uSMPECvzzi1/o+q9vGDdtBbl5ZXQ64LEXFz5ExFMq/CLlWIeUarw3tBtv3NiB2gkxjPxkCWf9+1vezVhLvuYBEJFjCOrCb2a9zGzMrl2hsU68hK4uTZL4aHg3XhrQnriocP767kLOf+w7vv55I6FwA6+IlFxQF37n3ETn3JCEhNBcK15Ci5lxbmpNvrz9bJ7s05rcvAJuem0elz89hbmrtSSwlE9r+l/Pmv7Xex0jqAR14RcJRWFhxmVt6zHpb+dyf69mrM06wJXPzWTQyzNZsTXb63ginjrRsry/Of300+nbt2/R859++ok2bdoUPR8/fjyxsbEcPHgQgIULF9K6dWsADh48yJ133kmTJk1o2bIlHTp04PPPP/fXj1MqKvwiQSoyPIwBXRox9c5zue3shsxctZPzH/+eu979ke3ZOV7HEymx3PWBW5Y3MzOTgoICpkyZwt69ewFo1aoVa9asYc+ePQDMmDGDZs2a8eOPPxY979KlCwD33HMPGzduZNGiRSxatIiJEycWva+sUOEXCXJx0RHcfkEqU//vHK5uV5u35m6g2yPf8fS3SzUEUMqFvA0bAvZZb775Jv3796dHjx58/PHHAISFhZGens6sWbMAmDt3LkOHDmXGjBlAYeHv3Lkz+/bt44UXXmDUqFFFS/XWqFGDa665JmD5S0KFXyREJFaM5pGr2/LViG60q5/Af75ezln//paJ89frBkAJGfv376dNmzZFj3vvvfeI19966y369OlD3759GT9+fFF7586dmTFjBnv37iUsLIyzzjrriMLfpUsXli9fTv369YmPjw/oz3SyVPhFQkzj6pV47abOvDG4AxWjIxg+YT5XjJ7KovUa/SJlx9ZRT5PZLJV9c+YAkNkslcxmqWwd9fQp7fdEy/LOmTOHpKQkkpOTOffcc5k3b17R4jxdunRhxowZzJ49m/T0dBo1asTy5cvZunUr2dnZNGzY8JRyBZIKv0iI6tI4iS9vP5v/d2lzVm3bR69R0xgxfi5b9+j6v3gvafgwUpdkEpueDkDqkkxSl2SSNHxYMe8svfHjx7NkyRIaNGhAo0aN2L17N++99x4AHTt2ZM6cOUybNo1OnToBULduXSZMmEDnzp0BaNy4MWvXri1z1/SPpsIvEsLCw4z+nVOYcue5DOhYj4kLN3HWo98x5vvlHMwvozMAivhBQUEB77zzDgsWLGD16tWsXr2ajz76qOh0f6VKlahXrx7jxo0rKvydOnXiv//9b1Hhj42N5cYbb+S2224jNzcXgI0bNx6x2l9ZoMIvIsTHRDKyd2u+HNGd0+vE88/Pl9LjsUlMW7bV62giRNSu7ffPmDJlCnXq1KFOnTpFbd26dWPx4sVs3LgRKDzdn5OTQ7169YDCwr9y5cqiwg/w4IMPkpSURPPmzWnZsiW9e/cmKSnJ7/lPhpblFZEjOOf4evEmRn60iA27c7kgNZH7e59OzYSYU9vxb/P0D/z01ENKuVSaZXl/m7wn+bVX/REpaGhZXhEpNTOjR4tadDutOs989wvPTVnF949+x5/Oaczg7o2JDNeJQgkcFXzf0zdYRI4pJjKc2y9I5Zvbz6JDcmUe+WoZPR6bxKyV272OJiKnIKgLvxbpETl19avF8upNnRnTvx37D+bTZ8wPjBg/l6y9uSe/s51rfB9QRE5KUBd+LdIj4js9WtTiuzvO4aYuyXy8YBNnP/odb81Zc3KT/+z61X8BpVwIhfvKAu1k/051jV9EjvTbTXjHEAv8A7gyuSr/2NSN/3svnw8+/ZyHa00hJep/Z9a2TssiqWvVI9+8aUGx+z8lummwzIuJiWH79u1Uq1YNM/M6TlBwzrF9+3ZiYkp+860Kv4ictGbRWbxT/0PG70zlX1s7ccGqPgyvNodbqv1EpBWwbfrO/xX+nWuO7OmvmVb4Z0I9qJwc+PDimbp167Ju3Tq2btUwUV+KiYmhbt26Jd5ew/lE5JRs2X2Af7z/E18v2UbjxAr8dfHHJM/4ktQlmUduOPbiwqI/UvfciATC8YbzBfU1fhHxv+rxMbxwwxlMiFrAqBeHkjzjS8B3c6uLiG+p8IuIT3R84C7qzF/AM+fdDMDgAaPY+enUI+dWT6jnUToR+Y0Kv4j4THxMJKOe/nPR874vzOL/3vmR7Jy8wgZd0xfxnAq/iPhc4tChfPOXs7mhYz3enruB8/4ziWl76xT/xlIY+MVABn4x0C/7FglGKvwi4nNJw4dRISqckb1b884tnYiKMP6w7jLu2njm/3r/IuIJFX4R8au0BlX56vazGVz5RybsbkmPxyZx+YQ71EsX8YgKv4j4XUxkOHfXmMnb9d7HzPhxfnd+WdaKAwfzvY4mEnJU+EUkYNJjN/P1X86idq0VrF/flAufmMyCdTu9jiUSUlT4RSSgYqMiaNp0Pqe3nsK+3HwuHz2dJ75aQl5+Qan3uT57vQ8TigQ3FX4R8UTVqlv4+i9nc2GL6jz53Qp6Pz2VVdv2lmpfG/du9HE6keClufpFxBPrs9eTUCGS0X9Ip8eP6/jHBwvp+d/vufviVPp1bFC0iEtxNwEuyVpSou0Axl449tSDi5RzQd3jN7NeZjZm1y7NDS5S1hzeS7+sbV2++cvZnF43nrs/WsyNY2eRtTf3hO9fn72ejM0ZZB/MBiBjcwYZmzN02l+kGFqkR0R86rg9700LC/+s2YolWUvIPphNWo0j1w9xDn79tQkrV7UkIiKH5qkZfHTdI8V+XsbmDBYOWOiL+EX5dXZAyjst0iMinltP3gl76WZQv/4y2rf7joiIg/y0oBsPfLyQ3LzS3/gnIkfSNX4R8anj9pTHXlz45w2flqiXvj83n5EfLeDlGWv5YcV2numfToPEuGNuWyuu1qnGFgkZ6vGLSJlUISqcR65uy7P92vLrzgP0fHIK7839tej1w5f7rVPRP+sAiAQjFX4R8URJe+k9W9XmixHdaVYjjr+8s4Db3sxgb04e20aP9nNCkeCkwi8injiZXnqdyhV459auDDsrhYkLNnPh39/2YzKR4KZr/CJSLkSEhzFg2Xdc/OH/evqZzVL5GzD9gjpwoXfZRMoTFX4RCYyBn57yLpKGDyNp+DC2ZeewNa0NPXv/h57Nk3i0TzsfBPwfzQUgwUyn+kWk3EmsGA3A7ec15svMrVz038ks2bTbZ/vXFMASzDSBj4iUS1tHPU3S8GHMXLGNYW/MZU9OPg9c2oJrz0g+7ntKMq3v8SYXOhZN8iNlmSbwEZGgkjR8GACdGiXyxYizaFM3njs/WMSI8XM5cDD/pPenKYAlVKjHLyJBIb/A8diXmTzz/SqaJFXgxRvOILnasSf8ORFfTwEs4hX1+EUkqIWHGX/r2ZyXBrRn0+5cLnpyCl/9rGv1IkdT4ReRoHJuak0++1M36leJYchr8/jnJ4vILzi5M5uaAliCmQq/iASdelVj+WBYN65qW4sx09bQb8wMdhSzzO/hNAWwBDMVfhEJSjGR4fynTzse6t2CuWt30fO/k1m0fpfXsUQ8p8IvIkGtX8cGvHNLJwocXPHMdN6Zs8brSCKeUuEXkaDXpn4VPvtzd06vU4k73lvEPR/8xMH8Aq9jiXhChV9EQkJixWjG39yFAR3r8dqsdVz3/HS2Z+d4HUsk4FT4RSRkRISHcX/v1jx6ZSvmr9/DxU9+z6L1O/36mQO/GFiiGQNFAkWFX0RCztXp9Xn30HX/K5+Zwcfz13kdSSRgtDqfiISk0+tV4dM/dWfw2B+4bcJPZG7YxR0XNi/z8+//dvagrOeUsiuoe/xm1svMxuzapSE8IvJ7SZWiefuPXbmiTU2enbKaG8f+QHZOHlC4CJBIMArqwu+cm+icG5KQkOB1FBEpo6IjwnmsTzvuvqgp3y/P4rKnvufXrH1sGz3a62gifhHUhV9EpCTMjMHdGjNuYDqb9+RyycOfeR1JxG90jV9E5JDUL9/mnbf/19PPbJYKQOLQoUXLAIuUd+rxi4gckjR8GKlLMqk17ycAevb+D+/++22qDR3qcTIR31HhFxE5SuXYKAD6tK/NSzPWMnjcLPbl5nmcSsQ3VPhFRI4hcehQ/nVVG/7RsymTftnOFaOnsnn3Aa9jiZwyFX4RkWNIGj4MM+Om7o15/g/tWL19P5c8OYWfN2h4sJRvKvwiIsXo0bIW793aBYCrnp3Ot4s3ndT712ev90cskVIpUeE3szZmNtzMHjaze83sCjPT4HgRCRkt6iTwyZ+6Ub9KBW56bS6vTF9R4vdu3LvRj8lETs4Jh/OZ2R+AEcA6YC6wBogBzgPuMbN5wH3OOU10LSJBr0Z8DO8PPZNbXp3NfROX8MKPH9K40SLMjv+eJVlLAEq0UI+m4ZVAKG4cfzWgm3Nu77FeNLM0IJXCXwxERIJeXHQEYwd15O73f2LCXDhwIJbmqRmEhxccsd367PVH9PQzNjOy6rsAACAASURBVGcAUCuuFnUq1gloZpHDnbDwO+eeLOb1DN/GEREp+yLCw3j4qjY0SIzjkS9h14rmjLuxE1Xion637cAvBpKxOYOFAxb67PN1z4CcimKv8ZtZlJn1NrPHzGy8mb1sZrebWbNABBQRKYvMjFvOPo0nrz2dxZv2ctmoKfyatS8gn617BuRUFHeN/27gSmAKhdf4v6bwGv9pwBNmZsBfnXOL/B1URKQsurRNXapXiuGmV+dw6agpvHJjR1rXrXzENrXiapVoXyW5D0D3DMipKu4a/0Ln3IPHee3fZlYLqOfjTCIi5UrHRom8/8eu9H9xJtc8N4Nn+rXjnNSaRa/74pq+7hkQXynuGv9HR7cd6uXHOuf2Ouc2AjrnJCIhr0mNSnw8vBv9XpjBTa/O5eHLW3JNh+ST2kdJeuj+uGdAQktJx/G/ambxZhYL/AysMrPb/RtNRKR8qX5ouF/7evH87f1FjPpmCc55nUrkSCWdua+Vc2430Bv4CqgL3OCvUCIi5VWlmEheH9KFns2TeOybFfyyrLXPi39J7xkQOZaSFv4oM4sALgM+dM7lAgXFvEdEJCRFRYQx+g/pDDijLhs2nMain9PJzfPdP5m6pi+noqSF/0VgLVAF+N7M6gPZfkslIlLOhYUZ919+Og1TfmLbtmT6vzCd7Bwt7SveK1Hhd8494Zyr7Zzr4ZxzFM7Ud45/o4mIlH/Jycto2nQ2c9bu5upnprI9O8frSBLiSrU6n3Ou4NDpfhERKUbtWmt5tl9bVmzbT++np7BuR2Am+hE5Fi3LKyISABe0rM3rg89gx748Ln96Kss27/E6koQoFX4RkQDpkFKNt2/pTF6B48pnpjN/7Q6vI0kIKslc/XFm9rtZKMyshX8iiYgEr+a1E3j/j12JjQqj75iZzL3/Ea8jSYg5YeE3syuB5cCnZrbQzNod9vJrfk0mIhIExl449ncz8qUkVeTDYd2oGR9F7PhxfL5Aq+1J4BTX478HSHPOtQRuBsab2aWHXjO/JhMRCWI1E2J4atEEAIaOn8+7c9Z4nEhCRXGL9IQ559YDOOdmmNk5wCdmVg/QRJQiIqWwddTTbBs9uugf4E8/+Ct8AF/36sv5j97raTYJfsX1+PeaWcpvTw79EnAWcDXQ3I+5RESCVtLwYaQuySR1SSYAKQsX8dDtL3JTeHtGfbPE43QS7Ior/EM56qyAc24X0IPCU/8iInKKYiLDeXlQJ85vWo3HvlnBw58swml1H/GTExZ+59w859yy356bWayZxQMxwAf+DiciEuwShw4FCuf3f27AGfRuXZ3np63hvg9/UvEXvyjuGj8AZjYY+H9APoWL8xiF1/jr+y+aiEjwSxo+rOi/w8OMx69NIzZqPq/OWs/+3HweubodYWG6l1p8p0SFH/g/4HTn3BZ/hhERCXVhYcZDV7YpPP0/81cOHMzgib7tiQgvPEF79NBAkZNV0pn7VgK7/RlEREQKmRn3XNqKYd0bMHHRFv742mwO5msldPGNkvb47wSmm9kPQNHSUs652/2SSkQkxJkZf+3ZgujIcB77ZgVDxv3A8wM6EhWhmdbl1JT0/6DngOnAfODnwx4iIuJHw89rxl0XNGbSsh0MenkmOXn5XkeScq6kPf4C59xtfk1yHGbWEPgHkOCcu+p4bSIiwerms5sSFRHO/Z8uZcALMxg3uDMxkeFex5JyqqQ9/m/NbJCZJZlZ/G+P4t5kZi+b2RYzW3RU+4VmttTMlpvZnSfah3NupXPuxuLaRESC2cAzG/PQpanMWrOb/mOmsz9XPX8pnZL2+Acc+vP+w9pKMpxvHPA08OpvDWYWDowGzgfWAXPM7GMgHHj4qPcP0kgCEZFC/To3JCLcuPODxfxhzHReG9KZ2KiS/jMuUqhE/8c45+qVZufOuSlm1uCo5g7AcufcSgAzmwBc5px7GLikNJ9zLGY2BBgCUL++phsQkeDQ54wUDLjzw8X0e346rw/pQly0ir+UXIlO9ZvZLWZW+bDnVQ4V1tKoA/x62PN1h9qO99nVzOw5oK2Z3XW8tqM558Y459Kcc2lJSUmljCoiUvZcc0YK/76iJT9tyOa656exNyfP60hSjpT0Gv8tzrmdvz1xzu0Abi3lZx5rCqrjzkvpnNvunLvFOdfo0FmBY7aJiISSq9KTeeyqlizcuFfFX05KSQv/EbePmlkYEFnKz1wHHH7poC6woZT7EhEJWZe3T+axq1qp+MtJKWnh/9rMxptZdzPrBrwBfFPKz5wDNDGzFDOLAq4FPi7lvkREQtrl7esX9fz7PT+dfbkq/nJiJS38d1A4gc8I4C/ANOCvxb3JzMYDM4GmZrbOzG50zuUBw4AvgUzgbeecJgMSESmly9sn8+8rWrBgY7aKvxTLgnnZRzPrBfRq3LjxTcuWLSt2exGR8uyd2av5vw9+pk2dirwxpCsVojTJTygzs7nOubSj20/Y4zezD82sp5n9bqyImSWb2b1mNsiXQX3JOTfROTckISHB6ygiIn53dYcG/Kt3c35cl82AF6dz4KAm+ZHfK+5U/1AKJ9r5xcxmmtnHZvaVmS0HxgI/O+de9ntKEREpkWvOSOHBS5sxe+0ebnhxhub2l985YeF3zq13zt3unGsI9AceBf4OtHXOneOcey8QIUVEpOT6dW7E/Zc05Yc1u7nxpZnk5hWwddTTXseSMqLE6zs655Y756Y65zKcc3v8GUpERE7NgK6NubtnE6at2sXAEc+zbfRoryNJGaF5HkVEgtTg7qeRl+/411eFz/MLHOFhx5pDTUKJCr+ISJDaOuppuo8eTfdDz39p3hyAxKFDSRo+zLtg4qni7up/xsx6mVlsoAL50qHsY3bt2uV1FBGRgEsaPozUJZnEpqcD0LP3f3h55GskDhvqcTLxUnHX+N8EzgC+OXQ3/1/MrEUAcvmEhvOJiEDya4Uro99wRh3emb+Ze97/kWCew0VO7ISn+p1z0yicpQ8zqw5cAPzDzJoBGcAXzrn3/Z5SREROSeLQodzX+3Ry8vJ5fc5GKkZHcOclrb2OJR4o8TV+59wW4DXgNTMzoANwIaDCLyJSxv12Tf+hK9uxN2cOz037lbjoCIaf39zjZBJopbq5zxWeI5p16CEiIuVEWJjxxHXp7H/lBx77dhUVoiIY3P00r2NJAJV4HL+IiASH8DDjmevPoGtKAg9+vozxM1d6HUkCqLi7+rXCg4hIEIoMD+PFQZ1oX7ci//g4k4k/rvU6kgRIcT3+DWb2rJl1C0gaEREJmJjIcF69qQvNkiow4p2FTMrc6HUkCYDiCn8rYBHwkJmtNbP/mFn7AOTyCY3jFxE5sbjoCN68pSv1K0dz6xs/MnvlVq8jiZ8Vt0jPFufcaOfcmUAXYCPwnJn9Ymb3ByThKdA4fhGR4lWOjWLCLV2pWiGCQeMyWLRuh9eRxI9OZpGeX4FngSeAvRQu2SsiIkGgenwME27pTEyE0f/FWazaqrXYglWxhd/MoszscjN7G1gFXAzcC9TydzgREQmc+tUq8sZNnch3jr7Pz2DTrv1eRxI/KO6u/leBdcAACifqSXHO9Tt0Cv1gIAKKiEjgNK2VwNgb0tl1IJ++z01n575cryOJjxXX4/8eOM0519s5N8E5ty8QoURExDvtUxJ5tl9bft2VQ7/np7EvN8/rSOJDxd3c95JzbqeZJZnZ82b2KYCZNTezGwKSUEREAu6s1Fo8dmVLMrfs54YXZ5CbV+B1JPGRkt7cN47C3n/dQ8+XAX/xRyARESkbLmufzL0XncbstXsY9tosCgq0ol8wKGnhr+6cexMoADh0fT/fb6lERKRMuOHMJgzrVp+vlmbxj/fmaTnfIFDSRXr2mllVwAGYWTpQ5sd6mFkvoFfjxo29jiIiUm79pWdLtmXnMH7uJpIq/sztPVt6HUlOQUl7/H8FJgINzex7YDww3G+pfEQT+IiInDoz459Xtef806rw1PdreGXqMq8jySkobjhfRwDnXAZwNtAd+BPQ3Dk33//xRESkLAgLM0Zf35EO9Spy/2e/8Nn8X72OJKVUXI//md/+wzmX65z7yTk33zmngZ0iIiEmKiKMsYO70CQxhj+/s5BZy7d4HUlKocRT9oqIiMRFR/D6kC4kxkYw+NW5LN24s+i1Nf2vZ9k553qYTkqiuMLf0Mw+Pt4jIAlFRKRMSaoUwxtDOhEeBv1fnMXGnf+b2y1vwwYPk0lJFHdX/1bgsUAEERGR8iMlqRJjb0jnuhdnc93z0/nwtu5eR5ISKq7wZzvnvg9IEhERKVfaNkhk1LWn8/3dj7Bh7G1F7ZnNUgFIHDqUpOHDvIonx1Fc4V8VkBQiIlIundeyDlv+dhs9P76As7NX8bdvRpO6JNPrWHICxV3jf/xEL5pZvJlpJgcRkRB2XefGDOtWn0kVU7yOIiVQXOG/0sxmmNm9ZnaxmXUws25mNsjMXgM+ASoEIGepmFkvMxuza9cur6OIiAS1v/RsSa/dS3m96fk8/516/GWZFTfvsplVAa4CugC1gP1AJvCpc26a3xP6QFpamsvIyPA6hohIUMsvcFw/ZiozVu/hyWtacmm7ZK8jhTQzm+ucS/tdeygsuKDCLyISGPtz87ni6cms2J7D64PS6NCouteRQtbxCr8m8BEREZ+pEBXOK4M7Uy02nJtencuKzbu9jiRHUeEXERGfqh5fgVcGdaTAwfUv/cD27ANeR5LDqPCLiISwgV8MZOAXA32+39NqJfBcv7ZsyT7IgBdmcOBgvs8/Q0qnRIXfzGLN7B4ze+HQ8yZmdol/o4mISHnWpWlNHu6dyqLN+xn66g8UFAT/PWXlQUl7/GOBHKDToefrgAf9kkhERILGVR0aclv3+ny7bCcPfPij13GEkhf+Rs65fwMHAZxz+wHzWyoREQkaIy5syWUtqzFu9kbGTlnqdZyQV9LCn2tmFQAHYGaNKDwDICIickJmxn/6diC9bhwPfr6cbxat8zpSSCtp4b8P+AKoZ2ZvAN8Cf/NbKhERCSqR4WG8dGNn6leO4ra3FvDzuiyvI4WsEhV+59zXwBXADcB4IM05N9l/sUREJNjEV4jilRs7EhMRxsCxc9i8a7/XkULSCQu/mbX77QEkAxuBDUD9Q20iIiIlVj+xEmP6t2fngXxueGkG+3PzvI4Ucorr8T926DEamAWMAV449N9P+TfaqdMiPSIiZU9awyT+3TuVJVsOaJifB05Y+J1zZzvnzgbWAO2cc2nOufZAW2B5IAKeCufcROfckISEBK+jiIjIYXqnp3Bb9/p8t3wXD36kYX6BVNKb+5o55xb+9sQ5twho459IIiISCv58YUsuaV6Vl2dt5PVpy7yOEzJKWvgzzexFMzvLzLofmsFPCy6LiEipmRlP9DuDNrViGfnZL0xfutHrSCGhpIV/IPAz8Cfgz8DiQ20iIiKlFhkexss3dqJGxUj++OZ8Vm3Ran7+VtLhfAecc0845y4/9HjCOafllkRE5JRVrRjDyzekk1fguOHlH9i1T/PD+VNJF+npYmZfm9kvZrbyt4e/w4mISGhoWrsKT/Vpza+7DjJk3A/k605/vynpqf6XgMeBrkD6YQ8RERGfOLdlXe48vyGz1mbz97czAFjT/3rW9L/e42TBJaKE2+1yzn3u1yQiIhLyhpyTyrLNe3hr/haa1lzCuV4HCkIlLfyTzOxR4H0OW5zHOTfPL6lERCRkPXxNGiu3fc8/v1pB1YhqtFizwOtIQaWkhf+MQ3+mHdbmgHN8G0dEREJdRHgYLw7qzCVPfs+91brw1swvvI4UVEpU+A/N3iciIhIQVeKiefmGdK54YhIAu/fnEl8hyuNUweGEhd/Mbj/R6865x30bR0REBLaOepqC0aN599Dz9W1PZz2QOHQoScOHeRmt3Cvurv5Khx5pwK1AnUOPW4Dm/o0mIiKhKmn4MFKXZBKbXjiArGfv//D2P19X0feBE/b4nXP3A5jZVxQu0rPn0PORwDt+TyciIgJcnFqFsbM3kVp7Odd0bOx1nHKtpOP46wO5hz3PBRr4PI2IiMhRImrX5rHrOtCiegx3T/yFuSu3eB2pXCtp4X8NmG1mI83sPmAW8Ir/YvmGmfUyszG7du3yOoqIiJRSVJ06xERG8NKgjiREh3HL6/PYtHOv17HKrZLO1f8QhYvy7AB2AgOdcw/7M5gvOOcmOueGJCQkeB1FRERKIfm1V0l+7VUAalaO47k/tGNXTgGDx/5AzsF8j9OVTyXt8eOcm+ece/LQ40d/hhIRETmW9g2r8/8uacKizQf424TZXscpl0pc+EVERMqCPp2a8If21fno5yxe+G6x13HKHRV+EREpd+6/Mo0OdeN45JtVTF2ywes45YoKv4iIlDvhYcbzN3SkelwEwyf8xNpte7yOVG6o8IuISLlUpWIML1yfRk6eY/C4WRzIzfM6Urmgwi8iIuVWi3rV+OelTfllWw4j3pyFc87rSGWeCr+IiJRrl3doxKAONfl8yU6e+fpnr+OUeSr8IiJS7t3dux0d68Xx+KQ1fL94nddxyjQVfhERKffCwoznbuhIzYoR3Pb2Qt3sdwIq/CIiEhQqx8XwXP/25OQ5bho3iwMHdbPfsajwi4hI0GhVP5EHLzmNpdtyuGP8HK/jlEkq/CIiElSu6tiYP7SrzsTFWbw4STP7HU2FX0REgs7IK9vTtnYFHvl6FT8s2+h1nDJFhV9ERIJORHgYzw/oSJUK4Qwb/xObd+3zOlKZocIvIiJBqXpCLKP7tmF3Tj43j/uBg3laxhdU+EVEJIilN67JneelMH/jfu5/f57XccoEFX4REQlqg85uTq/UKrw+bwvvz1rudRzPqfCLiEjQ+3ffdE6rFsXdn/xC5rosr+N4SoVfRESCXoWoSJ6/vgPhZtzyegZ79ud4HckzKvwiIhISUmok8O/LU1m78yB/fmN2yK7kp8IvIhLi1mev9zpCwPRs24DBHWvy7fLdPPP1Iq/jeEKFX0QkxG3cG1oT3NzZqy0d6sbyxOS1zFy6wes4AWfBfKrDzHoBvRo3bnzTsmXLvI4jIhJQA78YWOw2S7KWkH0wm7QaacVuO/bCsb6IVSZs3bWPi5+agnPw6Z/OpHpCnNeRfM7M5jrnfndgg7rH75yb6JwbkpCQ4HUUEZEyZX32ejI2Z5B9MBuAjM0ZZGzOCJnT/kkJsTzV53R2HMjnj6/OJj+/wOtIARPUPf7fpKWluYyMDK9jiIiUOQO/GEjG5gwWDljodRRPPPvVQh75bi03nlGTey5v73UcnwrJHr+IiMiJ3HJ+S85vHM/LszbxxfzVXscJCBV+EZEQVyuultcRPGNmPN6vA/USIvm/DzJZu3WX15H8ToVfRCTE1alYx+sInqpUIZrR17UlJ99x62sZ5Ab5Yj4q/CIiEvJaJSdxd4+G/LzlAPe+G9z3hKnwi4iIAH/o1oxLm1dmwvxtfBDEi/mo8IuIiBzyr2vSaVy1cDGf5Rt3eB3HL1T4RUREDomNieKZP7THAX98Yy4HcvO8juRzKvwiIiKHOa12Ve7r2ZhftuXwj3fmeB3H51T4RUREjtKn82lc0bIq7y3M4p2Zv3gdx6dU+EVERI7hwavTOK1aFPd9tpxfNmR5HcdnVPhFRESOITY6ktH92mPA0Dfmsj/3oNeRfEKFX0RE5Dia1K7KfT0bsWx7Lve8Exzj+1X4RURETuCazk25vEUV3l2Yxbs/lP/r/Sr8IiIixXjomnQaVY3ivk+Xs3xj+b7er8IvIiJyHGv6X8+a/tcTGx3J09e1pcDBsDfnkXOw/I7vV+EXEREpgdS6ify9RwpLtuYw8r3ye71fhV9ERKSE+ndL5eJmCYyfv51P5q70Ok6pqPCLiIichEf6pJOcEMnfP17K2q27vI5z0lT4RURETkLFCtE8de3pHMgrYPgbGeTlF3gd6aSo8IuIiJyk01Nq8Jez6vPTpgM88vE8r+OcFBV+ERGRUrjp3BaclVKRl2ZvZtKitV7HKTEVfhERkVIICwvj8b7pVI8L5473F7N1116vI5WICr+IiMgJ5K5ff9zXqsbH8p8rW7Bjfz4j3pxDQUHZv96vwi8iInICeRs2nPD1rqn1GNKxJtPW7OW5bxYFKFXpqfCLiIicor9c3IZ2tSrw38m/8uPKTV7HOSEVfhERkaNsHfU0mc1S2TdnDgCZzVLJbJbK1lFPH3P7iIhw/ntdOypEhvHntxeQvT8nkHFPigq/iIjIUZKGDyN1SSax6ekApC7JJHVJJknDhx33PfWTKvPAxY1Zs/Mgfy/DS/iq8IuIiPjIZR2acHWrqny8eCfvziybS/iq8IuIiJxARO3aJ7X9yCvb07BKJPd/voI1W3b6KVXpqfCLiIicQFSdOie1fVxMFP/tczq5+Y4/vTmXvLx8PyUrHRV+ERERH2vdoAZ/7laH+ZsO8Nin80+47Zr+17Om//UBSqbCLyIi4hc3n9+KLslxjPlhE9Mzf/U6ThEVfhERET8ICwvjsT7tqRwTxh3vL2Zn9n6vIwHlpPCbWUMze8nM3j2srbeZvWBmH5lZDy/ziYiIHEvNqpX456VN2bgnj7veLhtD/Pxe+M3sZTPbYmaLjmq/0MyWmtlyM7vzRPtwzq10zt14VNuHzrmbgBuAPj4PLiIi4gMXtG3ItW2q8fkvu5kwbYnXcQLS4x8HXHh4g5mFA6OBnkBzoK+ZNTezVmb2yVGP6sXs/+5D+xIRESmT7undjsZVI3nwq5Ws2rTD0yx+L/zOuSlA1lHNHYDlh3ryucAE4DLn3ELn3CVHPbYca79W6BHgc+fcvGO8PsTMMswsY+vWrb7+sUREREosNiaKx69uzcF8+POEeZ4O8fPqGn8d4PBbHNcdajsmM6tmZs8Bbc3srkPNw4HzgKvM7Jaj3+OcG+OcS3POpSUlJfkwuoiIhIrk114l+bVXfbKv1ik1+XO3Ovy06QCPFzPEz58iPPpcO0abO97GzrntwC1HtT0FPOXjXCIiIn5z8/mtmLoiizE/bKJbs3V0bFo34Bm86vGvA+od9rwucOIFj0VERMq5sLAw/tOnPZWiw7jjvZ/Zs+9A4DME/BMLzQGamFmKmUUB1wIfe5RFREQkYGpXi2fkRY35dXce9743N+CfH4jhfOOBmUBTM1tnZjc65/KAYcCXQCbwtnPuZ39nERERKQsu69CE3s0r88HPO5kenhjQzw7EXf19nXO1nHORzrm6zrmXDrV/5pw7zTnXyDn3kD8+28x6mdmYXbt2+WP3IiIipfb/rmpPvfgIHk3sxI4A3nJXLmbuKy3n3ETn3JCEhASvo4iIiByhUmwMj17Rgj3h0fynchrOHfced58K6sIvIiJSlnVsVpcBWT+SvHYJ05ZvC8hnejWcT0RERIDrDqwgd+kcmjV+MiCfpx6/iIiIh37rgZsda4ob31PhFxER8cDWUU+T2SyVfXPmAJDZLJXMZqlsHfW0Xz9Xp/pFREQ8kDR8GEnDh7Gm//XsmzOH1CWZAfncoO7xazifiIjIkYK68Gs4n4iIlAcRtWsH7LOCuvCLiIiUB1F1jrtArc+p8IuIiIQQFX4REZEQosIvIiISQlT4RUREQogKv4iISAgJ6sKvcfwiIiJHCurCr3H8IiIiRwrqwi8iIiJHUuEXEREJISr8IiIiIUSFX0REJISo8IuIiIQQFX4REZEQosIvIiISQoK68GsCHxERkSNFeB3An5xzE4GJaWlpN3mdRUSkLBp74VivI0iABXWPX0RERI6kwi8iIhJCVPhFRERCSFBf4xcRESnrkl97NaCfpx6/iIhICFHhFxERCSEq/CIiIiFEhV9ERCSEBHXh18x9IiIiRwrqwu+cm+icG5KQkOB1FBERkTIhqAu/iIiIHEmFX0REJISo8IuIiIQQFX4REZEQosIvIiISQlT4RUREQogKv4iISAhR4RcREQkh5pzzOoPfmdlWYM1RzQnAsab0O1b7sdoSgW0+CXhyjpfb3/soyXuK2+ZEr5f07/1Y7eX5WJRmPyXdvrTHQ98N379H343A7EffjSMlO+eSftfqnAvJBzCmpO3HacsoS7n9vY+SvKe4bU70ekn/3o/VXp6PRWn2U9LtS3s89N3w/Xv03QjMfvTdKNkjlE/1TzyJ9uNt6wVfZCnNPkrynuK2OdHrJ/P3XlaOh69ynOx+Srp9aY+Hvhu+f4++G4HZj74bJRASp/r9wcwynHNpXucQHYuyRsej7NCxKFvKyvEI5R7/qRrjdQApomNRtuh4lB06FmVLmTge6vGLiIiEEPX4RUREQogKv4iISAhR4RcREQkhKvw+Zma9zewFM/vIzHp4nSfUmVlDM3vJzN71OksoMrM4M3vl0Hein9d5Qp2+D2WHl7VChf8wZvaymW0xs0VHtV9oZkvNbLmZ3XmifTjnPnTO3QTcAPTxY9yg56PjsdI5d6N/k4aWkzwuVwDvHvpOXBrwsCHgZI6Hvg/+dZLHwrNaocJ/pHHAhYc3mFk4MBroCTQH+ppZczNrZWafHPWofthb7z70Pim9cfjueIjvjKOExwWoC/x6aLP8AGYMJeMo+fEQ/xrHyR+LgNeKiEB+WFnnnJtiZg2Oau4ALHfOrQQwswnAZc65h4FLjt6HmRnwL+Bz59w8/yYObr44HuJ7J3NcgHUUFv/5qKPhFyd5PBYHNl1oOZljYWaZeFQr9EUsXh3+12OBwn/I6pxg++HAecBVZnaLP4OFqJM6HmZWzcyeA9qa2V3+DhfCjndc3geuNLNnKcNTmAahYx4PfR88cbzvhme1Qj3+4tkx2o4765Fz7ingKf/FCXknezy2A/oFzP+OeVycc3uBgYEOI8c9Hvo+BN7xjoVntUI9/uKtA+od9rwusMGjLKLjUVbpuJQtOh7/v717C7GqiuM4/v1JDcoYBZkRQVik3aHSirLoIgW9FGZRdNOYHo0eiggKIkao7dCFPgAAA85JREFUqKfoQg81QcRJrcS0yCJFGBSTaWrUBp8SlcEYiwJjDJr+Pex1cHM4F2mm3J71+8Aw+6y119r/mQXnv/bZ++xVHZUbCyf+znYC8yWdL6kHeAD47ATHlDOPRzV5XKrF41EdlRsLJ/4SSTVgO3CRpIOS+iLiL2AlsAkYBdZExJ4TGWcuPB7V5HGpFo9HdZwsY+FFeszMzDLiM34zM7OMOPGbmZllxInfzMwsI078ZmZmGXHiNzMzy4gTv5mZWUac+M26nKRJSd+XftouZfx/kvRxWiN+R4ptv6TxUqzzWrRbJam/oWyRpJG0/Y2k0//7v8Ds5OPv8Zt1OUlHImL2NPd5SnowyVT6uAxYFRFLS2UrgEURsfI42q6LiAWlsteAXyLiJUl9wJyIeGUqMZp1I5/xm2VK0j5JL0r6TtIuSRen8l5J70naKWlY0t2pfIWktZI2AF9JmiHpLUl7JG2U9IWkeyUtkbSudJzbJX3aJISHgPXHEeedkranOFdL6k1PPjsqaWHaR8B9wEep2Xrgwan8f8y6lRO/Wfeb1fBR//2lusMRcTXwNvB0KnsO2BwR1wC3Aq9K6k111wPLI+I24B5gHnAF8HiqA9gMXCLprPT6MWCgSVyLgaF2gUuaCzwLLElxjgBPpuoaxXPP632NRcRPABFxGDhN0hnt+jfLkZflNet+ExFxZYu6+pn4EEUiB7gDuEtSfSIwEzgvbX8dEb+m7RuBtRHxN3BI0hYo1huV9AHwsKQBignBo02OfQ4w3iH2G4BLgW3FST09wGCqqwFbJT1DMQGoNbQdT8f4rcMxzLLixG+Wtz/T70mOvR8IWBYRe8s7SroO+KNc1KbfAWADcJRictDsfoAJiklFOwK+jIhHGisiYp+kMeAmYCmwsGGXmekYZlbij/rNrNEm4Il03RxJV7XYbxBYlq71nw3cUq+IiDGKNcefB95v0X4UuLBDLNuAmyVdkGLplTS/VF8DXgdGI+JQvVDSDGAOcKBD/2bZceI3636N1/hf7rB/P3AqMCJpd3rdzCfAQWA38A6wA/i9VP8hcCAifmzR/nNKk4VmIuJnoA9YLekHionAgtIua4DLOXZTX921wGBETLbr3yxH/jqfmf1rkmZHxBFJZwLfAovrZ96S3gCGI+LdFm1nAVtSm2lN0JLepFj3fOt09mvWDXyN38ymYmO6c74H6C8l/SGK+wGeatUwIiYkvQCcC+yf5riGnfTNmvMZv5mZWUZ8jd/MzCwjTvxmZmYZceI3MzPLiBO/mZlZRpz4zczMMuLEb2ZmlpF/AJCUau3UuUUGAAAAAElFTkSuQmCC\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();"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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
}