{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Likelihood Analysis with Python" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The python likelihood tools are a very powerful set of analysis tools that expand upon the command line tools provided with the Fermi Science Tools package. Not only can you perform all of the same likelihood analysis with the python tools that you can with the standard command line tools but you can directly access all of the model parameters. You can more easily script a standard analysis like light curve generation. There are also a few things built into the python tools that are not available from the command line like the calculation of upper limits.\n", "\n", "There are many [user contributed packages](https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/) built upon the python backbone of the Science Tools and we are going to highlight and use a few of them in this tutorial.\n", "\n", "This sample analysis will look at 2 years of data for the object 3C 279, which is a relatively birght example of a type of Active Galactic Nucleus (AGN) known as a blazar. At certain points we will refer to this article as well as the [Cicerone](https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/). This tutorial is based on the \"[Binned Likelihood Tutorial](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/binned_likelihood_tutorial.html)\" Analysis Thread, which is hosted on the Fermi Science Support Center website. \n", "\n", "This tutorial assumes you have the most recent [ScienceTools](https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/) installed. We will also make significant use of python, so you might want to familiarize yourself with python (there's a beginner's guide at http://wiki.python.org/moin/BeginnersGuide. This tutorial also assumes that you've gone through the non-python based [unbinned likelihood](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/likelihood_tutorial.html) thread. Running through this notebook should take several hours to complete (depending on your computer's speed) if you do everything but there are some steps you can skip along the way which would shorten that.\n", "\n", "> **Note:** This tutorial is generated from a [jupyter notebook](https://github.com/kialio/FSSC-Docs/blob/master/Threads/Python/python_tutorial.ipynb) which you can download and run yourself (the prefereed method). You can also run individual commands listed on this page. If you do that, be aware that some commands must be executed in an ipython/jupyter environment." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Get the Data\n", "\n", "For this thread the original data were extracted from the [LAT data server](https://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi) with the following selections (these selections are similar to those in the paper):\n", "\n", "* Search Center (RA,Dec) = (193.98, -5.82)\n", "* Radius = 15 degrees\n", "* Start Time (MET) = 239557417 seconds (2008-08-04 T15:43:37)\n", "* Stop Time (MET) = 302572802 seconds (2010-08-04 T00:00:00)\n", "* Minimum Energy = 100 MeV\n", "* Maximum Energy = 500000 MeV\n", "\n", "We've provided direct links to the event files as well as the spacecraft data file if you don't want to take the time to use the download server. For more information on how to download LAT data please see the [Extract LAT Data](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/explore_latdata.html) tutorial.\n", "\n", "* [L181126210218F4F0ED2738_PH00.fits](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH00.fits)\n", "* [L181126210218F4F0ED2738_PH01.fits](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH01.fits)\n", "* [L181126210218F4F0ED2738_PH02.fits](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH02.fits)\n", "* [L181126210218F4F0ED2738_PH03.fits](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH03.fits)\n", "* [L181126210218F4F0ED2738_PH04.fits](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH04.fits)\n", "* [L181126210218F4F0ED2738_PH05.fits](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH05.fits)\n", "* [L181126210218F4F0ED2738_PH06.fits](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH06.fits)\n", "* [L181126210218F4F0ED2738_SC00.fits](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_SC00.fits)\n", "\n", "Make a working directory and then download all of the files into that directory." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "!mkdir working" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import urllib.request" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "url_base = \"https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/\"\n", "datafiles = [\"L181126210218F4F0ED2738_PH00.fits\", \n", " \"L181126210218F4F0ED2738_PH01.fits\",\n", " \"L181126210218F4F0ED2738_PH02.fits\",\n", " \"L181126210218F4F0ED2738_PH03.fits\",\n", " \"L181126210218F4F0ED2738_PH04.fits\",\n", " \"L181126210218F4F0ED2738_PH05.fits\",\n", " \"L181126210218F4F0ED2738_PH06.fits\",\n", " \"L181126210218F4F0ED2738_SC00.fits\"]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "for datafile in datafiles:\n", " urllib.request.urlretrieve(url_base+datafile,\"working/\"+datafile)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "L181126210218F4F0ED2738_PH00.fits L181126210218F4F0ED2738_PH04.fits\r\n", "L181126210218F4F0ED2738_PH01.fits L181126210218F4F0ED2738_PH05.fits\r\n", "L181126210218F4F0ED2738_PH02.fits L181126210218F4F0ED2738_PH06.fits\r\n", "L181126210218F4F0ED2738_PH03.fits L181126210218F4F0ED2738_SC00.fits\r\n" ] } ], "source": [ "ls working/" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, you'll first need to make a file list with the names of your input event files:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "ls -1 working/*PH*.fits > binned_events.txt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can rename the spacecraft file to make it easier to work with:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "mv working/L181126210218F4F0ED2738_SC00.fits working/3C279_SC.fits" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the following analysis we've assumed that you've named your list of data files binned_events.list and the spacecraft file 3C279_SC.fits." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Tue Dec 15 22:20:21 EST 2020\r\n" ] } ], "source": [ "#Let's get a timestamp, so we can see how long this takes to run on your machine.\n", "!date" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Perform Event Selections" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The gt_apps module provides methods to call the various FermiTools from within python. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "import gt_apps as my_apps" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The python object for gtselect is called filter and we first need to set all of its options. It might not seem that convenient to do it this way (as opposed to entering the options via the command line) but it's really nice once you start building up scripts and reading back options and such. For example, generating a light curve entails running the likelihood analysis for each datapoint. It'll be much easier to do all of that within python and change the tmin and tmax in an iterative fashion. Note that these python objects are just wrappers for the standalone tools so if you want any information on their options, see the corresponding [documentation](/ssc/data/analysis/scitools/references.html) for the standalone tool." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "my_apps.filter['evclass'] = 128\n", "my_apps.filter['evtype'] = 3\n", "my_apps.filter['ra'] = 'INDEF'\n", "my_apps.filter['dec'] = 'INDEF'\n", "my_apps.filter['rad'] = 'INDEF'\n", "my_apps.filter['emin'] = 100\n", "my_apps.filter['emax'] = 500000\n", "my_apps.filter['zmax'] = 90\n", "my_apps.filter['tmin'] = 'INDEF'\n", "my_apps.filter['tmax'] = 'INDEF'\n", "my_apps.filter['infile'] = '@binned_events.txt'\n", "my_apps.filter['outfile'] = 'working/3C279_binned_filtered.fits'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, run gtselect:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time -p gtselect infile=@binned_events.txt outfile=working/3C279_binned_filtered.fits ra=\"INDEF\" dec=\"INDEF\" rad=\"INDEF\" tmin=\"INDEF\" tmax=\"INDEF\" emin=100.0 emax=500000.0 zmin=0.0 zmax=90.0 evclass=128 evtype=3 convtype=-1 phasemin=0.0 phasemax=1.0 evtable=\"EVENTS\" chatter=2 clobber=yes debug=no gui=no mode=\"ql\"\n", "Done.\n", "real 3.95\n", "user 3.42\n", "sys 0.28\n" ] } ], "source": [ "my_apps.filter.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, you need to run gtmktime. This is accessed within python via the maketime object:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time -p gtmktime scfile=working/3C279_SC.fits sctable=\"SC_DATA\" filter=\"(DATA_QUAL>0)&&(LAT_CONFIG==1)\" roicut=no evfile=working/3C279_binned_filtered.fits evtable=\"EVENTS\" outfile=\"working/3C279_binned_gti.fits\" apply_filter=yes overwrite=no header_obstimes=yes tstart=0.0 tstop=0.0 gtifile=\"default\" chatter=2 clobber=yes debug=no gui=no mode=\"ql\"\n", "real 21.44\n", "user 20.43\n", "sys 0.76\n" ] } ], "source": [ "my_apps.maketime['scfile'] = 'working/3C279_SC.fits'\n", "my_apps.maketime['filter'] = '(DATA_QUAL>0)&&(LAT_CONFIG==1)'\n", "my_apps.maketime['roicut'] = 'no'\n", "my_apps.maketime['evfile'] = 'working/3C279_binned_filtered.fits'\n", "my_apps.maketime['outfile'] = 'working/3C279_binned_gti.fits'\n", "my_apps.maketime.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We're using the most conservative and most commonly used cuts described in detail in the [Cicerone](http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Likelihood/Exposure.html)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create a 3-D (binned) counts map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For binned likelihood analysis, the data input is a three-dimensional counts map with an energy axis, called a counts cube. The gtbin tool performs this task as well, by using the CCUBE option." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time -p gtbin evfile=working/3C279_binned_gti.fits scfile=none outfile=working/3C279_binned_ccube.fits algorithm=\"CCUBE\" ebinalg=\"LOG\" emin=100.0 emax=500000.0 enumbins=37 denergy=0.0 ebinfile=NONE tbinalg=\"LIN\" tstart=239557517.0 tstop=255335817.0 dtime=86400.0 tbinfile=NONE snratio=0.0 lcemin=0.0 lcemax=0.0 nxpix=100 nypix=100 binsz=0.2 coordsys=\"CEL\" xref=193.98 yref=-5.82 axisrot=0.0 rafield=\"RA\" decfield=\"DEC\" proj=\"AIT\" hpx_ordering_scheme=\"RING\" hpx_order=3 hpx_ebin=yes hpx_region=\"\" evtable=\"EVENTS\" sctable=\"SC_DATA\" efield=\"ENERGY\" tfield=\"TIME\" chatter=2 clobber=yes debug=no gui=no mode=\"ql\"\n", "This is gtbin version HEAD\n", "gtbin: WARNING: No spacecraft file: EXPOSURE keyword will be set equal to ontime.\n", "real 0.51\n", "user 0.43\n", "sys 0.05\n" ] } ], "source": [ "my_apps.counts_map['algorithm'] = 'CCUBE'\n", "my_apps.counts_map['evfile'] = 'working/3C279_binned_gti.fits'\n", "my_apps.counts_map['scfile'] = 'none'\n", "my_apps.counts_map['outfile'] = 'working/3C279_binned_ccube.fits'\n", "my_apps.counts_map['nxpix'] = 100\n", "my_apps.counts_map['nypix'] = 100\n", "my_apps.counts_map['binsz'] = 0.2\n", "my_apps.counts_map['coordsys'] = 'CEL'\n", "my_apps.counts_map['xref'] = 193.98\n", "my_apps.counts_map['yref'] = -5.82\n", "my_apps.counts_map['axisrot'] = 0\n", "my_apps.counts_map['proj'] = 'AIT'\n", "my_apps.counts_map['ebinalg'] = 'LOG'\n", "my_apps.counts_map['emin'] = 100\n", "my_apps.counts_map['emax'] = 500000\n", "my_apps.counts_map['enumbins'] = 37\n", "my_apps.counts_map.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Livetime Cubes and Exposure Maps\n", "\n", "The next two steps pre-compute the livetime over the whole sky and the exposure map over your are of intere, respectively. Doing this in advance of running the likelihood tool will save a significant amount of time when running the gtlike task itself. This can be usefull if you wish to re-do the likelihood analysis multiple times (after, say, adjusting the model inputs or trying out different fitting algorithms). So, let's go ahead and do these \"extra\" steps.\n", "\n", "### Livetime Cube\n", "\n", "We can pre-compute the livetime as a function of sky position and off-axis angle, which will speed up the likelihood analysis. This step will take approximately 30 minutes to complete, so this may be a good time to respond to that email you've been thinking about." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time -p gtltcube evfile=\"working/3C279_binned_gti.fits\" evtable=\"EVENTS\" scfile=working/3C279_SC.fits sctable=\"SC_DATA\" outfile=working/3C279_ltCube.fits dcostheta=0.025 binsz=1.0 phibins=0 tmin=0.0 tmax=0.0 file_version=\"1\" zmin=0.0 zmax=90.0 chatter=2 clobber=yes debug=no gui=no mode=\"ql\"\n", "Working on file working/3C279_SC.fits\n", ".....................!\n", "real 1838.17\n", "user 1833.83\n", "sys 2.20\n" ] } ], "source": [ "my_apps.expCube['evfile'] = 'working/3C279_binned_gti.fits'\n", "my_apps.expCube['scfile'] = 'working/3C279_SC.fits'\n", "my_apps.expCube['outfile'] = 'working/3C279_ltCube.fits'\n", "my_apps.expCube['zmax'] = 90\n", "my_apps.expCube['dcostheta'] = 0.025\n", "my_apps.expCube['binsz'] = 1\n", "my_apps.expCube.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exposure Map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This next step will apply the livetime calculated in the previous step to your region of interest. This tool generates a binned exposure map, an accounting of the exposure at each position in the sky, that are a required input to the likelihood process. \n", "\n", "The LAT has a very large Point Spread Function (PSF), especially at low energies. In order to account for contributions from sources outside our ROI we need to make the exposure map cover a larger area on the sky, equal to the ROI plus 10 degrees." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "from GtApp import GtApp\n", "expCube2= GtApp('gtexpcube2','Likelihood')" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time -p gtexpcube2 infile=working/3C279_ltCube.fits cmap=none outfile=working/3C279_BinnedExpMap.fits irfs=\"P8R3_SOURCE_V3\" evtype=3 edisp_bins=0 nxpix=1800 nypix=900 binsz=0.2 coordsys=\"CEL\" xref=193.98 yref=-5.82 axisrot=0.0 proj=\"AIT\" ebinalg=\"LOG\" emin=100.0 emax=500000.0 enumbins=37 ebinfile=\"NONE\" hpx_ordering_scheme=\"RING\" hpx_order=6 bincalc=\"EDGE\" ignorephi=no thmax=180.0 thmin=0.0 table=\"EXPOSURE\" chatter=2 clobber=yes debug=no mode=\"ql\"\n", "Computing binned exposure map....................!\n", "Using evtype=3 (i.e., FRONT/BACK irfs)\n", "real 130.87\n", "user 128.80\n", "sys 1.74\n" ] } ], "source": [ "expCube2['infile'] = 'working/3C279_ltCube.fits'\n", "expCube2['cmap'] = 'none'\n", "expCube2['outfile'] = 'working/3C279_BinnedExpMap.fits'\n", "expCube2['irfs'] = 'P8R3_SOURCE_V3'\n", "expCube2['evtype'] = '3'\n", "expCube2['nxpix'] = 1800\n", "expCube2['nypix'] = 900\n", "expCube2['binsz'] = 0.2\n", "expCube2['coordsys'] = 'CEL'\n", "expCube2['xref'] = 193.98\n", "expCube2['yref'] = -5.82\n", "expCube2['axisrot'] = 0\n", "expCube2['proj'] = 'AIT'\n", "expCube2['ebinalg'] = 'LOG'\n", "expCube2['emin'] = 100\n", "expCube2['emax'] = 500000\n", "expCube2['enumbins'] = 37\n", "expCube2.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generate XML Model File\n", "\n", "We need to create an XML file with all of the sources of interest within the Region of Interest (ROI) of 3C 279 so we can correctly model the background. For more information on the format of the model file and how to create one, see the [likelihood analysis](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/likelihood_tutorial.html#createSourceModel) tutorial. We'll use the user contributed tool `make4FGLxml.py` to create a model file based on the LAT 4-year LAT catalog. You'll need to download the FITS version of this file at [http://fermi.gsfc.nasa.gov/ssc/data/access/lat/4yr_catalog/](http://fermi.gsfc.nasa.gov/ssc/data/access/lat/4yr_catalog/) and get the [make4FGLxml.py](https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/python3/make4FGLxml.py) tool from the [user contributed software](http://fermi.gsfc.nasa.gov/ssc/data/analysis/user/) page and put them both in your working directory. Also make sure you have the most recent galactic diffuse and isotropic model files which can be found at [http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html](http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html). In the following we assume that you have the galactic diffuse and isotropic files in your FermiTools install and we just make local symbolic links." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "('make4FGLxml.py', )" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "urllib.request.urlretrieve('https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/python3/make4FGLxml.py','make4FGLxml.py')" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "('working/gll_psc_v26.xml', )" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "urllib.request.urlretrieve('https://fermi.gsfc.nasa.gov/ssc/data/access/lat/10yr_catalog/gll_psc_v26.xml','working/gll_psc_v26.xml')" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "!ln -s $FERMI_DIR/refdata/fermi/galdiffuse/iso_P8R3_SOURCE_V3_v1.txt working/iso_P8R3_SOURCE_V3_v1.txt" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "!ln -s $FERMI_DIR/refdata/fermi/galdiffuse/gll_iem_v07.fits working/gll_iem_v07.fits" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3C279_BinnedExpMap.fits L181126210218F4F0ED2738_PH02.fits\r\n", "3C279_SC.fits L181126210218F4F0ED2738_PH03.fits\r\n", "3C279_binned_ccube.fits L181126210218F4F0ED2738_PH04.fits\r\n", "3C279_binned_filtered.fits L181126210218F4F0ED2738_PH05.fits\r\n", "3C279_binned_gti.fits L181126210218F4F0ED2738_PH06.fits\r\n", "3C279_ltCube.fits \u001b[1m\u001b[35mgll_iem_v07.fits\u001b[m\u001b[m@\r\n", "L181126210218F4F0ED2738_PH00.fits gll_psc_v26.xml\r\n", "L181126210218F4F0ED2738_PH01.fits \u001b[1m\u001b[35miso_P8R3_SOURCE_V3_v1.txt\u001b[m\u001b[m@\r\n" ] } ], "source": [ "ls working" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have all of the files we need, you can generate your model file:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "This is make4FGLxml_py3 version 01r06. This version is designed to be compatible with python3.\n", "The default diffuse model files and names are for pass 8 and 4FGL and assume you have v11r5p3 of the Fermi Science Tools or higher.\n", "The default isotropic template is now for P8R3_SOUCE_V3 IRFs.\n" ] } ], "source": [ "from make4FGLxml import *" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "mymodel = srcList('working/gll_psc_v26.xml','working/3C279_binned_gti.fits','working/3C279_model.xml')" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Creating file and adding sources from 4FGL\n", "Added 241 point sources and 0 extended sources\n" ] } ], "source": [ "mymodel.makeModel('working/gll_iem_v07.fits', \n", " 'working/gll_iem_v07.fits', \n", " 'working/iso_P8R3_SOURCE_V3_v1.txt', \n", " 'working/iso_P8R3_SOURCE_V3_v1')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can open the model file with your favorite text editor to have a look at it. In the 4FGL catalog, 3C 279 is referred to as \"4FGL J1256.1-0547\". The model for this source is reproduced below. Note the parameters that say \"free=\"1\"\" - the \"1\" indicates that these parameters will be left free during the likelihood analysis, which will also be true for several of the other sources in our ROI." ] }, { "cell_type": "raw", "metadata": {}, "source": [ " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Source Map" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The [gtsrcmaps](https://raw.githubusercontent.com/fermi-lat/fermitools-fhelp/master/gtsrcmaps.txt) tool creates model counts maps for use with the binned likelihood analysis. To do this, it takes each source spectrum in the XML model, multiplies it by the exposure at the source position, and convolves that exposure with the effective PSF." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time -p gtsrcmaps scfile= sctable=\"SC_DATA\" expcube=working/3C279_ltCube.fits cmap=working/3C279_binned_ccube.fits srcmdl=working/3C279_model.xml bexpmap=working/3C279_BinnedExpMap.fits wmap=none outfile=working/3C279_binned_srcmaps.fits irfs=\"CALDB\" evtype=\"INDEF\" convol=yes resample=yes rfactor=2 minbinsz=0.1 ptsrc=yes psfcorr=yes emapbnds=yes edisp_bins=0 copyall=no chatter=2 clobber=yes debug=no gui=no mode=\"ql\"\n", "Generating SourceMap for 4FGL J1118.2-0415 38....................!\n", "Generating SourceMap for 4FGL J1118.6-1235 38....................!\n", "Generating SourceMap for 4FGL J1119.9-1007 38....................!\n", "Generating SourceMap for 4FGL J1121.3-0011 38....................!\n", "Generating SourceMap for 4FGL J1121.4-0553 38....................!\n", "Generating SourceMap for 4FGL J1122.0-0231 38....................!\n", "Generating SourceMap for 4FGL J1122.5-1440 38....................!\n", "Generating SourceMap for 4FGL J1124.1-1203 38....................!\n", "Generating SourceMap for 4FGL J1124.6-0809 38....................!\n", "Generating SourceMap for 4FGL J1125.9-0742 38....................!\n", "Generating SourceMap for 4FGL J1129.2-0529 38....................!\n", "Generating SourceMap for 4FGL J1129.2-1014 38....................!\n", "Generating SourceMap for 4FGL J1129.8-1447 38....................!\n", "Generating SourceMap for 4FGL J1131.1-0944 38....................!\n", "Generating SourceMap for 4FGL J1131.4-0504 38....................!\n", "Generating SourceMap for 4FGL J1132.1-1448 38....................!\n", "Generating SourceMap for 4FGL J1132.7+0034 38....................!\n", "Generating SourceMap for 4FGL J1133.8-2048 38....................!\n", "Generating SourceMap for 4FGL J1134.8-1729 38....................!\n", "Generating SourceMap for 4FGL J1135.7-0427 38....................!\n", "Generating SourceMap for 4FGL J1136.3-0501 38....................!\n", "Generating SourceMap for 4FGL J1137.9-1708 38....................!\n", "Generating SourceMap for 4FGL J1141.5-1408 38....................!\n", "Generating SourceMap for 4FGL J1142.8+0120 38....................!\n", "Generating SourceMap for 4FGL J1145.5-0340 38....................!\n", "Generating SourceMap for 4FGL J1145.7+0453 38....................!\n", "Generating SourceMap for 4FGL J1146.0-0638 38....................!\n", "Generating SourceMap for 4FGL J1147.8-0724 38....................!\n", "Generating SourceMap for 4FGL J1148.2-0029 38....................!\n", "Generating SourceMap for 4FGL J1150.8-1905 38....................!\n", "Generating SourceMap for 4FGL J1151.3+0957 38....................!\n", "Generating SourceMap for 4FGL J1151.5-1347 38....................!\n", "Generating SourceMap for 4FGL J1151.6-2115 38....................!\n", "Generating SourceMap for 4FGL J1152.3-0839 38....................!\n", "Generating SourceMap for 4FGL J1153.3-1104 38....................!\n", "Generating SourceMap for 4FGL J1153.6-2553 38....................!\n", "Generating SourceMap for 4FGL J1154.0-0010 38....................!\n", "Generating SourceMap for 4FGL J1155.2-1111 38....................!\n", "Generating SourceMap for 4FGL J1156.6+0640 38....................!\n", "Generating SourceMap for 4FGL J1156.6-2248 38....................!\n", "Generating SourceMap for 4FGL J1158.8-1430 38....................!\n", "Generating SourceMap for 4FGL J1158.9+0818 38....................!\n", "Generating SourceMap for 4FGL J1159.0+0939 38....................!\n", "Generating SourceMap for 4FGL J1159.2-2227 38....................!\n", "Generating SourceMap for 4FGL J1159.3-2142 38....................!\n", "Generating SourceMap for 4FGL J1159.5-0723 38....................!\n", "Generating SourceMap for 4FGL J1200.2+0201 38....................!\n", "Generating SourceMap for 4FGL J1200.6+1229 38....................!\n", "Generating SourceMap for 4FGL J1200.8-1429 38....................!\n", "Generating SourceMap for 4FGL J1201.1-0332 38....................!\n", "Generating SourceMap for 4FGL J1201.7+1429 38....................!\n", "Generating SourceMap for 4FGL J1202.5-0528 38....................!\n", "Generating SourceMap for 4FGL J1203.3+1119 38....................!\n", "Generating SourceMap for 4FGL J1203.5-1748 38....................!\n", "Generating SourceMap for 4FGL J1204.0+1146 38....................!\n", "Generating SourceMap for 4FGL J1204.2-0709 38....................!\n", "Generating SourceMap for 4FGL J1204.8+0407 38....................!\n", "Generating SourceMap for 4FGL J1205.7-2635 38....................!\n", "Generating SourceMap for 4FGL J1207.2-0524 38....................!\n", "Generating SourceMap for 4FGL J1207.4-1840 38....................!\n", "Generating SourceMap for 4FGL J1207.7-0106 38....................!\n", "Generating SourceMap for 4FGL J1207.7-2229 38....................!\n", "Generating SourceMap for 4FGL J1208.2+1158 38....................!\n", "Generating SourceMap for 4FGL J1208.8-1751 38....................!\n", "Generating SourceMap for 4FGL J1208.8-2754 38....................!\n", "Generating SourceMap for 4FGL J1211.6-2735 38....................!\n", "Generating SourceMap for 4FGL J1212.0-2326 38....................!\n", "Generating SourceMap for 4FGL J1212.7-1402 38....................!\n", "Generating SourceMap for 4FGL J1213.3-2618 38....................!\n", "Generating SourceMap for 4FGL J1213.6+1306 38....................!\n", "Generating SourceMap for 4FGL J1214.5-2318 38....................!\n", "Generating SourceMap for 4FGL J1214.6-1926 38....................!\n", "Generating SourceMap for 4FGL J1215.0+1656 38....................!\n", "Generating SourceMap for 4FGL J1215.1+0731 38....................!\n", "Generating SourceMap for 4FGL J1215.8-1733 38....................!\n", "Generating SourceMap for 4FGL J1216.1+0930 38....................!\n", "Generating SourceMap for 4FGL J1216.1-0242 38....................!\n", "Generating SourceMap for 4FGL J1216.2+0537 38....................!\n", "Generating SourceMap for 4FGL J1216.2-0600 38....................!\n", "Generating SourceMap for 4FGL J1216.2-1550 38....................!\n", "Generating SourceMap for 4FGL J1216.4-1358 38....................!\n", "Generating SourceMap for 4FGL J1217.2-2500 38....................!\n", "Generating SourceMap for 4FGL J1218.0-0028 38....................!\n", "Generating SourceMap for 4FGL J1218.5-0119 38....................!\n", "Generating SourceMap for 4FGL J1218.6-1049 38....................!\n", "Generating SourceMap for 4FGL J1219.6+0550 38....................!\n", "Generating SourceMap for 4FGL J1219.7+0444 38....................!\n", "Generating SourceMap for 4FGL J1219.7-0313 38....................!\n", "Generating SourceMap for 4FGL J1220.1-2458 38....................!\n", "Generating SourceMap for 4FGL J1220.5-2745 38....................!\n", "Generating SourceMap for 4FGL J1221.4-0634 38....................!\n", "Generating SourceMap for 4FGL J1222.5+0414 38....................!\n", "Generating SourceMap for 4FGL J1223.0+1100 38....................!\n", "Generating SourceMap for 4FGL J1223.3+1213 38....................!\n", "Generating SourceMap for 4FGL J1223.5+0818 38....................!\n", "Generating SourceMap for 4FGL J1225.0+0330 38....................!\n", "Generating SourceMap for 4FGL J1225.4-1550 38....................!\n", "Generating SourceMap for 4FGL J1225.5-2851 38....................!\n", "Generating SourceMap for 4FGL J1226.7+0637 38....................!\n", "Generating SourceMap for 4FGL J1226.8-1329 38....................!\n", "Generating SourceMap for 4FGL J1228.7-0318 38....................!\n", "Generating SourceMap for 4FGL J1229.0+0202 38....................!\n", "Generating SourceMap for 4FGL J1230.8+1223 38....................!\n", "Generating SourceMap for 4FGL J1231.1-1412 38....................!\n", "Generating SourceMap for 4FGL J1231.5+1421 38....................!\n", "Generating SourceMap for 4FGL J1233.0+1333 38....................!\n", "Generating SourceMap for 4FGL J1233.1+1703 38....................!\n", "Generating SourceMap for 4FGL J1233.7-0144 38....................!\n", "Generating SourceMap for 4FGL J1234.7-0434 38....................!\n", "Generating SourceMap for 4FGL J1238.3-1959 38....................!\n", "Generating SourceMap for 4FGL J1238.5-1201 38....................!\n", "Generating SourceMap for 4FGL J1239.4+0728 38....................!\n", "Generating SourceMap for 4FGL J1239.5+0443 38....................!\n", "Generating SourceMap for 4FGL J1240.4-2606 38....................!\n", "Generating SourceMap for 4FGL J1241.8-1456 38....................!\n", "Generating SourceMap for 4FGL J1241.9+0636 38....................!\n", "Generating SourceMap for 4FGL J1242.4-2948 38....................!\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Generating SourceMap for 4FGL J1243.7+1727 38....................!\n", "Generating SourceMap for 4FGL J1243.9-0218 38....................!\n", "Generating SourceMap for 4FGL J1244.5+1616 38....................!\n", "Generating SourceMap for 4FGL J1245.8+0232 38....................!\n", "Generating SourceMap for 4FGL J1246.3+0112 38....................!\n", "Generating SourceMap for 4FGL J1246.7-2548 38....................!\n", "Generating SourceMap for 4FGL J1249.2-2809 38....................!\n", "Generating SourceMap for 4FGL J1249.3-0545 38....................!\n", "Generating SourceMap for 4FGL J1250.6+0217 38....................!\n", "Generating SourceMap for 4FGL J1251.2+1039 38....................!\n", "Generating SourceMap for 4FGL J1251.3-0201 38....................!\n", "Generating SourceMap for 4FGL J1251.3-1719 38....................!\n", "Generating SourceMap for 4FGL J1252.7-2520 38....................!\n", "Generating SourceMap for 4FGL J1253.8+0327 38....................!\n", "Generating SourceMap for 4FGL J1254.2-2205 38....................!\n", "Generating SourceMap for 4FGL J1254.9+1138 38....................!\n", "Generating SourceMap for 4FGL J1256.1-0547 38....................!\n", "Generating SourceMap for 4FGL J1256.2-1146 38....................!\n", "Generating SourceMap for 4FGL J1256.8-1333 38....................!\n", "Generating SourceMap for 4FGL J1258.6-1759 38....................!\n", "Generating SourceMap for 4FGL J1258.7-0452 38....................!\n", "Generating SourceMap for 4FGL J1258.8-2219 38....................!\n", "Generating SourceMap for 4FGL J1259.1-2311 38....................!\n", "Generating SourceMap for 4FGL J1300.0+1753 38....................!\n", "Generating SourceMap for 4FGL J1300.4+1416 38....................!\n", "Generating SourceMap for 4FGL J1301.6+0834 38....................!\n", "Generating SourceMap for 4FGL J1304.2-2412 38....................!\n", "Generating SourceMap for 4FGL J1304.4+1203 38....................!\n", "Generating SourceMap for 4FGL J1304.6+0539 38....................!\n", "Generating SourceMap for 4FGL J1304.6-0348 38....................!\n", "Generating SourceMap for 4FGL J1304.9-2107 38....................!\n", "Generating SourceMap for 4FGL J1306.3+1113 38....................!\n", "Generating SourceMap for 4FGL J1306.7-2148 38....................!\n", "Generating SourceMap for 4FGL J1308.7+0347 38....................!\n", "Generating SourceMap for 4FGL J1309.7+1153 38....................!\n", "Generating SourceMap for 4FGL J1310.2-1158 38....................!\n", "Generating SourceMap for 4FGL J1311.0+0034 38....................!\n", "Generating SourceMap for 4FGL J1312.4-2156 38....................!\n", "Generating SourceMap for 4FGL J1312.6-1900 38....................!\n", "Generating SourceMap for 4FGL J1312.7+0050 38....................!\n", "Generating SourceMap for 4FGL J1312.8-0425 38....................!\n", "Generating SourceMap for 4FGL J1312.8-2350 38....................!\n", "Generating SourceMap for 4FGL J1315.5+1135 38....................!\n", "Generating SourceMap for 4FGL J1315.9-0732 38....................!\n", "Generating SourceMap for 4FGL J1317.5-0153 38....................!\n", "Generating SourceMap for 4FGL J1318.1-1740 38....................!\n", "Generating SourceMap for 4FGL J1318.3-1148 38....................!\n", "Generating SourceMap for 4FGL J1318.7-1234 38....................!\n", "Generating SourceMap for 4FGL J1319.5+1404 38....................!\n", "Generating SourceMap for 4FGL J1319.5-0045 38....................!\n", "Generating SourceMap for 4FGL J1321.3-2641 38....................!\n", "Generating SourceMap for 4FGL J1322.2+0842 38....................!\n", "Generating SourceMap for 4FGL J1322.3-0606 38....................!\n", "Generating SourceMap for 4FGL J1322.6-0936 38....................!\n", "Generating SourceMap for 4FGL J1322.6-1418 38....................!\n", "Generating SourceMap for 4FGL J1322.6-1617 38....................!\n", "Generating SourceMap for 4FGL J1322.9+0437 38....................!\n", "Generating SourceMap for 4FGL J1323.9+1405 38....................!\n", "Generating SourceMap for 4FGL J1325.6-0227 38....................!\n", "Generating SourceMap for 4FGL J1326.1+1232 38....................!\n", "Generating SourceMap for 4FGL J1326.7-0503 38....................!\n", "Generating SourceMap for 4FGL J1328.6+1145 38....................!\n", "Generating SourceMap for 4FGL J1329.4-0530 38....................!\n", "Generating SourceMap for 4FGL J1331.2-1325 38....................!\n", "Generating SourceMap for 4FGL J1331.6+1711 38....................!\n", "Generating SourceMap for 4FGL J1331.7-0343 38....................!\n", "Generating SourceMap for 4FGL J1331.7-0647 38....................!\n", "Generating SourceMap for 4FGL J1332.0-0509 38....................!\n", "Generating SourceMap for 4FGL J1332.6-1256 38....................!\n", "Generating SourceMap for 4FGL J1337.6-1257 38....................!\n", "Generating SourceMap for 4FGL J1337.9-1956 38....................!\n", "Generating SourceMap for 4FGL J1338.9+1153 38....................!\n", "Generating SourceMap for 4FGL J1339.0-2400 38....................!\n", "Generating SourceMap for 4FGL J1339.1-2620 38....................!\n", "Generating SourceMap for 4FGL J1339.9-0138 38....................!\n", "Generating SourceMap for 4FGL J1340.0-1501 38....................!\n", "Generating SourceMap for 4FGL J1340.8-0409 38....................!\n", "Generating SourceMap for 4FGL J1341.8-2053 38....................!\n", "Generating SourceMap for 4FGL J1342.6+0944 38....................!\n", "Generating SourceMap for 4FGL J1342.7+0505 38....................!\n", "Generating SourceMap for 4FGL J1344.2-1723 38....................!\n", "Generating SourceMap for 4FGL J1345.8+0706 38....................!\n", "Generating SourceMap for 4FGL J1345.9-2612 38....................!\n", "Generating SourceMap for 4FGL J1348.9+0756 38....................!\n", "Generating SourceMap for 4FGL J1349.5-1131 38....................!\n", "Generating SourceMap for 4FGL J1351.0+0029 38....................!\n", "Generating SourceMap for 4FGL J1351.3+1115 38....................!\n", "Generating SourceMap for 4FGL J1351.4-1529 38....................!\n", "Generating SourceMap for 4FGL J1353.3+1434 38....................!\n", "Generating SourceMap for 4FGL J1354.2+0353 38....................!\n", "Generating SourceMap for 4FGL J1354.3-0206 38....................!\n", "Generating SourceMap for 4FGL J1354.7+0623 38....................!\n", "Generating SourceMap for 4FGL J1354.8-1041 38....................!\n", "Generating SourceMap for 4FGL J1356.2-1726 38....................!\n", "Generating SourceMap for 4FGL J1356.6+0234 38....................!\n", "Generating SourceMap for 4FGL J1357.5+0127 38....................!\n", "Generating SourceMap for 4FGL J1358.9-0703 38....................!\n", "Generating SourceMap for 4FGL J1359.1-1152 38....................!\n", "Generating SourceMap for 4FGL J1359.4+0202 38....................!\n", "Generating SourceMap for 4FGL J1400.0-2415 38....................!\n", "Generating SourceMap for 4FGL J1400.6-1432 38....................!\n", "Generating SourceMap for 4FGL J1401.2-0915 38....................!\n", "Generating SourceMap for 4FGL J1402.5-1827 38....................!\n", "Generating SourceMap for 4FGL J1404.8+0402 38....................!\n", "Generating SourceMap for 4FGL J1405.3-0640 38....................!\n", "Generating SourceMap for 4FGL J1405.9-1853 38....................!\n", "Generating SourceMap for 4FGL J1406.4-1654 38....................!\n", "Generating SourceMap for 4FGL J1407.4-0820 38....................!\n", "Generating SourceMap for 4FGL J1408.9-0751 38....................!\n", "Generating SourceMap for 4FGL J1410.1+0202 38....................!\n", "Generating SourceMap for 4FGL J1411.5-0723 38....................!\n", "Generating SourceMap for 4FGL J1415.9-1002 38....................!\n", "Generating SourceMap for 4FGL J1415.9-1504 38....................!\n", "Generating SourceMap for 4FGL J1418.4-0233 38....................!\n", "Generating SourceMap for 4FGL J1419.3+0444 38....................!\n", "Generating SourceMap for 4FGL J1419.4-0838 38....................!\n", "Generating SourceMap for 4FGL J1420.3+0612 38....................!\n", "Generating SourceMap for 4FGL J1421.1-1120 38....................!\n", "Generating SourceMap for 4FGL J1421.4-1655 38....................!\n", "Generating SourceMap for 4FGL J1424.1-1750 38....................!\n", "Generating SourceMap for 4FGL J1424.2+0433 38....................!\n", "Generating SourceMap for 4FGL J1425.4-0119 38....................!\n", "Generating SourceMap for 4FGL J1428.7-1017 38....................!\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Generating SourceMap for 4FGL J1429.8-0739 38....................!\n", "Generating SourceMap for working/gll_iem_v07.fits 38....................!\n", "Generating SourceMap for working/iso_P8R3_SOURCE_V3_v1 38....................!\n", "real 501.95\n", "user 482.46\n", "sys 17.98\n" ] } ], "source": [ "my_apps.srcMaps['expcube'] = 'working/3C279_ltCube.fits'\n", "my_apps.srcMaps['cmap'] = 'working/3C279_binned_ccube.fits'\n", "my_apps.srcMaps['srcmdl'] = 'working/3C279_model.xml'\n", "my_apps.srcMaps['outfile'] = 'working/3C279_binned_srcmaps.fits'\n", "my_apps.srcMaps['bexpmap'] = 'working/3C279_BinnedExpMap.fits'\n", "my_apps.srcMaps['irfs'] = 'CALDB'\n", "my_apps.srcMaps.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute the diffuse source responses.\n", "\n", "The diffuse source responses tell the likelihood fitter what the expected contribution would be for each diffuse source, given the livetime associated with each event. The source model XML file must contain all of the diffuse sources to be fit. The [gtdiffrsp](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtdiffrsp.txt) tool will add one column to the event data file for each diffuse source. The diffuse response depends on the instrument response function (IRF), which must be in agreement with the selection of events, i.e. the event class and event type we are using in our analysis. Since we are using SOURCE class, CALDB should use the P8R3_SOURCE_V3 IRF for this tool.\n", "\n", "If the diffuse responses are not precomputed using [gtdiffrsp](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtdiffrsp.txt), then the [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) tool will compute them at runtime (during the next step). However, as this step is very computationally intensive (often taking ~hours to complete), and it is very likely you will need to run [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) more than once, it is probably wise to precompute these quantities.\n" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "time -p gtdiffrsp evfile=working/3C279_binned_gti.fits evtable=\"EVENTS\" scfile=working/3C279_SC.fits sctable=\"SC_DATA\" srcmdl=working/3C279_model.xml irfs=\"CALDB\" evclsmin=0 evclass=\"INDEF\" evtype=\"INDEF\" convert=no chatter=2 clobber=no debug=no gui=no mode=\"ql\"\n", "adding source working/gll_iem_v07.fits\n", "adding source working/iso_P8R3_SOURCE_V3_v1\n", "Working on...\n", "working/3C279_binned_gti.fits.....................!\n", "real 8569.05\n", "user 8553.75\n", "sys 6.79\n" ] } ], "source": [ "my_apps.diffResps['evfile'] = 'working/3C279_binned_gti.fits'\n", "my_apps.diffResps['scfile'] = 'working/3C279_SC.fits'\n", "my_apps.diffResps['srcmdl'] = 'working/3C279_model.xml'\n", "my_apps.diffResps['irfs'] = 'CALDB'\n", "my_apps.diffResps.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run the Likelihood Analysis\n", "\n", "It's time to actually run the likelihood analysis now. First, you need to import the pyLikelihood module and then the UnbinnedAnalysis functions (there's also a binned analysis module that you can import to do [binned likelihood analysis](https://fermi.gsfc.nasa.gov//ssc/data/analysis/scitools/binned_likelihood_tutorial.html) which behaves almost exactly the same as the unbinned analysis module). For more details on the pyLikelihood module, check out the [pyLikelihood Usage Notes](https://fermi.gsfc.nasa.gov//ssc/data/analysis/scitools/python_usage_notes.html)." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3C279_BinnedExpMap.fits L181126210218F4F0ED2738_PH01.fits\r\n", "3C279_SC.fits L181126210218F4F0ED2738_PH02.fits\r\n", "3C279_binned_ccube.fits L181126210218F4F0ED2738_PH03.fits\r\n", "3C279_binned_filtered.fits L181126210218F4F0ED2738_PH04.fits\r\n", "3C279_binned_gti.fits L181126210218F4F0ED2738_PH05.fits\r\n", "3C279_binned_srcmaps.fits L181126210218F4F0ED2738_PH06.fits\r\n", "3C279_ltCube.fits \u001b[1m\u001b[35mgll_iem_v07.fits\u001b[m\u001b[m\r\n", "3C279_model.xml gll_psc_v26.xml\r\n", "L181126210218F4F0ED2738_PH00.fits \u001b[1m\u001b[35miso_P8R3_SOURCE_V3_v1.txt\u001b[m\u001b[m\r\n" ] } ], "source": [ "!ls working/" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "import pyLikelihood\n", "from BinnedAnalysis import *\n", "obs = BinnedObs(srcMaps='working/3C279_binned_srcmaps.fits',\n", " binnedExpMap='working/3C279_BinnedExpMap.fits',\n", " expCube='working/3C279_ltCube.fits',\n", " irfs='CALDB')\n", "like = BinnedAnalysis(obs,'working/3C279_model.xml',optimizer='NewMinuit')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By now, you'll have two objects, 'obs', an UnbinnedObs object and like, an UnbinnedAnalysis object. You can view these objects attributes and set them from the command line in various ways. For example:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Source maps: working/3C279_binned_srcmaps.fits\n", "Exposure cube: working/3C279_ltCube.fits\n", "Exposure map: working/3C279_BinnedExpMap.fits\n", "IRFs: CALDB\n" ] } ], "source": [ "print(obs)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Source maps: working/3C279_binned_srcmaps.fits\n", "Exposure cube: working/3C279_ltCube.fits\n", "Exposure map: working/3C279_BinnedExpMap.fits\n", "IRFs: CALDB\n", "Source model file: working/3C279_model.xml\n", "Optimizer: NewMinuit\n" ] } ], "source": [ "print(like)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are a lot of attributes and here you start to see the power of using pyLikelihood since you'll be able (once the fit is done) to access any of these attributes directly within python and use them in your own scripts." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we're ready to do the actual fit. This next step will take approximately 10 minutes to complete. We're doing something a bit fancy here. We're getting the minimizating object (and calling it `likeobj`) from the logLike object so that we can access it later. We pass this object to the fit routine so that it knows which fitting object to use. We're also telling the code to calculate the\n", "covariance matrix so we can get at the errors." ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "72593.89794181124" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "likeobj=pyLikelihood.NewMinuit(like.logLike)\n", "like.fit(verbosity=0,covar=True,optObject=likeobj)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The number that is printed out here is the -log(Likelihood) of the total fit to the data. You can print the results of the fit by accessing like.model. You can also access the fit for a particular source by doing the following (the source name must match that in the XML model file)." ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4FGL J1256.1-0547\n", " Spectrum: LogParabola\n", "422 norm: 1.386e+00 1.518e-02 0.000e+00 1.000e+02 ( 1.000e-10)\n", "423 alpha: 2.293e+00 1.028e-02 0.000e+00 5.000e+00 ( 1.000e+00)\n", "424 beta: 6.866e-01 6.540e-02 -5.000e+00 1.000e+01 ( 1.000e-01)\n", "425 Eb: 5.537e-01 0.000e+00 1.000e-02 1.000e+02 ( 1.000e+03) fixed" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "like.model['4FGL J1256.1-0547']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Test Statistic (TS) value is a useful quantity to record when evaluating likelihood results. A TS > 25 is generally considered a strong detection with a high likelihood of being real. " ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "29315.27972251101" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "like.Ts('4FGL J1256.1-0547')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can plot the results of the fit by executing the plot command. The results are shown below:\n" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Data 207751.0\n", "srcName 1.1504089856783581\n", "srcName 0.8788688981057069\n", "srcName 0.15820582237595157\n", "srcName 0.08109743641834478\n", "srcName 2.549122407017949\n", "srcName 0.2767968977519388\n", "srcName 0.4958909457972051\n", "srcName 0.20051919011898275\n", "srcName 0.17708226027712862\n", "srcName 0.3426429412482797\n", "srcName 2.754296571959106\n", "srcName 0.1824172436014473\n", "srcName 2.7826608455110713\n", "srcName 0.2888488572283987\n", "srcName 0.2764626981104433\n", "srcName 0.2542905757678804\n", "srcName 4.108657368700294\n", "srcName 0.07057590484918022\n", "srcName 0.03347189061549695\n", "srcName 2.2506870238087617\n", "srcName 1.7344154603658009\n", "srcName 0.003753050331835279\n", "srcName 0.5387637828022851\n", "srcName 0.700103020448729\n", "srcName 0.2992937680808332\n", "srcName 2.096307949009182\n", "srcName 0.6266209677939509\n", "srcName 9.050661290226904\n", "srcName 1.7915784920471234\n", "srcName 0.045931850792591\n", "srcName 1.6161343742928487\n", "srcName 2.6084321442857457\n", "srcName 0.19703599394953045\n", "srcName 22.535412245163425\n", "srcName 9.899551839764582\n", "srcName 0.10115531181125208\n", "srcName 1.9049553331982365\n", "srcName 1.4515575687182454\n", "srcName 1.0434020893916236\n", "srcName 0.44891801617440386\n", "srcName 1.2886786252852835\n", "srcName 0.46492401076803525\n", "srcName 1.9671520270671683\n", "srcName 4.826082295051852\n", "srcName 6.495554434116405\n", "srcName 518.9485981848627\n", "srcName 7.769305750106749\n", "srcName 0.8891681154632013\n", "srcName 3.211467402760401\n", "srcName 4.544102538491933e-05\n", "srcName 0.49040099918545055\n", "srcName 214.39199876138804\n", "srcName 0.7530134840782118\n", "srcName 5.635015762503352\n", "srcName 0.3210458381938306\n", "srcName 299.2966516440408\n", "srcName 8.820743888884317\n", "srcName 1.7277450527828846\n", "srcName 4.409944077625369\n", "srcName 1.7099383273721922\n", "srcName 227.22404072214985\n", "srcName 0.9658472048305758\n", "srcName 0.3659989407303855\n", "srcName 4.0907764918478815\n", "srcName 0.18174902772512985\n", "srcName 0.7241894413920065\n", "srcName 0.49428443657507826\n", "srcName 6.584345830806923e-05\n", "srcName 0.14388571587248353\n", "srcName 2.0036502715090148\n", "srcName 0.7087979907551332\n", "srcName 0.5463972379938437\n", "srcName 1.0764120803407025\n", "srcName 1.3105132090911287\n", "srcName 31.20002541103885\n", "srcName 2.293044584273701\n", "srcName 29.358131122319907\n", "srcName 4.272520599730501\n", "srcName 136.38896579887006\n", "srcName 86.25002970027987\n", "srcName 57.00987891623088\n", "srcName 0.019652170688910127\n", "srcName 245.07378725118267\n", "srcName 432.38024712765855\n", "srcName 142.76847033175372\n", "srcName 0.0037912570789370268\n", "srcName 165.09366675013024\n", "srcName 250.99736421097174\n", "srcName 0.4583677724166481\n", "srcName 0.7173553522256257\n", "srcName 270.8185544033446\n", "srcName 626.3937982257435\n", "srcName 0.756183420907142\n", "srcName 1.6400452042904314\n", "srcName 2.2524005186805067\n", "srcName 2.323464800713032e-05\n", "srcName 64.62572393679028\n", "srcName 0.37965068008166597\n", "srcName 826.2168612903702\n", "srcName 422.6544023967016\n", "srcName 294.83292203957626\n", "srcName 21345.569748897677\n", "srcName 5.717314658254754\n", "srcName 3717.365424971787\n", "srcName 0.34894190669246755\n", "srcName 2.5987165355157087\n", "srcName 0.09608645394797094\n", "srcName 561.2072081629782\n", "srcName 159.1405415456183\n", "srcName 4.002879167580427e-05\n", "srcName 628.673695728791\n", "srcName 0.003998432812931409\n", "srcName 967.2196831804422\n", "srcName 0.9671795867720775\n", "srcName 84.31688600676772\n", "srcName 0.06680054134944322\n", "srcName 0.7185092499104513\n", "srcName 0.07267075046905339\n", "srcName 84.58199413452941\n", "srcName 0.38941562075782843\n", "srcName 24.68851249775048\n", "srcName 61.7811914998598\n", "srcName 21.49043325276374\n", "srcName 0.5006922529553298\n", "srcName 908.8945627451053\n", "srcName 169.0927689860548\n", "srcName 2.0806372292025066\n", "srcName 2.9114407278603482\n", "srcName 0.012201996036989373\n", "srcName 2.1988467944551564\n", "srcName 163.4544129617077\n", "srcName 4.057097515626703\n", "srcName 9.01696059624693\n", "srcName 20764.834806147544\n", "srcName 475.5750663903168\n", "srcName 539.8278462125837\n", "srcName 36.259461573067675\n", "srcName 24.36827385184856\n", "srcName 26.864027639841566\n", "srcName 4.224329138406296\n", "srcName 0.4386048561841071\n", "srcName 6.0093063372576045\n", "srcName 0.17879104157720582\n", "srcName 1.6215378689226634\n", "srcName 0.08120584997947902\n", "srcName 45.956181585453734\n", "srcName 242.4758131403227\n", "srcName 5.651678151273757\n", "srcName 0.7240951857464616\n", "srcName 6.039491906398843\n", "srcName 79.65753128207245\n", "srcName 3.9158138808203344\n", "srcName 170.27680941847598\n", "srcName 202.08265138283173\n", "srcName 6.155970369444775\n", "srcName 0.4315806026239708\n", "srcName 754.8190315475374\n", "srcName 180.52262806288027\n", "srcName 1.2378912070087067\n", "srcName 0.5633957363372201\n", "srcName 259.3888866099624\n", "srcName 58.09605705199973\n", "srcName 91.19525182457193\n", "srcName 21.266570121862156\n", "srcName 476.5298111172637\n", "srcName 0.35889299472934105\n", "srcName 0.014943146319039315\n", "srcName 0.8310241294340641\n", "srcName 7.522070920525799\n", "srcName 38.33107270211055\n", "srcName 1004.5127216961654\n", "srcName 28.90801219583335\n", "srcName 340.6930790863642\n", "srcName 148.37894699012944\n", "srcName 0.5551241782288616\n", "srcName 29.89347400614889\n", "srcName 0.23756815052330466\n", "srcName 36.87668197053701\n", "srcName 0.12068389231460536\n", "srcName 59.74702552929831\n", "srcName 194.26963198822355\n", "srcName 0.4645042563458583\n", "srcName 273.9132906818952\n", "srcName 22.804565199690348\n", "srcName 4930.160237122009\n", "srcName 2003.6052598410477\n", "srcName 856.3178998255675\n", "srcName 4.349523630677995\n", "srcName 1.126078423841181\n", "srcName 3.182998634506742\n", "srcName 0.492623400755473\n", "srcName 29.521079033619856\n", "srcName 168.25021828208276\n", "srcName 28.365189950384874\n", "srcName 10.598289401592744\n", "srcName 1.6940345051205248\n", "srcName 12.6422907875694\n", "srcName 30.89668730867493\n", "srcName 17.519680318673583\n", "srcName 0.01028996524362737\n", "srcName 0.8289423729671689\n", "srcName 172.93770003615188\n", "srcName 14.233107704530584\n", "srcName 0.8568915040247556\n", "srcName 4.1309564803999645\n", "srcName 0.5725309390039445\n", "srcName 0.8840772770859929\n", "srcName 10.19604938827653\n", "srcName 0.7464054105377751\n", "srcName 19.927341005780388\n", "srcName 0.004607212508551199\n", "srcName 0.29651407728572965\n", "srcName 5.140964358554897\n", "srcName 0.028065053579295066\n", "srcName 0.09374963820605788\n", "srcName 3.9310424341375825\n", "srcName 0.03993497565506368\n", "srcName 2.0141650052449083\n", "srcName 2.2677289998609687\n", "srcName 3.2491116709481958\n", "srcName 1.7682005554888058\n", "srcName 1.471146461470995\n", "srcName 0.2323508365186142\n", "srcName 0.9829053505544236\n", "srcName 3.1452035942344976\n", "srcName 12.208680429270407\n", "srcName 2.794541914949428\n", "srcName 0.929506556188149\n", "srcName 0.15312241267622176\n", "srcName 0.42853856401857643\n", "srcName 1.0109281231913287\n", "srcName 0.14286276456733635\n", "srcName 5.178784324341854\n", "srcName 0.5285790288202092\n", "srcName 1.4779111853794438\n", "srcName 0.015758184349051928\n", "srcName 0.10231209779810721\n", "srcName 0.2690416639875203\n", "srcName 0.060737452903332\n", "srcName 0.01579278643778187\n", "srcName 0.008112506216967531\n", "srcName 80497.0901290422\n", "srcName 57778.47906938563\n", "Resid 207751.0 207750.6617432125 -2.9638954360465752\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/opt/miniconda3/envs/fermi/lib/python3.7/site-packages/fermitools/MPLPlot.py:53: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.\n", " self.fig.show()\n", "/opt/miniconda3/envs/fermi/lib/python3.7/site-packages/fermitools/MPLPlot.py:64: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.\n", " self.fig.show()\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAADCsElEQVR4nOxdZXgTWRc+E2/a1L1UKIVCgWKlFAoUd3f3xd3dXQvFfdFlYWEpFGexxRZ3K67FoUYleb8fZe6XtKmnAsz7PPdpmszcuZkk59xj7+EAkAABAgQIEEBEJMrpBQgQIECAgNwDQSkIECBAgAAGQSkIECBAgAAGQSkIECBAgAAGQSkIECBAgAAGQSkIECBAgAAGSU4vIDOwtraGm5tbTi9DgAABAn4oXLp06T0AG32v/dBKwc3NjS5evJjTyxAgQICAHwocxz1N7rUf0n3EcVx9juNWfvnyJaeXIkCAAAE/FX5IpQBgD4DuZmZmOb0UAQIECPip8EMqBQECBAgQkDUQlIIAAQIECGAQlIIAAQIECGD46ZTCxIkTs+0YAQIECPjZwP2I1Nkcx9UnovoeHh6/PXjwIPFrlPg9xcfHU3R0NEVHR1NsbCw5OzsnOUbPNVI9RoAAAQJ+RHAcdwmAj77Xfsg6BQB7iGiPj4/Pb/xz9+7dozZt2hBRQv0CrwSio6MpPj4+yRxWVlaUJ08eNpycnHT+FyBAgIBfET+kUkiMiRMn0qRJk9j/T58m1GWUKVOGqlSpQkZGRnTq1Ck6fPgwO+bjx4/08eNHevfuHV28eJHevn2bZF6O44iIqHXr1hQUFERWVlZJriu4mQQIEPAz4Yd0H/Hw8fFB4ormtLh99B0TExNDr1+/phcvXtDTp0+pXbt2VLp0abpy5QqzNAoWLEjlypVjw8vLS3AxCRAg4IdDSu4jQSmkckxUVBRdvHiRzpw5w8aHDx/Ycb1796aGDRtSpUqVSCaTZfzNCBAgQEA2ISWl8NNlH02YMMGgxyiVSqpYsSKNHDmSgoODqU+fPjrHLV26lGrWrEkqlYpatmxJmzdvpk+fPrHXBfeSAAECfiT8dJZCdoLjOIqKiqKjR49ScHAwBQcHU1hYGInFYqpYsSI1aNCABg0aJLiYBAgQkKvwS7mPshOJ3VAajYYuXLhAu3fvpuDgYLp16xYREVWsWJHatm1LzZo1I0tLy5xargABAgQQ0U/oPsotLKmJ3VAikYjKlClDMpmMKQQiopMnT1KPHj3IxsaGGjVqRH/++SdFR0ez1wUXkwABAnILBEshi8FxHGk0Grpy5Qpt3ryZtm7dSq9fvyaVSkWNGzemtm3bUs2aNQUXkwABArINgvsoB5HYxaRWq+nEiRO0efNm2rFjB339+pWIiMaNG0ddu3YlV1fXnFqqAAECfhH8dO6jHwmJXUxisZiqVKlCzs7OTCEQEU2ZMoXc3NzIw8OD/vrrL4qLi9M5T3AxCRAgIDsgWAq5ABzH0ZMnT2jt2rW0du1aevHiBdna2lLnzp2pW7du5OHhIXAxCRAgwGAQ3Ee5HNoCX61W04EDB2jVqlW0d+9eUqvVVLlyZTp27BjFxcWRRPJTMJMIECAgByG4j3I5tF1MYrGY6tatS3///TcNGDCAiIiOHTtGRERSqZQ4jqORI0fmyDoFCBDw80NQCrkAycUL5s2bRwAY91KlSpWIiGjx4sXUr18/SkwbLsQdBAgQkFkISuEHgFgsJqIEi+HKlSvUrFkzWrlyJXl6elL9+vXpn3/+IQA6TLECBAgQkBEISuEHAe9iKl68OK1fv56ePn1K48ePp/Pnz1PVqlWpePHiRERJspYECBAgID3IVUqB4zhjjuMucRxXL6fXktuQ2DVkb29PEydOpO7duxMR0fXr14mISCaTEcdxNHbs2OxeogABAn4CZKlS4DhuLcdxbzmOu5no+Vocx93jOC6U4zjtqOkIIvozK9f0s2Hq1KkEgDQaDRER+fn5ERHR2rVrad68eRQREZGTyxMgQMAPhqy2FNYTUS3tJziOExPREiKqTUReRNSa4zgvjuOqEdFtIgoz9CIA0NevX+nJkyd07do1unfvHr169YoiIiJ+mtx/vkvcmTNn6OjRo1SoUCEaOnQoubm50ZQpU+jz589EJASjBQgQkDKyvE6B4zg3ItoLoMj3/8sS0UQANb//P+r7oSZEZEwJiiKaiBoD0OiZrzsRdScicnFxKcW33nz+/DnNnj2bPn36xFpt8o8/ffpEarVa7/pEIhGZmJiQqakpqVQqMjU1JVNTU8qXLx95eXlR4cKFqXDhwmRra8sEb25F4vag586do2nTptHevXtJpVJRnz59aObMmT+NIhQgQEDGkKPFa3qUQjMiqgWg2/f/2xNRGQB9v//fiYjeA9ib2tzaxWt37tyhcuXKkaWlJVlaWpKFhUWSxxYWFmRmZkYxMTEUHh5OX79+1fv306dP9ODBA51mOZaWllS4cGEdReHj40OmpqaGvWFZgGvXrtH06dPpzz8TPHNDhgyhYcOGkZ2dXQ6vTIAAATmBlJRCTpTH6ttuM80EYH2qE3BcfSKq7+HhwZ4rVKiQjhDPLADQmzdv6Pbt23Tr1i26desW3b59m/788092HZFIRCVLlqSAgACqVKkSlS9fnszNzQ22BkNh165dTCEQJdQ/zJs3j/z8/GjXrl1kb2+fg6sTIEBAbkKucR8BmJHeudNKc6HRaOjVq1cUGhpKz549I5VKRQ4ODuTg4ED29vYkl8vTfE1eWdy4cYNOnz5Nx48fp3PnzlFsbCxxHEclSpSggIAACggIoAoVKuS6pjocx9G9e/do2rRptGnTJpLJZNSzZ08aPnw4OTg4EFFSN5QAAQJ+LuQ295GEiO4TUVUieklEF4ioDYBbyU6SDLSVQnx8PD179oxCQ0MpNDSUHj58yB4/evSIvn37luw8FhYWOkrC0dGRfHx8KCAgIE0ulujoaDp//jydOHGCTpw4QWfPnqVv374Rx3Hk6+tL9evXpwYNGlCRIkVyPC6hzbMUGhpK06ZNo40bN5JUKqUePXrQiBEjyNHRUYg7CBDwEyPHlALHcVuJqBIRWVNCVtEEAGs4jqtDRIFEJCaitQCmpXNe3n30G0/1cOHCBfL19WXHGBkZUb58+Shfvnzk4eHB/rq6ulJERAS9fv2aXr9+TW/evEny+NWrVxQbG0tECW6pSpUqUaVKldKsJGJiYui///6jY8eOUUhICP33339EROTq6koNGjSg+vXrU0BAAMlksvS8bYNAnxXw8OFDmjZtGm3YsIEkEgnFxMTQq1evmOUgQICAnwu/BEtqeHg47dixgwl/BweHDO/K4+Pj6fLly3T8+HE6fvw4nTp1iuX7ayuJGjVqpCmG8Pr1awoJCaHg4GA6cuQIRUdHk0qlolq1alH9+vWpXr16ZGFhkaG1GgoTJ07US5MxZMgQmjt3bg6sSIAAAVmFlJQCAfjhBhHVJ6KVHh4eyA7ExcXh/PnzmDVrFmrXrg0TExMQEWQyGRo2bIg//vgDERERaZorMjISwcHB6NatG+zt7UFEkEqlqF+/PjZv3ozw8PAsfjepg4jQqVMniEQiGBkZYejQoXj79q3OMRMmTMiZxQkQICDTIKKLSE6+JvfCjzBKlSpl0BuVVsTGxuL06dMYOHAgHBwcQEQwNjZG69atsXv3bnz79i1N86jVapw7dw5DhgxBnjx5QEQwMjJC8+bN8ddffyEqKiqL34l+JOwVgPv376N9+/YQiUQwNjbGyJEj8f79e51jBAgQ8ONBUApZiPj4eBw7dgzdu3eHpaUliAjm5ubo0qULDh8+jPj4+DTNo1arcerUKfTp0we2trYgIqhUKrRv3x4hISGIjY3N4nfyfyS2Au7cuYPWrVuD4ziYmJhg9OjRglIQIOAHxk+nFLLbfZRWxMTEICQkBO3bt2cuJmdnZ0ycOBHPnj1L8zxxcXE4fPgwunbtCnNzcxARbG1tMWTIENy6dSsL30HK6NWrFyihpkRnCK4kAQJ+LPx0SoEfucFSSA5RUVHYtm0batSoASKCSCRC3bp1sXv3bsTFxaV5npiYGOzevRuNGzeGRCIBEcHPzw8rV67Ely9fsvAdJI9r164xhVCgQAHs2bMHGo0mR9YiQICA9ENQCjmMR48eYcyYMSz+4OjoiLFjx+Lx48fpmicsLAzz5s2Dl5cXiz906NABJ06cyHahTETYu3cvPD09QUSoXr06bty4ka1rECBAQMYgKIVcgri4OPz999+oU6cOOI4Dx3GoWbMmgoODoVar0zyPRqPBuXPn0L17d6hUKhARPDw8MHPmTLx79y4L38H/wbuMYmNjsXDhQlhYWEAkEqFnz54sU0lwKwkQkDvx0ymF3BpTSA+ePn2KCRMmsKyj/PnzY8mSJWlObeURERGB33//HRUrVgQRQaFQoHPnzrh8+XIWrVw/Pnz4gP79+0MsFsPU1BRz5swRgtECBORS/HRKgR8/mqWgD3Fxcfjjjz/g6+sLIoKFhQVGjhyJFy9epHuuGzduoEePHlAqlSAilCtXDlu3bkVMTEwWrFw/bt++jTp16rCYw9GjR7Pt2gIECEgbBKXwA0Cj0eD06dNo1qwZRCIRJBIJ2rZti4sXL6Z7rk+fPmH+/PnIly8fiAgODg6YOHEiXr9+nQUr18WECRP0ZigNGzYsy68tQICAtEFQCj8YHj16hIEDB7J4QcWKFRESEpLuYLJarUZISAhq167NKqfbtm2LK1euZM3CE4GIMHr0aEgkElhZWeH3338XspQECMgF+OmUws8QU0gLPn/+jHnz5sHFxQVEhDJlyuDgwYMZEqz379/HgAEDWP1EtWrVMjxXWsHHFG7cuIFy5cqBiFC1alXcv38/y64pQICA1PHTKQV+/KyWQmLExsZi5cqVTDn4+/vj6NGjGRLonz59wsyZM1l6rLe3NzZs2JAlcQft7CO1Wo1ly5bBzMwMcrkcU6dORUxMjJChJEBADkBQCj8Jvn37hqVLl8LJyQlEhICAAJw4cSJDc8XExGDdunUoXLgwiAhOTk6YM2dOlhfEvXr1Cs2bNwcRsWsLECAge5GSUvhpqLOfPn1Ko0ePJrVaTfHx8aRWq9nQ/h8AOTo6Ur58+cjd3Z31XHBwcCCRSJTD7yht+PbtG61atYqmT59Ob968oapVq9LkyZOpXLly6Z4LAB04cIDmzJlDx44dI1NTU+rRowcNGjQoS/sphISEUNOmTSkmJob69etH06dPJxMTkyy7ngABAv6PX6Kfwt27d6levXokkUhILBazof2/RCIhAPTixQt69uwZqdVqNpdCoaC8efOSu7s7eXh4ULly5ahq1apkZWWVU28vVURHR9Py5ctp5syZ9PbtW6pbty7NmTOHChUqlKH5Ll26RHPnzqU///yTpFIpdenShYYNG0Z58+Y16LqT693Qvn172rBhg0GvJUCAgKT4JZQCD41GQ7GxsRQXF0exsbFJHqvVarK1tSULCwt6/vw5PXz4kB4+fEiPHj1ij0NDQykqKoo4jqNSpUpRjRo1qEaNGlS2bNkc6ZaWGiIjI2nx4sU0ffp0ioyMpF69etGECRPI2to6Q/OFhobSnDlzaP369aRWq6l169Y0cuRIKly4sIFXntAe9NSpU9S1a1e6f/8+devWjebOnUtmZmYGv5YAAQIS8Es02bl69SpEIpHeHHl9Qy6XI3/+/KhatSo6d+6MCRMmYM2aNThy5Ahu376Nf//9F5MmTYK/vz/EYjHrmVC3bl0sXLgQd+7cyXXplW/fvkWvXr0gEolgbm6OefPmZSqA/OLFCwwePJgVwzVq1Aj//fefAVcMFlOIiorC8OHDIRKJ4OjoiODgYHaMEIwWIMCwoF8hpvDmzRtavHgxyWQykkqlJJPJ9D4WiUQUFhZGz549o6dPn9KzZ8/o2bNn9Pr1a9K+F2ZmZlSuXDkqX748lShRgiIiIuj48eN06NAhCg0NJSKiIkWKUKdOnahdu3Zp6t2cXbh16xYNGTKEDh48SB4eHjR79mxq1KhRhtuTvn//noKCgmjRokX0+fNnqlatGo0dO5YCAgIyvdbEPaMvXLhAXbp0oZs3b1KbNm1o4cKFZGNjQz/y91SAgNyGn85S4Ichs49iYmLw6NEjHD9+HGvXrkX37t0ZGyl9b73p7++PESNGYPXq1ZgzZw7KlCkDIoJEIkGDBg2wa9eubKWUSA379+9n76FSpUqZ5kP6+vUrZs+ezdqI1q5dG9euXTPQav+PmJgYTJw4ERKJBDY2NkKGkgABBgYJKakZx7t377B7924MHToUfn5+rKcBx3EoW7Yshg4dit9++40JShsbGwwcODBLhGVGEBcXhyVLlsDa2hocx6FTp054+fJlpuaMiorC7NmzYW5uDo7j0KFDBzx58sRAK05AcnQZgitJgIDMQ1AKBkRkZCSOHz+OiRMnomTJkkxYFS1aFK1bt0a1atWY4ihZsiRWrVqF6OjobF9nYnz69AlDhw6FVCqFUqnEpEmT0s3ImhgfPnzAsGHDIJfLIZPJMHjwYNbD2VCIi4tjFB329vbYs2ePQecXIOBXhKAUshCPHz/GggULUKFCBXAcByKCm5sbqlSpwgjpbG1tMWXKFIMLzIwgNDQUzZo1YwVrv//+e7p6OejDs2fP0LlzZ3AcBzMzM8yYMQORkZEGWjFARLhy5QqKFi0KIkKXLl1yrOucAAE/AwSlkE148+YNVq5cidq1a0MqlYKIYGdnx5SDkZERevfujQcPHuT0UnHq1Cn4+Pgwi+b48eOZnvPGjRuoX78+6y63Zs2aTCsc4P/ZR9++fcOoUaMgEong4uIi0HILEJBB/HRKgX4AQrzPnz9j06ZNqFmzJkuVtbW1ZemtjRs3xunTp3N0jWq1Gps2bWKNfho3bmwQhXXixAn4+fmxng6Gjq+cPXsW+fPnBxGhX79+iIyMFGINAgSkAz+dUuBHbrMUksOLFy8wc+ZMFCpUiGUr8ZZEmTJlEBwcnKM1D5GRkZg6dSqMjY0hlUoxaNAgfPjwIVNzajQarFu3DtbW1hCLxRg8eDC+fv1qoBUnrLl///6sa52QoSRAQNohKIVcAo1Gg/Pnz6N3796wsLAAETHLoXjx4jh8+HCOKodXr16ha9eu4DgOFhYWWLBgQaZTbD98+IDu3buzGMb27dsN+h6PHj0KuVzOMpPi4uIMNrcAAT8rBKWQC/Ht2zds374dtWvXZgFq+k5lferUqRxd27Vr11C9enUQEfLly4cdO3ZkWpCfPXsWxYoVAxGhVq1aCA0NzfQ6k0tb7d+/f6bnFiDgZ4agFHI5Hj16hMGDB8PY2JgJNi8vL4MEfzMKjUaD/fv3M3rr8uXL4/z585maMy4uDoGBgVCpVJDL5Zg0aRK+fftmkPUSEbZs2QIzMzOoVCps2LAh19GQCBCQWyAohR8E0dHRWLFiBQv8EhHy5s2Lv//+O8fWFBcXhxUrVsDOzg5EhFatWuHx48eZmvPly5do2bIliAgFChTAv//+m+l18jGFJ0+eoHz58iAitG7dGp8+fcr03AIE/GwQlMIPiH/++Ye5W4gIDg4O2LZtW47tfr9+/YqxY8fCyMgIcrkco0aNQlRUVKbmPHjwIFxdXcFxHAYOHJip2gbt7KP4+HhMmTIFYrEYrq6uzB0nZCgJEJAAQSn8wLh37x4qVKjAlIOTkxM2btyYYwHV58+fo3379iAieHh44NixY5ma7+vXr+jduzeLX2S0k5w+nDt3Du7u7hCJRBg7dqyQoSRAwHcISuEnwPXr15l/n7cclixZkundekZx5MgRuLu7g4jQrVu3TLtpjh07hrx587Lag8xScPD4+vUrOnXqxO7bo0ePDDKvAAE/MgSl8JNAo9Fg586dsLW1ZULOysoKkydPznRdQUYQGRmJYcOGQSQSwcHBATt37szUfBEREaz2IG/evPjnn38yvUaBWE+AgKT4IZQCERUiouVEtIOIeqXlnF9NKfCIiorCxIkTIZPJWLW0sbExBg8ejDdv3mT7ei5evMjiH02aNMGrV68yNd/Jkyfh4eEBIkKvXr0MVvRGRChdujSICH369MkVRIUCBOQEckwpENFaInpLRDcTPV+LiO4RUSgRjUz0moiI1qRl/l9VKfB4/PgxmjRpAiKCiYkJRCIRlEolhg0bhrdv32brWmJjYzFjxgzI5XKYmZlh1apVmQqKR0ZGYtCgQeA4Dq6urgaJNRARYmJiMHjwYFYweP/+/UzPK0DAj4acVAoViaiktlIgIjERPSQidyKSEdE1IvL6/loDIjpDRG3SMv+vrhR4HDp0CAULFmQMrRzHwdjYGCNGjMC7d++ydS33799HQEAAiAiVK1fOdJ+Ff//9F/ny5QPHcRg+fHim6hq0XUZ79uyBpaUlTExMsHnzZr3HCBDwsyJH3UdE5JZIKZQlooNa/48iolGJzglJy9yCUvg/vn37hnHjxkEikcDCwgL+/v7gOA4mJiYYPXp0tsYc1Go1VqxYARMTE6hUKqxbty5TVkN4eDh+++03EBGKFSuGmzdvGmSdz549YzUNXbt2RWRkpJChJOCXQG5TCs2IaLXW/+2JaDERVSKiRUS0goj6pDBfdyK6SEQXXVxcsu6u/aC4du0ao8SuXLkyo7JWqVQYO3YsPn78mG1refToEUunbdSoEcLCwjI13+7du2FjYwO5XI4FCxYYhJY7Li4OY8aMAcdxLLtLgICfHblNKTTXoxSCMjK3YCnoR1xcHObMmQOFQgEzMzNMnDiRNdYxNzc3eBOclBAfH4+5c+dCJpPBxsYm09XZYWFhTNFVrVoVz58/z/QahQwlAb8acptSSNV9lIY5c30/hdyABw8eMP9+lSpVsGfPHtSrVw9EBHt7eyxZsiTTLKhpxY0bN1CiRAkQETp16oTPnz9neC6NRoNVq1bB2NgY5ubm2LJli0HW+Pr1a6YQOnbsaLBaCQECchtym1KQENEjIsqrFWgunJG5BUshdfD+fZVKBSMjIyxYsAAnT55kvnR3d3ds2rTJIK6Y1BATE4OxY8eyzmmZrUN48OABa+ZjKJ4j3kLgOA5eXl64detWpucUICC3ISezj7YS0WsiiiOiF0TU9fvzdYjo/vcspDEZmFewFNKJ58+fo27dusxqePr0KUJCQlh9gbe3N/bs2ZMt3ErandP69++P8PDwDM8VFxenw3OUWXI93mV05MgR2NraQqlUYv369ZmaU4CA3IYctRSycgiWQvqg7XYxMzPDpk2bEB8fjy1btrA+0v7+/jh58mSWryUiIgJ9+/YFEcHZ2RnBwcGZmk+b58hQzXZevXqFypUrg4jQuXNnFocRYg0CfnQISkGADkJDQ1GuXDkQEVq0aIEPHz4gNjYWy5Ytg4ODA4gI9evXzxbXyenTp1nWT5MmTfDixYsMz/Xlyxd06NAB9L03dGYpvoGEQPn48eNZdtKtW7eEDCUBPzxSUgpcwuvJg+M4GyL67XtsQMI/D6BLiidmITiOq09E9T08PH578OCBzmsAKCIigsLCwujNmzf05s0bncfv3r0jqVRKxsbGSYaJiQn76+HhQZ6eniSTyXLmTWYx1Go1zZ49m8aPH0+2tra0du1aqlmzJkVFRVFgYCDNmjWLIiIiqHPnzjRp0iRycnLKsrXExsbS/PnzadKkSSSVSmnGjBnUs2dPEovFGZpvy5Yt1KtXLyIiWrFiBbVq1SrTazxy5Ai1bduWIiIiKCoqilL73QgQkJvBcdwlAD56X0xOW/CDEiqMZxFRCyJqyo/UzsuOoW0p3LhxA+7u7lAqlXrTC0UiEezt7VG0aFEUKlQILi4usLKygkKh0Hs8EUEikaBIkSJo3bo1pk+fjj179uDJkyc/VUevy5cvw8vLC0SEvn37MhfJu3fvMHDgQEilUhgZGWHUqFGZyhhKC0JDQ1kbUF9fX1y9ejXDcz169Ahly5ZlmUSZ5U8S0lYF/EygzLiPiOhqasfk1NBWCi9fvkTbtm0xZMgQzJ49Gxs2bMChQ4dw7do1hIWFIT4+PtkbFB8fj69fv+LVq1d48OABLl26hC1btmDUqFGoV68eXF1ddQSBqakpypUrh+HDh+Pw4cM5Rl9tKERHR2PQoEEgInh6euLChQvstUePHqFNmzag74ysCxYsMFgLTX3QaDTYvHkzbG1tIRaLMWzYsAynhsbFxWH8+PEQiUTIly9fptuJ8nPym4x8+fLhv//+y/ScAgRkNzKrFKYSUZ3UjsvOQTmQffTlyxecPn0ay5cvR58+feDv7w+pVAoigkKhQLVq1TBr1ixcunQpW9I7swJHjhxBnjx5IJFIMHv2bJ33cenSJVSrVo3xK23dujVLLaYPHz4wagtXV1fs2rUrw9c7efIkXFxcIJFIMH/+/Eyvm4hw6tQpNufMmTN/2M9cwK+JzCqFcCLSEFE0EX39/v/X1M7LjpHTgebw8HCEhIRg0KBBKFKkCLMkrK2t0bJlS6xatSrTNNLZjY8fP6Jp06YgIlSvXh2vX7/Wef3QoUMoXrw4iAgVK1bEtWvXsnQ9J0+eZPe2Vq1aGWY1/fjxIxo1asQoNzJT08C7jD5+/MgqxatWrYqXL19meE4BArITmVIKuXnktFJIjFevXmHDhg3o0KEDy+LhOA4VK1ZEUFDQD6MgNBoNVq5cCSMjI9jY2GDfvn06r8fHx2PlypWwsrKCSCRCv379spRTKTY2FoGBgTA1NYVMJsPo0aMz5FLSaDRYsGABJBIJ8ubNq+Mmyyj4NF+lUglra2vs2bMn03MKEJDVyJBSIKKC3/+W1DeSOy87Rk64j9ILjUaD69evY+LEiSzl8kdTELdu3ULRokVBRBg0aFCSWMKHDx/Qp08fiEQiWFtbY9WqVVnqRnn9+jVLOXV2dsb27dsz5Ao6e/YsnJ2dIZPJsHjxYoO4we7cucMKAfv16yc08BGQq5FRpbDy+99jesY/yZ2XnSO3WQop4datW0kURIUKFRAUFJRp9tCsRFRUFPr06QMiQsmSJXHv3r0kx1y9epXRZpQuXdogAd2UcOrUKSaAq1Wrhjt37qR7jvfv36NOnTqsVuPLly+ZXld0dDQGDBgAIoKdnZ1AkSEg10JwH+Uy3L59G5MmTWK+cqlUihYtWuDIkSO5NmD5999/w9LSEsbGxnr7I2g0GmzatIm5zbp06ZKlyi4uLg5BQUEwMzODVCrF8OHD0+1SUqvVmDlzJsRiMfLnz5+pFFhthISEgIhgZGSEFStW/FQpzAJ+DmQ20KwgosFEtJOI/iKigUSkSO287Bg/qlLQxs2bNzFo0CBYWlqCiODh4YFZs2blSuvh+fPnjHW1bdu2eoXw169fMWzYMEilUlhYWGS6wU5qCAsLQ+fOnVmWUkZ8+idOnICDgwMUCkWm24gCYPeIz9Zq0qRJtjY5EiAgNWRWKfxJRGuIqPL3sZKItqd2XlaOHyGmkF5ER0dj06ZNqFixYq62HuLj4zFp0iRwHIeiRYviwYMHeo+7c+cOcylVr17dIJQTKeHkyZOsCK9Jkybp7rMQFhbGhHi3bt0yVIuRXIGbSCSCs7NztnBKCRCQFmRWKVxLy3M5MX4GS0Efbt++ncR6CAwMzFX8/gcOHICFhQXMzMywd+9evceo1WosWbIEJiYmUCqVCAwMTLGIMLOIiYnB9OnToVAoYGJiku7rxcfHY/To0SAilClTJlM8TKTFj/Tff/8hX758BiXrEyAgM8isUlhPRH5a/5choqWpnZcd42dVCjx468Hf3x9EBEtLS0ycOBHv37/P6aUBSKh25msWJk6cmKxF8/TpU9SuXRtEBD8/vywPwD58+BC1atViwfH0pp7u2LEDxsbGsLOzw6lTpzK0Bm2lACS41fjMKX9/fzx58kSgyBCQY8ho9tENIrpORHe+F689IaLH3x/fTO687Bw/u1LQxpkzZ9CgQQMQEZRKJQYOHIhnz57l9LIQGRmJ9u3bg4hQt27dZOsV+EC0lZUVpFIpJk2alKVd3zQaDbZt2wZ7e3tWS5GeDKObN28if/78kEgkWLJkSbrjDMkJ/I0bN8LExATm5uYC26qAHENKSiFZllSO41z1vvAdAJ6m9Hp2wMfHBxcvXszpZWQrbt26RbNmzaItW7aQSCSiGjVqUO3atUkmk1FkZCSp1WrSaDSkVqv1DplMRqampqRSqXT+aj82NzcniUSS+mK+AwAtXbqUBg4cSK6urrRz507y9vbWe+zbt29pwIAB9Mcff1CRIkVozZo15Ovra6jbkwRfvnyhMWPG0NKlS8nBwYECAwOpWbNmxHFcqud+/vyZ2rVrRyEhIdS5c2daunQpKRSKTK/p4cOHVKJECQoPD6eePXvS/PnzycjIKNPzChCQVmSKJTU3j8SWwrlz536K9D+NRoPnz5/j4MGDWLBgAX777TfUrl0bpUqVgrOzc4rMrskNkUjEuJrScqyzszP8/f3RunVrjBgxAkuWLMGePXtw9epVfPz4Ue99/vfff+Hg4AClUplq3+Tg4GA4OTlBJBJh7NixiI2NzarbCQA4f/486xFds2ZNhIaGpuk8tVqNcePGsRqM9AawEyO5YHSfPn0yNa8AAekBZaafQm6Evn4K//77L1WoUIFatmxJy5cvJ3Nz8xxdY1rx/Plzun79Ot2+fZtu375Nd+7codu3b1N4eDg7xsrKitzc3MjW1jbJkMlkdOTIEdqxYwd9+fKFGjduTOPHj6eCBQuSWCwmsVhMHMexnbFGo6GIiAgKDw+nr1+/0tevX9lj/u/bt2/p2bNnbLx48YJiY2N11m1tbU1+fn5UtmxZ8vPzo9KlS5NKpaLXr19TixYt6N9//6XBgwfT7Nmzk+2L8OXLF+rfvz9t2LCBSpYsSRs3biQvL68su9fx8fG0dOlSGjt2LMXGxtLo0aNpxIgRJJfLUz13165d1KFDB1IqlbR9+3aqWLFiptfDcRyFhIRQhw4d6Nu3b7R8+XJq165dpucVICA1/BKWQnx8PKZNm2awXr1ZAbVajRs3bmDp0qVo06YNnJ2ddXaL9vb2qFKlCvr27YulS5fi+PHjePv2bZrm/vz5M8aNGwdjY2OIxWJ07949U9kzidf9+vVrnD9/Htu3b8fcuXPRqVMnFCxYUMe68Pb2Rvfu3bFq1SpGt12/fv1UezD/9ddfsLKyglwux4IFC7I8Bffly5do0aIFiAgFChTAkSNH0nTe7du3UaBAAUgkEgQGBhqEbRVIqP+oUKECiBLafuamLDMBPycog4Hmg0Q0iL5zIOXGoS/QrN2rd9KkSTma/vft2zf8+++/mDlzJurWrcuCi7wCaN68ORYuXIh///3XYMVNb968Qb9+/SCVSqFQKDBixIgsJav7+PEj9u/fjwkTJqBmzZowMzNj71GlUjHBm5rb5fXr16hbty6ICFWqVMHTp0+zbM08Dh48yHpTt2nTJgkjrD58/vyZBfxbt26dKQGuHYyOi4vD2LFjwXEcvLy8cPPmzQzPK0BAasioUrAnok5E9AcRXSaiZUTUkIhMkjsnu0dy2UdfvnxBu3btQEQoX758tggYHm/evMHq1avRsGFDnS5wBQsWRLdu3bB+/Xo8fPgwy2Mfjx49Qvv27cFxHMzNzTFjxgzWVS0roVarcfv2baxZswZNmzaFTCYDEUEsFqNdu3Y4c+ZMsu+dZ2c1NjaGqakpNmzYkOX3KSoqCuPHj4dMJoOZmRmCgoJS3Uio1WpMnz4dIpEIRYoUyTCdtz4cPnwYdnZ2MDIywpo1a6DRaITUVQEGR0pKIU0xBY7jRJRQn1CbiKpSQm+FQwBmp82DlTXQzj6KiIigtm3bUpcuXahevXokFotp06ZN1Lt3bxKLxbRq1Spq0KABPXz4kD5//kwRERFsREZG6jyOiooilUpFNjY2ZG1tTTY2NmxYW1uTsbEx89EDoBs3btCePXtoz5499N9//xEAcnFxofr161O1atXI39+fbGxscuQeXb9+ncaMGUN79+4lBwcHmjRpEnXp0iXD/Y/Ti/DwcFq8eDFNnjyZvn37RkRELi4u1Lx5c2rZsiX5+PgkyQR6+PAhdezYkU6fPk1Nmzal5cuXk7W1dZau8/79+9SnTx86cuQIFSpUiGbNmkX16tVLMUvp8OHD1Lp1a4qLi6MNGzZQw4YNDbKWN2/eULt27ejo0aPUtm1b2rx5M6XldypAQFph8JgCEVkTUduMnGvIoW0pXLx4EXZ2dowDZ/To0di6dSuGDBnCKoNFIlGqmTdGRkawtLRMMVNHoVDAxsYGTk5OzEVCRChVqhSmTJmCq1ev5rosqFOnTrEiuJIlS2Z7zOXFixcoWrQoOI5DkSJF2P0tUKAAli1blqSlaXx8PGbOnAmpVAp7e/s0+/0zA41Gg127dqFAgQIgIgQEBKRa+PbkyRP4+PiAiDB69GiDVWzHx8djypQp7Lt15coVg8wrQACQQffRjzC0lcKBAwdARJBIJOA4LkkA18PDgz1esWIFTp06hStXruD+/ft49eoVvn79qvOD1mg0+Pz5Mx48eICzZ88iODgY8+fPR40aNWBqagr6Tn8tFot1ruXm5oZ69eph5MiR2Lx5Mx48eJBrFIRGo8Eff/yBPHnygL6T2hkqGJ0WhIeHs7hB7969sWrVKpQuXRpECd3qJkyYkIQI8OrVqyhYsCA4jsOYMWOyJUYUGxuLJUuWwMbGhsUOUuJuio6ORrdu3UDfeZ7evXuX6TUkl7o6fvz4TM8tQMAvoRQuXrwIIyMjVjGb+MdUsmRJTJw4EXZ2djAzM8OxY8fSdPPi4+OxZ88eNGjQgCmAKlWqYOvWrYiMjIRarcajR4+we/duTJs2Da1atUKRIkUgkUjYte3s7NC0aVMEBgbi4sWLOc59ExERgbFjx0Iul8PY2BjTp0/PEAFcRhAfH49+/fqBKKEtZnh4OE6cOMGCtwqFAj169MDdu3d11tulSxfQd4qI7IoRffnyBWPGjIGRkRFkMhmGDBmSYtB+1apVkMvlcHFxMUhXNx5ExGhCmjdvjs+fPxtsbgG/Jn46pUB6WFLj4+PRq1cvfPnyBbGxsTh58iTGjBmD4sWL61gORkZGsLCwgEQiwebNm5O9aU+fPsX48ePZrtrOzg4jRoxIlhU0MWJiYnD9+nWsWLEC7dq1g5ubG1uDsbExqlWrhkmTJuHo0aNJXCfZhYcPH7K+xfny5UNwcHC2WTULFy6ESCRC6dKlWdrtnTt30L17d8jlcnAchwYNGuDkyZNsTZs3b4aJiQksLCywa9eubFknkJAy2rlzZ3AcBwsLC8ybNy/ZzmoXLlyAi4sL5HI5Vq1aZZDrExHUajVmz54NsViMvHnz4r///tM5RghGC0gPMqUUiGgAEZkSEUcJFNqXiahGaudlx9C2FBYtWgQigkwmw9y5c3VcQW/fvsWGDRtQqVIltoPnFcXIkSN1BOG5c+dQt25dcBwHjuNQs2ZN/PXXXwapuH3+/Dn++OMP9OnTB8WKFWNrUCqVaNCgAVatWpUjbToPHTqEQoUKsWpf7V16VmL37t1QKBTw9PTU2f2HhYVh/PjxsLKyAhGhbNmyOHv2LADgwYMHKFmyJIiyv+3ltWvXULNmTRARHB0dsWTJEr0W1rt371C9enUQJdBwZ3aN2gL/zJkzcHFxgVQqxYIFC9h3l77XPAgQkBZkVilc+/63JhEFE1ExIrqc2nnZMbSVgkajwYgRI1gw2c3NDYcOHUpyM2JiYjB06FAYGRmxnbujoyOCgoJQr1495t8eO3ZslvcA+Pz5M0JCQtCnTx+4uLiw9fj4+GDSpEm4fPlytu3cY2NjsWDBApiamkImk2H8+PHZInBPnToFMzMzODk5JWFPjYyMxNKlS1k3t/bt2+Ply5f49u0bBg0aBCJC8eLF9bYIzUocO3aM9YpwcXHBypUrk2watGm4fXx8DOry+vDhA3O3NWjQAOXKlROUgoB0IbNK4fr3vwuJqPH3x1dSOy87hr46hVu3bsHV1ZUJ2Jo1a+qlao6JicGYMWN0AsUcx6F69epp5sUxJDQaDa5fv45p06bBz8+PWRFOTk7o0aMHDh8+nKW9CHi8efMGbdu2BREhf/78OHr0aJZf89q1a3BwcICFhQWzCLQRHh6OUaNGQSaTsRhIdHQ09uzZAysrKxgbG+P333/P8nVqQ6PR4NChQyhTpgyICHnz5sXatWuTxIv+/vtvmJqawsrKCocPHzbo9XmrJfEQXEkCUkNmlcI6IjpERA+ISElEKiK6lNp52TESKwWeHiE6Oho9evRggl4kEqF27dpYtGgRgoKCsHnzZmzduhVNmjTRm6oqlUrRtWtXPHnyxCAfQEYQFhaGdevWoUmTJjAxMQERIU+ePBg5cmSGGtWnF4cOHWLVvu3bt08z3UZG8ejRI+TLlw9KpRL79+/Xe0xoaCiLgbi7u2PXrl06LUIN4apJLzQaDfbt28fSUj08PLBhwwYdBX7v3j0ULlwYIpEIM2bMMKj1d+HCBbi7u4OIMGPGjFzVpU9A7kVmlYKIiEoSkfn3/62IyDu187JjaCuF48ePQyKRQKlUJklJTW2oVCrmoiBKqL7lWUV79eqVrWmb+hAdHY1t27ahbt26zLLx9fXF4sWLs7ThTlRUFMaOHQupVApLS0usXr06S4XOmzdvULx48VSTAA4fPsxab1arVg3Xrl1jrpoSJUrg4cOHWbbG5KDRaLB7924UK1YMRARPT09s376dKYDw8HC0bNmSZV0ZMoPo8+fP7LtbvXp1vHnzxmBzC/g5kVmlcDQtz+XE0FYKJ06cYMVrSqWSWQbaufDatBMBAQEIDAzEzJkz0b9/fzRt2pQVVyVWGmKxGLVq1cLVq1cN9ZlkGK9fv8a8efOY8JFKpWjcuDF27dqVZU1rbt26xQjbKlSokKWd0z5//sx2/gsXLkz2uLi4OCxatAjm5uYQi8UYMGAAduzYAXNzc5ibmyM4ODjL1pgS1Go1duzYwQL3Pj4+rPBOo9FgwYIFEIvFKFCggEH5jcaPH4+VK1dCoVDAzs7OoK4qAT8fMqQUiEhBRJZEdI2ILL4/tiQiNyK6k9x52Tm0lcKbN28wZ84cdO3alZnTMpkM3t7eCAgI0AksSyQS5lLasmULTp8+jefPnyM+Ph6PHj2Cu7s7ZDIZpFIpRCKRjnvJwsICbdq0wYYNGwxGYpdRXL16FYMHD2bK0N7eHhMnTkwTsVt6oVarsWbNGlhYWEAqlWL8+PFZ1gMhOjqauYnGjh2borvl3bt36NmzJ4gIhQoVQkhICMtOGjlyZI7VhMTHx2PdunWMCbdatWqsdoHfwBgbG2P79u0Gve7169dRqFAhcByH0aNH53hNjIDciYwqhQGU0H4zhogefX/8+LuS6Jvcedk5tJXC9evXmeCWyWQwMTGBXC7XCSK7uLgw5ZC4Epl/zsXFBaVLl4ZKpYJSqUSpUqVARChSpAgKFy6cJPbQrFkz7N27N0d/fHFxcdizZw8rcJJKpWjbti3Onz9v8GuFhYWxQHSpUqVw+/Ztg18DSHhPXbt2BRGhR48eqbqtDh8+DAcHB8hkMsyaNQu//fYbiAiVK1fOUXdKdHQ0FixYwNJrmzdvjrt37+Lly5fw8/MDkWHpMYCEYj/+3mkX+wkBaAE8Mus+6pfaMTk1EvdTCA8PZz+u8+fPo0iRIiBK4ELiOA5GRkbo0aMHyxhp2LAhUx5t2rTBqFGj0L59e/j7+zN2TyKCg4MDpFIp5HI5Bg0ahObNm4O+F8LxQWA7OzsMGTIE169fN8RnlmHcu3cP/fv3Z5xMvr6+2LRpk8FdSzt37oS1tTUUCgUWLlyYJbEGjUaDkSNHgojQoUOHVAXnu3fvWKpmzZo1sXDhQhgZGcHBwQEnT540+PrSgy9fvmD8+PGs38Vvv/2Ghw8fMnqMunXr4tOnTwa95pYtW1ixH7+hESAAyKRSSDifyhFRGyLqwI+0nJfeQUSNiGgVEe2mNBTI6UtJjYyMxJAhQyASieDk5IS9e/cCAO7evYtOnTpBLBbD3NycxRlatGiBatWqgYhQq1YtFlSOjo7Gxo0bYWxszFxJvJIwMTFBw4YN4e3tzTJOKlasyArjSpYsiUWLFhmEAyej+Pr1K4KCguDp6cmU1vjx4w3qWtLugVCtWrVMt6pMDjwxXKtWrVJ1WWk0GixbtgxGRkawsbFBUFAQ8ufPD7FYjLlz5+Y4D1Xifhdjx47FggULIJFIUKBAAYNbXtrFfkSUY9XzAnIXMmspbCSiM0S0lIiCvo9FqZ2ndf5aInpLRDcTPV+LiO4RUSgRjUz0mgURrUlt7sRKYc2aNSye0KNHD70ZHrdv30blypVB34vWiBK4jObMmcMoMLZu3cqOv3HjBqysrJAnTx5s3boVjRo10nE92drasl15ixYtMHnyZNYLWCqVol27dlkamE0NarUaBw4cYFXacrkcvXr1wqNHjwwyv0ajwYoVK6BUKmFubp5qb+aMYvbs2SAiNG7cOE1Wz61bt1gwvnv37mjYsCH7jHJDZ7NHjx6xbKQ8efJg/PjxsLGxgUqlMmiQPDliPcGV9Gsjs0rhDlFC34WMDCKqSAkprTe1nhMT0UMicici2fc4hZfW6/OIqGRqc+ujubC1tU2V7E6j0WDLli1wcHBgdQwFChTA0aNHmZ+3c+fOLOf98uXLMDc3h7u7O168eIFXr16hVq1aICJYWVmxjCVe6E6cOBEXL15E//79WcZT48aNDUqSlhHcv38f3bt3h0wmY01vDJUB8+DBA3bvWrVqlSVB+IULFzJXS1rqEbQrnwsXLoxBgwZBJBKhaNGiOVKgqA8nT55kmwgfHx/W4nTy5MkGdcnxGV12dnaQy+UICgrKcatJQM4hs0phOxE5pHZcKnO4JVIKZYnooNb/o74PjohmEVG1FObqTkQXieiii4sLe5NxcXGwtbWFq6trmhk/v3z5wgQFx3FQqVT4999/MXbsWPYjffbsGYCEGIVKpYKnpyfevHkDjUaDoKAgyOVy2NraomvXrqzYi3cxTZgwAa9evcK4ceNYK84aNWrgxIkTOfqDfPHiBQYPHgxjY2PQ97x5QwSl4+LiMHXqVEgkEjg6OmZJWuTy5cvZfUxrJ7n9+/fDzs6OuWssLCxgYWGBAwcOGHx9GUF8fDxWrVoFW1tbcBzHaN4bN26Mr1+/Guw6RISwsDDUqVMHRIR69epleVGigNyJzCqFY0T0iRJ6NgfzI7XzEs2RWCk0I6LVWv+3J6LFRNSfiC4R0XIi6pnavInrFHiBPGvWrHTdoOvXrzO/K8dxmDt3Lv7++2+oVCrY2Ngwy+PUqVNQKpUoUqQIixdo8/2PGDEC586d06HZNjIywtixY/Hs2TPMnDkTtra2LCskJCQkR5XD+/fvMWHCBFhYWICIULVqVRw5ciTTa7p48SJLi5w0aZLBg9Dr1q0Dx3GoVKkSwsPD03TOmzdvULZsWXAch5EjR8Lb2xscxxm8wjgz+Pz5M4YMGQKJRMKYYgsVKpRmZt7UwLuMNBoNFi1aBLlcDnt7exw8eDDJMQJ+bmRWKQToG6mdl2iOxEqhuR6lEJSO+ZJQZy9ZskQn/TS95jH/Q+GDxbVq1cK1a9dQsGBBiMVixkh59OhRKBQKlChRgmWLREREsBRIX19fPHz4EBEREWjfvj2rcVAoFBg4cCDu3buHoKAglr9evHjxbKWs1oevX79i7ty5rKrbz88Phw4dytSaIiMj0b59exAl9AIwtDtp8+bNEIvFKFeuXJqrg6OioljmWJcuXZhPv1mzZmlWLtmBu3fvsgC+SCSCsbEx9u3bZ/DrXLt2jVWGDx48mBUoCvj5kSmlYIiRVvdReufVthRiY2PRqlUrnWCag4MDlixZkq4f/MuXL1kPhbx58+LChQuskKpNmzaIjIzE/v37IZVKERAQoOOq2r59O8zNzaFSqVjA9eHDhywFlr7XQrRq1QpnzpzBunXrkD9/frZLv3btWprXmRWIjo7G8uXLmcIqX748/vnnnwzPp9FosHTpUkilUri5ueHSpUsGXG3C/ZZIJChdunSKzW+0oVarWZprzZo1MXXqVIhEIhQuXBj379836Poyi3379un04UhM824IREVFoU+fPszlKSiFXwOZtRTCiejr9/GNiNRE9DW18xLNkVgpSCihIC4v/T/QXDgd8yWxFF69eoUGDRqgdevWSTItlEolWrZsiV27duHWrVv48uVLijcsMjKS0RTI5XIsWbIEU6ZMAcdxKFasGB49eoQtW7aAKKFVo7Z75MmTJ6wXcseOHREeHg6NRoNNmzbBwsKCBaOJCBUrVsRff/2FwMBAWFpaQiQSoXv37klaUmY3vn37hiVLlrDsrEqVKmUqz//8+fNwdnaGXC7H6tWrDbhSIDg4GDKZDCVKlEiXNbJy5UqIxWJ4e3tj8+bNsLKygpmZGUJCQgy6vswiOjoaY8eOZRZn8eLFDd55TchQ+vVgUEuBEmoJpqfj+K1E9JqI4ojoBRF1/f58HSK6TwlZSGPSuw4kshQ2bdqk94ud3LC1tYW/vz86deqEkSNHYuLEiRg/fjwGDBiAZs2aoXbt2izVlL6nr5YtWxYymQwymQw1atRApUqVmGJ49OgRq2qOi4vD+PHjIRKJULBgQZaS+v79e3Tq1IllLfHxBW9vb2zduhUDBgyARCKBSqXCrFmzsq1FZnKIjo7GwoULYW9vD6IEsrUzZ85kaC7txjNdunQxaL78/v37IZPJULJkyTRbDEACE6ypqSkcHR2xd+9e1qVv+vTpuSbOwOPevXsoUKAAc0VmRee5smXLsu97uXLlcp3lJMBwMLj7iIjOZeQ8Qw9tpXDjxg0mZB0dHRkfkEgkgkKhgEgkgouLC6tUNjY2TpZRVSwWw9TUFHZ2djo1CWKxGNbW1qyQLfG5fKvEqlWr4rfffkPXrl1hbm4OpVKpw/p59OhRlmFSvnx55iIoX748Nm/ejPr16zP3lTbTZk4hMjIS8+bNY43sa9eujYsXL6Z7nvj4eJbZZWg205CQEEilUpQuXTpdO+kbN27A2dmZ8RDxlmbLli1zRT2DNjQaDUaPHs2+dw0aNDCo1cCnrW7atAnm5uYwMjJCYGCgQMf9EyKz7qMmWqMZEc0korOpnZeVQ5/7aO/evVAoFIxjhv+b2rCyskK9evUwbNgwzJ8/H9OmTUPHjh1ZERyvWCQSCaytrUFE6NmzJyuGat++PcqUKQOO49CiRQu0bt0afn5+TEFpj4IFC2LZsmW4desWIiIiMGbMGEgkElhZWaFNmzZM6JYtWxb9+/dn3dg8PDwwbtw4bNy4Edu3b8fevXtx9OhRnDlzBpcvX8adO3fw5MmTLA+WRkREYObMmbC0tGTusZcvX6Z7nr179zI2U77i3BAIDg6GRCKBn59fqi5Cbbx69QqlSpWCSCTCokWLMGvWLHAch+LFi+doT43kcPXqVfb9VqlU+OOPPwy2ceBdRi9fvmSpqwEBAUyBBwQEwNXV1SDXEpBzyKxSWKc1VhHRGCKyTe287BjalsLHjx9RsGBBWFhYoEiRIhCLxUkUQ548eSAWi2FnZ4fGjRvrFdz58+dHs2bNsGTJEvz777/YsGEDmjZtmoRAz9TUlLVBrFu3LooXLw6lUqnTUP3Tp0/4888/8dtvvyW5llgs1iHsM9QwNTWFl5cXqlevjk6dOmHs2LFYtmwZgoODcfnyZYPsLL98+YLhw4ezTmjTpk1Ltzvo4cOHzF0ze/Zsgwm1nTt3QiwWw9/fP11KMiIigvEmzZ49GyEhITAzM4ONjQ1OnDhhkLUZEhEREcwdR5RQt2HIlp9AgmWydu1amJqawtjYGMuWLUPFihVBQjD6h4fB3Ue5ZWgrhSdPnqBdu3Y6NNdSqRQqlQre3t7o2LEj84tbWlrCysoKJ06cwOXLl5lwt7Kygp+fn06/ZC8vLwwdOhSrV69mdQsVKlRI0q3N0dERVlZWMDExQfPmzVGqVCkoFAq9gpvjOFhZWbGCNiJieelKpRIDBw5Eu3btoFAoIJFI0LJlS1ZB7eXlha1bt+L48eM4cOAAdu3ahS1btmDNmjWYNWsW+vXrhyZNmsDX1xdOTk5J1skrx5o1a2Lw4MFYu3Ytzp8/nyErIzQ0FI0bNwZRAungn3/+mS7hHhkZyVJEO3fubDDSvj///BNisRgVK1ZMlwsoLi6OZbDNnz8fd+/eRYECBSCRSLBs2TKDrM2Q0Gg0zKrhOA7GxsYICgoyuLvn6dOnTAHx31kBPzYyaynkIaJdlMBfFEZEfxFRntTOy8qhz330+PFjGBkZoUKFCpBKpfDy8oJSqUTevHlBRNi2bRvGjRsHooRsmgIFCkAqleL333+HRqPBnj17GMVApUqVsHfvXsyfPx/VqlVjMQSlUgmxWAwHBwdMmjSJ7f4TxxY4joOfnx8GDx6MpUuXYs+ePbh69So+fPiABw8eoHjx4iAijBkzBu/evcOmTZvQtGlTHSXi4uKCuXPnomvXrozEr0ePHrC2tmb00GmhW46Li8OLFy9w7tw57NixAzNmzEC7du1QokSJJJaKq6sr6tati7Fjx2L//v1ptiqOHj3KyAErVqyIy5cvp+k8ICFFdPz48excQ5EIbt26FSKRCJUrV05z5TOQcL+aNWsGIkJQUBA+ffrE3Cg9evTIskZGmcHhw4dhZmbGamzKlStncGI9/jNKPIQMpR8TmVUKh4moMyWkkUqIqBMRHU7tvOwYiS2FDh06oGnTpowIje+doC+YrJ0a6uLigsaNG6NLly6oXr06lEolJBIJhg0bhpiYGHz+/Bnz58+Hj48P++Hxg+c2srGxgUwm07lW2bJl9bbMjIqKYpTJVapUYXz/UVFR+OuvvxgXDm/tVK5cmRUZ+fr6MlZXPz8/3Lt3L8NfjPj4eNy/fx+7du3ClClT0KpVKxQpUoRZFxzHwdvbG7169cLmzZvx5MmTZC2B+Ph4LF++HNbW1uA4Dl27dk1XH4MtW7ZALpfD3d3dYAJt48aN4DgO1apVS5d7KzY2ltWmLFu2DPHx8ay2oXz58jmeMqwPjx8/1vney2QyTJo0yeBKzNfXl303q1evnitjLgJSR2aVwtW0PJcTQ1sp3L17F87OzvDy8oKfnx8rCtP25Tds2BBBQUGsWtTNzQ2urq7MF29nZ6dDkc0P7XhCvnz5EBAQwMx1XnCLxWKYmZnBxcWFHW9mZsZeb9KkCXbv3q1D/bxu3TooFAo4OTkl4R66fPkyew+8cjM1NYWRkRGkUikaNmwIc3NzKBQKg2eIhIeH48iRI5g4cSKqV6/Oipp411PLli2xcuVKvbv6T58+YfDgwZBIJDAzM8Pq1avT7FI6d+4c7OzsYGZmpkO9kBmsX78eHMehVq1aaSLR4xETE4N69eqBiLBq1SoACYrLyMgIzs7OGcq+ympERkaiXbt2ICI4OTmBKKE51Llz5wx6HV5ZmpiYQKVSYcWKFTmeIScgfcisUjhCRO0ogdlU/P1xruvRfP/+fdSqVQutWrVCjx49MGzYMLar4ekbFAoFRo8ejbVr16J3794gIrRt2xZz5swBx3EoWbIk5s+fDz8/vyTWRWJXC1/D4O3tjRYtWujEB7Qfd+rUCYMGDWKZRTY2Nhg4cCCuXLkCALhy5Qrc3Nwgk8mSFHbFxsZi6tSpkEqlsLCw0NmlESUwXvLPVaxY0WB02IkRFxeHy5cvIygoCC1btmQCRywWo3r16li9enUSa+ju3bssxbFKlSppZiV9+vQpihUrBrFYjKCgIIOsf/Xq1SBKSAhIz87527dvqFWrFjiOw/r16wEAly5dgouLC+RyOdatW2eQ9RkSGo0GCxcuhFgsRp48eWBvbw+O4zBgwACDZafxLqPHjx+jatWqIEropyFYDT8OMqsUXCiBBO/d97jC30Tkmtp5WTn0xRT++OMPvT7PtAxzc3Mm6IgSgsZjxozB2bNn8eLFC5aVUrp0aaxfvx4LFixAp06ddIS/iYkJihUrplPwxo9Vq1YhNjYWwcHBaNq0KauV8Pb2xrx583D37l3mEurZs2cSwXXt2jXWRa5Dhw4YOnQos0LoeyzA2NgYZmZm2L17d+a/MalAo9HgypUrGDVqFGOGlUgkqFWrFtauXcsKyNRqNZYvXw6VSgUjIyPMnTs3TXGQ8PBwds979+5tkF7Qy5YtA1ECz1F6WqdGRUWhWrVq4DgOmzZtAgC8ffsWVapUYevLjXGG48ePw8bGBiYmJqxNq5ubG44cOWLQ62g0GixfvhwmJiYwMTHB8uXLBavhB0CmlEJuHtqWQkxMDJ48eYLLly8jJCQEa9aswdSpU9GtWzeYmpomycLhOA6WlpZ63UW8kO/QoQN27dqFiIgIbNy4MUlBj1qtZj+4WrVqMWoM7Wvwj+vWrcvSVd+/f4/Fixez7m9SqRQtW7ZEmzZtQN8Dha9evdL5EKOjozFw4EAQJTSoP3fuHNavX894ivg1ExFGjBiRbT2jNRoNLl26hBEjRrCgvlQqRZ06dbB582bExsbi+fPnzBVTunTpNLUsjY+Px7Bhw5jvOj11B8lh3rx5IEqor0iPuy0yMhKVK1eGSCTCH3/8ASDBeho6dCiIEhhvE39euQHPnj1j37EOHTowd2S3bt0MTpXx5MkTtrHhN1kCci8yayn8TkTmWv9bENHa1M7LjqGvHac2NBoN3r17h507dzL/O+/2SFx3wAtx7UAyL9TFYjHc3NxQo0YNRjVQpEgRHDlyBJ8/f0bp0qWhVCpx6dIlhIWFYfv27Tr02dojX758WLJkCcuIuXnzJvr37892/k5OTpDJZLC3t8fp06eTvKdDhw6xntGzZs2CWq3Gzp07WQU3P3x9fbNdUGk0Gly4cAHDhg1jsRpnZ2fMmzcPX758wdatW2FjYwOJRIJx48alicZj7dq1EIvFKFGiRLoC18lh4sSJICL07ds3XTvaiIgIVKhQAWKxGDt27GDPb926FUqlEg4ODhmmAMlKREdHo0uXLky59uvXDyKRCI6Ojgbt8AYkfP48pxQRCVZDLkZmlcKVtDyXE0NbKXz69Al//fUXpk+fjo4dO8LPz4/1CdA3ePrrunXrokqVKsznn9iiMDIygkKh0JvBxB/v5OQEuVwOpVKJadOm4ezZs/j06RMePXoET09PvedKJBJUqFABe/bsgUajQWRkJNauXctiBHxHuBEjRiTZ1b5//x5NmjQBUUL67LNnz1iGjPb6jYyMsqw9ZmpQq9XYs2cPiyuYmZlhxIgRuHHjBguGFipUKE2CdN++fVAqlXB3d890bwGNRoMhQ4aAiDBq1Kh0nfv161eUK1cOEolEp4HQtWvX4O7uDqlUimXLluU6Qci7eGQyGVxcXLB+/XoULVoURAnMv4buJa7NCly9enWDF9UJyDwyqxSuEZGF1v+WRHQjtfOycpCemMLVq1d1YgKVK1dGz549MWvWLHTu3BkmJiZJBH7FihVRrFgxndiAIQdfX6AdZ5DJZKhYsSJMTU11hHezZs1w+vRpqNVqXLp0CR07dmQ7LktLSyxbtkwne4avNjUxMYG5uTnrK33x4kWdDnBEhFKlSmUqdTWzOH/+PJo1awaRSASpVIrOnTszim6RSIQxY8akGjc4d+4cIxHMbOaPRqNB9+7dQUSYPn16us79/PkzihQpAjMzM53e2x8/fmSuxC5duqQr0ym7cOHCBbi6ukImk2Hx4sWYMGECpFIpbGxsDEKVkRzbqkwmw6pVq3KdsvyVkVml0IES+jRPIaLJRHSXiNqndl52jMQxhUuXLrH2hfHx8Zg7dy7j6eHTOrWHVCpFrVq10KdPH8yaNQtbtmxhFbqjR4/GkydPmE+WTz81NTVF69atMX/+fLRu3ZoJbrlcziqj7e3tUaZMGbi5ubHz9A2JRAKlUqmjrExMTNCpUyccPnwYYWFhrJKZKKHiesKECTp58qGhoaw3cvv27fH161dERkay7Co+qE2U0Ds5J3sTh4aGok+fPuyzqFGjBisM8/Hxwd27d1M8/+7du3B1dYWJiQkOHTqUqbXEx8ezGE56s5yePn0Ke3t7uLm56bi04uPjMWbMGBAlxE74Vq65Ce/fv2ffqQ4dOuD8+fPw8fEBUULKdka4rBKDtw4B4NGjR6hcuTKIEvpXPH/+PNPzC8g8MqUUEs4nLyLqS0T9iMgrLedkx9BWClevXoWpqSksLS2TsJ9KpVLkzZsXlStXRo8ePVghGBGhXbt2OjsYtVqNtm3bgiihCjomJoZVuHbo0AEdO3Zku3wXFxdMmjQJLVq0YNfjA9etWrVic0ZFReH06dNo2rRpssohOUujYMGCqF69OqRSKYuJyOVydO3aFTdu3ACQEPScMGECRCIR8ufPzyqKQ0JCYGtrC7FYzNbHF+UZOtCYHrx79w6TJk1iBINly5ZlQfylS5emuKN8+fIlvL29IZVKM+0ai42NZcSG6U0vvXDhAoyMjFCmTJkkhXG7du2CSqWCra1tpvpQZBXUajUmTZoEjuNQtGhR3L59G7Nnz4ZCoYCpqSmWLFmSqboXbaXAX2/JkiVQKpUwNTXF2rVrBashh5FppZBbh7ZSWLRoURKhqlQqYWdnBwcHB1hYWCRLQOfn54cZM2Zgz549ePLkCaKiolCuXDkoFAqcO3cOcXFxjDtp+PDhiIqKwvbt29kOSC6Xo379+rCxsWGxAKIE6uu///5b5wewb98+9rpEIkGPHj3QsmVL1rOAKGlcI/FwdHRkFkD16tVx4MABaDQaHD9+nAWqFy5cCI1Gg7CwMJbeKZPJ2Nw2NjZYsWJFmlJEswqRkZGYPHkyjIyMIJfLmdurTp06eP36dbLnffr0iQme+fPnZ2oN3759Q/Xq1SESibB9+/Z0nbtz505wHIdmzZolEaK3b99mvEmLFy/OlULwwIEDsLS0hKmpKXbu3IkHDx6wugM/P780ZYklB330F6GhoYxQr06dOnjx4kUmVi8gM/gllMJff/3Fsoc8PDzg5uamN2hcuHBhNG7cmPl/9QWBVSoVqlWrxojzHj9+DLVajV69eoEoITedFwI3b95Er169mJuID25rxwzKlCmDo0ePsrXu3buXXZfjOGzYsAEajQYPHjzAmjVr0LFjR526icQCXd+wtLREr169cOPGDZb+2aBBA7x//55lhRgZGTGrhCf98/b2zlTLTUPg2bNnrI+Bubk5pFIprK2t8ffffyd7TnR0NLO8hg0blqmdbUREBPz9/SGVStPdeW3u3LkgSkgDTozPnz+zz6Jz5865Ms6g7SIdNmwYYmNjsWHDBlhZWUEikWDUqFEGbYikVquxcOFCGBkZwcLCgqX4Cshe/HRKgfQEmj98+IA2bdpg/fr12L59OzZv3ozVq1dj6tSp6NOnD5o1a4YKFSrA09MzSXxBJBLBwcEBv//+O4KCgtCtWzfWBId3CXXp0gU7d+5E//79QZSQ665dC/Dp0ycEBgbqnCcWiyGVSpkVUK1aNUZnsWLFCp01LFy4MMkH9+zZMyxduhQlSpTQUV7aqa5mZmYsbsIPe3t7+Pv7QyKRwNHRkVE/3717F4UKFWJzlS1bltU5NGrUKNOZPZnFqVOnGO8Tr2S7deuWbCVufHw8i5106NAhU7UZnz9/RsmSJaFQKNLl8tFoNOjRoweI/k+HoQ21Ws2IGH19fXOlT/3bt29sw1OuXDk8efIE7969Y9Zxvnz5dLKtDIH79++zLKU2bdqkq2OegMzjp1MK/EgcU0jJ5ZLeIZVK4eLiwlwa/C5dLBazHPy6desmcb+o1Wrs379fJy3P2NgYQ4YMYT70xo0b49atW+jXr5/ONVNinPz48SP69eunN2DOD2dn52SbC5UrVw7BwcF4/fo1o4cmIhQrVgwjR46EsbExpFJpjscb4uPjsXr1apYizLvhkmNe1Wg0mDx5MruvmWlh+u7dOxQsWBCmpqaMhiQtiIuLQ82aNSEWi5MVnjt37oSJiQns7Oxw6tSpDK8xK7FlyxaoVCqYm5vjr7/+AqDbJbBdu3Z4+/atwa4XFxeHyZMnQyKRwMnJyeCKR0Dy+CWUwvv37zFkyBAMHz4cLVq0QIkSJZJk/uTLlw/169dH8+bN0aRJE9SrVw958uTROYZv3Zmc4BWJRDqxCVNTU0ybNk0v78uxY8d00lGbNGmCnj17QqVSQSQSoV27dihdurSOa6h3794p+p+/fv2KmTNnMuGfXA0Ez/ckFot1FIlUKkXVqlXRqFEjiEQicBwHW1tbHDhwAJ07dwbHcbC3t093bwRD4/PnzxgyZAgLkkul0hSLoQIDA0GUUFmeGXfHs2fP4OzsDFtb23T1KP7y5QuKFCkCU1NT3Lx5U+8xt27dQv78+SGRSFINqOcUQkNDWTZSz549ERUVhaioKNYl0NLS0uCB4gsXLjDa+v79+xvUXSVAP34JpbBy5Uom+MzNzVGxYkX07dsXq1atwvnz55NttvLy5UvIZDLWRczY2Bh37tzBwYMH0bhxYyZQUwv+8jvaESNG6FQSf/nyhe20eLdPpUqV0KJFC8jlckgkEhgZGcHd3Z0pm2bNmiXJiErsM4+IiMC8efN0lEO5cuVYv2d+8HOKRCKIxWIoFApm6Wi/L4lEgo0bN+LChQsoWbIkiAj16tXL8cKjO3fusP4TRIQWLVok605atWoVOI5DpUqVWGpyRnD37l1YW1vDzc0tXcHQJ0+ewM7ODq6urslWX2v3Z+jatWumLJusQkxMDKMYKVKkCFNyN2/eZGnXAQEBuHPnjsGuGRUVxVyzBQsWzJUstD8TfgmlEB0djeLFi0OlUumld4iLi8ODBw9w8OBB7NixAxs3bsTKlSuxcOFClCtXju24RSIRLC0t0a9fP0yaNAnz5s1D//79dZhTLSwsEBAQAB8fH71UFkQJNQUtW7bE1atXERsby+gx8uTJw2IAzs7OOkKcD+4RJQTF7e3tdawdqVQKExMTWFlZwdHREXnz5oWnpyccHBzY2kxMTNCkSRPUr19fb7tRfrRs2RKLFy9GjRo1dBRevnz5sHDhQowbNw5KpRLGxsYIDAzM0Syl2NhYjBo1iq3R1dU12cyYTZs2QSwWw8/PD58+fcrwNS9evAgTExMULlwYHz58SPN5//33H4yMjODr65vsjle7niG39mcAErKTbG1tYWRkhJUrV0Kj0UCtVmPlypUsIWDcuHEZCqAn1+v50KFDcHJygkQiwZQpU7KNw+tXwy+hFIAE81wqlaJSpUoICgpC//79Ubt2bWayp7bT13bFaBd9pTTMzc1ZINnc3FyvkpBKpXBwcEjW2uDdOLxQ5+sRLCws0L9/f0yYMAETJ07EyJEjMWjQIPTu3Rtdu3ZF27Zt0axZM9SvXx+VKlViDX/SOoyNjbF7926Eh4cz2gztUbRoUabMfHx80uVnzwocP36cxWWkUmkSqnEeO3fuhFQqRfHixTPlA//nn38gl8vh5+eXrraef/31F3MFpoRt27Yxyy0z6Z9ZidevXzOiu+bNmzNF++bNG1bP4+Hhke54QOJaBm18/PiRZaOVK1cOjx8/zuS7EJAYv4RSuHjxYhLBb2xsjOLFi6N58+YYNWoU1q5di+PHj+PatWu4f/8+nj9/jvfv3yMyMhIDBgyAWCxGuXLlmGA/cuQI3r59i3v37uHcuXPYt28fFi1aBDMzM5iYmMDLyytNyiY5SyI4OBgbN25k/lR+LjMzM+biyZMnT5oZQjUaDWubyL8HPz8/DB8+HIMGDUKdOnX08kGJRCJUqFABHTt2ZMV3Tk5OKFWqFDuG9+23bt06U66ZzOLDhw8szZN3tekT2Pv374dCoYCXl1emqnR37doFkUiEGjVqpIsie/DgwSAi7Ny5M8XjLly4AEdHR5iYmBicoM5QUKvVmDlzJiQSCdzc3HD27Fn22uHDh5l7tG3btmm2elJSCjwKFSrEmldt27YtU+9BgC5+CaXw+vVrlC5dGjNnzoSTkxPy5s3LgmS3bt1CcHAwFixYgIEDB6J3797o2bMnunfvjt9++w1du3ZFy5YtWWaRkZERRCIRzM3NsXr1ahw7dgz37t1DREQEPn36hGnTprGAs75Ar1wuh6urK1xcXGBlZZViPMLR0RGdO3dG586dk/R85jmZlEolVq5ciatXr6YpM2jbtm2QyWQwNzdnNBqdOnViO66YmBhs3bqVUSknN1QqFQIDAzFv3jxUrVqVvQ+RSIQ6derg33//zZFgqUajwapVq5gSzZMnj94WnseOHYOxsTHy5cuXqQYwa9euBVFCPCOtbrSYmBiUKlUKFhYWqcZlXrx4AR8fH3Ach1mzZuXKADSQwD/l5uYGsViMmTNnsjhXdHQ0xo0bB6lUCnNzc6xYsSLZupHk+JH0Zd4FBATA0dGR0bh069YtXRabgOTx0ykF0lOncODAAUilUshkMrbz1tfwRi6Xs4Yg2sPY2JgJmeR6LKR3yGQyWFhYwMrKirXSzKhloT0UCgVcXFzg7++Pzp07Y968edi7d6+Oq+Ts2bOwsbGBqakpmjZtCrlcDqlUir59++oEQSMjI1GpUqVUr1myZEn07dsXrVu31mnP6e7ujkmTJmVZ17eUcO/ePabYJBIJIwXUxtmzZ2FmZgZnZ+d0ZRMlBl+k1qNHjzQL7QcPHsDExAT+/v6p+sYjIyPRokULECXUwOTGADSQEChv3rw5iBKq6bW/S3fu3GEWQLly5RgNiz6kxVLgj4mNjcXo0aPBcRwKFiyIq1evGuz9/Kr46ZQCP7QthTNnzugIMX63bWpqynz0+gbHcRCLxZDJZJlSBhKJBC4uLqhTpw6jJeYzf6pUqYIjR44wYfL161eW9pfekRLBHlFCpXLHjh3x999/48SJEyhUqBAkEgnmzJmD7t27QywWw9jYGOPGjdNxS23ZsoXVKuizfpydndn94TiO9W/QVhAVK1bE6tWrDdIQJ62IiYlhhVdEhMGDByfZzV++fBnW1tawt7fPVMbMyJEjQUQYM2ZMms/ZtGkTiAjjxo1L9ViNRsP6Pfj7++faADRfIa9QKGBnZ6dDTqjRaLBu3bpUK6LToxR4HD16FA4ODjo0LgIyhl9CKUyaNMkgu/uUhL69vT2KFSuGGjVqoG3btmjXrh3MzMxYHr2+86RSKXvNysqKmd/aiiu9yqhAgQL47bffEBgYiPXr12PkyJFo0KABo6LWPlYulzOLqWLFili8eDHq1q0LIoK1tTUWLFjAdqX37t1DsWLFQJRQ7Zw4cM3HHlq2bMmsMaKEeEPZsmWZb1mhUKB169Y4fPhwtv1w161bx957QEBAksyjW7duwc7OLlOKQZtye968eWk+r2PHjuA4Ls10ItoB6GvXrmVordmBGzduoHDhwiAijBw5Uof+/N27d+jUqRN4a/LgwYM65yaXfQSk7GJ69+4diynVrVvXoMV0vxJ+CaUwY8aMFAWpIawBfYMnwFMoFGjTpg1atWqFokWLppoJZGRkxAjTLC0tmXWRmOFVm2AvuWFsbIwCBQqgTZs22LhxIwIDA1G+fHl2nqmpaRIyQDMzMxZ0tra2xowZMxATE4Po6Gi28/bz88P8+fP1BqdNTU3Rpk0bNG7cmN1TmUyGRo0aoWnTpiwe4uXlhRUrVrBOc1mJ48ePs/vu7Oys0+8AMIxiiI+PZ+6TNWvWpOmc8PBwFChQAI6OjmluaKMdgN6zZ0+G1podiIyMZIrSz88vSabQP//8wzLY2rRpw9xNAQEBCAgISHHu5KwJjUaDRYsWQSaTwcHBweB9p38F/BJK4f79+6hQoQLatGmDsWPHYtOmTbh06RI8PDx0hKxCoUCtWrUwePBg1K5dmxWnWVhYMK5+bWHr6OioI0j1KQV9gloul+u4ekxMTFJMS3VxcWGpf35+fhgxYoSO28vV1ZX9SLR37nny5NFZMz/Mzc1RtWpV1KtXD+7u7kwx8tdq27YtSpUqlaT9aIECBdCrVy/07t0bxsbGsLS0RHBwMGbMmJHsey1VqhQLBvKjYMGCaNWqFby9vUGUQNg3cuTILOf+uXfvHiMTVCgUSbJ/DKEYYmJiULNmTYhEIp3WnCnhypUrkMlkqFu3bpqtpxcvXqBUqVIQiUQIDAzM1e6SP/74A6ampjAzM0vCNhsdHY0JEyaw5IeVK1emicAwNRfT1atXUbBgQXAcl8RSSTxPagroV8MvoRQS49q1a0yw8TtzsVgMpVKJQoUKsZ1tu3btcPr0afaDU6vVOHLkCEQiEQoXLpyEBsPDw0OnU1ti9lI+RqFPeHp6emLChAnYtm0bmjRpkoTHKHGdRNOmTRlhHT+cnJwY/bD2sW3atMGoUaNQu3ZtvUpCJBKxHH/eOrh9+zZiY2Nx5coV9O7dmyk97ffDP/b19UVQUFCSYjrtY3hFyHEcW7e1tTV69OiBBg0asKrqVq1a4dy5cxn4KqcN796904nZjB07VkcIGUIxREREoFy5cpDJZGlu+MPTuy9YsCBd12nUqBGIEuoecnMx18OHD1k72V69eiUpatMORPv7+yex5BIjJRcTj/LlyzNKF19fXzx8+FDvPIJS0MUvoxSioqKwfv16lC1blgmpdu3aMVoBbTeIhYUF/P394evri0KFCiFPnjwwNTVNspPnaxL0CXkLCwtUr14d06ZNw/r165O0weQFJi8otZ83MjKCt7c36tevr2ONEBHry5DYjaTtkkqOGM/d3R0TJkxASEgIZs2ahQoVKjBrQN+c/fr1Y4Hh2NhYLFu2jP3I+BoPbXI6bYuDL/DTfn/aCjFPnjwoX748Uw7Dhw9H3759Ga14mTJlsG3btiyplo6OjmbZPEQJ/P3a6byGUAwfP36Et7c3lEplmnpNazQaNGjQAFKpNF00Dmq1mtFO1KpVK1sD+elFTEwMhg4dCqKEjLXEnf60A9EymQwzZ85MVtGl1cUUEBCA7du3s9a3mzdvTnJMasrlV8MvoRR27NjBBJNCoUCePHng7u6exLcvk8mYUFIqlShfvjyaNWuGzp07o3///hg7diyGDx+uw3JKpNsfwdbWFubm5sm6UziOQ+HChVG2bFmUKlUqyW5fpVIlEerasQ5t/79EIoG/vz9br77rJVY2/OOKFStizZo1eP36Nfbv34+BAwfqUHtrr7dy5cpYtGgRHj9+jIiICEybNo0pyU6dOmHBggVQKpVQKpUICAjQWS+vdFQqFbMktF+vVq0aUw5WVlaYOHEi5syZw9ZSoEABrF+/PtU+zemFWq1mdBIcx8HDw0OnV7UhFMObN2+Y9ZiWoPD79+/h5OQEDw+PdBcBrly5EmKxGEWKFMlU3UV2IDg4GBYWFjA1NWWMq9oICwtj/TB8fX1TtRqSg7bAf/LkCfutdOjQgd3ftGQ6/Wr4JZTCli1bUhSW/A62WLFimDFjBmbPng1ra2sYGxtjw4YNABLM3w4dOjBrgbcSeJO0Z8+eTOk8fvwYGo0GmzZtYm6ZxC1A8+TJk6ogl0qlMDU11Uuroe2a4tc/ePBgLFu2DAEBASnWPBQsWJBVRRsZGaFNmzY4cOAA4uPj8eDBAwQGBiZRfPzw9vbGokWLcO/ePQwZMgQymQxyuRxdunRhAfFBgwaxeIGxsbHO+5RKpXrjL35+fsytY2VlhenTp2PDhg0s48nNzQ3Lly83eI7+unXrIBaLIRaLYW5urtPwyBCK4fHjx3BycoKdnV2aelKcOHECIpEI7du3T/e1Dh8+DDMzM9jZ2bHeHLkVjx8/Zg18BgwYkKQiXKPR4I8//kiT1ZAcEgt87da0Hh4euHDhgqAU9OCXUAoRERFYu3Yt9u3bhx49erCexs2aNcO0adPQu3dvWFtbJ2lWw6drOjg4QCKRQKFQYPDgwXj06BH++ecfEBGWL18OICHTIl++fOA4Dr6+vizrws3NDQcOHEBkZCT+/vtvnYCwkZERGjRogOnTp2P16tUYMGAAlEolJBJJqjUHyQ0/Pz/s2rULMTEx2LBhQ5Iubdrvr1GjRujcuTOLg+TJkwdjx45lxWbv379ngl7f+TVr1sSaNWuYslSpVEyZBAQEMNbMPHny4J9//kGzZs10AuTa7iTeenBxcWF9sm1tbbF582YEBwczf7STkxMWLlxo0IylI0eOQKFQMGbaFStWsNcMoRhu377NUo7TwqzK1yOkNVCd+Fp58+aFQqFIdwvR7EZMTAxjP/X19dVr4WTGakhO4J88eVLHutceKfUt+VXwQygFInInojVEtCOt52grhfDwcNStW5f1KujWrVuSH2doaCikUikaN26M33//HR06dEjiLzczM4O/vz+GDx+OM2fOoGTJkvDw8GB+70uXLukIusGDB+stvR8wYABznfCWhKmpKTp06ICVK1fCxsaGpU2ePXsW48aNg5+fX7LxC33D0dERs2bNwocPH5J0ctMn3KdMmYJatWoxxVi1alVs3bqV8eXzx5uYmOjQa/PWQJ06dVj1s7m5OWQyGRwdHVkMx87ODh8+fEBMTAwCAwOZtZDYouF/rMbGxqwIrmbNmnj48CEOHz7MAum2traYNWuWwbiWjh07BoVCwTYCAwYMYDtTbcWg7WJKDy5cuACVSgUvLy+8f/8+xWPj4uJQvHhxODo6Zuj9hYWFsfs+Y8aMXJ2ZBCS4d01NTWFhYaGX4ym9VkNa6DI+fPigQ/Qo9IT+P3JMKRDRWiJ6S0Q3Ez1fi4juEVEoEY1M9FqGlAJPBCcSieDp6YlKlSqhZs2aqFKlCvz8/FCwYMEkCiA9w9bWFpUqVUL16tWZUhCLxcn6kWNjY+Hv7w9TU1Pcv38fBw8eROfOnZmg5F1GTk5OSb6sU6ZMAVFCplPi+gJ9QyQSoVKlSli6dKnezCft4Lm5uTlatmyJ1q1bs1gHz8Y6depUSCQS1uNh9erVCAwMRLFixXQsLDMzM8bTxO+8+Z4HlpaWjIDu69evGDNmDKPYSHz/jYyMWPBbIpFALpdj9uzZiI2NxYkTJ1CjRg0QJbia5s+fb5Aex0eOHIFcLmeKunbt2iwAfevWLdjY2CBPnjwZpu04duwY5HI5SpcunaqwP3fuHDiOw8CBAzN0rejoaNZFr2vXrgaPyRgaoaGhrN0q3w86Md68ecMEua+vr15OK22k5hrSaDSMCkXoCf1/5KRSqEhEJbWVAhGJiejhd8tARkTXiMhL6/UMKYWdO3dmSNjL5XLY2dmhePHiqFOnDtsh+/j4oE6dOjr9j5Ojy8iXLx9atmyJOXPm4NixYyw75MmTJzAzM4Ofnx/b9Xz79g3BwcFo06YNm08ikaB///6MPlmtVqNGjRowMjLC9evXcfjwYfTq1StJL2Z9ge7EsQltJcE/5nfucrkcZcqUgY+PD3Pt8D2s+WPmzZsHjUaDz58/Y9WqVSxvnp9T2wrglYxKpdJJDXz27BnatWvHXhOLxZBIJGz9iQv08uXLx/zl58+fR/Xq1UGU4HZav359prOVDh48CLlcDhcXF0gkEnh5ebH1Xr16FRYWFnBzc8OzZ88yNH9wcDDEYjGqVKmSanykZ8+eEIlEybYbTQ0ajQZjx44FEaFGjRq5OjMJgE5xZPny5fUy2GpbDUqlEuvXr092vrTSZWgT67Vu3fqX7wmdY0oh4drklkgplCWig1r/jyKiUVr/Z0gpHD16NIlwrFu3LiZMmIDNmzdj8eLFrIiL71ZmZ2eXYiBYKpWibt26TChpj/z58zNB5uPjk8TdwlcY803dx44dm+SDiYyMxJgxY3SEe7FixTBv3jxcuXIFlpaW8PHxYTuq+Ph47N69O0kQ187OLln/qfbQPk8kEsHHx0enyQ/fhjKx0unevbvOri4mJgbTp09nrh99SkmpVCYhRLtw4QJzDUmlUojFYgQEBOgoMu17Ubt2bUZXcfjwYUblXaRIEQQHB2fKZRISEgKZTAZPT0+Ym5vD2toaJ0+eZOs0NTVF/vz59TZsSgt+//13EBGaNm2aohL7+PEjbG1t4evrmyllt3r1aojFYnh7e2d5gaAhsGXLFiiVStjZ2eH48eN6j3n58iUT+p06ddLrpk1LuimfthoXF4cpU6YIPaGR+5RCMyJarfV/eyJaTERWRLScEqyIUSnM152ILhLRRRcXF/YmlyxZwgTf0qVLYWtrCysrK5w9exarVq2CUqmEpaUldu3ahdjYWHh6eqJAgQKIjY3F58+fcebMGSxbtgxdunRB8eLFk6089vLygr+/v05NglgsxsqVK/Hs2TPs378fU6dORaNGjVi+Py/o+vTpo9ctwWdOeXt7s2wNkUjEsnJGjBihc/z79++TtPjkBW3BggWTrWHgj+GVIn+/pk2bhmnTpul0l0usLAsUKKC32cnhw4d12mUmVhCJ+xVrNBrs2rWLKWiihL68fKe3xHPw8aFv375BrVZj27Zt7L37+/vj33//TfnbnwKCg4MhlUpRrFgx5M+fH1KpFOvWrQMAnD59GsbGxvDy8sowv868efOYUk1JgW3cuBFEhGXLlmXoOjwOHjwIlUoFJyenH4JJ9ObNm/D09IRYLMbs2bP13qO4uDiMGzcOHMfBy8sryfcpPbUMPLR7Qg8YMABRUVG/XIFbblMKzfUohaCMzK1tKURGRjIfNH03TZ2dnZnQrFKlio7vPjg4GESEoKAgvTctMjISc+bMSRIklUgkOoJY+3VTU1N0794dZ86cgUajgUajweXLlzFy5Eid3bC3tzfGjx+PK1eusB/CwoULQZTgG759+zbGjBmjY33UrVsXp06dYse/f/+eFb2VK1cOxsbGEIlETLAaGxunyXrgh729PY4fP443b95g9erVaNCggV7lUqdOHZw8eTLJD/jChQs6zW+0d/5Tp05N0lc5JiYG8+bNY24rHx8fPH36FNOmTWMWTeJaiBEjRiA6OpoV2fEd7xo0aJBEWKQVu3btgkQiga+vLypXrgwiwujRo6FWq3H8+HEYGRmhWLFi6WrJqQ2+jWhKzKoajQaVK1eGubl5sr2d04pr167ByckJKpUKBw4cyNRc2YGvX7+iWbNmICI0atQo2X4hhw8fZq1B165dy75/GRXmkZGR6NevH4gIhQoVQsmSJQWlkINKIUX3URrnTNJPYc2aNSAiVKhQgQkL7d1mYr+kRqNBlSpVYGVlpbeXb2xsLDp06ACihACVSCSCVCqFm5sbc5twHKeTLcTnwhMR8ubNi6lTpzJT/sKFC5BIJPD29kbFihWZJeLq6oqRI0ciNDSU+YZHjRoFICG2EBISAhMTE7aDz58/P2bMmIGXL1/i7du3LJNm9uzZ6NixI4j+HxRPLeVVLBYniVOYmppi8ODBOH/+PL5+/Ypdu3YxiyWxEpk9e3YSWuTr16+z9ELtoVAo0KNHD1y6dEnn+NevX7P0VIVCgVWrVuHz58+YMWMGW5t2LEcqlaJjx474/PlzkiK7Xr16pZlwThvbt2+HWCxm/SmICK1atUJ0dDQOHjwImUyG0qVLp6nBUWJoNBr89ttvICIEBgYme9ydO3cglUrRrl27dF8jMV68eIFixYpBLBYn27I0N0Gj0WDBggWQSCTw8PBINnnj1atXTHG3b98e4eHhmd7hHzp0CI6OjuA4Dm5ubrmaRsSQyG1KQUJEj4goL/0/0Fw4I3NrWwofP35k+fqJs1z4Xfrs2bN1bsyVK1fAcRyGDh2q83xkZCSjl548eTI+fvzI4gocx+Hhw4e4cuUKxo8fzwQa/1piQcxXC+/ZswczZ84EEWHVqlV4+/Yt1q5di7p16zJFUrlyZVStWhVEhPnz57P1nDlzBiKRCGXKlGE+eZFIhLp16yIwMJAFbXfv3o2TJ0+yuoOaNWti4cKFOm01kxs2NjZJgtSmpqZo1KgRFi9ejClTpkAkEiVxq0kkEtSuXTtJ/ODOnTuMHVNbCREl0B8sX76cBUW1KamJCKVLl8Z///2H8PBwzJo1i7m7tOMNIpEI1atXR2hoKN6/f49+/fqx4rQFCxakOxPnjz/+YFlcU6dOBVGCe+rdu3cIDg6GRCJBuXLlklg8aUF8fDzLqNm4cWOyx/GbAu3iuoziy5cvqFmzJrNScnvKKgCcOnUKDg4OMDIyYgWliREfH4+JEyeyhjuG6G394cMHJjN8fHwybHUCPw7PUo4pBSLaSkSviSiOiF4QUdfvz9chovuUED8Yk4F5k1gKQEJKIe8+MTIywpYtW3QCp0QJLRW1ydE6deoEmUzGfP2RkZHw9/cHx3E6Pt64uDi0b98eRAn8QtpBrxs3biS5jqOjI4spaAvZPHnyQKFQ6BRJvXz5ElOnTmVBX144z5o1ix3Dp6muX78e9+/fx+jRo5n7iLcWxGIxTpw4gdjYWMyfPx8mJiaQy+Xo1asXGjdunKpi0De0XTh84xReMZQoUUJnF29ubo6GDRtiwYIFCAoKwvr163Ua4NB3C8rT0xNECS6u3377je0Mg4KC2DU5jkPXrl0RFhaGL1++YMKECclyPxUtWhQHDx7ErVu3mAvR09MTISEhafl9MPANcdq1a4dt27ZBLpfDw8MD9+/fx44dOyAWi1GpUqUMFdVFR0ejcuXKkEgkya4rKioK7u7u8PT0NEhVd2xsLLNS2rRpk2u7uWnj9evXrBamZ8+eya756NGjsLe3h0KhwOrVqzOt9AICAmBtbQ1ra+sMV1fz8whKIQdH4uwjvkpYuzfBwYMH2S6NF2YFCxZkroAXL15AqVSiRYsW0Gg0aNmyJTiO09soXKPRsKrbokWL6sQobt++zXbCvXr1YgFjogQXUWKlIZfLsWjRIp0vnlqtxsGDB9G4cWO2K/by8sLvv//OTGVjY2NWXBUfH499+/ahefPm7HiRSISBAwdi6dKlqF+/vk6dA08IyAvUxCmtyXE5pfX15EbiVF6RSISuXbuiY8eObC2VKlXCzp07sXXrVkgkEqaAzMzM2H16+/Yt+vbtyz5HOzs7HRZbZ2dnbNy4EX///TfLTa9Vq1aque7a4K2EqVOn4vTp07C2toalpSVOnTqFzZs3g+M41KxZM0MC9suXLyhZsiSMjIySDZDv378fRIQpU6ake3590Gg0mD59OogSKtB/hFTMuLg4jBgxgu3ck+N5evPmDaOb79ChQ6b6N/NZTmFhYUxelClTJt0V7oJSyEVK4fTp0zpBX97Xz3EcJk6ciGXLlkEulzPBrVAo0LFjR/Ts2ZNRafNCk8+bl8lkUCqVMDU1hZWVFezs7HSEu0wmQ/ny5dGmTRv069eP+TutrKywb98+hISEYPTo0aw7Fe8z1Y5DSCQSVK1aFRcuXND50EJDQ5mvkyiBZbR///4wNzdHqVKldHhkNBoN1q5dm6SBkLGxMRo0aIBhw4axbKmSJUuCiJjQ9PX1TeJus7GxSTaDSdtS0B58j4ZatWrpuNBkMhmKFi3K7oH24OnBR4wYARcXF6ZAu3fvDhMTEzg6OjKCs5IlS+K///4DkFD3wCtnjuNQr1499hny92revHmYOXMmzMzMIBaL0b9//zQFizUaDaup2L59O0JDQ1GgQAHIZDJs2bKFxa7q16+fhMsnLQgLC0OBAgVgbm6erOujefPmkMvlSRhGM4MtW7aw/uU50U87I9i1axdMTU1haWmJ/fv36z0mPj4eEyZMYCSUGaUp0a530Gg02Lp1KywtLSGXyzF37lzEx8dnKNMpo8dkNX46pUDJuI+ePHmCS5cuYeHChahTp06aWEWVSiXrgEaUUDg1atQojB07FiNHjsTQoUMxYMAA9OnTBz169EDXrl2RJ08e1maTz3nWRwCnLUgdHR3h5OTEhC0vWLUFuZ2dHcaNG8esGJ6BU6VSoXLlyqzyl98ZaTQa7Nu3j5HM8cKYrw7mr1G5cmVs3rwZU6ZMgUKhYDt3XlA3aNAAq1evTpJpxXEcE6r6lID2/XVycmLn29jYoFKlSklScr28vJJVNlWrVsWoUaNYzEQul0OhUMDc3BwzZsyAg4MDOI5Dnz592P3ZuHEjs4SUSiV69+7NSPr4jcHgwYPRsWNHiEQiWFlZpanBS3R0NMqVKwcjIyNcuHABHz58QIUKFZgFsXjxYhARmjVrliEXw5MnT+Dk5AQHBwe9ab4vXryASqVCrVq1DBoLOHHiBCwsLGBra8sUbG7H/fv34e3tzTZ3yX12Bw8ehLW1NUxMTLB169Y0z58SXcbr16/RoEEDECVk+JUuXTpNAj+tdRM5iZ9OKfBD21L4+++/IRKJdASvtkDiBY0+F4iLiwsTfN26dcOKFSuwdetWHDlyBNeuXcPLly91doUnT54EEWH8+PFwcXGBSqXCP//8g9jYWBw7doxRNyxcuBCLFy/GyJEj0a5dO1SuXBkeHh56GVETj3z58mH69Ok4ffo07O3t4ezsjFOnTqFv375M+PJC2dXVFWvXrkVsbCy6d+8OjuNgamrKMpv4WIWtrS26d+/Odt9ExILBvHLSl8bKcRyKFCmiU+imvdvnH4tEIvj6+qJatWpsHlNTU0YpkdiSSdyciH8v/fv3R5s2bdjcYrEY06dPZ64je3t7bN26FRqNhmXa8PPkyZMHo0aN0qmdkEqlaN68OVOevr6+qfYzCAsLg5ubGxwcHPD8+XN8+/YNbdu2BRGhS5cumD17NogSfPUZKTq7efMmzM3N4enpqdeC4VOU//zzz3TPnRLu3LkDNzc3KJVK7N6926BzZxUiIyNZJmCtWrWS5ZV6/vw5+2737t07XS6+5CqjNRoNNm7cCHNzc4hEIuTLly/FzzutFdaCUsgGpXDjxg0mRBwcHNCnTx+MHTsWAwcO1PHvZ2aYmZkhX758KFOmDMzNzWFubo7p06fD1dUVUqmU8caPHDkSRAmxC307G41Ggw8fPmDr1q2MqVVfkxptS0MkEsHS0hI9evTQyXjiBfuyZcsQGRmJL1++wNnZGW5ubjAxMYGXlxfCwsKwf/9+NGzYkM3v7u6exA1UuHBhfPz4Edu3b9eJA/AWA1FCPMLZ2VmvchWLxUwZlCpVCqNHj0bLli3ZXPzfxNaCtrtJW4E0bdoUbdu2Zev08PDA1KlTmQusRo0aePDgAb59+8bSYPnge9GiRTFjxgzGtcO/j7Jly8LS0hIcx6F3794p+thv3LgBlUqFEiVKICIiAhqNhnFsValShT3u0qVLmtpLJsbJkychk8ng7++fJLU3Li4OJUqUgIODg8GpK968eYPSpUtDJBIlW6eT26DRaLBixQrIZDK4uLgka+nExsayJj8+Pj5pdpWlJsxfvnzJUqT9/PySZXIVlEIODdLjPrp58yZKlSqVROAoFAoULFgQxYsXZwLRyMgITZs21fFDaw9fX19MnDgR06dPx8yZMzFt2jRMmDAB/fr1Q6tWrVCtWrUkjXP44e7ujg4dOjCq6pEjR6boAuCrXoOCgjB37lwm8KVSaZqotVUqFavLsLOzw/r167Fv3z62i5XL5fDx8cGXL18QFxeH5cuX66xduzGRWCxGiRIl8PbtWzx69AhFihTRuZa5uTlq1KihoyD0rSl//vwsANyqVSs8efIE69atQ5UqVdIUzFYqlahSpQpTMEWKFGGCnCiBqrx58+ZQqVSQy+WYPHkyIiIi0KVLF7ab5GMolSpVQlBQEFMOvBWXP39+cBwHGxsbrFu3Llmhvm/fPohEIjRq1Igd8/vvv7PqcZ4WulevXhly9fz555/gOA5NmjRJsgM9f/48iAjjxo1L97ypISIiAg0bNgRRQm+MjCi1nMCFCxfg6uoKmUyG5cuXJ3vP//77b5iZmcHc3FwvK2tipMXtU7FiRVhZWcHKygpSqRSTJk1iHoS0sLam51pZjZ9OKfBD21LYtWsXiBJcKnzmiiGsA+1haWkJLy8vVKlSBa1atYJKpUKBAgUQGBiIiRMn6hDCJT6vRo0aGD16NP766y88ffqUfZnj4+Ph7+8Pc3NzvHz5EhqNBufPn0fPnj2Z8E1sPXAcBx8fH1aqn1jQmpqawsPDAyKRiHG9eHt7M1eRu7s75s+fjz/++IO1KuUH35zkxYsXiI6OZo2FtBVAhQoVsGDBAuZnT6xgtC0A3qV35MgRAAkm/qxZs5A3b94k9zexe8nPzw+TJ0+Gk5MTs1bkcjlTnBYWFkxxeXp64vjx4xg8eDCICG3btkVgYCALordo0QIrVqxg5/LWB7/78/f3T5YaIjAwEES6dCMnTpyApaUlrK2tmVtpwIABGVIMCxYsAFEC3Ufi81u0aAFjY2OEhYWle97UEB8fz5RakyZNklgruRXv379HrVq1QJQQW0suRfjhw4fMqhwyZEiK7qS0BoiJEjKUeHbaIkWK6DQ7SqulkNoxWY1fQikcO3aMNYbnBRS/0/bw8MDmzZtx7do1HDt2TGfX2aRJE4SEhDAlUr9+fR2hNmjQIKxYsQJTpkxBnz590LRpU/j7+yNv3rxJhKClpSUrtPLz82MZMg4ODihevLiOorK2tkbNmjUxZswYLF++HAqFAnXr1tURConXmpwCEIlEWL16NVasWIHKlSsnEa6820YqlWLatGlJfvyPHz9OUuCmVCoZQdzGjRthZGQEU1NTphykUilWrFiBW7duoW/fvskyyGqvuUuXLjr011euXGHFevzgWzjy/4vFYvTo0QN9+/YFx3EQi8WQy+VYuHAhKzCUy+XsnG7dumH06NEgIjRs2BBhYWGMV0kqlaJ///5YuXIlU0ra/aY5jkP//v2TuGs0Gg1Tjjw3EpAQBOXpzXkhNWLEiAwphkGDBoGIMHfuXJ3n7969y9KMswoLFiwAx3Hw8/PLMM9TdkOtVmPSpEngOA5FixZNtgdGdHQ0evfuDaIEssnERZY80qMUeAQHB8PJyQkikQiDBw9GZGSkoBRyapAe99Hnz5/h6+uLs2fP4sKFCxgxYoQO6ZpUKsWAAQMQFhaGrVu36gjMY8eO4eLFiywmUbhwYR3hW6ZMGb0+zFevXkGhUKBGjRqYN28eevTokSy9hJubG1q3bo2+ffuid+/eaNmyJby9vXWUGFGCr3rz5s2YMGECJBIJ8ubNi0uXLuHBgwcYNmxYEloKPvuG4zi0b98ely9fRmRkJMaNG6djYWjXKygUCgQEBGDMmDHYt28fPn36pEPHoK1UfHx8sGPHDly5cgUFCxZkXef49VaqVAnh4eGIjo7GzJkzU+VbUqlUWLdunY7gPHDgQJLjnJycdNavUqkwbNgwZhmIxWL88ccfuHHjBjp06ACxWAyRSMRcQrwrqUqVKvj69StevXqF7t27QyQSwczMDDNnzsSSJUtYhlRihb19+3adNcbGxqJq1aqQSqU4ceIEe/79+/fMYuJjVxMnTkzL71IHarUazZs3BxElyZ7p0qULZDIZnj59mu5504q//voLCoUCefPmTVddR07jwIEDsLKygkqlSrGL3Z49e2Brawu5XI7AwMAk7rKUlEJKrqHPnz8zJmR3d3d4e3sn6xpKj4spq/HTKQV+6Guyw5vh4eHhjJCue/fuTOBzHAeFQgEHBwc4OjqyeoSQkBBcu3ZNR8jyRXD8LrJLly5JCMt69OgBuVyus8PSaDRst8rXA0il0iRtMy0sLODv749mzZqhadOmSQSqkZERmjdvjhUrVuD+/fvQaDT49u0by4ZIbkdesmRJlpHDvx8zMzNG1VGiRAn4+PjoKCRvb2/07NmTpYTyVg4/bGxsMHDgQLY7DwgIYLEMY2NjRr6mvYPT3u0nXqu9vT02bdrE6Ch4qunEllBiN6ClpSVq167NjuvZsydiY2Px+PFj9OzZU+d4Ly8vlhHFZ6zcvHkTtWvXZop6w4YNmDdvHlO22gqxXLlyOimjHz9+hKenJ6ysrHQKqr59+8aq3flYxowZM1L6TepFdHQ0KlSoAJlMhmPHjrHnnz59CplMhq5du6Z7zvTg3LlzsLW1hZmZGXP3/Qh4+vQpaxE7aNCgZClOwsLCmCegevXqOsWnGbEUtHH8+HHG3mtvb5+ixSVYCtmkFL58+aJDhJcnTx4dlsjw8HAd15BYLEbLli2RJ08e5nbatm0brly5whSISCTC5MmTIRKJYG1tDalUClNTU8ybN48FmG7fvg0i/RWofFoh71IaOHAgPn78iBMnTmDx4sXo0aMHypYtqzdgmz9/fpQtW5Zl4xAlUET07NkTu3btwvPnz1GwYEFIpVKd8xPvrocOHQoLCwtWMMcX6M2YMQMRERE4evQoJk+ejBo1aiRpBapdv8G3OSUi9uUvXbo0OnfuzAS0dt7+vXv3knSN0xdktrKyYv2Y+Z2UQqHQa21pC3zeHcQLYt6v+/z5c/Tr148dy5MUenl56TR0OXToEKtpKFOmDA4ePIjJkyfrUIbw34H+/fszQXP//n2oVCqULVtWR/hoNBpGRcLHMRYuXKj/F5kCPnz4gEKFCsHMzEyHg2fAgAEQi8W4e/duuudMDx4/fozChQtDIpFg1apVWXotQyImJgZ9+/YFUUJ8KLnWm3wWk1KphIWFBetxnVmlACTQlPCZeebm5ggKCtJbxyIohWxSCvfu3UPhwoWZYuCFQvv27dkucdy4cUw42draQiaTQSqVwtzcnGWlrFmzhmXvECUUVR0+fBhmZmawtLRk3Zs8PT2ZG6FWrVqwt7fXW+HK907mq6T1pcjFxsZiwoQJrCCOF0bau+oSJUqgSJEiTAHwdM98tTVPA8wPuVzOMmx4AVetWjVmuehTZHFxcbh8+TLmz5+fxKrhh4uLi46VYm5ujjVr1jDlZWFhgVOnTgFICMZqKwNeIekbpqammDJlClq0aAGiBBeOp6cnK9jTvh+8m0j7fJFIhCVLlrD38ubNG3Tv3l1HSdra2uoohvj4eKxbt46lsTZp0gTnz5/H8OHDoVAodAgATU1NGUkb734cOXJkks9yy5YtkMvlTKmtWbMmyTGp4cmTJ3BwcECePHmYcAsLC4OxsTFatmyZ7vnSi8+fPzMyvWHDhv0wmUlAwmdjbGwMW1vbFIkF7927x9x9HTt2TFPab1oL03x8fFiszNvbm8Xm0jNPVuOnUwqkJ6bw6tUr1jmMz1aRSqUQiUSwsbHBzJkzIRaL0b59e+zZswdKpRJOTk5o0aIFI2DjhceCBQtYcRJRAt3BnTt34O7uDplMhqFDh8Ld3R0cx2H48OHYs2cPiChZZsfff/+dzV+4cGGd1yIjIxnldbVq1fD69WtUrFgRKpUKf/31F+bPn4+WLVvqCGkTExN4eHgkEdy8QrG3t2fXk8vlOkFuT09PtGjRgv1fq1YtncyW6OhoXLx4EZs2bWLEdVZWVqmmkspkMh33l4+PD0JCQjBmzBimPHhhX6VKlWSbGCkUCva+ChYsyAju+FhC4nUktm7q1Kmjo5zfvXvHUi/5dSbuuBUREYHJkyfD2NiYtUa9ceMG+vTpk6R/Rt68eXH69GnG6qqvZwHPmcR/r/TxaKWGK1euQKVSoWjRoqyCm2dRvXLlSrrnSy/i4uIYmWHjxo0zRAKYU7h9+zYKFSoEkUiEadOmJavUYmNjMX78eIhEIri5uSUR3omRHgoLjUaDHTt2sIzEtm3bsg2JUKeQTZaCRqNh3Pqenp5wcHBgwUf+g5HJZCzt8OLFi7C3t4eZmRk2b97M0gr5MWPGDBb4k0gkePjwId69e8eCiqNGjWKBWW9vb7i7u6NEiRLJZp788ccfbG5+R/vu3TsUL14cHMdhwoQJLE89NDQUSqVSh+ZAo9Hg0aNHWL9+Pbp27aqz4+eHtsD08vJCyZIlddwoPP0FUYJvnq+IViqVGDhwoE4arLbFQZTgX8+TJw84jmPWir7K7MQxAJVKBVtbWygUClZlzHEcBg0ahA4dOugoB31B6rJly6Jbt25M4GtnSWkrBO15jI2NsWTJEh3l8PDhQ53EA39//yR9kV+/fq0TjJ49ezbu3r2LTp06Jekj7ePjAw8PD9jY2OjtM/zgwQPky5ePnbd3716934uUcOjQIUgkElSvXh2xsbH49OkTLCwsUKdOnXTPlRHwfQ74FOiMtibNCYSHh6NNmzbse5MS79WZM2fYJm/gwIHJKsCM8BpFRkZi7NixkMvlMDExwezZs1GhQgVBKWTV0FYKvNC1srKCWCxG2bJlmTDS3iWamJhg06ZNABLM9MKFC0MqleL333/H5cuXdXz4gYGBzLVgb2+P9+/f49u3b2xn36lTJ+zatQs2NjZMGCbXbxb4PzW0SCTCtWvXUKxYMSgUCuzbty/JsYsWLWLuh9jYWNy7dw979+7FggUL0KdPH1StWjVN3E76dtcqlSpZpVK+fHn8+eefuH79OiIiIhAVFYWqVatCJBJhzZo1KFeuHIj+T6zn7u6Otm3bptgCVFsxaVNlFC1aFCdPnkSjRo10LIXEhYFubm4YMGAAiBKsqdWrVzMFolKpdHpFa79XW1tbzJ07V8c1wM/Dj7p16yZRDjdv3mT1G66urtiyZQtu3rzJNgmJ+zr4+PjopT549+4dC4CKxeIM9UlYu3YtiBJSbTUaDevJwbvosgPBwcEwNjaGs7PzD9Hmk4dGo8GSJUsglUrh6uqaIrVJeHg4S1318PDQe38zQ3YXGhrKYppGRkYoWrRoet+OQfFLKIWnT58mKaaqUKGCzs63SJEiLLjYs2dPxMTE4NOnT6hSpQqIEriMHjx4oCO4+E5cRAlZO5GRkdBoNCwo2qRJEzx79oxltNjY2KTYOJ331YpEIsjlchw6dCjJMWq1GhcuXIC7uzvEYnESgatSqZjPunLlyvj9999x9uxZxu6pLSR5gZQW5aG90y5evDgqV64MPz8/5M+fP9ViQN4qI6Ik7jjtnbz2e+F99pMmTcLt27d1PitTU9MkXEulS5dmhXgPHjxgBUS8ktLOHNKOB5iZmWHEiBFsRz9r1qwkSqRRo0ZJ3DJHjhxh1k3p0qVx4sQJXLp0iWVgaY/8+fMnYboFEoKPvPtKIpFkqKc074KbMWMGIiMjYW9vjwoVKmRr45zLly/DyckJJiYmejcxuRnnz5+Hs7Mz5HI5Vq1aleJ9++eff+Dm5qbXajCE2yckJIT9BmrUqJEtrkB9+CWUwr179+Dn56eX0yelIZFIWOeytJ4jlUqhUqlYGqOdnR169uzJKoaVSiVmz56Njx8/JvkCvn37lq3P2toaz58/h0ajwb1797B06VI0bdpUpxaB4zi4uLhg3bp1OH36NE6ePAk3NzcoFAqWOcFDo9Ewk5mnlOjbty+qVKmS4QpviUQCBwcHlCxZElZWVhCJROjbty8TlnXq1GFsq3Z2diwekDdvXri6uuq9r9bW1knour29vXHlyhW2s+aHvqI4fk03b97E0aNHWcaQSqXSSz/CM6zKZDJ0794doaGhLGWYj3XwCiWxclCr1fj999/Z+2rYsCHu3r2Ly5cv6+1JXadOHZw7d07nc1Gr1cw/L5FIcObMmRR/sImhVquZAty2bRuWLFkCIv3xjKzEixcvUKJECYjFYqxYsSJbr51ZvHv3jsWnOnfunGL1dnJWg6FiARUqVEC+fPlYYWq7du2S7RmRVfgllMLFixchlUpRunRpNG7cOMluOfHQTju1tbWFo6Mj232npiA8PT3Ru3dvtGnThgVjk9uNm5qaolSpUmjZsiXra6CttFQqlU6HNmdnZ3Tq1AkbN27Ey5cvGU3zqlWrcODAAZiamsLe3j5ZQrDo6GiUL18ecrkcFhYWcHJyQtOmTVN072gPc3NzlmHF3wt9WUM1atRA+fLlQZTQOnTLli2wsLCAkZERq3Vo0aIFwsPD8e+//+oQ0yU3RCIRBg0axHbiiYPIiT8XqVSKrVu34sOHDzqU2XwGmkgkYgpDIpGgWLFijFK8VatWzG1VvXp1dizvkmvcuLGOqyQyMhLTpk2DiYkJxGIxunXrhmfPnuH69eusj4b2GmvUqJFE+PMpqxKJJN3un+joaPj7+0Mul+P48eNwc3NDyZIlsz0zKDw8nFnFI0eO/KEyk+Lj41k9U7FixVLtV/HPP/8w5oIBAwagfPnyBlEKvHL59OkTRowYAYVCAZlMhiFDhqSp54ch8NMpBdKTfcRXFSY3+B+7WCyGhYUFpFIpOnXqBDs7O4hEItSvXx+jRo1CtWrVWMqqvsYw2qNw4cLo2bMny4svUqQIatSoAaVSiYYNG7KAZbly5fRy/fBDLBajZs2aWLNmDeM/4qFWq1GpUiVGWOft7a23slWj0eDWrVtYvnw5mjZtmsQykEqlaNiwIfbv34+HDx+iRIkSMDIySlZYu7q6onTp0kzZ8eePGTNGx83DK9LRo0fj5cuXjO6Br2UoVqwYS8P99OlTmhhr+R4VRAlFdHzxIC90ec4nfvj7++Pu3bvMTZPYutCOn0ilUhQvXpytm2+aNHPmTKZYChcuzBRSkyZNdBrJh4WFoV+/fpDJZJDJZOjfvz/evHmDHTt26Cgt7ewubW6c5cuXs888cRZUanj37h08PDxgbW2NOXPmgIiSWIvZgbi4OPZ7a9Wq1Q/R5lMbISEhsLCwgJmZWar04eHh4ejTpw/7XhUrVizT109scTx79owlM5ibm2P27Nk6dDBZgZ9OKfAjcUwh8W6Y4zid1papDZ50LaV8en3naP/PWyD9+vXDsWPHYGNjA1NTU+Za4gWWdge3xHNYWVmhcuXKGDBgAFavXs0afdja2uoETTUaDc6ePYsBAwboWBt2dnYoUaKEzrzaFbJAAq++TCZDixYt8O7dO6xatSoJJTe/Njs7OybkpFIp2rdvD3d3dygUCtajgCgh6Lxu3TrMmDEDCoUCcrmcFde1bt0agwcPRs+ePXV29Ynfv77PysLCgt0/fqhUKmal8ed17NiRJQEkDqTb2NhALBbD1taWxTy8vb2Z+4jjOEybNg2TJ0+GTCaDhYUFGjduzCyNpk2b6iiHJ0+eoGvXrhCLxVAqlRg1ahTmzp3LrqWthIgSAtp8oHPr1q0sK2nXrl1p/R0DSCies7KyQv78+eHp6YmCBQtmqNFPZqEd9K5QoUK27XANBW2+r1GjRqV6D3mrgYjw22+/ZUlb02vXrjErzNnZGevXr89Qr4604JdQCo8ePULDhg31ZuRYWloylw0vPLV3k/369cO4cePAcRyKFSuGhw8fYuXKlSBK6Erm6uqapDqXKME9UatWLRQuXDhZmuvEDX8aN26MEydOQKPRYMaMGSAiFC9eHGZmZrCzs8OECRPQrVs3+Pr66rwXXljWrFkTM2fORM+ePeHq6gqihLTRxo0bY+3atTh69CjLnMmbNy/EYjGkUilKlCiR5IvP9yPW3i29fv2aUTWkpARFIhGMjY1hbm6OhQsXsiytlIZMJoONjQ1cXV2TuIbSonQTp50SJbiYtFNZRSIRypYtC6KESmX+HvEFbxzHoXTp0hgwYABbg6enJ7OIChYsiHnz5jEXWrVq1TBgwAB2jWbNmum00bx37x5at24NogRXIV8NPGPGDFbJrq0cGjRogMuXL2PXrl1sPVu2bEnXD/rUqVOQyWSM+n3t2rXpOt+Q2Lp1K2QyGTw9PfHw4cMcW0dGEB0dzepNqlSpkioTbUREBIYOHco2F1u2bMmSYP/Ro0eZwipcuDB2795t8Ov8Ekph6dKlTEgQJZjv2k3dvb292Q+T/7HyLS6JCBUrVsTu3bthbm4OCwsLHD9+nGWptGjRAgqFgu0atXsnOzo6sl3D58+fsXnzZp14hT5BZ21tjerVq2PMmDFsLUuWLIG5uTnc3NwY3860adNAlOCf5knRtOeRy+UoWbIkJk6ciMOHD2PcuHFQKBQwMTHBvHnzEBsbyziFiAjTp0/X+WLExMSgaNGicHJyYgVSQILvlQ/K8fcyLcF7XnE6OjoyhaadOsrfy4iICLx9+xb29vYoUKAAqx/QVjg8L5X2/GKxWOceSCQS5kpKzAzLJw94e3vrzRYqU6YMXrx4galTp7I0ZO0MKg8PD7Rs2RJGRkYwMTHBnDlzMG7cOPZ+GjRooBMzuHbtGrPoRCIRLCws8OzZM0yaNElnbfzjxo0bMxcQEWH58uVp+zV/B19VbWlpCRcXlxx14Zw8eZK1+dR2lf0oWLduHSuaTEsSwJUrV5gbtGbNmlnS81qtVmPbtm3M4i1btqwOEWNm8UsohcOHDycRInXr1kX79u2TPM8XYHEch8DAQNY03szMDJs2bUKhQoVgZGSEI0eOYPjw4SAitvvmhRIvAIgSUl1DQ0Px+vVrXLp0ie1U+TFjxgzMmjVLbztL7TWVL18eCoUCdnZ2LDvG399fp0k9v5udOXMm2rdvz3z32gqnf//+OHLkCEun4+k9xGJxEgbM8+fPQyQSoVevXjrPR0ZGwtfXlzXpIUqwrlJyxXl7e+tQSG/cuFFvELZQoUJ48OAB9u7dy5RV4voBooSAcY0aNVJ1L8nlcp3CtMRuRAsLC530VX4ULFgQUVFRiIyMRFBQELMiZTKZTgtT/odZoUIFnDt3DhMnTmTKvGLFiggJCWE7uXPnzjGBoVQqsXTpUjx+/FinqppXWrxy4p+bNm1aunaE/KYhI0rF0Lh79y7y5s0LIyMj/P333zm6lozg6tWrcHd3h1QqRVBQUKqfQ3x8PBYtWgQTExMYGRlh5syZyRLxZQaxsbFYuXIls8Rr165tkDTWX0Ip8JkdaR28O4HvC6AtyE1NTXV86OmZV99wdXVF2bJlUa9ePVSqVIlVFru6uuoEhBNfixeANjY2GDp0KJ4/f85cPrwv+vbt26wvddOmTeHr66sTHC5fvjzGjBnDetcWKFAgScYIz+WfuMz/8ePHzAefWFgnblajPXjX3PDhw6HRaBAaGor+/fvruOD4zlmNGzeGkZERHj9+jDNnzrCdOJ9Cyiu6lPo18EMsFrN0XrFYzJQ9Px/vbtF+3tLSkgXuY2NjGaka/z75uIVSqYRcLodMJsOMGTPw6dMnLFy4kKXAent7Y/PmzcxFx3M4ESUU361duxb79+9nyov/3BPzOA34X3tXHlZV1fXXvhN3YBa4igioMSg4K6aiIs5DvSqUIyU5lYr5mKZmGppNpOWQpZJo6WtOqWnOH4nT61SWSpqJZWZaDg2IAwr39/1x2btz7sCkAuL5Pc9+4J5hn3333WetveYXXyy2R4/FYhEpwqtWrVrunkC///67SKv+wQcflOtYSoM///xTuBkPGDAAOTk5Rd7z66+/olevXmIN2Loj3y/cvHkTKSkpwt7Zt29fnDlzptT9PXJMwdY9tFevXvjyyy8RHx9fqBpE+nLq9Xp4eHhArVajb9++wjjq5+cn0z/bpphu2LAhevfuLfrTaDTQ6XRo3rw5GjVqhMDAwGK7h9oyCl9fXzz11FMICAiAj48Ptm3bBm9vb5jNZllU7j///IMtW7bg5ZdfljEJ3lq3bo39+/eLnU1OTg6Cg4MRFhaGmzdv4siRI3j++ecFgeZZYo8dO4ZPP/1UENV27doJt1Rn8xoZGSl0tbm5uUhISJDNc0BAAPR6PZ544gkA1jTUHTp0EOdNJpOInpbORZcuXcQ88nTojn7HKlWqyIzwXEU0evRo0Z9arZbp9TmT5FKYt7e3qLHB5zIkJAR79+4VKjrOJIODgzF//nxcuXIFNWvWRPXq1cX4Q0JCsGTJEiQnJwvvJUdroU+fPg6TKzrCnTt3hJfc22+/Xax7HiRu3LghpKLJkyeXaYDd/UB+fj5ef/11MMYQGRmJH3/8sVj3bdiwQaSBGT58+ANLCfLXX3/hlVdegdFohEajwfDhwwsNlnWGR4IpTJo0SfZi2ebRcZbf35aYcV0091Dx9/eHVqvFsmXL0KFDB6hUKri7u4s+eM0F3uLj4+Hq6oqoqCgMHToUOp0OZrMZfn5+OHXqlBjvxYsXBUHt0aOHTP2hUqmcSii2BN5oNCI5ObnQAizZ2dnYunUrRo4cKZsHV1dX9OjRA++//77I5srFVIPBgISEBGRkZGD37t3Q6XSIjo7GrVu3cPPmTSQnJ4t8LiNGjBC1GhwxMyJrrqD09HRYLBbs3btXBJRJr5FmOf34449lEs/cuXPt3Gfj4uJkbsNSom37u1etWlUmlTHGMGfOHNnv17FjR+Tk5ODOnTuIjo6GwWBAWlqasK+YTCZZYSTOnC5fvoz8/Hx88cUXQuL09fUVnlCTJ0/Ghg0bxLMiIiIwf/58wfyCgoLsHBXq1q0rs/MUhsuXL4sgzIpg7L17967IVzV48OBy8Y66V2zfvh1VqlSBu7t7sT3EsrOzMWbMGGg0GhgMBrz88ssiQ/P9xqVLlzBixAhotVq4uLhgzJgxdrVeCkOlYwrkIE7h6tWrdlGyjlpISAiSk5MLjRsoabPNIqrRaDBz5kx89NFHILIWdff19YW/v78sYObkyZNiHFLmtG7dOgBWNzij0QiDwSALrHOUiI4ziCZNmmDSpEk4ceKEQ3UCz+hqMpkwZMgQhzmQunfvju+//152H88t1adPH9FvVlaWMOLWrVsX8+fPlxFpHx8fOy8jb29vpKSk4Mcff5TlPOJt2LBhQrK4cOGCcBZgjGH58uVYu3atrE8+r0SEJk2a2CX0s2VStsyqUaNG6NOnj2wON2zYgIsXL8JsNiM0NBT//PMPvv76a8THx4v5r1+/vpBO1Go1Bg0aJFKg7NmzR9iguJF+x44dyM/Px8qVK4VKqnHjxnjppZdE3ExsbKyM4eh0OnzyySfF2m1PnDgRRFZJpbjM5EHCYrGIrK5PPvnkQ1P/WYpz584Je9qECROKzdzOnj0rJGJ3d3dMmzYN2dnZJXp2caOnf/75Zzz33HPCNXrChAnFYkSVjinwJpUUOAF21JxFGztyMzUYDOjRo4cwJNruZmvUqOHQK4Z71zhrPD+/Xq9Hr1690KJFC+FNwwkzJ/y+vr6CEBw5cgSenp6oXbs21qxZY2fEVqvVaN26NRo2bGgnHanVajz22GNITEzEp59+iszMTOTl5YkcSQ0bNhT6UFdXV9SvX1/2fcPDw5GUlISNGzciOztbeGNJC9gD1oRp3DDbr18/EewjZQS286FSqdCzZ0+MHDlSNp+MMbi5uWHGjBmCyEqJ9tChQ3H37l27mAouBaxatQqpqakOmR0ntI6OS2MeiKxSw44dO6BWqxEXFyd+j9OnT2PYsGFC7VOrVi2xjrg4f+nSJQDA8ePHZbaFgQMH4sSJE7h79y6WLl0q5qxFixYi/35YWJhdNH5ERITDvEpS/PPPPzAajWCMoXPnzhVmd/7BBx8IJ4oH4dv/oFFSt1UpMjMzhSrZx8cHM2fOLDZzLGlKjdOnT6N///7i/Zk6dWqhm4NHgilMnz69UKIsLZYi3WU621kW1UqSX4kT6OIYS/k13bt3x/79+3H16lUcPHgQJpNJGJkGDBiAAwcO4JlnnpHd27x5cxGAFRsb6zAVNQ8444xSo9EgOTlZ7GS4q2O3bt3QuXNnWVGf6OhosXOyNSTevHkTU6dOFSolvlPmrqlS119bZhwQECBTn/AYAX9/fyxevBh5eXn4+OOPZQT8999/t2NifJw///wzLBYLMjIyRI4m2+Yo/oEf45KIwWAQKdVnzZol+75Xr17Fm2++KaQULy8vWbW2fv36iWA1nvFWGsi2Z88e3L59G/PmzROBjM2bNxf/c+lBOubu3bsXqj8eO3as+C6jRo1yel1ZY/Xq1dDpdIiIiHBaEa2iIy0tTbitltSYfPjwYaGC9Pf3x4IFC4r0VCptnqXMzEzExcWJNfTGG2/g+vXrdtc9EkwhNTXVjggXRnw5AdBqtTJVkjT+gBMeo9Fop/OdNGmSIB785a1Vqxbq1KkjPms0Gvj7+8Pd3R2+vr6yimUlScBnMpnE893c3LBnzx6xqLgr57Bhw2T1Clq3bo158+Zh3759mDlzJlq2bOk0KZ6Xlxfi4uKQkpKC3bt3o2fPntDpdMjMzMStW7eQnp6OiRMnomnTprJxt2/fHqtXr5ZFWWdlZQmGwKWtbt26yRhiYmIi3njjDTvpSkqcZ8+eLdw169evj8OHD+Obb74Rc+vi4oLPPvsM1apVszPW6nQ6me/4gQMHilQXenh4yOaH59fn30OlUjkswnLnzh2sWLFCuKHy8fF7W7ZsiZUrV6JLly4wmUwYO3asiIto0aKFkMJmzJgBd3d3MMZk+bS4y6N0XQ8bNszhi/7LL79ArVaLwKeK5AHEExcGBgbKbGsPE44ePYrg4GDodDosWrSoxPdnZGQIL8Dg4GDMnj3bqVrpXpPvSbP5+vr62tmaHgmmsHHjxmITWZVKhcDAQFmWS6nhkxOnxo0bIzMzUxA0adAaEQlxTUoMFi9eDA8PDxiNRhiNRuh0OrEzrlatGlavXo3du3eLoChey2HDhg2YO3cuYmJiisXYVCoVqlWrhjZt2sDDwwP+/v64du0aTp8+jenTpwvdvkqlQmxsLIYMGSLKjoaEhMiibTmz4f9zo6Wfnx+WL1+OX375RahPrl69iqVLl4q++Hfv3Lkz5s+fj/Pnz8NisWD9+vUwm83img8//FBWyMjPzw/79+/HzJkznUY3v/TSS1i1apUIZBs/fjwyMzNlElB8fDwMBoP4PlIJ7uWXX8Zff/0FwKrjTk1NLTJbrFSNpNPpZIzczc3NYUEd3v/+/fvx1FNPycbAf3uz2QydTof27dsjJycH8+fPF+qj+vXrY+XKlfjjjz8wbtw4uLi4QKPRiN/EZDLBbDbDxcVFrAm9Xo/p06fbpUHo27cv3N3d0aVLF6jV6jLPpFoYjh49CrPZDG9vbxw4cKC8h1MqXL16Vez6hwwZUuIcRRaLBVu2bBFp/t3d3TFu3Dg7R5H7lZH1wIEDGDlypJ1d6pFgClwna7tjdPbyu7i4wNvbG++++67YiTlS8fAXU6VSoX379jhw4ICMcEgJWq1ateDp6SlywtiqmAYMGCDGyz1TatSoIQgXB4/O5rtxIsLYsWMxZ84csZgKy8rasGFDDB8+HHPnzsXIkSNlY4yOjsbSpUvx559/YseOHUIPzc/7+/ujefPmdoF21atXR3x8PN577z0cPHgQp06dgpeXF4KCgpCUlCTT4Tdq1AjJycnIyMgQhWmIrN5Fp0+fljGg9u3b46effsJrr73mUL3m5eWFBQsWCH/8xx57DJ999plM786zovLoZql6ytPTE7NmzRIRv5cuXRLxCo4YBDckS89J+zMYDEUWmjl37hzGjh0r2+FLnSDatGmDb7/9Vriz8nGHhIRg8eLFOHv2LIYOHQq1Wi1T81WrVg1GoxGdOnUSx3mNbI4jR46AyBoQ2KBBA7i7u9s5DZQnsrKyULt2bRgMhoeuLgNHXl6eCC5t1qwZzp8/X6p+Dh8+jH79+onfuW/fviL78YMu2flIMAXbwvXFaW5ubtDr9Vi2bJnM3dFWHREYGChEck9Pz1LbITjR4cSPExtXV1fExsbi6aefxvDhw0W2Vn5Pnz59cPLkSaGmWbZsGRhj6NKlC7799lukpqYKn3pn1di4rz3fUWu1WvTs2ROjR48WO89q1aohNjZWSD08O2nv3r0RHx8vYxRGoxFRUVFQqVRo3bo17t69i1OnTuGdd95Bq1atZNHLvXv3Fgxy5MiR+Pvvv2WGYq1WixkzZuC3337Diy++6JBY85oVfAwDBw5EQECATI9PRKIynDQ4jMjq9vnpp58KzynOZDgzlz5L+vtKGYKUeY4ePbrIYu85OTl46623ZAF5UsbXpk0brF+/Hnfu3MGaNWvEGgwMDMS8efPw3XffyZgqX386nQ5Lly6VJXusXr26KNjUpk0bBAYG4uzZszCbzahZsyYuX75cNKUoI/z+++9o1KgRNBpNqepXVxSsW7cObm5u8PX1xVdffVXqfn755ReMGzdOrJPo6GhERESgTZs293G0cjwSTIEbSB21OnXqOFUbcF1948aNZaoC2+br6yvcI9u1aye8dhwRX07oqlWrJggs3w26ubnBxcUFVapUQUJCgiBi3t7eCA0NFdk8C2Nk4eHhor86depg3rx5WL58Odzd3dGqVStcuHBBSBR6vd7hDlzq088Ts6lUKnTt2hV//vknVq9eLVwwOYPo2rUr3nrrLSxYsACjRo2SSUxGoxH9+/dHWloazp8/jz/++AMfffSRrBoeJ9CRkZE4evSoYCr8GQEBAdixYwd+/vlnsZu3ba6urmjWrJmIITGbzbL0G0ajETExMXZSGpcsmjRpIkqmLly4UNwXGBgoYwbSOhI6nU6WR4s3Ly8vpKWlFRlJnJ+fj7fffltml/Ly8hJBdbVr18a8efNw/fp1bNmyReidzWYz3nnnHWRkZMjSmfAAvCVLluDXX39F27ZtxbnIyEjMnj0bRISVK1fi4MGD0Ov1iI6OrlAprv/++29ER0eDMYbU1NTyHk6pcerUKYSHh0OlUmHmzJn3FKyXnZ2N2bNni42PXq/H3LlzHdqP7hWPBFMAIIq7OGodO3YsNPcQJ7DSl8/WS2bu3Lno3LmziFJu164dtm/fbkeAuEfK+++/L3YSfDfHX2ppGuZZs2aByJriwGKx4OjRozCZTIiMjBR916hRA8nJyRg9ejTi4uLQvHlzpxILJ3QNGjTAlClTsGzZMmzcuBELFixA//79ReZQZ81kMiEuLg5Lly4VqrDo6GjZjrpZs2aYNWsWDh06JPIbSXX9oaGhSEpKQnp6OrKyspCSkiKLYVCpVIiPj0dERAR0Op0sy2q3bt3w3Xffwc/Pz6EHFZFVEuDf39XVFSaTSailuIcVnzv+O4aEhAhVU69evfDjjz9i69at0Ov1YtMgrdFt28xms8PiTU2aNCmWR0pubi7Gjh0rkzjq1q0rmKuXlxcmTpyIX3/9FRkZGSIg0MvLC8nJycKAL20zZswAYE3IJ2WkBoMBkZGRsFgsWLVqFYgIzz77bIWKML5x44bIlWXr3fUwITs7W3j88MJS94K8vDzUrVtXrH0PDw+MHz++1GoqR3gkmIKjOIWS5C2SvqiJiYmCoNgS/LS0NJFM7+DBg8jPz8fYsWNl16jVajRt2hRGo1EkLatatapsZ96uXTtZ/dcxY8aAyFonukaNGggICMDFixfx/vvvi3saNmxoF5jCa0U/+eSTYsfIGEOjRo1Qo0YNOy8nd3d3NG3aFHFxcejXrx+6d++OsLCwQr2huCTxxBNPYNy4cUhMTJTVRGjVqhVCQkKg0+nw3//+F++99x66du0qJBQvLy8kJCRg3bp1WLFihayGM2NMeNgMGjRISA4ajUYEt40YMUJGrLmtSNqHSqWC0WgUKTgYY7JyoDVr1oSrqyu8vb2RkJAAk8kErVaLMWPGiHQhfCfvTKrkz3VW1e/ZZ58VMQqF4fTp0+JZfE0EBwejSZMmUKlU0Gg0SEhIwNGjR3Ho0CGRfNHV1RWhoaFgjMmkjtjYWEHsv/jiC5n9Ijo6GmfPnkVycjKICCkpKUWOryyRm5srVGRTpkypUEyrJOD1JVQqFSIiIoqs6lYUuE3h4MGD6NOnj8zucD8y0T4UTIGITET0CRGlEtGA4twjZQqceBbVvLy87LyIbFv16tVF4jln8QhGoxENGjQQ3kItW7aUXevh4QE/Pz+EhYUhLi5OxmR4Hp2OHTsK74X8/HwR6OLi4iLLhMiTrvFI2itXrohzFosFgwYNEoRwyJAhYIwJP/Vbt27h5MmT2LRpE2bPno1Ro0ahU6dOoji59DvxYvfOmKmjQL569eqJ3TcfY0pKCq5du4acnBysW7cOzzzzjFDFGAwGtGrVSnhPqdVqmddXr169ZAnrdDod9Ho9Ll26JOoWEFndRHmfUnWbVquVjUdqY+nTpw8aNGgAIqsf/+DBg0Wa64kTJyIgIACurq4yScpR/W4p07B1Vdbr9Xj33XeLzF2UkZEh1gnRv2pHHx8ftGjRQvTbrl07bNq0CUePHhWBcCqVCmq1WqbeM5lMIrVyXl6eLBcYYwwDBgxA9+7dwRgrstpYWSMvL0/YeJKSkso9sd+9YMeOHfDy8oKXl1eJK+tJ0bZtWwQFBYnPtnaHVq1aYe3ataUuwlNuTIGI0ojoMhFl2hzvQkSniSiLiCYWHEsgoicK/l9VnP6lTIEnMXPUbNVGer1e7BRtr+UERq1WCz95Rzp+HmWr0WiwaNEiWCwWbN26VXYNd3PlRmOpNxSXILp164bbt28jPz9fiKBqtVpmuNqzZw+IrHYNvV6PevXqCcPh5s2bZbWFU1JSkJSUBJVKJYKnnOHWrVv4/vvvsX79eqSkpIhEZkUl7GOMwcfHB/Xr15epfWzzS8XGxuLzzz9Hbm4u7ty5g/T0dCQlJQkmR2Q1nPI+pEyVq9r4bxQQEIBTp05h+fLlMgmhbdu2ggnYBqRxwirtd/HixSJCNTo6Gjt37hRqmqCgIFSvXh16vV4waCrYSEi/vyPDtO2cBQQEYPPmzYXOf69eveDq6iqiut3d3UV+JIPBgJYtW8o8qxYuXIhvv/1WFiXdu3dvWcK/du3aiXocL730kmxeGGNCIpKqLysCLBaLeIefffbZChORXRpkZWUJ1e97771XKumHS/22yM7Oxpw5c0TcTc2aNUtldyhPptCGiBpLmQIRqYnoLBHVIiIdER0jorpENImIGhZcs6I4/UuZwrx585zu5mwJsvS8I8Zgmy1VrVY7TInBxXTpzkvqi0/0r47aZDIJFRAVcHqehK5v376i5kFycjIiIiLg7u4uq/DFDdIvvvgiDAYDIiIiMH36dKhUKjRq1Ajnzp0TjGrevHkwm82Iiooq8a6Lv5gfffSRiH94+eWXERcX59Szyc3NDY0bN0aPHj3s6jvwOezSpYvYOVksFnz99ddCn2w7525ubmLR2+aUmjJlCrZt2yb7Pby9vTFgwACZhMCbbZ1tXlPiv//9L0wmE3x8fLB161Zs2bJF6OTd3NygVqsxfvx4wVD8/PxkY2nQoIGMOPPvaStZRkVF4YcffnA41z/88APUajVGjRqFY8eOiUjuZs2aoXfv3mJtNmvWTNgdfH198dprr2Ht2rVCPadSqewCMF944QWcOHFClFzlfUsN8nv37i3R2njQsFgsmDZtmmB2FckwXlJcv35dbCwSEhJKnPvJGVPgyMvLw+effy6cEry8vPDKK68US30JlCNTsD6bgm2YQgsi2i75PKmgJRBRj4JjKwvpbxgRfU1EXwcGBsomKTIy0ilxL6zZBnJJ7+VSgrM8SQ0aNICXl5cIPvn7778d7rQ1Go1Q83BGsWHDBlmhlMTERFgsFpw/fx7+/v6oXr26MC5duHBB7H43b94smFrXrl1F3vfc3Fx06tQJKpVKRDovXLiwWIuE4+bNmwgPD0dAQAA2bdoEjUaDLl264O7du/j999/h7e2Nhg0bIi0tDVFRUQ49m6R68tDQUBnjdXFxQXR0NFJTU3Hx4kVRxKht27ayutU6nQ7jxo1D//79ZbYYIutO/YMPPoCnp6fMe6lBgwYYOXKkXaS0n5+fbFNgNpuRnZ2NU6dOITIyEowxvPrqq7h9+zY+/PBDWa6mYcOGCUKv1Wpllf14ym8ikqmcHKkc4+PjHeaiGT58ODQaDc6cOYP8/HwsWLAAHh4e0Ol0GDt2LMaPHy+84urVqycip11cXPDkk0/CxcVFBMbZPtNoNKJx48bQ6/W4cuUKVq5caVc2NT4+HidOnCjRGnnQ4N5TPGvtw4r8/HyRfqdZs2bFSvHBbYS27bXXXnN6z//+9z/07t1bxNgMHjy4yNiUisYU4onoY8nnBCL6gKw2hSVE9BGVwqZw7do1NG3atNCI1cKYhLRCmJS4aTQapx4wRFZPIzc3N7Rs2VKknpg6darddTqdDjVr1hQBRfylXbNmjSAimzZtEt/n2LFjouYvD27jTIW/2FqtFuHh4bLc7Tk5OXj88ceh0+kEwyqpj/qRI0eg0WgwcOBAUat6zJgxAP51/X3nnXfE9WfPnkViYmKh88Szmdoy16CgIISGhoKI8PHHH2Pjxo0y5sCLB0mJrjSPktlsFkyDE+ynnnoKiYmJdmOQEuvHH38cFosFN27cEPrsmJgYXLx4EX/++SeSkpLEc7hRnX/m34FvJjw8PODm5oa9e/eiT58+Tm0ynFlLpbeLFy/CaDTi6aefFsek9pOQkBBhD+JzEBwcjNatW8sYclRUFMaOHSurESF9do8ePZCbm4ucnBy88sordnEcvXv3ltXkKG+kpaVBpVKhVatWFSLr671gw4YNcHV1RdWqVbF///5i3VOUpOAIZ86cwYgRI8SmtHv37ti1a5dD9VVFYwpPOWAK80rYp13q7MmTJzslSET/qouc2RKISJaNU8oYpNdzw6iUwHAVxaRJkwBYGZQzvXz79u1FHWVuMAwODkbdunXh6ekpq6aUnp4OrVaLtm3b4vbt2zh27JggTLNnz8aePXtgMpkQFhYmS79w7do1REREwGg0Qq1W47nnnivR4gIgxPg1a9YIzyhuO+nVqxf0er1dbVqLxYJvvvkGzz77rJ0B1sPDQ+ZCa1txjLfmzZtjzpw5Ig5EpVLJdv62nkc8Qy0ngjVr1oSLiwuMRiOef/552SbBlmn16NFDjP2TTz6B0WiEn58f/u///g+ANbU5t0dxRqDVamXxLAEBAUKtVrduXeTn5yM3NxdTpkxx6jLMq7dxcNUhj2bl2L59u7BfcM+mFStWiCA3X19fxMbGirk2m81YuHChXXQ/n6vg4GAsXboUeXl5yMrKEsyYp2QhskqezghJWWPNmjXQaDRo2rQprl27Vt7DuSdkZmaidu3a0Gq1xYrLKA1T4Lhy5QqmT58u1NtNmjSxS6NR0ZiCQ/VRafqWSgrp6emFMgVHxMfRy+qIMUh3VFJiJyX8nDFs27YNt2/fFioiafAWb1IjJpHVZ/7s2bMi6lhqNFqxYgWIrJKMl5eXeD6Pdty3bx9cXV0REhIiE09/++03BAcHizHu27evRAvrzp07aNasGapUqYLz58+jS5cu0Gg02LVrF3799VeYTCZRLc3Z/YsWLXLIgMPCwhATEyO8sJz9RlK13sCBAwURk/5etv2r1WoEBQUJ47G/v7+d9Chl6LVr1xYEJzMzE3Xq1AFjDMnJycjLy4PFYkFCQoLsfl6cSMpsuMdUZGSkTHKzTckhbQaDAampqcjOzoavry9iYmLsiPHNmzfFzr5KlSpYsmQJ8vPzkZ6eLmwyPHsn7zcoKEg4DUi/K187derUwZo1a5CXlyeL7QkPDxeqs6ioqHvybrlf2LRpE1xcXFCvXr0SFZGpiLh27ZrImzRy5MhCM6Xaeh+VBjdv3sTChQvRoUMHO8N9RWMKGiL6iYhq0r+G5ogS9mknKdjGKRSV2lqlUgkXPelxo9FoZ5zkTVqGkh+TqjqqV68OX19fkTrCxcUF/fr1c9gf3wHy3eTOnTuxc+dOqFQq9O7dW6ZikOZkP3PmjDCoclF0//79cHNzQ926dWV5lM6cOSMipMPDw0vs0fHDDz9Ar9ejW7du+Ouvv1CnTh14e3vjzJkzePfdd0FERbo37t69GxqNBiEhIWLnIlVdxMbGYty4cXjqqafsIopdXV3tSncajUbUqlVLZtCWpq3m/RuNRrz++ut2diZna2HChAm4fv06cnJyREry9u3bC0I0d+5cEa/BCezQoUNl4+Xfz9PTEytXrpTNw/bt250WgfLw8BCVypx5LGVmZgpng5iYGGG8Pn78OJ555hnZ95emPOGR3/yci4uLWD+NGjXCunXr0KhRI2i1WhgMBuh0OnTq1EkYrkNCQrBw4cISJ367n9i5cyeMRiPCwsJKVXqyIuHu3bsYN26c+B2l7uVSVNrcR0T0GRFdIqK7RHSBiAYXHO9GRD+S1Qtpcmn7l0oKgwcPlr1ojtIS2DZnaqTCXDLd3d0RGhqKlJQUGTPg/3P1xjPPPIPRo0dDo9Hg8OHDhRqzAwMDERwcjOvXr4vo5tdffx2AtUiJ9PssWbJEuL4GBQWJneWuXbug1WrRrl07mY/8t99+K1QDvM+SgNcCWLRoEbKysuDt7Y3w8HBcuXIFERERCAoKkgXhOQL3DHv11VfxxRdfiJS+0rl2d3fHsGHDkJKSItI4FPbb8apltllNba+Li4vDe++9J57D+5W6cfLm5uaGN998E3///TcWL14MvV6PqlWrYteuXQCAVatWQavVyiQPHkwmjUPhto3nnnvObm5SU1Odri+VSoXg4GCnu/P8/HwsXLhQ5D9KTk4WHjrnz5+XFSyqV68e2rZtK/M2kq45Dw8PwaSaNm2KKlWqwN/fX3hUVa1aFSNGjBA5v8xmM954441yK5Szd+9e4ZVmq7Z8GPHpp5/CxcUFwcHBMi9DjkrLFB50kzIFWxGfiMSCLqw9//zzduocR4Tbto0bN06WzdR2FzhhwgScO3cOarUaY8aMwdGjR+2YEP/MI4qTkpJgsVgwcOBAMMbw/PPPg8garZyTk4MOHTpAo9EgPT1d6JUXLFgg5mDZsmUgsrrASdUQu3fvFoFpp0+fLtHiyc/PR4cOHWAymZCVlYWMjAxoNBp07txZqOwmT55caB8Wi0UYfnmp0V9++QVTp04Vu1Y3NzdB1Hkkdnh4uAgOtP0dpJ9r1aoFnU5nx0ikhuHRo0fbuaw6S/fh6emJadOmYd++fQgLC4NKpcKMGTOQn5+PHTt2CHuQ9Dlubm6y+It69eqBMWvxd9v6ATk5OaL6nbONB7drOILUEB0WFiaYFmB1UJCOr169emjcuLHou2rVqrKMtl5eXsJGolKpEBkZiT179ggX1qZNm2Lu3LlCVeXq6oqxY8eWy4798OHD8PLyQkBAQInXcUXEoUOH4O/vD5PJJN4LDoUplFwCsVMfbdmyxeEL5qgUpG3bs2cPrl696tQP39Zoyl8gnqucH7MtDP/VV1+JlApXr17F+vXrnTKdnj17gjGGvXv34ubNm0JVFRMTI3b+PLuop6cnNm/eLMYm9YHmkaxTp06VLQLuRWQ2m0vs/33+/Hl4eHigVatWyMvLE32NHz8eCQkJ0Ol0Tn3xOW7duoWoqCi4urrK3OVyc3Px2WefCX9rnU4n28VXq1YNs2bNwpAhQ+wyovJm62lTpUoVhwGHjiQQR78tXwfu7u4YP368CCrs3LkzLl++jKNHj8rUhtxmoNfrZfUw6tWrJ4LFli9fbjcnu3fvLlSi9fPzE5lPHWHbtm1CzZOYmChSoBw+fBgGgwG1atVC7dq1RV+2fUslK56okX+fr776CsuXLxdScL9+/bB9+3YMGDBAOF8MGjQIJ0+eLNFauld899138PX1hdlsrnCutKXBb7/9hqioKBBZY5S42lhhCvdBUti2bVuRxN+2SYn4unXrROZRR81RDeY6derg9u3bsjxA0ujRqlWrinQG06ZNw08//SR7EaW7XR8fH9SoUQOhoaFCZaTT6RAaGipzyfv5559FOmSeE2fEiBHivMViES6WaWlpsoXAjY/t2rUrcVDb8uXLQUR4++23AQAvvPACiKzFczw8PNChQ4ciPVYuXLgAs9mMxx57zK6GBGCtFpWYmCiIk9TLx8/PD6NHj7YzNgcEBKBv374O9fW26iTu7WU799LfTHqMq6ZcXV3RuXNn6HQ6VK9eHXv37sW5c+dkkgZ/sYnk6kS9Xi8I99ChQ+2CmK5fvy4kQmd1vr29vfH55587nNMbN25g4sSJ0Gg08PHxwbJly2CxWLB27VoQWRO0rV27Vuz8+Zh4ug9bd14pM23Tpg02btyIV199FXq9HgaDAcnJyTh58iSSkpKEGuw///kP/ve//5VoPd0LTp48CX9/f1SpUgXffPNNmT33QeHWrVvCjhUXF4fr168rTKG0TcoUOJEqqtnuFKVuiiqVCmlpaU7vdWSDePfdd3HixAmnaqaYmBh069YN3t7eaNWqFVxdXWWRvNImPd61a1fs3LkTGo0GPXr0kBHxQ4cOwWAwoGnTpkJtInU5u3PnDjp27AiNRiPLv3Ljxg1BaLmqqriwWCyIj4+HVqvFd999h9zcXLRp0wZ6vV4EoBUnN/6+ffug1WrRtWtXp7rzK1eu4K233pLtonn1NcaYnX1GpVJh+PDh2LFjB6ZMmWKXhsJRc8QIiKxqQClR9/HxQceOHUUtBB4wN336dFy5ckWmiomMjIRGo5GpkTiR5WOuX7++Q6lq+/btoj6EI+mFM425c+c6ZOjHjx8XhL9jx47IysrCO++8A6J/1Xt79uyRreE2bdqgW7duMsO/1IDPr3388cexZMkSYW8ICgrC+vXrcfnyZbz22mtCGm/dujW+/PLLMnFnzcrKQlBQENzd3cuUIT0oWCwWzJo1CyqVCvXr10fz5s3LjSkw6/mHC4yxJ8iqQkogou/L6LEeRPRPGd5f1PWFnS/pOdtjjq7xIaKrhYznfqIs57o415Zmroszz46OVdZ5Ls71pT2vrOmSXxsCwMPhGWfc4mFoVAi3ewDPWlSW9xd1fWHnS3rO9piTayrlXBfn2tLMdXHm2cncV8p5Ls71pT2vrOmSX1vYNaoiuImCf7GpjO8v6vrCzpf0nO2xe/2u94qynOviXFuauS7OPBf3+Q8KD9OaLuy8sqZLfq3Tax5K9REHY+xrAE3LexyPApS5Lhso81x2UObaMR52SWFReQ/gEYIy12UDZZ7LDspcO8BDLSkoUKBAgYL7i4ddUlCgQIECBfcRClNQoECBAgUCClNQoECBAgUClYopMMZ6MsZSGWNfMMY6lfd4KisYY3UYYwsYY2sZYy+U93gqOxhjJsbYN4yxHuU9lsoKxlgMY2xvwbqOKe/xlCcqPFNgjKUxxi4zxjJtjndhjJ1mjGUxxiYSEQHYAGAoEQ0ioj7lMNyHFiWc51MAnieip4lIcekrIUoy1wWYQESry3aUDz9KOM8gohwi0pM1zf8jiwrPFIhoKRF1kR5gjKmJaD4RdSWiukTUjzFWV3LJqwXnFRQfS6kE88wYe5KI9hFRetkOs1JgKRVzrhljHYjoJBH9UdaDrARYSsVf03sBdCUrA55WxuOsUKjwTAHAHiL60+ZwFBFlAfgJwB0iWklE/2FWvENEWwEcLeuxPswoyTwXXL8RQEsiGlC2I334UcK5bkdEjxNRfyIayhir8O9sRUFJ5hmApeD8X0TkUobDrHDQlPcASonqRPSr5PMFImpORElE1IGIPBhjjwFYUB6Dq0RwOM8FOtfeZH15tpT9sColHM41gFFERIyxQUR0VUK8FJQOztZ0byLqTESeRPRBOYyrwuBhZQrMwTEAmEtEc8t6MJUYzuY5g4gyynYolR4O51r8Aywtu6FUajhb0+uIaF1ZD6Yi4mEVRS8QUQ3J5wAiulhOY6nMUOa57KDMddlAmeci8LAyhSNEFMIYq8kY0xFRXyLaWM5jqoxQ5rnsoMx12UCZ5yJQ4ZkCY+wzIjpARGGMsQuMscEA8ohoFBFtJ6JTRLQaQFkV26mUUOa57KDMddlAmefSQUmIp0CBAgUKBCq8pKBAgQIFCsoOClNQoECBAgUCClNQoECBAgUCClNQoECBAgUCClNQoECBAgUCClNQoECBAgUCClNQUCnBGMtnjH0naROLvuvBoyBp41eMMfeCz2CMLZOc1zDGrjDGviykDxNj7BpjzMPm+AbG2NOMsR6MsUc606eC0kNhCgoqK24BaChpb99rh4yx+5ErrBsRHQOQXfD5BhFFMsYMBZ87EtFvhXUA4AYR7SCinpKxeRBRNBF9SUSbiehJxpjxPoxXwSMGhSkoeKTAGDvHGJvGGDvKGDvBGAsvOG4qKMpyhDH2LWPsPwXHBzHG1jDGNhHRDsaYkTG2mjF2nDG2ijF2iDHWlDE2mDH2vuQ5Qxlj7zkYwgAi+sLm2FYi6l7wfz8i+kzSj8NxFVzTV9JHLyLaBuAmrBGpGUSkVGpTUGIoTEFBZYXBRn0krcR3FUBjIvqIiMYVHJtMRF8BaEbWGgbvMsZMBedaENGzAGKJaAQR/QWgPhG9TkRNCq5ZSdbdubbgcyIRLXEwrlZE9I3NsZVE1Jcxpiei+kR0SHLO2bi2EVETxliVguv6koSZENHXRNTa6ewoUOAED2vqbAUKisItAA2dnOMpkr8ha10IIqJOZCXqnEnoiSiw4P+dAHixlmgimkNEBCCTMXa84P8bjLGviKgHY+wUEWkBnHDwbG8A16UHABxnjAWTVUqwrU/hcFwATjHGNhJRPGPscyJqSFaVEsdlIvJ38v0VKHAKhSkoeBSRW/A3n/59BxgRxQE4Lb2QMdacrHp/klznDB8T0StE9AM5lhKIiPIYYyoHxXI2EtFMIoohoiqS4w7HVYDPyFp6lhHRFwDuSs7piehWIWNVoMAhFPWRAgVWbCeiJMYYIyJijDVyct0+Inq64Jq6RFSPnwBwiKy5+vuTXJUjxWkiquXgeBoRTXcgXRQ2rl1EFEJEIx08L5SIMkmBghJCYQoKKitsbQpFeR+9TkRaIjrOGMss+OwIHxKRb4HaaAIRHSeifyTnVxPRfgB/Obl/M1mlARkAXAAwpyTjKpA2PierZLHH5r52Bc9SoKBEUFJnK1BQAjDG1GS1F9xmjNUmonQiCi0oAk8F8QXvA0h3cn81IvoUQMcHOEYzEa0A0P5BPUNB5YViU1CgoGQwEtGuAi8jRkQvALjDGPMkosNkjUFwyBCIiABcYoylMsbcJbEK9xuBRPTSA+pbQSWHIikoUKBAgQIBxaagQIECBQoEFKagQIECBQoEFKagQIECBQoEFKagQIECBQoEFKagQIECBQoEFKagQIECBQoE/h+VTRLxw0hVtgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdEUlEQVR4nO3df5QdZZ3n8fdHfgjmAFGJoIg0sKiDu+qQFhRwUzBBEgZBQdgwgwNZ1z5odODMsEhWl3RkHNBBEByjQxwIDjP8mBEhKCIS7UYXlXQwYgTiZjOJiYQhEQ0oIAG/+8et29xu7q2u2/feuj/68zrnnnRVPbfq20+q61v11FP1KCIwMzOr5SXtDsDMzDqbE4WZmWVyojAzs0xOFGZmlsmJwszMMjlRmJlZpp3bHUAr7L333tHX19fuMMzMusaqVau2RcSMast6MlH09fUxMjLS7jDMzLqGpI21lrnpyczMMjlRmJlZJicKMzPL5ERhZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmHSBJEpIkaXcYVTlRmJlZJicKMzPL5ERhZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmlqmtiULSNZIek7SmxnJJukrSOkkPSDqs6BjNzKa6dl9RLAPmZCyfCxySfgaALxYQk5mZVWhrooiIe4DHM4qcDHwlSn4ITJf06mKiMzMzaP8VxUT2AzZVTG9O55mZWUE6PVGoyryoWlAakDQiaWTr1q0tDsvMbOro9ESxGdi/Yvq1wCPVCkbE1RHRHxH9M2bMKCQ4M7OpoNMTxXLgL9LeT28HtkfElnYHZWY2lezczo1LugFIgL0lbQYWAbsARMSXgDuAE4B1wFPA/PZEamY2dbU1UUTEGRMsD2BBQeGYmVkVnd70ZGZmbeZEYWZmmZwozMwskxOFmZllcqIwM7NMThRmZpbJicLMrMWSJCFJknaHMWlOFGZmlsmJwszMMjlRmJlZJicKMzPL5ERhZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmlsmJwszMMjlRmJlZJicKMzPL5ERhZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmlsmJwszMMjlRmJlZpp2zFkr6q6zlEXF5c8MxM7NOk5kogD0KicLMzDpWZqKIiMVFBWJmZp0p1z0KSa+XtELSmnT6zZI+0drQzMysE+S9mb0UWAjsAIiIB4B5rQrKzMw6R95E8bKIuG/cvOeaHYyZmXWevIlim6SDgQCQ9D5gS8uiMjOzjjFRr6eyBcDVwBsl/RL4d+DMlkVlZmYdI1eiiIj1wGxJ04CXRMSTrQ3LzMw6xaQeuJME+IE7M7OpIO8Dd28A3gYsT6ffDdzTqqDMzKxz5HrgTtJdwGHlJidJg8C/tjw6MzNru7y9nl4HPFsx/SzQ1+jGJc2RtFbSOkkXVlmeSNouaXX6uajRbZqZWX3y9nr6J+A+SV9Lp98DXNfIhiXtBHwBOA7YDKyUtDwiHhxX9HsRcWIj2zIzs8nL2+vpU5K+CbyT0rMU8yPixw1u+3BgXdqjCkk3AicD4xOFmZm1UT3jUTwP/KHi06j9gE0V05vTeeO9Q9JPJH1T0ptqrUzSgKQRSSNbt25tQnhmZgb5Xwp4LvDPwN7Aq4DrJX20wW2ryrwYN30/cEBEvAX4PHBrrZVFxNUR0R8R/TNmzGgwNDOz4m3YsKHdIVSV9x7FB4AjIuJ3AJI+DfyA0sF7sjYD+1dMvxZ4pLJARDxR8fMdkpZI2jsitjWwXTOzjrRx48Z2h1BV3qYnUWp6Knue6lcE9VgJHCLpQEm7Unob7fLKApL2Vfp0n6TD03h/1eB2zcysDnkTxbXAjyQNps9Q/BD4x0Y2HBHPAR8BvgU8BNwcET+TdI6kc9Ji7wPWSPoJcBUwLyLGN0+ZmXWtwcFBJDE8PAyU3nwhicHBwfYGVkF5j7uSDgOOpnQlcU8Tej21TH9/f4yMjLQ7DDMzAJIkAWBoaCizzPDwMO06F5a0KiL6qy3Le48CSm+MfS79jiQdFhH3NyNAM7Oi5Tl4W0muRCHpYuBs4P/xQs+kAI5tTVhmZlPPAQcc0O4Qqsp7RXE6cHBEPDthSTMzm5S+vr52h1BV3pvZa4DpLYzDzMw6VN4rikuAH0taA/y+PDMiTmpJVGZm1jHyJorrgE8DP6U5r+8wM7MukTdRbIuIq1oaiZmZdaS8iWKVpEsoPTld2fTk7rFmZjl06nuc8sibKP44/fftFfPcPdbMLKdOfY9THpmJQtI7gB9GxDEFxWNmZh1mou6xZ1FqdrpR0tmS9i0iKDOzXtAN73HKIzNRRMQ5EXEYMAi8HFgm6QeS/lbSf02HMzUzsyoGBweJCGbNmgVARBARvZUoyiLi4Yi4IiLmULov8X3gNOBHrQzOzMzaL/dLASW9HHgN8DRwZ0Tc0bKozMx6TKe+xymPiW5m7wUsAM4AdgW2ArsB+0j6IbAkIr7b8ijNzLpcp77HKY+Jrij+DfgK8M6I+E3lAkkzgfdLOigiGhrEyMzMOldmooiI4zKWrQJWNT0iMzPrKBM1PR2WtdxPZpuZ9b6Jmp4+m7HMT2abmU0BEzU9+YlsM7MpLtdzFJJeJukTkq5Opw+RdGJrQzMzs06Qd4S7a4FngSPT6c3A37QkIjMz6yh5E8XBEfEZYAdARDwNqGVRmZlZx8ibKJ6VtDulG9hIOpiKcSnMzKx35X2FxyLgTmB/Sf8MHAWc3aqgzMy6RZIkAAwNDbU1jlbKlSgi4tuS7qc0cJGAcyNiW0sjMzOzjlDvA3db0n9fJ+l1fuDOzLr5jLqbhyctUt4H7nYD+oGfULqieDOlV4wf3brQzMxaa6LhSbs5CTbTRAMXHZM+dLcROCwi+iNiJqUxtNcVEaCZ1S9JktGDnFmj8vZ6emNE/LQ8ERFrgLe2JCIzsxbqleFJi5Q3UTwk6cuSEkmzJC0FHmplYGZmrdArw5MWKW/32PnAh4Bz0+l7gC+2JCIz6xhuozfI3z32GUlfAO6m9NDd2ojY0dLIzKwh7tEzsW4enrRIuRKFpAS4DthAqdfT/pLOioh7WhaZmTXEPXom1s3DkxYpb9PTZ4F3RcRaAEmvB24AZrYqMDPrHU5K3S3vzexdykkCICJ+DuzSmpCsm/mGYHu1q0ePm7l6W95EMSLpH9NeT0na68njZdsYSZKwePHidofR1Rp9/qFdPXomauay7pY3UXwI+Bnwl5R6Pj0InNOqoDpVN54td2LMzYqpqN8tSZIJ27L9gJv1slyJIiJ+HxGXR8QpEfHeiLgiIhp+zbikOZLWSlon6cIqyyXpqnT5A1XePZVLngNKnjJ5zpaLPAgWHXOtcvU2dzTjqqPoq5c8Z8zNan6ZaD15k1KzevTUiscPrk0deYdCPVHSjyU9LukJSU9KeqKRDUvaCfgCMBc4FDhD0qHjis0FDkk/A0zy2Y08B5SJyuQ9W2zWwasTY65VrtzcERFAdnNH3pi68WDTrOaXZq0nT4+ePMmtVjx+cG3qyNv09DngLOCVEbFnROwREXs2uO3DgXURsT4ingVuBE4eV+Zk4CtR8kNguqRX17OR8oFp2bJlAOzYsYMkSbj++usBeOqpp5g+ffpo+e3bt5MkCbfccgsAF1xwQdWzpjPPPBOATZs2kSQJd99995htlsuvXbuWJEm49957AVizZg19fX2sXLkSgNWrV5MkCatXrwZg5cqVYw6m9957L0mSsHZtqS/B8PAwSZJwxBFHAHD33XeTJAmbNm0C4M4776Svr6/mmd4tt9xCkiRs374dgEMPPXS0HgCuv/56kiRhx44do/VWeQa7dOlSZs+ePRrfkiVLmDt37pg6P+mkk0Z/vuyyyzj11FNrnn0ec8wxo3UJcNFFF7HvvvuOJqWFCxcyMDAwuvzII4/MPItdsGAB559//mj5mTNnsnDhwtHp+fPnc9FFF41O77PPPmP+/+fNm8ell14K1D5jnjNnzmj5uXPnsmTJktHp2bNns3Tp0tHpJEky970kSbjpppuA0r5X3g8Atm3bRpIk3H777QA8+uijJEnC448/Drx431u/fv2Yfa+8/sp9L0mS0X3vt7/9LRs3bnzRvrdmzRqgtO9VxlPe99avXw+8sO8988wzQGnfS5KERx99FIDbb7+dJEnYtq00IsHWrVvH7Hs33XQTSZJMuO+VTbTvXXnllVX3vbJLL72UefPmjU5v3LiRhx564QUTF110EfPnzx+dXrhwIT//+c9Hp88//3wWLFgwOn3eeedx3nnnjU6P3/cGBgbG7HsPP/zwmH3vzDPP5OKLLx6dfvDBB/nFL34xOn3qqady2WWXjU6fdNJJXHnllaPTje579cibKDYBa6J8ytgc+6XrLduczqu3DACSBiSNSBrZunXri/7I58+fjyQ++clPjn5ncHCQadOmje64kpg+ffqYs6wLLriAWbNmsXz5cgC2bNnCrFmzxhzcNmzYwHHHHYdUGh22/AdV68wq6wzt8MMPH3NgOuqoo8bEc+211zI8PMx9990HwHHHHcfw8DCXX375aJm+vj62bNkyeoY/a9YsynVSdskllyBp9A9l2rRpSBpNkGW33norw8PDozENDAywYsWKmr/b8ccfX/N3iwgOPPBA4IWzz2OPPbZq+VqOPPJIPvzhD4+exZ577rmce+65Na9e7r9/4jfhl///q8V8yimncNBBB43G/O53v3vMwWndunUsWLBgtH5WrFjBwMBA3WfVg4ODTJ8+fcy+OGPGjJpn/K3uZTQ4OMhRRx01Jp4kSWpu1w+u9bjyH2zWB3gbpRHuFgJ/Vf7k+W7GOk8Dvlwx/X7g8+PKfAM4umJ6BTBzonXPnDkzKpV+zWxFlZk1a1ZHxZO3TD3lJrueRYsWBaUn/8d8Fi1aNKl48tR1kWXyaGY8s2bNqrqsnnpudFv1lClas+IuskwrASNR45ia94riU8BTlMal2KPi04jNwP4V068FHplEmcIsWrRo0t9t142/RmJulVox1XOvY6L1TFTX9fx/DA0NdWQ9ZhkaGqr5cFs331vIeyPfmitvonhFlHo8LYqIxeVPg9teCRwi6UBJuwLzgOXjyiwH/iLt/fR2YHtEbBm/oonk+SPPUybPH1KrD4L1lmkk5smWm0ire4blORDWe7DME/NEzS95D3K11tOukw03K1neRHG3pHc1c8MR8RzwEeBblF5ZfnNE/EzSOZLKz2jcAaynNEjSUuDDk9lWs7qaNmtbzVpP0TEXecbZjKRU5OsihoaGct03yFOmVm+lVlwJ5EkCzXofUlFPb/uqo/nyvutpAXCBpN8DOyi9GDCiwZ5PEXEHpWRQOe9LFT9Huu2e0m3NGO3QrKQ0UV0PDQ0VmgAn6vpa9LuQinwpnp/e7l55XzPe6P0Iq9ANbcG9osgrsyK5OciKlLfpycwmqRX3Frrl9dh+ers3OFGYtVg39zJq1FT+3XtJ3nsUZtYhOm1Mh06Lx5qv7kQhaSAirm5FMGa9rpPuLeQ5wDczCXTS7271mUzT05R7vbhZM+TtQturuuW+ir3YZBKFmh6FmVmHmsrJvWwyieLdTY/CzKxD+fmPSSSKiNjcikDMzLpVr191uHusmdk49T7/0etXHU4UZmbj+PmPsfIOhXqapD3Snz8h6ZbJjl9tZtaIZjXzZL2KPY+p9NR53iuK/x0RT0o6GjgeuI5Jjl9tZlZN3re+Ft3Mk/Xa92ZedTSauFopb6J4Pv33T4EvRsRtwK6tCcnMrHP4+Y/8ieKXkv4BOB24Q9JL6/iumVlDOr2ZZ6Knzjv5aiGPvAf70ykNMDQnIn4DvAL4n60Kysx6SzPuB9TTzFN0d9Vev+rImyj+ISJuiYj/C5AOR/r+1oVlZjZ5vd5dtWh5E8WbKick7QTMbH44ZmbZpvLLBds1zGtmopC0UNKTwJslPZF+ngQeA24rJEIzswpZY4p38n2MbpaZKCLiknQY1L+LiD3Tzx4R8cqIWFhQjGZmE/JDcq2Td8zshZL2Aw6o/E5E3NOqwMzMrDPkShSSLgXmAQ/ywjMVAThRmFnHmcr3MVoh7wh37wXeEBG/b2UwZja1Natba693Vy1a3l5P64FdWhmImZm7tXamvFcUTwGrJa0ARq8qIuIvWxKVmZl1jLxXFMuBi4F7gVUVHzOzhrhba+fL2+vpulYHYmZT0+DgIIODgyRJwvDwMBHR7pBsnLy9nv6dUi+nMSLioKZHZGZmHSXvPYr+ip93A06j9GJAM7OmcbfWzpTrHkVE/Kri88uI+BxwbGtDM7Opxt1aO1PepqfKYU9fQukKY4+WRGRmZh0lb9PTZyt+fg7YQGmMCjOzwnTz4D/dLG+vp2NaHYiZmXWmXPcoJO0l6XJJI+nns5L2anVwZmbWfnkfuLsGeJJSc9PpwBPAta0KyszMOkfeexQHR8SpFdOLJa1uQTxmZtZh8l5RPC3p6PKEpKOAp1sTkpmZdZK8VxQfAq6ruC/xa+DsyW5U0iuAm4A+0h5UEfHrKuU2UGryeh54LiL6x5cxM7PWytvraTXwFkl7ptNPNLjdC4EVEXGppAvT6Y/VKHtMRGxrcHtmNkW4C23z5X3g7m+Bz0TEb9LplwN/HRGfmOR2TwaS9OfrgCFqJwozmwI68QDfiTG1Q957FHPLSQIgbSY6oYHt7hMRW9J1bQFeVaNcAHdJWiVpIGuFkgbK3Xe3bt3aQGhmZlYp7z2KnSS9tDwUqqTdgZdmfUHS3cC+VRZ9vI74joqIRyS9Cvi2pIcjouo43RFxNXA1QH9/v99TbGbWJHkTxfXACknXUjrL/++UmoxqiojZtZZJ+g9Jr46ILZJeDTxWYx2PpP8+JulrwOFA1URhZmatkfftsZ8B/gb4I+BNwMXpvMlaDpyV/nwWcNv4ApKmSdqj/DPwLmBNA9s0M7NJyLyikKRIh5uKiDuBO7PK1OFS4GZJHwB+QWl8CyS9BvhyRJwA7AN8TVI5zn9JYzAzswJN1PT0XUlfBW6LiF+UZ0raFTia0tXAd4Fl9Ww0In4F/EmV+Y+Q3iSPiPXAW+pZr5lZ0aZCz6iJEsUcSvcjbpB0IPAbYHdKTVZ3AVekz1iYmVmPykwUEfEMsARYImkXYG/g6cqusmZm1tvy9noiInYAW1oYi5mZdaC8D9yZmVkH2LBhQ+HbdKIwM+siGzduLHybThRmZpbJicLMrMMNDg4iieHhYQAkIYnBwcFCtq/6n5XrfP39/TEyMtLuMMzMmipJEoaHh2nFcVvSqlpj/viKwszMMjlRmJl1kQMOOKDwbTpRmJl1kb6+vsK36URhZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmlsmJwszMMjlRmJlZJicKMzPL5ERhZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmlsmJwszMMjlRmJlZJicKMzPL5ERhZmaZnCjMzCyTE4WZmWVyojAzs0xOFGZmlsmJwszMMjlRmJlZJicKMzPL5ERhZmaZnCjMzCyTE4WZmWVqS6KQdJqkn0n6g6T+jHJzJK2VtE7ShUXGaGbWrQYHB5tSpqxdVxRrgFOAe2oVkLQT8AVgLnAocIakQ4sJz8ysey1evLgpZcp2biSYyYqIhwAkZRU7HFgXEevTsjcCJwMPtjxAMzMb1cn3KPYDNlVMb07nVSVpQNKIpJGtW7e2PDgzs6INDQ0xNDRUddng4CCSRk/Ayz9XNjHlKVONIqIZ8b94xdLdwL5VFn08Im5LywwB50fESJXvnwYcHxH/I51+P3B4RHx0om339/fHyMiLVmlmNiVIYqJj+/gyklZFRNV7xi1reoqI2Q2uYjOwf8X0a4FHGlynmZnVqZObnlYCh0g6UNKuwDxgeZtjMjPreIsWLWpKmbKWNT1lblR6L/B5YAbwG2B1RBwv6TXAlyPihLTcCcDngJ2AayLiU3nW76YnM7P6tKXpKUtEfA34WpX5jwAnVEzfAdxRYGhmZjZOJzc9mZlZB3CiMDOzTE4UZmaWyYnCzMwytaXXU6tJehJYW9Dm9gK2F/T9PGWzytRaVm1+nnl7A9smiKdZiqznPOUnuzxPvVYr06t17X26uO9PVP6QiNir6pKI6LkPMFLgtq4u6vt5ymaVqbWs2vw883q1nvOUn+zynPVarUxP1rX36eK+38g+7aanxt1e4PfzlM0qU2tZtfl55xWlyHrOU36yy/PUazvruRnb9z6dT9fs073a9DQSNR4cseZxPRfHdV0M13N1vXpFcXW7A5giXM/FcV0Xw/VcRU9eUZiZWfP06hWFmZk1iROFmZllcqIwM7NMPZ8oJL1H0lJJt0l6V7vj6WWS/kjSlyT9m6QPtTueXiZpmqRVkk5sdyy9TFIi6Xvpfp20O5526cpEIekaSY9JWjNu/hxJayWtk3QhQETcGhEfBM4G/lsbwu1qddb1QxFxDnA64C6GdainnlMfA24uNsreUGddB/BbYDdKo25OSV2ZKIBlwJzKGZJ2Ar4AzAUOBc6QdGhFkU+ky60+y6ijriWdBHwfWFFsmF1vGTnrWdJs4EHgP4oOskcsI/8+/b2ImEspMS8uOM6O0ZWJIiLuAR4fN/twYF1ErI+IZ4EbgZNV8mngmxFxf9Gxdrt66jotvzwijgT+vNhIu1ud9XwM8Hbgz4APSurKv+N2qaeuI+IP6fJfAy8tMMyO0pYR7lpkP2BTxfRm4Ajgo8BsYC9J/ykivtSO4HpM1bpO23BPofQH5ZEJG1e1niPiIwCSzga2VRzMbPJq7dOnAMcD04G/b0NcHaGXEoWqzIuIuAq4quhgelytuh4ChooNpadVrefRHyKWFRdKz6u1T98C3FJ0MJ2mly5ZNwP7V0y/FnikTbH0Otd1MVzPxXFdZ+ilRLESOETSgZJ2BeYBy9scU69yXRfD9Vwc13WGrkwUkm4AfgC8QdJmSR+IiOeAjwDfAh4Cbo6In7Uzzl7gui6G67k4ruv6+aWAZmaWqSuvKMzMrDhOFGZmlsmJwszMMjlRmJlZJicKMzPL5ERhZmaZnChsypD0vKTVFZ8LJ/5W66UvrvyOpD3T6ZD0TxXLd5a0VdLXM9YxTdKvJO01bv6tkk6XdKKkKfv2U2uME4VNJU9HxFsrPpc2ukJJzXhf2gnATyLiiXT6d8B/lrR7On0c8MusFUTE74C7gPdUxLYXcDTwdeAbwEmSXtaEeG2KcaKwKU/SBkmLJd0v6aeS3pjOn5YOcrNS0o8lnZzOP1vSv0q6HbhL0ssk3SzpAUk3SfqRpH5JH5B0RcV2Pijp8ioh/Dlw27h53wT+NP35DOCGivVUjSstM69iHe8F7oyIp6L0ZO0Q4BHxrG5OFDaV7D6u6alyxMNtEXEY8EXg/HTex4HvRMTbKI0B8XeSpqXL3gGcFRHHAh8Gfh0RbwYuBmamZW6kdBa/Szo9H7i2SlxHAavGzbsRmCdpN+DNwI8qltWK605gpqRXpuXmUZFggBHgnTVrx6yGXnrNuNlEno6It9ZYVn6V9CpKY2oAvIvSgb6cOHYDXpf+/O2IKA9+czRwJUBErJH0QPrz7yR9BzhR0kPALhHx0yrbfkVEPFk5IyIekNRH6Wpi/NgeVeOKiIckLQfeJ+mrwFspNUeVPQa8psbvb1aTE4VZye/Tf5/nhb8LAadGxNrKgpKOoHQfgYpytXwZ+F/Aw1S/mgB4TtJLqgxAtBy4DEiAV1bMrxpX6gZKw/4KuC0idlQs2w14OiNWs6rc9GRW27eAj0oSgKQ/rlHu+8DpaZlDgf9SXhARP6I0zsGfMbYZqNJa4KAq868BPlnlKiQrru8ChwALqmzv9cCaGjGY1eREYVPJ+HsUE/V6uhjYBXhA0pp0upolwIy0yeljwAPA9orlNwP/JyJ+XeP736B01TBGRGyOiCvriSu9KvkqpSuQe8Z975h0W2Z18WvGzRokaSdK9x+ekXQwsAJ4fUQ8my7/OnBFRKyo8f1XA1+JiONaGOM+wL9ExJ+0ahvWu3yPwqxxLwO+m/ZuEvChiHhW0nTgPkrPSFRNEgARsUXSUkl7VjxL0WyvA/66Reu2HucrCjMzy+R7FGZmlsmJwszMMjlRmJlZJicKMzPL5ERhZmaZnCjMzCzT/wfyo3dJZcR/awAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "like.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The output of the plot function of the like BinnedAnalysis object shows:\n", "\n", "* Top: the contribution of each of the objects in the model to the total model, and plots the data points on top.\n", "* Bottom: the residuals of the likelihood fit to the data.Notice that the fit is poor in the second to last bin.\n", "\n", "Now, check for convergence." ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [], "source": [ "#print(likeobj.getRetCode())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you get anything other than '0' here, than NewMinuit didn't converge. There are several reasons that we might not get convergance:\n", "\n", "* The culprit is usually a parameter (or parameters) in the model that reach the limits set in the XML file. If this happens the minimizer cannot reach the formal minimum, and hence cannot calculate the curvature.\n", "* Often the problem is with spectral shape parameters (PL index etc..), so simply freezing the shape of all spectral parameters to their values from the full time period (and certainly for weaker background sources) when fitting a shorter time period may solve the problem. Remember that the 3FGL catalog used a full 4 years of data and we're using a much shorter time period here\n", "* Weak background sources are more likely to cause problems, so you could consider just freezing them completely (or removing them from the model). For example a background source from the catalog that is detected at TS~=25 in 2 years could cause convergence problems in a 1-month light curve, where it will often not be detectable at all.\n", "* If there are no parameters at their limits, then increasing the overall convergence tolerance may help - try using a value of 1E-8 for the absolute tolerance.\n", "* If that doesn't help then try to systematically simplify the model. Progressively freeze all sources, starting with those at the edge of the ROI in and moving in until you get a model simple enough for the minimizer to work reliably. For example if you are using a 10 degree ROI, you could start by freezing all background sources further than 7 degrees from the source of interest, and move to 5 degrees if that doesn't solve the problem." ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Wed Dec 16 02:19:54 EST 2020\r\n" ] } ], "source": [ "!date" ] }, { "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.8" } }, "nbformat": 4, "nbformat_minor": 2 }