{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"<div class=\"alert alert-info\">\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.18?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/docs/tutorials).\n",
"- **Source files:**\n",
"[analysis_mwl.ipynb](../_static/notebooks/analysis_mwl.ipynb) |\n",
"[analysis_mwl.py](../_static/notebooks/analysis_mwl.py)\n",
"</div>\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.modeling.models import Models\n",
"from gammapy.datasets import Datasets, FluxPointsDataset, SpectrumDatasetOnOff\n",
"from gammapy.estimators import FluxPoints, FluxPointsEstimator\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:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"datasets:\r\n",
"- name: Fermi-LAT\r\n",
" type: MapDataset\r\n",
" filename: Fermi-LAT-3FHL_data_Fermi-LAT.fits\r\n"
]
}
],
"source": [
"!cat $GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_datasets.yaml"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"components:\r\n",
"- name: Crab Nebula\r\n",
" type: SkyModel\r\n",
" spectral:\r\n",
" type: LogParabolaSpectralModel\r\n",
" parameters:\r\n",
" - name: amplitude\r\n",
" value: 0.018182745349064267\r\n",
" unit: cm-2 s-1 TeV-1\r\n",
" min: .nan\r\n",
" max: .nan\r\n",
" frozen: false\r\n",
" error: 0.003026327991562108\r\n",
" - name: reference\r\n",
" value: 5.054833602905273e-05\r\n",
" unit: TeV\r\n",
" min: .nan\r\n",
" max: .nan\r\n",
" frozen: true\r\n",
" error: 0.0\r\n",
" - name: alpha\r\n",
" value: 1.652368617859867\r\n",
" unit: ''\r\n",
" min: .nan\r\n",
" max: .nan\r\n",
" frozen: false\r\n",
" error: 0.05762513693893088\r\n",
" - name: beta\r\n",
" value: 0.03921700077803329\r\n",
" unit: ''\r\n",
" min: .nan\r\n",
" max: .nan\r\n",
" frozen: false\r\n",
" error: 0.00521472221220211\r\n",
" spatial:\r\n",
" type: PointSpatialModel\r\n",
" frame: icrs\r\n",
" parameters:\r\n",
" - name: lon_0\r\n",
" value: 83.63310241699219\r\n",
" unit: deg\r\n",
" min: .nan\r\n",
" max: .nan\r\n",
" frozen: true\r\n",
" error: 0.0\r\n",
" - name: lat_0\r\n",
" value: 22.019899368286133\r\n",
" unit: deg\r\n",
" min: -90.0\r\n",
" max: 90.0\r\n",
" frozen: true\r\n",
" error: 0.0\r\n",
"- type: FoVBackgroundModel\r\n",
" datasets_names:\r\n",
" - Fermi-LAT\r\n",
" spectral:\r\n",
" type: PowerLawNormSpectralModel\r\n",
" parameters:\r\n",
" - name: norm\r\n",
" value: 1.3004625872247901\r\n",
" unit: ''\r\n",
" min: 0.0\r\n",
" max: .nan\r\n",
" frozen: false\r\n",
" error: 0.07512322002655547\r\n",
" - name: tilt\r\n",
" value: 0.0\r\n",
" unit: ''\r\n",
" min: .nan\r\n",
" max: .nan\r\n",
" frozen: true\r\n",
" error: 0.0\r\n",
" - name: reference\r\n",
" value: 1.0\r\n",
" unit: TeV\r\n",
" min: .nan\r\n",
" max: .nan\r\n",
" frozen: true\r\n",
" error: 0.0\r\n"
]
}
],
"source": [
"!cat $GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_models.yaml"
]
},
{
"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": 4,
"metadata": {},
"outputs": [],
"source": [
"path = Path(\"$GAMMAPY_DATA/fermi-3fhl-crab\")\n",
"filename = path / \"Fermi-LAT-3FHL_datasets.yaml\"\n",
"\n",
"datasets = Datasets.read(filename=filename)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"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 : \n",
" Parameters:\n",
" amplitude : 1.82e-02 1 / (cm2 s TeV)\n",
" reference (frozen) : 0.000 TeV \n",
" alpha : 1.652 \n",
" beta : 0.039 \n",
" lon_0 (frozen) : 83.633 deg \n",
" lat_0 (frozen) : 22.020 deg \n",
"\n",
"Component 1: FoVBackgroundModel\n",
"\n",
" Name : Fermi-LAT-bkg\n",
" Datasets names : ['Fermi-LAT']\n",
" Spectral model type : PowerLawNormSpectralModel\n",
" Parameters:\n",
" norm : 1.300 \n",
" tilt (frozen) : 0.000 \n",
" reference (frozen) : 1.000 TeV \n",
"\n",
"\n"
]
}
],
"source": [
"models = Models.read(path / \"Fermi-LAT-3FHL_models.yaml\")\n",
"print(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": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SkyModel\n",
"\n",
" Name : Crab Nebula\n",
" Datasets names : None\n",
" Spectral model type : LogParabolaSpectralModel\n",
" Spatial model type : PointSpatialModel\n",
" Temporal model type : \n",
" Parameters:\n",
" amplitude : 1.82e-02 1 / (cm2 s TeV)\n",
" reference (frozen) : 0.000 TeV \n",
" alpha : 1.652 \n",
" beta : 0.039 \n",
" lon_0 (frozen) : 83.633 deg \n",
" lat_0 (frozen) : 22.020 deg \n",
"\n",
"\n"
]
}
],
"source": [
"print(models[\"Crab Nebula\"])"
]
},
{
"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": 7,
"metadata": {},
"outputs": [],
"source": [
"datasets_hess = 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_hess.append(dataset)\n",
"\n",
"dataset_hess = datasets_hess.stack_reduce(name=\"HESS\")\n",
"\n",
"datasets.append(dataset_hess)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Datasets\n",
"--------\n",
"\n",
"Dataset 0: \n",
"\n",
" Type : MapDataset\n",
" Name : Fermi-LAT\n",
" Instrument : \n",
" Models : \n",
"\n",
"Dataset 1: \n",
"\n",
" Type : SpectrumDatasetOnOff\n",
" Name : HESS\n",
" Instrument : \n",
" Models : \n",
"\n",
"\n"
]
}
],
"source": [
"print(datasets)"
]
},
{
"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": 9,
"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(data=flux_points_hawc, name=\"HAWC\")\n",
"\n",
"datasets.append(dataset_hawc)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Datasets\n",
"--------\n",
"\n",
"Dataset 0: \n",
"\n",
" Type : MapDataset\n",
" Name : Fermi-LAT\n",
" Instrument : \n",
" Models : \n",
"\n",
"Dataset 1: \n",
"\n",
" Type : SpectrumDatasetOnOff\n",
" Name : HESS\n",
" Instrument : \n",
" Models : \n",
"\n",
"Dataset 2: \n",
"\n",
" Type : FluxPointsDataset\n",
" Name : HAWC\n",
" Instrument : \n",
" Models : \n",
"\n",
"\n"
]
}
],
"source": [
"print(datasets)"
]
},
{
"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": 11,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"HDU 'MASK_FIT' not found\n"
]
}
],
"source": [
"path = Path(\"crab-3datasets\")\n",
"path.mkdir(exist_ok=True)\n",
"\n",
"filename = path / \"crab_10GeV_100TeV_datasets.yaml\"\n",
"\n",
"datasets.write(filename, overwrite=True)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"datasets:\r\n",
"- name: Fermi-LAT\r\n",
" type: MapDataset\r\n",
" filename: Fermi-LAT.fits\r\n",
"- name: HESS\r\n",
" type: SpectrumDatasetOnOff\r\n",
" filename: pha_obsHESS.fits\r\n",
"- name: HAWC\r\n",
" type: FluxPointsDataset\r\n",
" filename: HAWC.fits\r\n"
]
}
],
"source": [
"!cat crab-3datasets/crab_10GeV_100TeV_datasets.yaml"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"datasets = Datasets.read(filename)\n",
"datasets.models = models"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Datasets\n",
"--------\n",
"\n",
"Dataset 0: \n",
"\n",
" Type : MapDataset\n",
" Name : Fermi-LAT\n",
" Instrument : \n",
" Models : ['Crab Nebula', 'Fermi-LAT-bkg']\n",
"\n",
"Dataset 1: \n",
"\n",
" Type : SpectrumDatasetOnOff\n",
" Name : HESS\n",
" Instrument : \n",
" Models : ['Crab Nebula']\n",
"\n",
"Dataset 2: \n",
"\n",
" Type : FluxPointsDataset\n",
" Name : HAWC\n",
" Instrument : \n",
" Models : ['Crab Nebula']\n",
"\n",
"\n"
]
}
],
"source": [
"print(datasets)"
]
},
{
"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": 15,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"HDU 'MASK_FIT' not found\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"OptimizeResult\n",
"\n",
"\tbackend : minuit\n",
"\tmethod : minuit\n",
"\tsuccess : True\n",
"\tmessage : Optimization terminated successfully.\n",
"\tnfev : 305\n",
"\ttotal stat : -12695.04\n",
"\n",
"CPU times: user 3.42 s, sys: 22.2 ms, total: 3.44 s\n",
"Wall time: 3.45 s\n"
]
}
],
"source": [
"%%time\n",
"fit_joint = Fit(datasets)\n",
"results_joint = fit_joint.run()\n",
"print(results_joint)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's display only the parameters of the Crab spectral model"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"LogParabolaSpectralModel\n",
"\n",
" name value unit min max frozen error \n",
"--------- ---------- -------------- --- --- ------ ---------\n",
"amplitude 3.8966e-03 cm-2 s-1 TeV-1 nan nan False 3.119e-04\n",
"reference 5.0548e-05 TeV nan nan True 0.000e+00\n",
" alpha 1.2560e+00 nan nan False 1.747e-02\n",
" beta 6.1986e-02 nan nan False 9.778e-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": 17,
"metadata": {},
"outputs": [],
"source": [
"# compute Fermi-LAT and HESS flux points\n",
"energy_edges = MapAxis.from_energy_bounds(\"10 GeV\", \"2 TeV\", nbin=5).edges\n",
"\n",
"flux_points_fermi = FluxPointsEstimator(\n",
" energy_edges=energy_edges, source=\"Crab Nebula\"\n",
").run([datasets[\"Fermi-LAT\"]])\n",
"\n",
"\n",
"energy_edges = MapAxis.from_bounds(\n",
" 1, 15, nbin=6, interp=\"log\", unit=\"TeV\"\n",
").edges\n",
"flux_points_hess = FluxPointsEstimator(\n",
" energy_edges=energy_edges, source=\"Crab Nebula\"\n",
").run([datasets[\"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": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(array([0.00698646, 0.02015874, 0.05816603, 0.16783227, 0.48426319]), array([0.01186753, 0.03424258, 0.09880353, 0.2850877 , 0.82259195]))\n",
"(array([0.25892541, 0.41036912, 0.65039123, 1.03080063, 1.6337089 ,\n",
" 2.58925412]), array([0.32596778, 0.51662412, 0.81879405, 1.29770111, 2.05671765,\n",
" 3.25967781]))\n",
"None\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAF3CAYAAABE0Ck1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3yV5f3/8dcnGxIII2HvHSIIGpCloCiiAtaJaBFBRSvgoLbf+lMrWqxdtg6wFqvixIFWRBxVZA8loLICyJRNCDMQyLp+fyRQQEhOwjm5k3Pez8fjPOC+zj3eySF8ct/3dV+XOecQERGR0BDmdQAREREpOyr8IiIiIUSFX0REJISo8IuIiIQQFX4REZEQosIvIiISQiK8DlAWEhISXJMmTbyOISIiUmYWL1682zmXeGp7SBT+Jk2akJqa6nUMERGRMmNmm07Xrkv9IiIiIUSFX0REJISo8IuIiISQkLjHLyIi5UNOTg5btmzhyJEjXkcJGjExMTRo0IDIyEif1lfhFxGRMrNlyxaqVKlCkyZNMDOv41R4zjkyMjLYsmULTZs29WkbXeoXEZEyc+TIEWrWrKmi7ydmRs2aNUt0BUWFX0REypSKvn+V9Pupwi8iIiHFzBg8ePDx5dzcXBITE+nXr1+J9tOkSRN279591uuUtXJf+M2smZm9bGaTT2iLNbPXzOwlM7vFy3wiIlKxxMbGsnz5crKysgD48ssvqV+/vsepyk5AC7+ZvWJmu8xs+Sntfc1stZmtNbPfFbUP59x659ztpzRfC0x2zt0JDPBzbBERCXJXXHEF06ZNA2DSpEkMGjTo+Ht79uzhF7/4Be3bt6dLly4sXboUgIyMDPr06UPHjh256667cM4d3+bNN9+kc+fOdOjQgbvuuou8vLyy/YJKINC9+icC44DXjzWYWTgwHrgM2AIsMrOPgXDgqVO2H+ac23Wa/TYAlhX+vfx+d0VE5Iwen7qCldsO+HWfbetV5bH+ycWud9NNN/HEE0/Qr18/li5dyrBhw5gzZw4Ajz32GB07duSjjz7i66+/5tZbb+X777/n8ccfp0ePHvz+979n2rRpTJgwAYC0tDTeffdd5s2bR2RkJPfccw9vvfUWt956q1+/Nn8JaOF3zs02syanNHcG1jrn1gOY2TvA1c65pwBfb7BsoaD4f08FuF1RYjOegosf8jqFiEjQat++PRs3bmTSpElceeWVJ703d+5cPvjgAwAuueQSMjIy2L9/P7Nnz+bDDz8E4KqrrqJ69eoATJ8+ncWLF9OpUycAsrKyqFWrVhl+NSXjxXP89YHNJyxvAS4408pmVhN4EuhoZg8V/oLwITDOzK4Cpp5hu+HAcIBGjRr5KXoZmfUnFX4RCXq+nJkH0oABA3jwwQeZOXMmGRkZx9tPvIR/zLGe86frQe+cY8iQITz11KkXrcsnLwr/6Z47+Pl3+dgbzmUAd5/SdggYWtRBnHMTgAkAKSkpZ9x/qb16ld93CcCOpYHb/9Bp/t+niEgFNWzYMOLj42nXrh0zZ8483n7RRRfx1ltv8eijjzJz5kwSEhKoWrXq8fZHHnmEzz77jL179wLQu3dvrr76ah544AFq1arFnj17OHjwII0bN/boKyuaF4V/C9DwhOUGwDYPcpQv+zbB/hMuhGyaW/BnfEOoVj7/8YiIVGQNGjTgvvvu+1n7mDFjGDp0KO3bt6dy5cq89tprQMG9/0GDBnHeeefRs2fP41eT27Zty9ixY+nTpw/5+flERkYyfvz4clv47XSXNPx6gIJ7/J84584pXI4A1gC9ga3AIuBm59yKQGVISUlxqampgdr9GeXn55Obm0teXt7xP/Pz88nOySXzSA6HjuZw6GgumUdzOZKTx5GcfI7m5jF4dk/+ecHX5OU7cgtf+fknfE5mGBAeZoSHGRFhRkR4GNERYURFhBMTGU50ZDixUeHERkcQGxNJXEwU8ZWjiIyIIDw8nIjCP0VEylJaWhpJSUlexwg6p/u+mtli51zKqesG9IzfzCYBvYAEM9sCPOace9nMRgJfUNCT/5VAFv1AyMvLIysri5ycHHJycth36Ag79mex68ARdmdmszszm31HcjlwJI8DR/M5mJ1PZrbjUE4+h7Mdh3OL/mVrcAz8edYOv+c2oFKEUTnKiI0Mo0pUGFVjwomvFEG1ShHUjI0isUo0iVViqB1fmfrVK1O9SmUiIyMJCwu+PpQiIqEo0L36B52h/VPg00AeO1DSDx7liSk/sGnXfjKy8tiTlU/WaQp5GBAXHUbV6DCqRBm1Y8OJjYokNtKIjQonLrrgjLxydASVoyKoHBVOdEQ4MVHhbPpxOO+2SyIiPIyIcCMyPJwwg7Cw/3UucQ5y8x15+Y68/HyycwtfefkczcnjSE4eh7MLXlk5eRw6mkdm4dWFg9n5ZB7N4+DRPDbty+bgrqNkHs0n/zRfb6UIo0alMBIqR1C7ShT14qOpX60SDWvG0jSxCo0SqhITE60hOEVEKoignp3PzPoD/Vu0aOG3fUaEGd9u2k98pKNh1Qg61AmnZuUIEuMiSYwrOFuuVTWGmlViiIqMJDIykoiIiOOX1iMifPiWn/9XyuLOkHPu+G2I7Jwcdh88ws79h0k/cIRdB46w8+BRdh3MJj0zm/RDuSzacoi9aw+etI+IMKgdG07dKlE0qhFDs4TKtKhVhVZ1q9EwoarP00SKiEjZCOrC75ybCkxNSUm501/7rB4bxewHLyIrK4uoqCiioqJ8K+blkJkd/6UkJiaGqlWq0KzemdfPycnh4OEsftqdyU8ZmWzKOMTmPVls2X+UrfuzWbozi+y8vcfXrxxpNIqPpGmNGFrUiqVtvWq0b1SDejXjdYVARMQjFbNieSw6Opro6GivY5S5yMhIasRHUiO+Kh2a//z9rKwj/LT7AD/u2M/aXZmsTT/Exj1HmLfxIJ+tOQBsB6BGpTCa1Yimda1Y2tWP57ymCTSvU12dDUVEyoAKv/hNpUoxtG4YQ+uGJ49YlZeXx5bd+1m2eQ8rt+1n9c5DrN2dxZJtu3nru93AOqpEGa0SYjinXhU6NqpO5+a1qFezqjdfiIiUL8fGNdFYJH6hwi8BFx4eTuPaNWhcu8ZJYzIfyDzMD5t288PmvSzfdoBVu7J4I3UXr6XuAlZTJy6cc+rEcn6jalzQIpFzG9ciPFxPF4jI2QkPD6ddu3bHlz/66COaNGni12N069aN+fPn/6x9zJgxxMXF8eCDD/7svfT0dOrVq8e4ceO46667GDFiBPPmzSM7O5sNGzbQunVrAB555BGuv/76UmdT4RfPVI2rzIXJjbgw+X9DKu89eJjFG3axZOMefthygEWbM/lq7QH4+idiI412dSqT0jiebi1qkdK8NlGR+icsEhL2bfLbripVqsT3339f4u1yc3N97tN1uqJfnPfff58uXbowadIk7rrrLsaPHw/Axo0b6devX6kyn47+15RypXqVylzavgmXtm8CFAyCtHrrHhau3cmijXv5ftshFm7exri526gUYbSvW5kuTavTs00dzm2cqCsCIsHqxJFNA2Dx4sWMHj2azMxMEhISmDhxInXr1qVXr15069aNefPmMWDAAKZOnUrHjh1ZvHgx6enpvP766zz11FMsW7aMgQMHMnbsWADi4uLIzMwsUYZJkybx9NNPc/PNN7N161bq168fiC81uAt/IB7nk7IVFhZGUsMEkhomHJ+c4af0/cxdvYMF6zJYvOUg38zewrOztxAfHUanhnH0aJnApec0oEHNKkXvXPcNRcqnU+cqOdMcJqX82c3KyqJDhw4ANG3alPfee49Ro0YxZcoUEhMTeffdd3n44Yd55ZVXANi3bx+zZs0CYOrUqURFRTF79myeffZZrr76ahYvXkyNGjVo3rw5DzzwADVr1ixxps2bN7Njxw46d+7MjTfeyLvvvsvo0aNL9fUVJ6gLfyAe5xPvNUqM5+bEeG7uUbC8fsc+ZqZtY966jOO3BsZ8tp6m1SPp3rQ6l51Tl26t6hIZcZqnBvx4+VBE/CxAc5iceql/+fLlLF++nMsuuwwo6JBct27d4+8PHDjwpO0HDBgAQLt27UhOTj6+brNmzdi8eXOpCv8777zDjTfeCMBNN93E7bffrsIvcibN6lSjWZ1qDLsY8vLy+W7jLr5euZ256/Yw6btdvLlkF7GRS+ncqAqXJtXiyg6NqR4XU7BxgC8fikgpnHom/+pVBUV/zP6AHM45R3JyMgsWLDjt+7GxsSctH3ucOyws7KRHu8PCwsjNzT1p3Ycffphp0wq+nqLu0U+aNImdO3fy1ltvAbBt2zZ+/PFHWrZsWfIvqBgq/FL+lWCK4nAgpfD123jYHxvJrMx6zMhszKwNTZmx7gCPfbKa86O30Iea3G4l23+J6BaCSIXQunVr0tPTWbBgAV27diUnJ4c1a9aQnJx81vt+8sknefLJJ4tcZ/Xq1Rw6dIitW7ceb3vsscd45513ePTRR886w6lU+CWoxUfkMKDaJgZU20S+m83iQ4kc2beNC/O++d9KhZcP98S2oEZCHY+SikiR4hsWv04pRUVFMXnyZO699172799Pbm4u999/v18K/+mMHTuWZ5555vjyHXfcwTXXXHPSOtdddx033XRTQAp/wKflLQ+8mpZXyre1Lw6ixY5PuTzuA1bvPgpA0+pRXJ6UwLWdmtCqbnWPE4oEn1JNy6uOuMUqN9PySoh49aqCTjgPLPc6SYm0iD4AwBcPXspPuw/y0eKNfL5iFy/O38aL87fRokY0V55Ti+s7N6VRQjFPCIhI4Kjg+1VQF349zleGKmonucLLh40SqnDv5e2493L4afdB/pO6kU+X7+S52Zt5bvZm2tepRL/2dbm+c1NqHOsYKCJSAelSfygJVCe2HUvh6AFo3CMw+/fwt/21O/fz/jfr+XRFOpv35xARBt2bVOX6lIZc3r4RUREaMEikJEp1qV+KpUv9UjYC9IxtedKidjwPDejI7/o7Fm9I571vN/JFWgaz1q+g6pQ0rkhKYHD35pzTsIbXUUVEfKLCH0oCdeYc4GdsywMzI6VZLVKa1eLJvHz+u3Qz7y36iQ+W7uLd73fRKiGG686rx01dmxNfKcrruCIiZ6TCL2dv6DQYE+91ijITGR7GVR0bc1XHxuw+eIR3F67jgyXbeOq/63l6+gYubVWdW7s344LmtTAzr+OKVHhDPy8YsPvVvq96nCQ46Aal+EfP33mdwBMJVWIYcVky0397KR/c1Zkrk2ow48e93PTvVC7+y1e8OD2NA1nZXscUkRPExcWdtDxx4kRGjhwJFEybW79+fTp06HD8tW/fPg4fPswtt9xCu3btOOecc+jRo8fxSXiefPJJkpOTad++PR06dOCbb7752THLE53xi39c/JDXCTxlZpzfNJHzmyaSeSSH979Zz6RFW/jTl+t5ZsYG+rapye09W9JOfQFESmVr5tbiV/KTBx54gAcffPCktqeeeoratWuzbNkyoGC0vcjISBYsWMAnn3zCkiVLiI6OZvfu3WRnl+9f9lX4RfwsLiaSoT1bM7RnaxZvSOfVOWv5NC2Dj5bvpl2dSgzt3pT+5zUmUlMIi/hs+6Ht3h5/+3YaN/5fp+XWrVsfb09ISDg+Zn9CQoIn+UpCj/OJlIG9h47yxty1TErdwvaDuSRUDmfg+fUY1rM1NeOii9+BSJDw5XG+Y/f0j1m1ZxWZOZmk1D75ybTS3vMPDw+nXbt2x5f37NnDgAEDGDduHGPGjOGll14iMTERgOrVqzNjxgy+//57+vTpQ/PmzenduzdDhgyhZcuWZGZm0qNHDw4fPsyll17KwIED6dmzZ6lynY2SPM4X1KccZtbfzCbs3x+8vc2lYqgeG829lycz96E+jB/YjiY1Yhg/ZzNdn5rO/W99w+rt+jcqcqqtmVtJ3ZlKZk7BvfTUnamk7kw968v+x6blPfZ64oknTnr/gQceOP7ejBkzAOjQoQPr16/nN7/5DXv27KFTp06kpaURFxfH4sWLmTBhAomJiQwcOJCJEyeeVb5AC+pL/c65qcDUlJSUO73OIgIQHmZc1bERV3VsRNq2fbz49WqmrdjNR8vm0qVRFe6+uCU929TR0wAS0k49kx/6+VBSd6aybMgyjxIViIuL49prr+Xaa68lLCyMTz/9lKSkJMLDw+nVqxe9evWiXbt2vPbaa9x2222eZi1KUJ/xi5RnSfWq8ewvL2DBQ735VY+GrN51mNteW0Kfp79m8rcbyM3L9zqiiBSaN28ee/fuBSA7O5uVK1fSuHFjVq9ezY8//nh8ve+///6kvgDlkQq/iMcSqsTwf/3as+DhyxhzZUuycvJ48MOV9HjqSybMWM2RnDyvI4p4rm5s3TI71j/+8Y+THufbuHEj69ato2fPnrRr146OHTuSkpLCddddR2ZmJkOGDKFt27a0b9+elStXMmbMmDLLWhrq3CdSzuTnO6b9sJkXZ65lxc4sqsWEc2uXBtzRqzVVYyK9jidyVkozVr8G8CmexuoXqcDCwoz+HRvRr0ND5q3ZyXPT1/DczE38e95P3JxSnxGXJlE9VsMCS+hQwfcvFX6RcsrM6NG6Dj1a1+H7TRk8899VvLxgC28u2srA8+oyqk9bEvQooIiUkO7xi1QAHRrXZOKd3fnsvh5cFp3GG99uoftT03l48hLSDx71Op6IVCAq/CIVSJu68Tzf4Gu+bPo2vVvVYNLi7fT403Qe/WAJuzP1C4CIFE+X+kUqoOZRB3hhSBfW7TrI3z5bwVup23nvux3cnFKP+/q0pVpl9QEQkdPTGb9IBda8VhX+OaQL/73/Qnq1qM7EhVvp9tRX/GXaMg4eyfE63slevargJVJCmwbfyqbBt3odI2gEdeHXkL0SKlrUrsq/buvKp/f2oFOjeF6Y8xPdn/qK8V+t0jgAIqcoalreY84991wGDRp0fPmaa67ho48+Or7cunVrxo4de3z5uuuu48MPPwTg22+/5aKLLqJ169a0adOGO+64g8OHDwfiSymVoC78zrmpzrnh8fHxXkcRKRNJ9eJ57c7u/OdXXWlTO5a/frWOi/70FW8tWE9efvCP2SHBK3tr2U3Lm5aWRn5+PrNnz+bQoUMAdOvWjfnz5wOQkZFBXFwcCxYsOL7NggUL6NatGzt37uSGG27gz3/+M6tXryYtLY2+ffty8ODBMstfnKAu/CJBa9+mIt/u2LgG791zEW8MS6FmbCQPT0mj91+n8/nSrYTCoF0SfHK3bSuzY7399tsMHjyYPn368PHHHwPQvXv344V//vz59OvXj/T0dJxzbNiwgUqVKlGnTh3Gjx/PkCFD6Nq1K1DwWO71119P7dq1yyx/cVT4RSqi/Zt9Wu3CVrX59IGLeW5gO/LyHXe//T3XPD+b737aE+CAIuVXVlbWSUPy/v73vz/p/XfffZeBAwcyaNAgJk2aBMD555/P8uXLyc7OZv78+XTt2pXWrVuTlpbG/Pnz6d69OwDLly/n/PPPL/OvqSTUq18kUALVkW3H0hLt34ABwJW1jUlRbXhmR2eueSGTq2JX83+1vqFRVObJGwyd5te4IqWV/vw4do8ff3w5rU3BkLQJI0aQOGrkmTYr1rFpeY+ZOHEix4Z1X7RoEYmJiTRu3JgGDRowbNgw9u7dS/Xq1UlOTmbJkiUsXLiQ3/72t6xfv5758+fz3Xff0a1bt1LnKWs64xepKPZtgk1z4eiBguVNcwtexVz2PybCHINrpDGr+duMqP4t0w83o/eGm/njzs4cyNPjf1L+JI4aSdKqNCp36gRA0qo0klalnVXRL86kSZNYtWoVTZo0oXnz5hw4cIAPPvgAKLjPP3v2bA4ePEj16tXp0qUL8+fPP+mMPzk5mcWLFwcsnz/ojF8kUAJ15vzqVQUFf0zpnlaJA34DDN5/hLEf/8CEFSm8d+QC7u/dgl92a6b/FCRk5efn8/7777N06VLq168PwIwZMxg7dix33HEH3bt359e//jW9evUCoH379ixcuJCdO3eSnJwMwMiRI+ncuTNXXXUVF1xwAQBvvvkml156KXXq1PHk6zqVzvhFQlSd+BjGDb6AKSO60bhGDGOmrabP0zOYvWaX19FEfiaiXr2AH2P27NnUr1//eNEHuOiii1i5ciXbt2+nW7durF+//njHvYiICGrVqkVKSgphYQXltHbt2rzzzjs8+OCDtG7dmqSkJObMmUPVqlUDnt9X+uVepCKKb+i3XZ3bsDofjerJJz9s5clpK7n1lUVc3LI6T1zTgYY1KvvtOCJnI+qEYny2MjNP7tdy2223cdtttwGwcOHCk94LDw9n+/btx5dPfSpm5syZP9t/165dmTNnjn/CBoAKv0hFVK2xX3dnZvTv0IDLkuvywvTV/GvuRi7520yGdWvI/X3aUikq3K/HEymJxm+87nWEoKJL/SJyXExkOKP7tmXmby6md+sa/GvuT/T6y3Q++X6Lnv8XCRIq/CLyM3XjK/HikC5MuqMzcdHhjHznB256cS5rd5Wf0cdEpHSCuvBrrH6Rs9O1RSJfjL6Y/9e3Fcu2ZdL3mdn8ceoysrI1/r+Unq4e+VdJv59BXfg1Vr/I2YsID2N4r5bM+u0lXJ6UyIR5P3HxX6fz1YrtxW8scoqYmBgyMjJU/P3EOUdGRgYxMTE+b6POfSLik8Qq0Ywf3Jmbf9zF//twKXe8sYSLW1bnj9d3pG58Jd924uNgQxK8GjRowJYtW0hPT/c6StCIiYmhQYMGPq9vofBbV0pKijs2HKOInL3s3HxenLGa8TM3EGYw+rKWDLuwBeFhBhQMtfqz0dXOcuAhESkZM1vsnEs5tV1n/CJyMh/mAIgC7gV+0agqD+/owZOfOT6ePps/1Z1FckwGu8evJzHus5M3KuEcAyWmOQZEfBLU9/hFJLAaRR3g9Yaf8mydL9ieW4UBm27goVnnnLzSWc4xICL+pUv9IuIXm//+LJkTXvxZ+/GZ1HSpX6RMnelSv874RcQvGo6+j6RVaeyZOhuAK37xN14a8yZRd9zlcTIROZEKv4j4VfeWiQDc3q0hH/2wg95/+5ovjz3658c5BkSkdFT4RcTvEkaM4NEB7floRHfiYyK4840l3LftEvZWaeV1NJGQp8IvIn537FG+9g2q8dkDFzOiZxOmHWzJZetv4rNl2/x6rKGfD2Xo50P9uk+RYKbCLyIBFRURxm+uSGZK4/epFZHJr976js5Pj+eXU4d7HU0kJKnwi0iZSI7Zw5Qm/2Fkz6ak727EN99epmF/RTygwi8iZSbS8nnwiracf97XREYe5Y43lvDApFQOHsnxOppIyFDhF5EyV6XKPjqlzODO7o2YsnQnfZ6eycJ1u72OJRISVPhFxBPbD2/m4f7teO+uroSFwaCXvmHMR0s5klPyKX+3Zm4NQEKR4KTCLyKe2H6o4P5+SpMa/Hf0xdxwXl0mLtzMVc/OYtWOA6Xal4gUL6iH7DWz/kD/Fi1a3Pnjjz96HUckJJzx0bodywr+rNOOVXtWkZmTSUrtk0cT3b27NqtWp5CbG0WzZsto2GAtZkUf70z7Op1X+77qy5cgEhRCcnY+59xUYGpKSsqdXmcREdhKLtt3/m/ejNTCv9eNrUv9uPokJOykc9WvSFvVkXXrOpCRUYe2SalERx/5+b4yt550pn/qvkTk9IL6jP8YTdIjUg4cm4536DSGfj6U1J2pLBuy7LSrOud4a+FGxk5LIyoijL9efy6Xn1P3tOsWty+RUKVJekSkwjAzftm1KdPuu4i6VaO5680l/G7ydyd1/Et/flxAjq2RACXYqfCLiCfqxp7+DP5EzRPjmHpvT4Z2bcg7qdu44h8zj3f82z1+fIn2JSIFVPhFxBO+3oePigjjsavb89rQTuw/ksuA5+fy7F1jSrUvEQnyzn0iEjx6tq7Ff2LTOPz2v463pbVJ4rfAvMvrQ1/vsolUJDrjF5EKo/Gv76f1ypXMeGEKAEMHP8vDj1/JvCsaepxMpOLQGb+IlI2h0/yym7Aw455LWpEGGMZ33/eiWdOlOOew4h76FxGd8YtI2Xu176tnPZhOwogRfD66Fxe3TGTd+g4MfWUh+w5n+ymhSPBS4ReRCilx1EjiK0Xy8tALeOTK1sxdt4fL/zGTJZv2nPW+Nfa/BDMVfhGp0MyMOy5qweS7uxFmxo3/WsDLc9ZxNoOTaex/CWYauU9Egsb+wznc+3Yqs9bu4fKkBJ6+6Xziov/XlcmXgXk09r8EC43cJyJBL75yJK8O68KDl7Xgy1W7ufKZmazecdCnbbdmbiV1ZyqZOZlAwdj/qTtTddlfgo7O+EUkKC1Yt5sRby3mcHYeT13bjmvO8+2RP439L8FCZ/wiElK6Nk/g8/t7klQnjgfeW8ojH/5ATl6+17FEPKfCLyJBq1bVGN77VQ+GXNCAN7/dwvUvzGXXgZ9P8Xsqjf0vwUyFX0SCWmR4GI9fcy7PDjyX1Tsz6fvMLBZtyChyG439L8FMhV9EQsLVHRvw8agLqRwZzqCXFjJx3vqzeuRPpKJS4ReRkNGqdhWm3d+TLk2rM2ZqGqPfWcKRnDyvY4mUKRV+EQkp8ZUief32rtxzURP+88MOrhs/hx37i7/vLxIsVPhFJOSEhRm/vTKZf97SkfUZh7ny2Vks3nj2Q/2KVAQq/CISsq5oV48pIy+kUmQYN01YwKRvNh5/L/35cd4FEwmgoC78ZtbfzCbs37/f6ygiUk61ql2Faff15PxG8Tz0nxWsWXMu+fnG7vHj/bL/oZ8P9WmoYJGyEtSF3zk31Tk3PD4+3usoIlKOVascxZt3dmNIl4Zs3daSg1908zqSSMBEFL+KiEjwiwgPY+Tm2dz00f/O9NPaJAGQMGIEiaNGehXtJMeuHmiCICmtoD7jFxEpicRRI0lalUbmZ3MBuP7Gv7Pzo5nlpuiL+IMKv4jIKTo1rQlA3arRDHstlZfnrPU4kYj/qPCLiJxGwogRTBl1ERe1qMEfpq3mocnfkatJfiQIqPCLiJxG4qiRxEZH8PLQLgzr1ohJqdsY/NICDhzJ8TqayFlR4RcRKUJ4mPH7Ae344y+S+XbTPgY8N5vNew57HUuk1FT4RUR8cKKmBiUAACAASURBVHOXJrxxe2cyDmXT/7nZLNmkkf6kYiq28JtZAzN70MymmNkiM5ttZi+Y2VVmpl8cRCRkdGuRyEcjLyQ2OpybJizk4++3eB1JpMSKLNxm9irwCpAN/BkYBNwDfAX0Beaa2UWBDikiUl40T4xj6r09SaoTy73v/MBzX64qdnrfrZlbyyidSPGKG8Dnaefc8tO0Lwc+NLMooJH/Y4mIlF81YqN49+4e3D9pMX+fvo6f9hzmT9d3ICL89OdS2w9tL+OEImdWZOE/Q9E/8f1sQA+4ikjIiYkM50jiP2nUsC2Tv4PpGxdzTvK3RETknrTeqj2rAHwar1+j8UlZKHbIXjPrCvwSuBCoC2RRcMY/DXjTOacZcEQkJJlB8+YriYk5zJofz2PJdxdxbvv5REcfYWvm1pPO9FN3pgJQN7Yu9ePqexVZpOjCb2afAduAKcCTwC4gBmgFXAxMMbO/O+c+DnRQEZHy5sQz9OkrdzDi7SWsXXYNb97ZlVa1qwAFZ/qpO1NZNmSZVzFFTlJcr/zBzrnbnXMfO+e2OedynXOZzrklzrmnnXO9gPllkFNEpFzr3bYOk3/Vjdx8x7Xj57Jw3e6AHUudBeVsFFn4nXO7Acws9tije2bWyswGmFnkieuIiIS6c+pXY8rIC6lROZLBL3/L1MLH/erG1vXrcdRZUM6Gr9PyzgYuNLPqwHQgFRgI3BKoYCIiFVHDGpWZMuoihry8gHvf+YFmzZvTqKFv2/rSAVCdBeVs+Vr4zTl32MxuB553zv3FzL4LZDARkYqqemwU7/2qB/e88S1fr+lIdnYMzjnMrNT7VGdB8RefC39h7/5bgNtLuK2ISMiJiQznpdu60OP559i8OYn73k7l6ZvOJ/IMz/qDb2fo6iwoZ8vXIXfvAx4C/uOcW2FmzYAZgYslIlLxhYcZrVt9T5PGy/l42S6GvryAw9m5xW8oEkA+FX7n3Gzn3ADn3J8Ll9c75+4NbDQRkYrPDJo2XcUTA5KYv2EfN7wwl72Hss9qn/7uLCihRZPsiIiUgVu7NWP8oA6s2XWYa8bNZvv+rFLvS/f05Wyo8IuIlJEr2tfntWGd2JWZwy+en8PaXQe9jiQhSIVfRKQMdWuRyHt3dyU7L5/rXpjHsrF/9TqShJhiC7+ZXW5mt5tZk1PahwUqlIhIMDunfjU+HHEhlSLDiHjzFeas2eV1JAkhRRZ+M/sj8DDQDphuZqNOeHtkIIOJiASzpgmxjF8zGYBhE1P5bNk2jxNJqCjujL8/cIlz7n7gfOAKM/tH4XulH4lCRCSEpT8/jrQ2SVT+Zg4AH3/4a5rc0JuZvxvrcTIJBcUV/gjnXC6Ac24fBb8IVDWz94GoQIcTEanoXu376s8G5kkcNZKkVWkkrUoDoOEPy3j03gncRkf+PetHL2JKCCmu8K8zs57HFpxzec6524HVQFJAk4mIhIi46AjevLMbPVtUZ+xna/jHFyu9jiRBrLjCfwPw7amNzrlHAB+nnRARkTNJGDECKBji9+WhXeiXnMizMzbwh4+X4pzzOJ0EoyLH23fOHR9hwszaA01O2ebDwMQSEQkNiaP+1086IjyM527pRKX3l/Dy/M0cOprLH6/rSFiYulSJ//g00Y6ZvQK0B1YA+YXNDhV+ERG/Cgsz/nLjecRGL2Xiwi0czs7j7zedT0QRk/uIlISvM+x1cc61DWgSEREBwMx47Or2xEZHMH7WRrKyv+WFWzsXObOfiK98/Ve0wMxU+EVEyoiZ8ZsrkvntZS34cnUGd7y6kKO5eV7HkiDg6xn/axQU/x3AUQqe4XfOufYBSyYiItzTuzVREWGM/WwNt/17AZENwwgPzy9+Q5Ez8LXwvwIMBpbxv3v8ZcLMmlEwemC8c+76M7WJiASrO3q2pFJkOI98nMZ5uUN4/c5uXkeSCszXS/0/Oec+ds5tcM5tOvYqbiMze8XMdpnZ8lPa+5rZajNba2a/K2ofzrn1hWMHFNkmIhLMbunWjL9cl8x3Ww8y6MW5HDyS43UkqaB8PeNfZWZvA1MpuNQPgHOuuF79E4FxwOvHGswsHBgPXAZsARaZ2cdAOPDUKdsPc85p9goREeCGTk2IjgjngfeWctOLc5l0dw+qxkR6HUsqGF8LfyUKCn6fE9qKfZzPOTf71Fn9gM7AWufcegAzewe42jn3FNDPxzzFMrPhwHCARo0a+Wu3IiKeGtCxIWFm3P/eD9z4whzevftC4iur+IvvfLrU75wbeppXaaflrQ9sPmF5S2HbaZlZTTN7EehoZg+dqe00mSc451KccymJiYmljCoiUv7069CAcTd1YO3uLK5/YTZ7D2V7HUkqEJ8Kv5m9ZmbVTliuXjioT2mcbgiqM45L6ZzLcM7d7ZxrXnhV4LRtIiKhpG/7+vzzlo5s3HOEG16Yo+IvPvO1c1/7wtn5AHDO7QU6lvKYWzh5nP8GgCaiFhEpocuS6/GvX57Hpr0FxX+Pir/4wNfCH2Zm1Y8tmFkNfO8fcKpFQEsza2pmUcBNwMel3JeISEi7pG1dFX8pEV8L/9PAfDP7g5k9AcwH/lLcRmY2CVgAtDazLWZ2u3MuFxgJfAGkAe8551aULr6IiFzSti7/Gnwem/cdVfGXYpmv0z4WDtl7CQX36Kc758r9hNFm1h/o36JFizt//PFHr+OIiATUjLQd3PXmEhpWi+b9ey6kRmyU15HEQ2a22DmX8rP2ogq/mcU55zKL2XGx63gtJSXFpaameh1DRCTgZqTtYPgbi2lcI4bJ91xItcoq/qHqTIW/uEv9U8zsaTO7yMxiT9hZMzO73cy+APr6O6yIiJTOxUl1ePGWjmwq7O2//7BG+JOTFVn4nXO9genAXcAKMztgZhnAm0AdYIhzbnLgY4qIiK96J9fjhZs7smHPEW7852z2Z+WQ/vw4r2NJOeHzPf6KTJf6RSQUfb50CyPf+YHmh9N55vM/k7QqzetIUoZKe6lfREQqqL7tG/DcwHNZG10TgENHcz1OJOVBaZ/FFxGRci79+XE0HT+eTwqXfzq3HQAJI0aQOGqkd8HEU0F9xm9m/c1swv79+72OIiJS5hJHjSRpVRqVO3UC4Kpr/sajo16kyt2/8jiZeKnIwm9mqWb2rJn1NbOYsgrlL865qc654fHx8V5HERHxTOM3CmZGHzugDYs3ZzL05flk5+Z7nEq8UtwZfxfgP0AvYJaZfWpm95lZq4AnExERv0kYMYKbuzbnsatasWDjAYZPXEBunop/KCpRr34zqwtcQcGz+y2Ahc65ewKUzW/Uq19E5H/+OX0Vf/5yHX3b1GD8rV0IDzvdpKlS0fmlV79zbrtz7hXn3I1ACvCWvwKKiEjZ+FXvNtzbszGfr9rD6EmLCIXHuuV/St25zzmX75yb588wIiJSNkZfcQ53dG3AlGXpPPT+EhX/EKLH+UREQtTDA9qTlZPHW6nbiYteyiNXn+t1JCkDRRZ+M0txzunmuIhIEDIzxl7XkUNHc/n3gi3ExURw/+XJXseSACvuUv9LZvajmT1ROC1vhaLn+EVEimZmPD2oE5e1qs4zMzby0ozVXkeSACtukp6OQD8gD5hsZt+b2f+ZWeMySXeW9By/iEjxwsOMF4Z0oXuTqjz5xVrenr/O60gSQMV27nPOrXbOPe6cawsMAaoBX5uZOvaJiASJyPAwXr69G+fVj+WRqav45LufvI4kAeJzr34zCwNqAbWBWCA9UKFERKTsxUSG8/rwHrRKiOGB95cxa9V2ryNJABRb+M3sQjN7AdgC/AaYC7R2zv0i0OFERKRsxUVH8PZdPahfNYq73/yO1A06xws2xY3Vvxn4E5AGdHTO9SkcwEe95UREglSNuGgm3d2d+Jhwhk1MJW3bPq8jiR8Vd8bfwznX3Tn3vHNup5nFlkkqERHxVN1qlXnrzq6Ehxm//PdCftqd6XUk8ZPievVvAjCzrma2koIzf8zs3MLL/yIiEqSa167KG8M6cyQ3n5snzCf9QJbXkcQPfO3c9wxwOZAB4Jz7AbgoUKH8Rc/xi4icnXMa1uClX57HrkO53DJhHgeP5HgdSc6Sz736nXObT2nK83MWv9Nz/CIiZ69bqzo8e8M5rM04yq0T5nIkp9z/9y9FKK5zX5fCv242s26AM7MoM3uQwsv+IiIS/K7o0Iix/Vvz3bbD3PXqfPLyNalPRVXcGf+x+/h3AyOA+hQ81tehcFlERELEzd1a8OAlTZi1/gC/fvtbzehXQfk0O59zbjdwS4CziIhIOTeyTzLpB4/w2qId1Jr6A/9vQAevI0kJFVf4m5nZx2d60zk3wM95RESknBtz7XnszlzIhPlbSawSzZ0XJwGwafCtZG/dSsuvp3ucUIpSXOFPB54uiyAiIlIxmBnP/vIC9kyYwx+/WE/NuGiu7dQMgNxt2zxOJ8UprvBnOudmlUkSERGpMCLCw3h5WHeuHz+L3/4njZpxMTTxOpT4pLjOfRvKJIWIiFQ4laMjeHN4d4av+4pa1/Tm8KJFAKS1SSKtTRLpz4/zOKGcTnFn/H8v6k0zqwo0cs4t918kERGpKGrExTBw3BNc88Ll2IEDvDntMZJW6Wnv8qy4M/7rzGy+mf3ezK4ys85mdpGZDTOzN4BPgEplkFNERMqpxolVePW2FLLCIgHYe+iox4mkKEWe8TvnHjCz6sD1wA1AXSCLgsF7/uWcmxv4iKVnZv2B/i1atPA6iohIUDu3cQLjbuvE2z9exvcvzeP9ey4iJsqnJ8aljFkoDMCQkpLiUlNTvY4hIhL03l3wI/83ZQ29mlflldt7EBZmXkcKWWa22DmXcmq7z2P1i4iIFGdg15bc17MhM9cd4P+9v8jrOHIaug4jIiJ+dX/fdmzbl8U736VTv9pyRl1+jteR5AQ64xcREb8yM/40sDM9mlTh7zM28cE367yOJCfwqfCbWWUze9TMXipcbmlm/QIbTUREKqrwMOOlYd1okxjD76asYu4qjehXXvh6xv8qcBToWri8BRgbkEQiIhIUKkVF8Pqd3UiMjeCeST+wZtteryMJvhf+5s65vwA5AM65LEBdNUVEpEiJVSsxcWhnAG579Vt27T/scSLxtfBnm1klwAGYWXMKrgCIiIgUqVW96vzz5nNJP5TLbS8vIOtojteRQpqvhf8x4HOgoZm9BUwHfhuwVCIiUiaGfj6UoZ8PDfhxureux5P9W5G26wh3v7aA/PzgH0OmvPKp8DvnvgSuBW4DJgEpzrmZgYslIiLB5sauLRl5YQNmrT/Io5P1jL9XinyO38zOO6Vpe+GfjcyskXNuSWBiiYhIMBp9ZXt+2nOYt5ak07jmCob3TvY6UsgpbgCfpwv/jAFSgB8o6NTXHvgG6BG4aGdPY/WLiJQvZsbfBnVm6z9n86evNtKoZix9OzTxOlZIKfJSv3PuYufcxcAm4DznXIpz7nygI7C2LAKeDefcVOfc8Pj4eK+jiIhIociIcP49rBuN4iMZ/cFKftiY7nWkkOJr5742zrllxxacc8uBDoGJJCIiwa5abDQTh11AdIRx5xuL2bYn0+tIIcPXwp9mZv82s15m1rNwBL+0QAYTEZHg1qRWPC/e3IH9R/IZ+spCDusxvzLha+EfCqwA7gPuB1YWtomIiJTaBS3r8scBrViz+yi/mqjH/MqCr4/zHXHO/cM5d03h6x/OuSOBDiciIsHvugtacE/3eszacJDHP0z1Ok7Q83WSnu5m9qWZrTGz9cdegQ4nIiKh4cF+HbiidTyvpe7itVm6kxxIvl7qfxn4OwWP73U64SUiInLWzIx/3NKFc+vE8Icv1jNzxWYANg2+lU2Db/U4XXDxtfDvd8595pzb5ZzLOPYKaDIREQkpMVER/HtoV2rFRnDve8tZu32P15GCkq+Ff4aZ/dXMuprZecdeAU0mIiIhJzG+Mi8NPp+8fMftr6WSSTjZW7d6HSuoFDdy3zEXFP6ZckKbAy7xbxwREQl1yY0S+Nu1SYx4byWPV7uAxxct9DpSUPGp8BeO3iciIlImrujYlNE7D/C3mV4nCT7FTdIzuqj3nXN/928cERERSH9+HL3Hj6d34XJamyQAEkaMIHHUSO+CBYHi7vFXKXylAL8C6he+7gbaBjaaiIiEqsRRI0lalUZkp84A9L/mb+x4/wsVfT8o8ozfOfc4gJn9l4JJeg4WLo8B3g94OhERCWmROHKAWnER3Dd5BR/UjKNlvRpex6rQfO3V3wjIPmE5G2ji9zQiIiKniKhXj3/9siO5+Y7hb6SyLzPL60gVmq+F/w3gWzMbY2aPAd8ArwUuloiISIGo+vVp17gWfx7Qik17c7jn9W/Iy8v3OlaF5Wuv/ifN7DPgwsKmoc657wIXyz/MrD/Qv0WLFl5HERGRUmj8xuvH/96/UwtW7zjAuHnbefzDVJ64obOHySouX8/4cc4tcc49W/gq90UfwDk31Tk3PD4+3usoIiLiB6Ov6sDlLavy+uJ03p6jMf1Lw+fCLyIi4rWwsDD+fssFtE6I4vHPN/DtGo3qV1Iq/CIiUqHExkQx4dZOxEUZI95ZxraMA15HqlBU+EVEpMJpXKsaz96QzL4jeQx/7VuOZOd4HanCUOEXEZEKqUfbRjx0aSOW7zrKbyZ9g3PO60gVggq/iIhUWMMuaccN7aozNW0/L/53qddxKgQVfhERqdDG3tCJjnVjeHrWFmYs2+R1nHJPhV9ERCq06KhIXhjciZqVwhn9wUo27tzrdaRyTYVfREQqvLo1qvL8wHPIynHc9XoqmVlHvY5Ubqnwi4hIUOjcqgGP9GnM6oxsHpz0rTr7nYEKv4iIBI1f9kxm0Lk1+HzNAcZ//oPXccolFX4REQkqj1/fifPqxvDM7K18vXSD13HKHRV+EREJKlGREYwf3ImalcP59Yer2Lhjj9eRyhUVfhERCTp1a1Tl2RuSOZzj+NVbizmkzn7HqfCLiEhQ6tKmIb/r3ZC09Gx+9+4idfYrpMIvIiJB67aLz+G65GpMXbWfl77UyH6gwi8iIkHMzBh7QwrtakXzt5lbmLfyJ68jeU6FX0REglqlmGjG3XIecdFhjP5gJdsy9nsdyVMq/CIiEvQa167BX69uRUZWHqPeSiU7hKfxVeEXEZGQ0LtDc0Z1q8PibUd44sPUkO3sp8IvIiIhY+QVHbiseRxvfr+HyfNXeR3HEyr8IiISMsLDw/nrTSk0rRbB459vYPmG7V5HKnNBXfjNrL+ZTdi/P7Q7coiIyP9UqxLLcwPb44B731vGvoOHvI5UpoK68DvnpjrnhsfHx3sdRUREypF2TevyaJ/GrN+bw2/eSSU/P9/rSGUmqAu/iIjImQzs0ZZB7avz5bpMxn/+nddxyowKv4hIiNuaudXrCJ4wM35/7fl0qBPNc3N3MGvpeq8jlQkVfhGRELf9UOh1cDumUkw0zwzsQNXoMH47ZQ1b0/d6HSngLBSeY0xJSXGpqalexxARKVNDPx9a7Dqr9qwiMyeTlNopxa77at9X/RGrXPpi8Y/cM3kN59eL4Y3hFxIdHeV1pLNmZoudcz/7YHXGLyISgrZmbiV1ZyqZOZkApO5MJXVnashe9r/8/JYM75zIt1uP8NepS4J6cJ8IrwOIiEhg+HKGPvTzoaTuTGXZkGVlkKh8G92vI0u3z+GV1Aw6NFpDv86tvY4UEDrjFxERASIjI3l64HnUqRLOI9PWs3bLLq8jBYQKv4hIiKsbW9frCOVGnZrV+MvVrTmcm88D7/3A4awjXkfyOxV+EZEQVz+uvtcRypUe5zTlvu51WLYrmz/8Z3HQ3e9X4RcRETnFXX3O5dJmsbyzdB+T567wOo5fqfCLiIicIiIigqdu6EjD+Aj+8OVPrNy4zetIfqPCLyIichqJ1eP569Wtyc6H0ZNXcPDQYa8j+YUKv4iIyBlc0LYJoy+sw6rd2Tz+4eKgmMxHhV9ERKQIw3q3p2+LOCavOMC7c5Z7HeesqfCLiIgUISIigrHXdaBJtQj+OH0Ly9dX7NENVfhFRETOYNPgW9k0+FYSqsfzl6tbk5MPD36wggOZh7yOVmoq/CIiIj7onNSEB3rUZlVGDo//Z0mFvd+vwi8iIuKj2y89l8tbxPLBigO8P6diPt+vwi8iIuKjiIgIxl7bgUbxETw5fTNpFfD5fhV+ERGREkisUY0/9W/J0Tz4zYcryKxgz/er8IuIiJRQ1+SmjOqayPJd2Tw15bsKdb9fhV9ERKSEzIzhfc6lV5NKTFq6jykLVnodyWcq/CIiIqUQFRXFH689lzpx4Tzx5WbWbt7hdSSfqPCLiIgUIXvrmQfsqVerJmOvbEpmdj7/9+FysrKyyjBZ6ajwi4iIFCF3W9E99y/u2Io7UmqyePtR/j7t+3J/v1+FX0RE5CyYGfdecS6d68fw6uI9fLVkjdeRiqTCLyIicor058eR1iaJw4sWAZDWJom0NkmkPz/utOtXqlSJp36RTHxMGL//bCObd6SXZdwSUeEXERE5ReKokSStSqNyp04AJK1KI2lVGomjRp5xm+YN6/DopQ3ZdSiPh/+zjKNHj5ZV3BJR4RcREfGT/l2SuKl9PLM3ZfHSl0txznkd6WdU+EVERIoQUa+ez+uGh4fzf/3OpW1iFOMWpPPNyg0BTFY6KvwiIiJFiKpfv0Trx1etwpP9WhIRBg9/8iMZe/cFKFnpqPCLiIj4WYdWjRl9YW3W7c1l7MdLycnJOeO6mwbfyqbBt5ZZNhV+ERERPzMzftmrHX1bxPKftINMnlt+pvBV4RcREQmAqKgoft8/mYZVI/jzzO2kbdjidSSgghR+M2tmZi+b2eQT2n5hZi+Z2RQz6+NlPhERkdOpVzuRMX0acSgnn4c/XsWhQ4e8jhT4wm9mr5jZLjNbfkp7XzNbbWZrzex3Re3DObfeOXf7KW0fOefuBG4DBvo9uIiIiB/06tiKYedXZ8n2ozz3+Q+eD+lbFmf8E4G+JzaYWTgwHrgCaAsMMrO2ZtbOzD455VWrmP0/UrgvERGRcic8PJxRl7fn/LrRvLJ4LzO/83ZI34AXfufcbGDPKc2dgbWFZ/LZwDvA1c65Zc65fqe8dp1uv1bgz8Bnzrklp3l/uJmlmllqenr5HTpRRESCX1xcHI9f1ZLYyDDGfLGJbTu9q0te3eOvD2w+YXlLYdtpmVlNM3sR6GhmDxU2jwIuBa43s7tP3cY5N8E5l+KcS0lMTPRjdBERCRWN33idxm+87pd9JTdvxIMX1uanA7mM/WQF2dnZftlvSUV4clSw07SdcVxD51wGcPcpbc8Bz/k5l4iISECYGTdcmMw3Px1g6ppDdJm9nMG9O5Z5Dq/O+LcADU9YbgAUPeGxiIhIBRcdHc1DVyXTOD6Cp+fsZOW6zcVv5GdeFf5FQEsza2pmUcBNwMceZRERESkz9Won8vvLGnE4J59Hp60hNy+vTI9fFo/zTQIWAK3NbIuZ3e6cywVGAl8AacB7zrnyM6yRiIhIAPXs0JKhHauxZPtRPohugjvz3W6/C/g9fufcoDO0fwp8Gshjm1l/oH+LFi0CeRgREZESiYiIYESfc1i09Vv+nd+Z87Zn0KSMjl0hRu4rLefcVOfc8Pj4eK+jiIiInCQ+Pp7H+jYnJj+HpxK6czS3bC75B3XhFxERKc/at2rCA7vm0X39t3z3U9lM36vCLyIi4pGwsDB65qUzaPVXdGlWs2yOWSZHERERkdMKDyvbUqzCLyIi4oH058eR1iaJw4sWAZDWJom0NkmkPz8uoMf1auQ+ERGRkJY4aiSJo0ayafCtHF60iKRVaWVy3KA+4zez/mY2Yf/+/V5HERERKReCuvDrcT4REakIIurVK7NjBXXhFxERqQii6p9xglq/U+EXEREJISr8IiIiIUSFX0REJISo8IuIiIQQFX4REZEQEtSFX8/xi4iInCyoC7+e4xcRETlZUBd+EREROZkKv4iISAhR4RcREQkhKvwiIiIhRIVfREQkhKjwi4iIhBAVfhERkRAS1IVfA/iIiIicLKgLvwbwEREROVlQF34RERE5mQq/iIhICFHhFxERCSHmnPM6Q8ClpKS41NRUr2OIiIiUGTNb7JxLObVdZ/wiIiIhRIVfREQkhKjwi4iIhBAVfhERkRAS1IVfI/eJiIicLKgLv0buExEROVlQF34RERE5mQq/iIhICFHhFxERCSEq/CIiIiFEhV9ERCSEqPCLiIiEEBV+ERGREKLCLyIiEkJCYlpeM0sHNp3SHA+cbki/07Wfri0B2O2XgCVzptyB3ocv2xS3TlHv+/p9P117Rf4sSrMfX9cv7eehnw3/b6OfjbLZj342TtbYOZf4s1bnXEi+gAm+tp+hLbU85Q70PnzZprh1inrf1+/76dor8mdRmv34un5pPw/9bPh/G/1slM1+9LPh2yuUL/VPLUH7mdb1gj+ylGYfvmxT3DpFvV+S73t5+Tz8laOk+/F1/dJ+HvrZ8P82+tkom/3oZ8MHIXGpPxDMLNU5l+J1DtFnUd7o8yg/9FmUL+Xl8wjlM/6zNcHrAHKcPovyRZ9H+aHPonwpF5+HzvhFRERCiM74RUREQogKv4iISAhR4RcREQkhKvx+Zma/MLOXzGyKmfXxOk+oM7NmZvaymU32OksoMrNYM3ut8GfiFq/zhDr9PJQfXtYKFf4TmNkrZrbLzJaf0t7XzFab2Voz+11R+3DOfeScuxO4DRgYwLhBz0+fx3rn3O2BTRpaSvi5XAtMLvyZGFDmYUNAST4P/TwEVgk/C89qhQr/ySYCfU9sMLNwYDxwBdAWGGRmbc2snZl9csqr1gmbPlK4nZTeRPz3eYj/TMTHzwVoAGwuXC2vDDOGkon4/nlIYE2k5J9FmdeKiLI8WHnnnJttZk1Oae4MrHXOrQcws//f3r2GWFHGcRz//rI105DogqhlVKjR3UzCCrpgvcmSUkoqSJHAIl8YQVRvComMeJwJBQAABJtJREFUIslAg8qgXngphMrICyrdKBNNV82EiF5YRjfLNEu0fy/OIw6nmd02z+457vP7wHJm5nlmzn/3z5z/PDM7ZxYBEyLiKWB8/TYkCZgNvBcRG7s34t6tEfmwxutKXoCd1Ir/JjzQ6BZdzMcXPRtdXrqSC0nbaVKt8I7YuaEcGbFA7YNsaAf9ZwDjgEmSpndnYJnqUj4knSrpRWCUpEe6O7iMVeVlKTBR0nxa+CtMe6HSfHh/aIqqfaNptcIj/s6pZFnltx5FxFxgbveFk72u5uNnwAdg3a80LxGxD5ja08FYZT68P/S8qlw0rVZ4xN+5ncCZhfkzgO+aFIs5H63KeWktzkfraLlcuPB3bj0wXNLZkvoCk4G3mxxTzpyP1uS8tBbno3W0XC5c+AskLQQ+AUZK2ilpWkQcBB4AVgDbgSURsa2ZcebC+WhNzktrcT5ax7GSCz+kx8zMLCMe8ZuZmWXEhd/MzCwjLvxmZmYZceE3MzPLiAu/mZlZRlz4zczMMuLCb5YBSYckbSr8dPg4456imjWSzirE9r2kbwvzfevWmZLuly4uO03Sj5JOkLRI0vCe/U3Mjh2+j98sA5L2RsRJDd7m8enLSY5mGzcB4yJiZmHZ48DeiHi2Yp2BwNfAsIj4Iy2bDoyJiGmSrgHuTs86N7M6HvGbZUzSN5KekLRR0hZJ56XlAyQtkLRe0ueSJqTlUyS9IekdYKWk/pKWSGqXtFjSOkmXS5omaU7hfe6V9FxJCHcBb3UQ32hJ70vaIGmFpMERsQf4ALi50HUycPgswIfAOEl+CJlZCRd+szycWHeq/45C208RcRkwH3goLXsMWBMRY4DrgGckDUhtY4F7IuJ64H5gd0RcDMwCRqc+i4BbJLWl+anAqyVxXQVsKAs4rfsCMCkiRgMLgCdT80JqxR5JQ4ARwFqAiPgb+Aq45D/8Xcyy4yNiszzsj4hLK9qWptcNwG1p+kZqhfvwgUA/YFiaXhURv6Tpq4HnASJiq6T2NL1P0hpgvKTtQFtEbCl571Mi4veKuEYCFwKrJAH0AXaltmXAvHTa/3bgzYg4VFj3B2AIFQcVZjlz4Tezv9LrIY58JgiYGBE7ih0lXQHsKy7qYLsvA48CX1I+2gc4KOm4NEqvJ2BbRIytb4iI/ZKWA7dSG/nPrOvSD9jfQWxm2fKpfjMrswKYoTTUljSqot9H1EbcSDofuOhwQ0Sso/Yc8js5cv293g7gnA7aTpc0Nm2/TdIFhfaFwIPAIODTunVHAH4anVkJF36zPNRf45/dSf9ZQBvQLmlrmi8zj1pxbgceBtqB3wrtS4CPI2J3xfrvAteWNUTEAWAS8LSkzcAm4MpCl5XUTucvjsLtSZIGUbu0sQsz+xffzmdm/5ukPtSu3/8p6VxgNTAiFW0kLQPmRMTqivUHA69FxA0NjGkmsCciXmnUNs16E1/jN7Oj0R9Ym/4DX8B9EXFA0snAZ8DmqqIPEBG7JL0kaWC6Ta8RfgVeb9C2zHodj/jNzMwy4mv8ZmZmGXHhNzMzy4gLv5mZWUZc+M3MzDLiwm9mZpYRF34zM7OM/AMhZOAUe3Q4zwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 576x432 with 1 Axes>"
]
},
"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
}