stack_intro.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:c0100d3b8676a720fbccd9f46689cd2f171e110c9f615bdf5580114dfc42ab1a"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Simple notebook to illustrate some basic stack functionalities"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Prerequisite"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "<UL>\n",
      "<LI> Install CernVM FS as described in : https://github.com/airnandez/lsst-cvmfs\n",
      "</UL>"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Getting some sample data and setting up the environment"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "```\n",
      "# The following should be executed in a shell script outside of the ipython notebook environment\n",
      "# The ipython notebook should be (re)launched after this initialization step\n",
      "#\n",
      "mkdir stack_intro\n",
      "# setup the LSST stack\n",
      "source /cvmfs/lsst.in2p3.fr/software/linux-x86_64/lsst-v10.1/loadLSST.bash\n",
      "# get obs_cfht package\n",
      "git clone https://github.com/lsst/obs_cfht\n",
      "# build obs_cfht\n",
      "cd obs_cfht\n",
      "setup -k -r .\n",
      "scons opt=3\n",
      "cd ..\n",
      "#\n",
      "# setup obs_cfht\n",
      "setup -j -r obs_cfht obs_cfht\n",
      "# get some sample CFHT data\n",
      "curl -o stack_intro.tar.gz  -O 'https://cernbox.cern.ch/public.php?service=files&t=2a449a45074a4498927f10874bf6296e&download'\n",
      "# untar data \n",
      "tar -xvzf stack_intro.tar.gz\n",
      "# create _parent link in output directory\n",
      "cd output\n",
      "rm _parent\n",
      "ln -s ../input _parent\n",
      "cd ..\n",
      "# setup standard stack environment and astrometry_net_data\n",
      "setup pipe_tasks\n",
      "setup -j -r astrometry_net_data astrometry_net_data\n",
      "```"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import lsst.daf.persistence as dafPersist\n",
      "\n",
      "butler = dafPersist.Butler(\"output\")\n",
      "dataid = {'visit':849375, 'filter':'r', 'ccd':12}\n",
      "\n",
      "# Read a calibrated exposure and the source catalog\n",
      "exp = butler.get('calexp', dataid, immediate=True)\n",
      "sources = butler.get('src', dataid, immediate=True)\n",
      "\n",
      "# Retrieve the metadata associated to the calibrated exposure\n",
      "md = butler.get(\"calexp_md\", dataid, immediate=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print exp"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "<lsst.afw.image.imageLib.ExposureF; proxy of <Swig Object of type 'boost::shared_ptr< lsst::afw::image::Exposure< float,lsst::afw::image::MaskPixel,lsst::afw::image::VariancePixel > > *' at 0x7fb302aff750> >\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#List the methods associated to the exposure object\u00b6\n",
      "dir(exp)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "['Factory',\n",
        " '__class__',\n",
        " '__del__',\n",
        " '__delattr__',\n",
        " '__dict__',\n",
        " '__doc__',\n",
        " '__float__',\n",
        " '__format__',\n",
        " '__getattr__',\n",
        " '__getattribute__',\n",
        " '__getitem__',\n",
        " '__hash__',\n",
        " '__init__',\n",
        " '__int__',\n",
        " '__module__',\n",
        " '__new__',\n",
        " '__reduce__',\n",
        " '__reduce_ex__',\n",
        " '__repr__',\n",
        " '__setattr__',\n",
        " '__setitem__',\n",
        " '__sizeof__',\n",
        " '__str__',\n",
        " '__subclasshook__',\n",
        " '__swig_destroy__',\n",
        " '__swig_getmethods__',\n",
        " '__swig_setmethods__',\n",
        " '__weakref__',\n",
        " '_s',\n",
        " 'census',\n",
        " 'clone',\n",
        " 'getBBox',\n",
        " 'getCalib',\n",
        " 'getDetector',\n",
        " 'getDimensions',\n",
        " 'getFilter',\n",
        " 'getHeight',\n",
        " 'getId',\n",
        " 'getInfo',\n",
        " 'getMaskedImage',\n",
        " 'getMetadata',\n",
        " 'getNextMemId',\n",
        " 'getPsf',\n",
        " 'getWcs',\n",
        " 'getWidth',\n",
        " 'getX0',\n",
        " 'getXY0',\n",
        " 'getY0',\n",
        " 'hasBeenCorrupted',\n",
        " 'hasPsf',\n",
        " 'hasWcs',\n",
        " 'init',\n",
        " 'magicSentinel',\n",
        " 'markPersistent',\n",
        " 'readFits',\n",
        " 'repr',\n",
        " 'setCalib',\n",
        " 'setCorruptionCallback',\n",
        " 'setDeleteCallback',\n",
        " 'setDeleteCallbackId',\n",
        " 'setDetector',\n",
        " 'setFilter',\n",
        " 'setMaskedImage',\n",
        " 'setMetadata',\n",
        " 'setNewCallback',\n",
        " 'setNewCallbackId',\n",
        " 'setPsf',\n",
        " 'setWcs',\n",
        " 'setXY0',\n",
        " 'swigConvert',\n",
        " 'this',\n",
        " 'writeFits']"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Look at the schema of the sources' table\n",
      "schema = sources.getSchema()\n",
      "print schema.getOrderedNames()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "('id', 'coord', 'parent', 'calib_detected', 'calib_psfCandidate', 'calib_psfUsed', 'flags_negative', 'base_GaussianCentroid_x', 'base_GaussianCentroid_y', 'base_GaussianCentroid_flag', 'base_GaussianCentroid_flag_noPeak', 'base_NaiveCentroid_x', 'base_NaiveCentroid_y', 'base_NaiveCentroid_flag', 'base_NaiveCentroid_flag_noCounts', 'base_NaiveCentroid_flag_edge', 'base_SdssCentroid_x', 'base_SdssCentroid_y', 'base_SdssCentroid_xSigma', 'base_SdssCentroid_ySigma', 'base_SdssCentroid_flag', 'base_SdssCentroid_flag_edge', 'base_SdssCentroid_flag_noSecondDerivative', 'base_SdssCentroid_flag_almostNoSecondDerivative', 'base_SdssCentroid_flag_notAtMaximum', 'base_SdssShape_xx', 'base_SdssShape_yy', 'base_SdssShape_xy', 'base_SdssShape_xxSigma', 'base_SdssShape_yySigma', 'base_SdssShape_xySigma', 'base_SdssShape_x', 'base_SdssShape_y', 'base_SdssShape_flux', 'base_SdssShape_fluxSigma', 'base_SdssShape_flux_xx_Cov', 'base_SdssShape_flux_yy_Cov', 'base_SdssShape_flux_xy_Cov', 'base_SdssShape_flag', 'base_SdssShape_flag_unweightedBad', 'base_SdssShape_flag_unweighted', 'base_SdssShape_flag_shift', 'base_SdssShape_flag_maxIter', 'base_CircularApertureFlux_0_flux', 'base_CircularApertureFlux_0_fluxSigma', 'base_CircularApertureFlux_0_flag', 'base_CircularApertureFlux_0_flag_apertureTruncated', 'base_CircularApertureFlux_0_flag_sincCoeffsTruncated', 'base_CircularApertureFlux_1_flux', 'base_CircularApertureFlux_1_fluxSigma', 'base_CircularApertureFlux_1_flag', 'base_CircularApertureFlux_1_flag_apertureTruncated', 'base_CircularApertureFlux_1_flag_sincCoeffsTruncated', 'base_CircularApertureFlux_2_flux', 'base_CircularApertureFlux_2_fluxSigma', 'base_CircularApertureFlux_2_flag', 'base_CircularApertureFlux_2_flag_apertureTruncated', 'base_CircularApertureFlux_2_flag_sincCoeffsTruncated', 'base_CircularApertureFlux_3_flux', 'base_CircularApertureFlux_3_fluxSigma', 'base_CircularApertureFlux_3_flag', 'base_CircularApertureFlux_3_flag_apertureTruncated', 'base_CircularApertureFlux_3_flag_sincCoeffsTruncated', 'base_CircularApertureFlux_4_flux', 'base_CircularApertureFlux_4_fluxSigma', 'base_CircularApertureFlux_4_flag', 'base_CircularApertureFlux_4_flag_apertureTruncated', 'base_CircularApertureFlux_5_flux', 'base_CircularApertureFlux_5_fluxSigma', 'base_CircularApertureFlux_5_flag', 'base_CircularApertureFlux_5_flag_apertureTruncated', 'base_CircularApertureFlux_6_flux', 'base_CircularApertureFlux_6_fluxSigma', 'base_CircularApertureFlux_6_flag', 'base_CircularApertureFlux_6_flag_apertureTruncated', 'base_CircularApertureFlux_7_flux', 'base_CircularApertureFlux_7_fluxSigma', 'base_CircularApertureFlux_7_flag', 'base_CircularApertureFlux_7_flag_apertureTruncated', 'base_CircularApertureFlux_8_flux', 'base_CircularApertureFlux_8_fluxSigma', 'base_CircularApertureFlux_8_flag', 'base_CircularApertureFlux_8_flag_apertureTruncated', 'base_CircularApertureFlux_9_flux', 'base_CircularApertureFlux_9_fluxSigma', 'base_CircularApertureFlux_9_flag', 'base_CircularApertureFlux_9_flag_apertureTruncated', 'base_GaussianFlux_flux', 'base_GaussianFlux_fluxSigma', 'base_GaussianFlux_flag', 'base_PixelFlags_flag', 'base_PixelFlags_flag_edge', 'base_PixelFlags_flag_interpolated', 'base_PixelFlags_flag_interpolatedCenter', 'base_PixelFlags_flag_saturated', 'base_PixelFlags_flag_saturatedCenter', 'base_PixelFlags_flag_cr', 'base_PixelFlags_flag_crCenter', 'base_PixelFlags_flag_bad', 'base_PsfFlux_flux', 'base_PsfFlux_fluxSigma', 'base_PsfFlux_flag', 'base_PsfFlux_flag_noGoodPixels', 'base_PsfFlux_flag_edge', 'base_ClassificationExtendedness_value', 'base_ClassificationExtendedness_flag', 'footprint')\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Loop over the source table and make  a simple plot"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import lsst.afw.image as afwImage\n",
      "\n",
      "# Get calibration object from the metadata\n",
      "calib = afwImage.Calib(md)\n",
      "\n",
      "mag = []\n",
      "# Loop over sources\n",
      "for src in sources :\n",
      "    # apply calibration to PSF flux to get magnitude\n",
      "    magPsf, errPsf = calib.getMagnitude(src.get('base_PsfFlux_flux'),src.get('base_PsfFlux_fluxSigma'))\n",
      "    mag.append(magPsf)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import matplotlib.pylab as plt\n",
      "\n",
      "fig, (ax1) = plt.subplots(ncols=1, figsize=(10, 10))\n",
      "\n",
      "ax1.hist(mag, bins=40)\n",
      "ax1.set_xlabel(\"PSF mag\", fontsize=20)\n",
      "ax1.tick_params(labelsize=20)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAJ0CAYAAABeGyReAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuUbGdd5+HvT05uhJCICAEdCBEIoMASRCQw5CRcRFAS\nMXIbJIKyBkREcSmjoAkMDDpyM8hldASBjAQhoiKIQpIm3BRJJCIiJBACiFxi7gmBwHnnj70bmkr1\n6a4+3f1Wdz/PWrX26V17V791dvr0J/tW1VoLAAB9fEfvAQAA7GRiDACgIzEGANCRGAMA6EiMAQB0\nJMYAADoSYwAAHe01xqrqZlX181X1lqq6sKqurarLq+o9VfXEqqpl1ju6qt5eVZeO65xfVU+vqmW/\nX1WdVFUfrKqrxu9xdlU9bF/fIADAPKu93fS1qp6c5BVJPp/k7CSfSXJ4kkckOTTJGa21n55Y5/gk\nZyS5Nskbk1ya5OFJjkry5tbaI6d8nxcmeUaSzyZ5c5IDkjw6yc2SPK219vJ9epcAAHNqpRg7NsmN\nW2tvm5h/yyQfTPJfkpzYWvvzcf5Nk1yY5JAk922tnTfOPyDJWUnuk+QxrbU3Lnmto5O8d1zvXq21\nK8b5t01ybpKDk9yptXbxurxjAIA5stfDlK21sydDbJz/xSSvGr88ZslTJya5eZLTF0NsXP6rSZ49\nfvmUiZd78jh9/mKIjetcnOTlGfaSPWHltwIAsPXsywn8X5+YJslx4/QdU5Y/J8lXktynqvafWKct\ns87fjNNj92GcAABza00xVlW7kjx+/HJpRB01Tj8xuU5r7RtJLkqyK8mR4+scnOTWSa4e97ZNunCc\n3nEt4wQAmHdr3TP2O0m+P8nbWmvvXDL/0Ax7ua6YutYwv8blsmS6t+WT5LA1jhMAYK7NHGNV9UsZ\nrnz8WJKfWfcRAQDsILtmWbiqfjHJS5N8NMkDWmuXTywyuedr0uL8xfWumJi/0vKT41n+UlAAgDnT\nWrvBPVpXvWesqn45yalJPpLk2Nbal6Ys9vFxetTkE+N5ZrdLcn2ST40DuibDPcxuUlWHT3m9O4zT\nG5yDtqi1tm0fJ598cvcxeNh2O/Fh+23dh223tR/bffstZ1UxVlXPTPLiJP+UIcQuWWbRM8fpQ6Y8\nd/8kByV5f2vt+ol1apl1fmycnrWacQIAbDUrxlhV/VaSFyT5UIZDk5fuZfE3J7kkyaOr6p5LXuPA\nJM8bv3zlxDqL9yt7VlUdtmSdI5I8Ncl1SV6z0jgBALaivZ4zVlUnJXlOkm9kuEv+L0/5OMqLWmuv\nTZLW2lVV9aQMUbZQVacnuSzDxyHdMcmbWmt/tnTl1toHqurFGS4K+OeqOiPJ/kkeleEqyqe11j6z\nb29za9q9e3fvIbBGtt3WZvttXbbd1rZTt99KH4d0cpKTM9yuYuqHgidZaK0dt3TG+BFHz8rw8UcH\nJrkgyauTnNqW+YZj+D01yV0yxN95SX6vtfb2vYxvuZcDAJgrVZU25QT+vcbYvBNjAMBWsVyM7cvH\nIQEAsI/EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA\n0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJj\nAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAj\nMQYA0JEYAwDoaFfvAQAA20NVrWm91to6j2RrEWMAwDqaNazWFnDbicOUAAAdiTEAgI7EGABAR2IM\nAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQk\nxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABA\nR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0NGu3gMAAHa2\nqppp+dbaBo2kDzEGAHQ2S1zNFm5bgcOUAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQk\nxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABA\nR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwB\nAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7E\nGABAR2IMAKAjMQYA0JEYAwDoSIwBAHS0YoxV1YlV9bKqek9VXVlVe6rq9csse8T4/HKPN+zl+5xU\nVR+sqquq6vKqOruqHrYvbw4AYN7tWsUyz05ytyRXJflckjslaSus8+EkfzFl/r9MW7iqXpjkGUk+\nm+QPkxyQ5NFJ3lpVT2utvXwV4wQA2HKqtb13VVXtTvLZ1tonq+qYJGcnOa219vgpyx6R5FNJ/qS1\n9sRVDaDq6CTvTXJhknu11q4Y5982yblJDk5yp9baxVPWbSuNHwDYHFWVlffX3GCtGdepbNXf/VWV\n1lpNzl/xMGVrbaG19snF11n3kSVPHqfPXwyx8ftenOTlGfaSPWEDvi8AQHcbdQL/91TVf6+q3xyn\nd93LssdlSOJ3THnub8bpses+QgCAObCac8bW4kHj45uqaiHJSa21zy6Zd3CSWye5qrX2xSmvc+E4\nveMGjRMAoKv13jN2TZLnJrlHksPGx+J5ZruTnFlVN16y/KHj9IpMtzj/sHUeJwDAXFjXPWOttS8n\nOWVi9nuq6sEZTtK/d5KfT3Lqen3PU0751rfbvXt3du/evV4vDQCwZgsLC1lYWFhxuRWvpvy2hYcr\nK8/KMldTrrDuzyX5oyRntNZ+epx3cIZbZlzVWjt0yjo3T/KlJF9srd1qyvOupgSAOeFqyr1b89WU\n6+iScXrw4ozW2jVJPp/kJlV1+JR17jBOP7HBYwMA6GIzY+xHxumnJuafmSGLHzJlnR8bp2dt1KAA\nAHpa1xirqnvUsI9ycv4DkvxKhv2Qp008/apx+qyqOmzJOkckeWqS65K8Zj3HCQAwL1Y8gb+qTkhy\nwvjl4qHEo6vqT8Y/f7m19mvjn1+c5PZV9f4k/z7Ou1uG+4S1JL/VWvv7pa/fWvtAVb04w8ch/XNV\nnZFk/ySPynAV5dNaa59Zy5sDAJh3q/k4pJOTnJwbnl23uAfs0621I8dln5jkJ5P8QJKbJ9kvyReS\nfCDJH7TW3reX73NShj1hd0nyjSTnJfm91trb97KOE/gBYE44gX/vljuBf6arKeeNGAOA+SHG9m4e\nrqYEAGCCGAMA6EiMAQB0JMYAADoSYwAAHYkxAICOxBgAQEdiDACgIzEGANCRGAMA6EiMAQB0JMYA\nADoSYwAAHYkxAICOxBgAQEdiDACgIzEGANCRGAMA6EiMAQB0JMYAADoSYwAAHYkxAICOxBgAQEdi\nDACgIzEGANCRGAMA6EiMAQB0JMYAADoSYwAAHYkxAICOxBgAQEdiDACgIzEGANCRGAMA6EiMAQB0\nJMYAADoSYwAAHYkxAICOxBgAQEdiDACgIzEGANCRGAMA6EiMAQB0tKv3AACA+VRVvYewI4gxAGAv\n2gzLire1cJgSAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEY\nAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAd\niTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA\n0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJj\nAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAj\nMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAA\nOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHS0q/cAAIDNUVW9h8AUYgwAdpQ2w7Li\nbTM4TAkA0NGKMVZVJ1bVy6rqPVV1ZVXtqarXr7DO0VX19qq6tKqurarzq+rpVbXs96uqk6rqg1V1\nVVVdXlVnV9XD1vKmAAC2itXsGXt2kqcmuVuSz43zlt3HWVXHJzknyf2SnJHkZUn2T/KSJKcvs84L\nk7wmyS2T/GGS05LcNclbq+qpq3kjAABbUbW292PHVbU7yWdba5+sqmOSnJ3ktNba46cse9MkFyY5\nJMl9W2vnjfMPSHJWkvskeUxr7Y1L1jk6yXvH9e7VWrtinH/bJOcmOTjJnVprF0/5fm2l8QMAg+EE\n/lnPGdvI5df2Pbbq7/6qSmvtBifirbhnrLW20Fr75OLrrLD4iUlunuT0xRAbX+OrGfawJclTJtZ5\n8jh9/mKIjetcnOTlSQ5I8oSVxgkAsBWt9wn8x43Td0x57pwkX0lyn6raf2Kdtsw6fzNOj123EQIA\nzJH1jrGjxuknJp9orX0jyUUZbqdxZJJU1cFJbp3k6tbaF6e83oXj9I7rPE4AgLmw3jF2aIa9XFcs\n8/wVGQ51Hrpk+cX5yy2fJIety+gAAOaM+4wBAHS03nfgn9zzNWlx/uVLll86f6Xlb+CUU0755p93\n796d3bt3r2KYAAAba2FhIQsLCysut+KtLb5t4eE2F2dl+VtbnJbksUke21o7feK5XRnia1eSm7TW\nrh/nfy7JrZJ8T2vtCxPr3CfJ+5K8p7V2zJTv59YWALBKbm3R15pvbTGjM8fpQ6Y8d/8kByV5/2KI\nLVmnllnnx8bpWes2QgCAObLeMfbmJJckeXRV3XNxZlUdmOR545evnFjnVeP0WVV12JJ1jshw5//r\nMtydHwBg21nNHfhPSHLC+OXhSR6c5FMZ7pqfJF9urf3akuWPzxBl12X4+KPLkjw8w+0p3tRae9SU\n7/HCJM/I8HFLZ2T4+KRHJfnOJE9rrb1imbE5TAkAq+QwZV/LHaZcTYydnOTk3PBvavHFPt1aO3Ji\nnaOTPCvDxx8dmOSCJK9Ocupy9VRVJ2XYE3aXJN9Icl6S32utvX0vYxNjALBKYqyvNcfYPBNjALB6\nYqyvzTqBHwCAGYgxAICOxBgAQEdiDACgIzEGANCRGAMA6EiMAQB0JMYAADoSYwAAHYkxAICOxBgA\nQEdiDACgo129BwAAzG740G+2AzEGAFtWm3F5ATePHKYEAOhIjAEAdCTGAAA6EmMAAB2JMQCAjsQY\nAEBHYgwAoCMxBgDQkRgDAOhIjAEAdCTGAAA6EmMAAB2JMQCAjsQYAEBHYgwAoCMxBgDQkRgDAOhI\njAEAdCTGAAA6EmMAAB2JMQCAjsQYAEBHYgwAoCMxBgDQkRgDAOhIjAEAdCTGAAA6EmMAAB2JMQCA\njsQYAEBHYgwAoCMxBgDQkRgDAOhIjAEAdCTGAAA6EmMAAB2JMQCAjsQYAEBHYgwAoCMxBgDQkRgD\nAOhIjAEAdCTGAAA6EmMAAB2JMQCAjsQYAEBHYgwAoCMxBgDQkRgDAOhIjAEAdCTGAAA6EmMAAB2J\nMQCAjsQYAEBHYgwAoCMxBgDQkRgDAOhIjAEAdCTGAAA6EmMAAB2JMQCAjsQYAEBHYgwAoCMxBgDQ\nkRgDAOhIjAEAdCTGAAA6EmMAAB2JMQCAjsQYAEBHYgwAoCMxBgDQkRgDAOhIjAEAdCTGAAA6EmMA\nAB2JMQCAjsQYAEBHYgwAoCMxBgDQkRgDAOhIjAEAdCTGAAA6EmMAAB2JMQCAjsQYAEBHYgwAoCMx\nBgDQkRgDAOhIjAEAdCTGAAA6EmMAAB2JMQCAjsQYAEBHYgwAoCMxBgDQkRgDAOhIjAEAdCTGAAA6\nEmMAAB2JMQCAjsQYAEBHYgwAoCMxBgDQkRgDAOhIjAEAdCTGAAA62tV7AABAUlW9h0AnYgwA5kab\nYVnxtl04TAkA0JEYAwDoaENirKo+XVV7lnn8xzLrHF1Vb6+qS6vq2qo6v6qeXlWCEQDYtjbynLHL\nk7x0yvyrJ2dU1fFJzkhybZI3Jrk0ycOTvCTJfZM8cuOGCQDQT7U2y8mCq3zRqk8n2dNaO3IVy940\nyYVJDkly39baeeP8A5KcleQ+SR7TWnvjlHXbRowfADbbcDXlrCfwz/o7cKO/x+aMaav+7q+qtNZu\ncOXFPBwCPDHJzZOcvhhiSdJa+2qSZ49fPqXHwAAANtpGHqY8sKoel+Q2Sa5Jcn6Sc1preyaWO26c\nvmPKa5yT5CtJ7lNV+7XWrt+w0QIAdLBRMdaSHJ7kdRPzL6qqJ7TWzlky76hx+okbvEhr36iqi5Lc\nOcmRST6+EYMFAOhlow5TvibDHq9bJrlxkrsm+T9JjkjyN1V1tyXLHpoh3q5Y5rWuyHBA+bANGisA\nQDcbsmestfbciVkfTfKUqro6ya8mOSXJI9bje51yyinf/PPu3buze/fu9XhZAIB9srCwkIWFhRWX\n25CrKZf9ZlXfl+SCJP/ZWvvucd4/Jrlnknu21v5pyjr/kuQuSe7cWvv4xHOupgRgW3A15eqX36q/\n++flaspLxunBS+YtBtZRE8umqnYluV2S65N8amOHBgCw+TY7xn5knC4NqzPH6UOmLH//JAcleb8r\nKQGA7WjdY6yq7lRVB0+Zf0SSPxi/PG3JU2/OsMfs0VV1zyXLH5jkeeOXr1zvcQIAzIN1P2esqk7J\ncJL+u5N8JslVSb4vycOSHJDkbUl+srX29SXrHJ8hyq5LcnqSyzJ8HNIdk7yptfaoZb6Xc8YA2Bac\nM7b65bfq7/7lzhnbiBi7f5InJ/nBDPcaOzhDXH04yetba6cts97RSZ6V4eOPDsxwov+rk5y6XHGJ\nMQC2CzG2+uW36u/+TYuxzSTGANguxNjql9+qv/vn5WpKAACWEGMAAB2JMQCAjsQYAEBHYgwAoCMx\nBgDQkRgDAOhIjAEAdCTGAAA6EmMAAB2JMQCAjsQYAEBHYgwAoCMxBgDQkRgDAOhoV+8BAMBWUFUz\nr9Na24CRsN2IMQBYtVniavZ4Y2dymBIAoCMxBgDQkRgDAOjIOWMAsEHWctI/O48YA4AN44R/VuYw\nJQBAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDo\nSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKCjXb0HAAD7qqp6DwHWTIwBsE20GZatGZdfXAfW\nn8OUAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IM\nAKAjMQYA0JEYAwDoaFfvAQCwtVTVzOu01jZgJLA9iDEA1mCWuJo93mAncZgSAKAjMQYA0JEYAwDo\nSIwBAHQkxgAAOhJjAAAdiTEAgI7cZwyAubOWG8vCViXGAJhTbizLzuAwJQBAR/aMAexwDglCX2IM\ngDgkCP04TAkA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwB\nAHQkxgAAOvJB4QBsuCofLg7LEWMAbII24/LijZ3DYUoAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwB\nAHQkxgAAOhJjAAAdiTEAgI7EGABAR2IMAKAjMQYA0JEYAwDoSIwBAHQkxgAAOhJjAAAdiTEAgI52\n9R4AAMAsqmqm5VtrGzSS9SHGAIAtZpa4mi3cenCYEgCgIzEGANCRGAMA6EiMAQB0JMYAADoSYwAA\nHbm1BcAcm/V+SsDWI8YA5t6s91Sa9QaXgg96EmNz6PLLL89HP/rRmde7973vnV27bFIA2Er85p5D\n5557bh760EfkoIN+YNXrXHnl3+eyyy7NoYceuoEjAwDWmxibUwcd9EO54oozV738/vuLMADYilxN\nCQDQkRgDAOhIjAEAdCTGAAA6cgL/NnLYYYfNvE5rs96PCFi01huy+rkDlhJj286sN4cE9o0brAL7\nxmFKAICO7BmDFazlUNROPAy1GZ+hOOvf67x+ruO8jgvoQ4zBqjj8uzob+RmKa/173Q5jArYzMbbD\nzfp/6Dtxj89ONK97buZ1XAD7Yq7OGauq762qV1fV56vquqq6qKpeUlWzXybIKrUZHuwss/y3sVn/\nfczjmAD2zdzsGauq70vy/iTfneQvkvxbknsneXqSh1TVfVtrl3YcIgDAupubGEvyigwh9rTW2ssX\nZ1bVi5L8SpLnJ3lKp7GxSebxJPDtwOE9gNXb7Au35uIw5bhX7EFJLloaYqOTk1yb5HFVdeNNHxwd\nzHoYymGr1fH3BLB6m/dv5lzEWJJjx+nfTT7RWrs6yfuSHJzkRzZzUP0t9B4Aa7SwsNB7COyThd4D\nYM0Weg+AfbLQewBdzEuMHTVOP7HM8xeM0ztswljmyELvAeyzqprpsV0ce+yxM7/3jXwwq4XeA2DN\nFnoPgH2ysCGvOu//Zs5LjB06Tq9Y5vnF+a6q3JJ24uGxk7Nxh1odngWYzXxfiT1PJ/CzxHXXnZ/9\n9/9sDjzw3FUtf/XV12zwiACAjTAvMba45+vQZZ5fnH/55BPb+xDMf+ZrX7tg5cW+zax/H7Mtv7a/\n740d03y+h+eMj416/bWss9HLb8b32MwxrXb77fS/p41afl++x07adtvhPUyus5rtN39/T/vSI/MS\nY/82To9a5vnFc8W+7Zyy1tp2LjEAYAeoebjnUlUdmeTCJBcluX1bMqiqOiTJf2Q4iHuL1tpX+owS\nAGD9zcUJ/K21T2W4rcXtkjx14unnJLlxktcLMQBgu5mLPWPJN/eOvT/JLZL8Zb71cUi7k3w8ydGt\ntcu6DRAAYAPMTYwlwweFJ3lukock+a4kn0/yliTPaa0td9sLAIAtay4OUy5qrX2utfbE1tqtW2sH\ntNZu11p7xnYJsao6sapeVlXvqaorq2pPVb1+hXVuVFU/X1XnVNVlVXVtVX2yqk6vqh12E9y+Zt1+\nVXVIVf1aVZ1bVZdU1VVV9a9V9ftVdZvNHPtOV1U3G3+O3lJVF44/R5eP2/KJtcxlUFV1dFW9vaou\nHdc5v6qeXlVz9W/ndjfr9quqO1TVM6vqrKr6bFV9taq+UFV/UVW7O72NHWmtP3sTr/F/x39v94xH\n0badudoztt1V1YeT3C3JVUn+PcmdkpzWWnv8MsvfJMMh22OT/FOSdye5Lsn3Jrlfkl9srb19E4ZO\nZtt+VXVQkn9I8gNJPpbkXUm+muSHk9w/w+1cjm6tfWxzRr+zVdWTk7wiw972s5N8JsnhSR6R4dY5\nZ7TWfnpineOTnJHhs3HfmOTSJA/PcNX3m1trj9y0N7DDzbr9qur0JI9M8tEk782w7e6UYfvdKMnT\nW2sv28z3sFOt5WdvYv2fyPB78OoMH4t4h/E88+2lteaxSY8M57993/jnY5LsSfK6vSz//8ZlnrTM\n87t6v6ed9Jhl+yU5aXz+76Y8d8r43B/3fk875ZHhf2geNmX+LZNcPG6PRyyZf9MkX0rylST3WDL/\ngAyflbsnyaN6v6+d8ljD9jspyd2nLH//DP9TdF2Sw3u/r53wmHXbTSzz3Um+kORPM4TcniRH9n5P\nG/Gwq30TtdYWWmufHL/c667ZqrpHksckOb219kfLvN7X13mI7MUs2y/DPyJJ8rYpz/3VOL35ugyM\nFbXWzm6t3WBbtNa+mORV45fHLHnqxAzb5/TW2nlLlv9qkmePXz5lg4bLhFm3X2vtta2186csf06G\nIwz7Jzl6g4bLEmv42VvqDzME2FOztru2bhnzctNXbuix4/QNVXVokp9I8l+S/GeSM5dEAfPprAz3\nxntoVZ3axv/NG/34OH3X5g+LKb4+MU2S48bpO6Ysf06GPWb3qar9WmvXb+TgWNG07bc3109M6WfZ\nbVdVP5vk+CTHt9Yu296ftiPG5tm9xukRSV6T5GZLnmtV9cokv9Ra27PZA2NlrbXzqupJSV6Y5CNV\ndWaSryW5Z5L7Jjk1ycs7DpEkVbUryeI5f0vDa/HTQD6RCa21b1TVRUnunOTIDLfeoYO9bL/llr9t\nkgckuSZDVNPJ3rbduJ1+P8P9Rd+62WPrwWHK+XWLcfriDHtZ7pTkJkkemOSTSX4hyW/1GRqr9HdJ\n3pThl/bTkvxqhvPOzknyBiE9F34nyfcneVtr7Z1L5h+aYc/mcldyX5HhsMlhGzs8VrDc9ruBqjog\nw3m4+yc5pW2Tq/S3sKnbbrxS+bVJrkzyS53GtunE2Pxa3DYfy3Ci8Cdaa9e21s7KcD7LniTPqKr9\nuo2QZVXV7ZJ8KMN5f0/OcPXQTZM8NMltk5xTVQ/vN0Kq6peSPCPDz9jPdB4OM5pl+1XVjZK8PsN5\nYqe31l608SNkOStsu1/JcKHFk3ZSMIux+XX5OH3rxPlGaa39c5JPZ9hTdudNHherc3KGvZu/2Vr7\no9bal1prV7fW3pEhpvfLsBueDqrqF5O8NMOtD45trV0+scjinq9Dl3mJxfmT67EJVrH9li57oySn\nZfi5e2OSx23KIJlqb9uuqu6Y5PlJXj3+Wzn1JTZ+lJtPjM2vfxuny/0jc1mG/ygP3JzhMKMfynCY\n6+zJJ8aYvjzJbarqOzd7YDtdVf1yhnP2PpLhl8GXpiy2eB7YUZNPjOe63C7DCeDb735Hc26V229x\n2f2SvCHJozIconys0wP6WcW2u0uGw8hPXHKT1z1VtSfD3rIkuWCcd/zmjXzjOYF/fr0rw+7bu04+\nMZ77cIcMv+w/vbnDYpW+liGWb5Hh/wC/adx+hyxZjk1SVc9M8oIMN1F+UGvt0mUWPTPDFc0PSXL6\nxHP3T3Jh9f3zAAAHfElEQVRQkne7knJzzbD9UlX7J/mzDDd6fW1r7QmbM0qmWeW2uyjJH2f43Tbp\nxzOc7vFnGc4nu2iDhtpH7xud7dRHhhO593bT0Bsn+VyGGxTea+K5543rvqv3+9ipj1Vsv/89Pv/O\nJPtPPPeC8bm/7/0+dtIjwwUve5J8MMlhKyx7SIabvl6X5J5L5h+Y5P3j6zyy93vaSY8Zt98BGe7x\ntyfDvaqq9/h38mOWbbeX11jINr7pq49D2kRVdUKSE8YvD0/y4AyHOd47zvtya+3Xliz/wCR/PX75\n5xk+TuLeGW6N8MUk92vuN7ZpZtl+VfVdST6Q5PYZ9l7+bYZ7U903w21Lrk3ygNbaP2zW+Heyqjop\nwy1ivpHkZRn+z3rSRa211y5Z5/gkb84QZKdnODXg4UnumORNrbVHbfS4Gcy6/arqNRnuwn9Jho/i\nmebs1tq7N2C4LLGWn71lXmchw17p27dt+HFIDlNurrtnuK/KYgG3DOeeLH7w6aeTfDPGWmvvqqof\nzvB/FQ/McNLwfyR5ZZL/2Vr7wuYMm9Gqt19r7T+r6l5JnpnhF/hJGT4T7/MZ/mH63dbaDe5hxYY5\nYpx+R5JfXmaZhQyX1CdJWmt/WVXHJHlWkp/KsFfsggxXe526UQNlqiPG6Wq33xEZfj6/K8lvT1m2\nZdjLIsY23hHjdNU/e8tomX74cluwZwwAoCNXUwIAdCTGAAA6EmMAAB2JMQCAjsQYAEBHYgwAoCMx\nBgDQkRgDAOhIjAEAdCTGgE1VVXsmHl+vqi9X1ZlV9Zgpy9+oqp5UVe+uqkur6mtV9cWqOr+q/qiq\nfmJi+d1Tvsfk4zab944B9s5nUwI9tCTPGf+8X5I7Jzk+ybFV9UOttV9NhhBL8tdJfjTDB3X/dZLP\nJdk/yQ8keWySo5K8dcr3+HSSP1nm+1+xHm8CYD34bEpgU1XVniSttXajifnHJXnn+OWRrbWLq+px\nSV6X5MNJjmmtXTWxzkFJfri19u4l83YnOSvJQmvtuI17JwDrw2FKYC601s5K8vEkleSHxtlHj9M/\nmQyxcZ2vLA2x9VRVp4yHNI+pqsdU1Yeq6pqq+nxVvaiq9h+Xe+B4CPXKqrqsql5fVTeb8nrHVtUf\nVtW/VtUVVXVtVX2kqn67qg5YZgy3qqrXVNWXxuX/qaoev+RQ7Mkb8d6BzeUwJTBPapwu7rK/ZJwe\n1WEsi56W5MeSvCXJ2RkOmf5Kku+uqr9KclqGw6evSnLfJP8tyXcleejE6/x6hvfx/gyHVQ9Mcr8k\npyTZXVUPbK3tWVy4qm6R5ANJbpPk3eN6t0ryinxrD6JDG7ANiDFgLlTVAzPEyp4k/zjO/vMkz0zy\n5Ko6JEMQnddau3gVL3m7qjplyvyzZ9yb9oAk92itfXwc57OSnJchun48yYNaa+8Zn6skf5vkIVV1\n99ba+Ute5xdaaxdNvnhVPTfJs5OcmOTPljz1ggwh9ruttd9YsvxLk3xwhvEDc06MAT3UeIitMpzA\nf1SSEzLs6Xlpa+2zSdJa+/B43tjvJ3nc+EhVXZZhb9GrW2t/vcz3uG2S356Y1zLE3iwxdupiiI1j\n+lpVvTHDBQh/tRhi43Otqk5L8sAkd0ty/pLnbhBio5dmiLEHZ4yx8RDoY5JcnuR53/YGWvvnqnpd\nkp+f4T0Ac0yMAb0snu/UMlwp+e4kf9xa+9OlC7XW3lRVb0lybIbDgD+Y4fDeCUlOqKrXtdZ+dsrr\nr9cJ/B+aMu8/xum5U577/Dj93qUzq+rgJE9P8pNJ7pjkJvnWYdkk+Z4lfz4qw2HMD7bWrpnyPd4X\nMQbbhhgDerjB1ZQrLPz1DOdJvTNJquo7kvxUklcneXxVvaW19pcbMtLpt8H4+iqe229xRlXtl+EK\nz3sl+UiSNyT5cpLrMwTZyUmWnsR/6Dj94jJjWm4+sAWJMWDLGU90f1NV3TXDIb5jk2xUjK2H4zOE\n2Gtaaz+39ImqulW+tZdw0ZXj9JbLvN5y84EtyK0tgK3s6nFae12qv9uP0z+f8twxU+Z9LMl1Se5W\nVTeZ8vz91mtgQH9iDJhb4/29HjhepTj53OFJnjR+ec7mjmxmiyfvH7t0ZlUdmeR3JxdurV2f5PQk\nh2XY87d0nbsnefzGDBPowWFKYJ79cIaT3r9QVe/N8BFHSXK7JA/LcJL7X7TWzugzvFV7a5ILkzxj\nPLT64Qy3rXhYhnuUPXrKOv8jyXFJfr2q7p3hnmO3SvLTSd6W4QKGPVPWA7YYMQbMsxcluSDfulXE\nj2YIsEsynBD/p5NXX2b9boTa9vJae3vuhgu3du34cU+/k2R3kv+a5JNJnpvkJZkSY621L1XV0Un+\nV4YbyN47yb8leUqSazPE2JWT6wFbj8+mBNhiqur5SX4jyY+21t650vLAfBNjAHOqqm7dWvv8xLy7\nZvhopOuSfE9r7WtdBgesG4cpAebXh6rqgiQfTXJNkjtkOM8sSZ4kxGB7sGcMYE5V1W9nODfsiCSH\nZPikgr9P8sLW2rxfQQqskhgDAOjIfcYAADoSYwAAHYkxAICOxBgAQEdiDACgo/8PB1stGXMClDEA\nAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f705242c9d0>"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Display image in ds9 and indicate source location"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from lsst.afw.table import Point2DKey\n",
      "import lsst.afw.display.ds9 as ds9\n",
      "\n",
      "# Get wcs from exposure and pass it to ds9\n",
      "wcs = exp.getWcs() \n",
      "ds9.mtv(exp.getMaskedImage().getImage(),frame=1,wcs=wcs)\n",
      "with ds9.Buffering() :\n",
      "    sourceCentroidKey = Point2DKey(sources.schema[\"slot_Centroid\"])\n",
      "    for source in sources:\n",
      "        sx, sy = source.get(sourceCentroidKey)\n",
      "        ds9.dot(\"+\", sx,  sy, size=5, frame=1, ctype=ds9.GREEN)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Match sources to the external astrometric reference catalog"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import lsst.afw.cameraGeom as cg\n",
      "from lsst.meas.astrom import AstrometryTask\n",
      "import lsst.afw.table as afwTable\n",
      "import lsst.afw.geom as afwGeom\n",
      "\n",
      "cam = butler.get('camera')\n",
      "ccd = cam[\"ccd12\"]\n",
      "wcs = exp.getWcs() \n",
      "\n",
      "#setup astrometry solver and load external reference catalog\n",
      "aConfig = AstrometryTask.ConfigClass()\n",
      "astrom = AstrometryTask(aConfig)\n",
      "astromRet = astrom.run(exp, sources)\n",
      "refCat = astromRet.refCat\n",
      "\n",
      "#match source and reference catalog\n",
      "match = afwTable.matchRaDec(refCat, sources, afwGeom.Angle(1./3600., afwGeom.degrees))\n",
      "print(\"Found %d matches between source and reference catalogs\"%len(match))\n",
      "\n",
      "# store matching distance in a list\n",
      "distList = []\n",
      "for m in match :\n",
      "    distList.append(afwGeom.radToMas(m.distance))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Found 55 matches between source and reference catalogs\n"
       ]
      }
     ],
     "prompt_number": 37
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Plot matching distance"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import matplotlib.pylab as plt\n",
      "\n",
      "fig, (ax2) = plt.subplots(ncols=1, figsize=(10, 10))\n",
      "\n",
      "ax2.hist(distList, bins=20)\n",
      "ax2.set_xlabel(\"Matching distance in mas\", fontsize=20)\n",
      "ax2.tick_params(labelsize=20)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAmgAAAJ0CAYAAACm+9+pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYbVddJ/zvTxIZAgYhhCmSAMZgI4EX7BdIlFTECLxo\nQAyD/SBDBJRWAtICaivc+KqoraAoSAsmQFBwYHp5kcmQgkBa02kSGhAMQ5B5CEkgkARCsvqPvSs5\nOffUvVX3nlO1qurzeZ7z7Ko9rXXqrHPOt/Zea+9qrQUAgH5812ZXAACAGxLQAAA6I6ABAHRGQAMA\n6IyABgDQGQENAKAzAhoAQGf2GNCq6qSq+rOqOruqvl5V11bVGXvZ5piq+sequqSqrqiqD1TV06tK\nGAQAWIMD9rL8N5McneTyJJ9Ncrckq17ZtqoeluR1Sa5I8rdJLklyYpIXJjk2yaP2v8oAANtb7elO\nAlW1lOQzrbVPVNVxSc5K8urW2uNmrPs9ST6e5BZJjm2tvX+cf+Mk70py/yQ/21r727k/CwCAbWSP\npx1ba8uttU+Mv9Ze9nVSkkOSvHYlnI37+FaGI3FJ8tR9rSgAwE4xz35hPzZO3zZj2XuSXJnk/lV1\n4BzLBADYduYZ0I4apxdOL2itXZPkogx93u4yxzIBALadeQa0gzMMIPjaKsu/luE06S3nWCYAwLbj\n0hcAAJ3Z22U21mPlCNnBqyxfmX/Z9IKqWn0oKQBAZ1prexs8uV/mGdD+Lcl9MvRFO39yQVUdkOTO\nSa5O8slZG+/pch/0p6qyh0vizauUme1i165d2bVr14LLZjvQVlgP7YW1Gr4DF2uepzjPHKcPnrHs\nAUlumuSc1trVcywTAGDbmWdA+4ckFyd5TFXdZ2VmVd0kye+Mv/7FHMsDANiW9niKs6oenuTh46+3\nG6fHVNUrxp+/0lp7VpK01i6vqidnCGrLVfXaJJdmuNXTDyT5+9ba3825/uxAS0tLm10FtghthfXQ\nXujJ3m719Lwkz8vunY1WTr5+qrV2l6ltjknyXzPc2ukmST6W5LQkL2qrFFZVqy2iU5vZBw0ANlNV\nLXyQwB4D2kYR0LYeAQ2AnWojAprroAEAdEZAAwDojIAGANAZAQ0AoDMCGgBAZwQ0AIDOCGgAAJ0R\n0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBnBDQAgM4IaAAAnRHQAAA6I6ABAHRGQAMA6IyA\nBgDQGQENAKAzAhoAQGcENACAzghoAACdEdAAADojoAEAdEZAAwDojIAGANAZAQ0AoDMCGgBAZwQ0\nAIDOCGgAAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBnBDQAgM4IaAAAnRHQAAA6I6AB\nAHRGQAMA6IyABgDQGQENAKAzAhoAQGcENACAzghoAACdEdAAADojoAEAdEZAAwDojIAGANAZAQ0A\noDMCGgBAZwQ0AIDOCGgAAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBnBDQAgM4IaAAA\nnRHQAAA6I6ABAHRGQAMA6IyABgDQGQENAKAzAhoAQGcENACAzghoAACdEdAAADojoAEAdEZAAwDo\njIAGANAZAQ0AoDMCGgBAZwQ0AIDOCGgAAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBn\nBDQAgM4IaAAAnRHQAAA6I6ABAHRGQAMA6IyABgDQGQENAKAzAhoAQGcENACAzghoAACdEdAAADoj\noAEAdGbuAa0Gj66qs6rqc1V1RVV9oqr+rqruN+/yAAC2m2qtzXeHVS9PcnKSi5O8cZwemeTEJAck\neVxr7a+ntmnzrgeLVVVJFv2aVbQLAHpTVWmt1ULLmOcXYFUdnuSiJF9McnRr7eKJZUtJ3pXkotba\nXae2E9C2GAENgJ1qIwLavE9x3mac/stkOEuS1tpykm8kOWTOZQIAbCvzDmgfynD07L5VdevJBVX1\ngCQ3T/JPcy4TAGBbWUQftP87yauTHJzkTUm+muSuSX4qybuTPHb66JpTnFuPU5wA7FQbcYrzgAXs\n838neUWS5yR50sT8jyd55XQ4AwDghuYa0KrqgCRnJrl/khck+fMMpzx/MMnzk/x1Vd2rtfac6W13\n7dp13c9LS0tZWlqaZ9UAAPbJ8vJylpeXN7TMeY/ifEKS05K8vrV20tSymya5MMntkxzZWrtoYplT\nnFuMU5wA7FRbcRTnD4/Ts6YXtNauTPI/xzLvNedyAQC2jXkHtG+P00NXWX6bqfUAAJgy74C2cgmN\np1TVHSYXVNVDkhyb5Mok58y5XACAbWOugwRaa/9YVW9M8vAkH6mqNyT5UoZBAj+ZodPSr7XWLp1n\nuQAA28kiroP2XUmekuTnkvxQkptluBbauUle1Frb7UK1BglsPQYJALBTbbl7ce5zJQS0LUdAA2Cn\n2oqjOAEA2E8CGgBAZwQ0AIDOCGgAAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBnBDQA\ngM4IaAAAnRHQAAA6I6ABAHRGQAMA6IyABgDQGQENAKAzAhoAQGcENACAzghoAACdEdAAADojoAEA\ndEZAAwDojIAGANAZAQ0AoDMCGgBAZwQ0AIDOCGgAAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCg\nMwIaAEBnBDQAgM4IaAAAnRHQAAA6I6ABAHRGQAMA6IyABgDQGQENAKAzAhoAQGcENACAzghoAACd\nEdAAADojoAEAdEZAAwDojIAGANCZAza7Aj2pqg0pp7W2IeUAAFuTgLabRYenjQmBAMDW5RQnAEBn\nBDQAgM4IaAAAnRHQAAA6I6ABAHRGQAMA6IyABgDQGQENAKAzAhoAQGcENACAzghoAACdEdAAADoj\noAEAdEZAAwDojIAGANAZAQ0AoDMCGgBAZwQ0AIDOCGgAAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkB\nDQCgMwIaAEBnBDQAgM4IaAAAnRHQAAA6I6ABAHRGQAMA6IyABgDQGQENAKAzAhoAQGcENACAzgho\nAACdEdAAADojoAEAdEZAAwDojIAGANAZAQ0AoDMCGgBAZwQ0AIDOCGgAAJ0R0AAAOiOgAQB0RkAD\nAOjMwgJaVT2wqt5QVV+sqquq6nNV9baqesiiygQA2A4OWMROq+oPk/xqks8keWOSi5McmuTeSY5L\n8tZFlAsAsB3MPaBV1ZMzhLNXJHlKa+07U8sXEgoBALaLaq3Nb2dVN85w1OybSY6cDmd72K7Nsx77\nqqqSLLoelR6e6/7ytwJgp6qqtNZqkWXM+2jWCUkOSXJGklZVD03yQ0muSvIvrbV/nnN5AADbzrwD\n2n8cp99KckGSu08urKr3JDmptXbxnMsFANg25j2K89Bx+qwk1yT5kSQ3T3J0knckeUCSv59zmQAA\n28q8A9rK/q5OcmJr7ZzW2hWttQ8l+ekkn01yXFXdb87lAgBsG/M+xXnZOD2/tfbpyQWttSur6u1J\nfj7DqdAb9EfbtWvXdT8vLS1laWlpzlXbOYYO/NvDdnkuBjsAbF3Ly8tZXl7e0DLnPYrziUn+Kslb\nW2sPnbH8vyX5L0l+rbX2hxPzjeKcZwkb9DyUsfYyemjfAMzHRozinPcpzjMzfNv9h5p96OOHxulF\ncy4XAGDbmGtAG09rvjnJ4UmePrmsqn4iyYOSXJrkbfMsFwBgO5nrKc4kqao7JjknyfdlOKJ2QZI7\nJ3l4hpGdj2mtvWFqG6c451mCU5zdldFD+wZgPjbiFOfcA1qSVNUhSZ6b5MQkt0/ytSRnJ3l+a+28\nGesLaPMsQUDrrowe2jcA87FlA9q6KyGgzbcEAa27Mnpo3wDMx1YcJAAAwH4S0AAAOiOgAQB0RkAD\nAOiMgAYA0BkBDQCgMwIaAEBnBDQAgM4IaAAAnRHQAAA6I6ABAHRGQAMA6IyABgDQGQENAKAzAhoA\nQGcENACAzghoAACdEdAAADojoAEAdEZAAwDojIAGANAZAQ0AoDMCGgBAZwQ0AIDOCGgAAJ0R0AAA\nOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBnBDQAgM4IaAAAnRHQAAA6I6ABAHRGQAMA6IyABgDQ\nGQENAKAzAhoAQGcENACAzghoAACdEdAAADojoAEAdEZAAwDojIAGANAZAQ0AoDMCGgBAZwQ0AIDO\nCGgAAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBnBDQAgM4IaAAAnRHQAAA6I6ABAHRG\nQAMA6IyABgDQGQENAKAzAhoAQGcENACAzghoAACdEdAAADojoAEAdEZAAwDojIAGANAZAQ0AoDMC\nGgBAZwQ0AIDOCGgAAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBnBDQAgM4IaAAAnRHQ\nAAA6I6ABAHRGQAMA6IyABgDQGQENAKAzAhoAQGcENACAzghoAACdEdAAADojoAEAdEZAAwDojIAG\nANAZAQ0AoDMCGgBAZwQ0AIDOCGgAAJ3ZkIBWVY+tqmvHx89vRJkAAFvVwgNaVX1fkj9P8o1xVlt0\nmQAAW9lCA1pVVZLTk3wlyUsXWRYAwHax6CNopyQ5PskTk1yx4LIAALaFhQW0qvrBJL+f5E9aa+9d\nVDkAANvNQgJaVR2Q5Iwkn0ryG4soAwBguzpgQft9bpJ7JTm2tfatBZUBALAtzf0IWlXdN8mvJ/lv\nrbV/mff+AQC2u7keQRtPbb4qyb8led5qq82auWvXrut+XlpaytLS0jyr1pVhcCsAsBUsLy9neXl5\nQ8us1uZ3WbKqumWSS9a4+p+21n5l3K7Nsx77aghOi66HMnZiGT20bwDmo6rSWlvo0ZZ590G7Kslf\nZfY33n2S/F9Jzs5whO2cOZcNALAtzDWgtdauSvLkWcuqaleGgPbK1tpp8ywXAGA7cbN0AIDObGRA\na3EfTgCAvZrrIIF9roRBAsrY5mX00L4BmI+NGCTgFCcAQGcENACAzghoAACdEdAAADojoAEAdEZA\nAwDojIAGANAZAQ0AoDMCGgBAZwQ0AIDOCGgAAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIa\nAEBnBDQAgM4IaAAAnRHQAAA6I6ABAHRGQAMA6IyABgDQGQENAKAzAhoAQGcENACAzghoAACdEdAA\nADojoAEAdEZAAwDojIAGANAZAQ0AoDMCGgBAZwQ0AIDOCGgAAJ0R0AAAOiOgAQB0RkADAOiMgAYA\n0BkBDQCgMwIaAEBnBDQAgM4IaAAAnRHQAAA6I6ABAHRGQAMA6IyABgDQGQENAKAzB2x2BWAnqKqF\nl9FaW+j+N+I5JIt/HgBbgYAGG2LRoWNjwtP2eR4AfXOKEwCgMwIaAEBnBDQAgM4IaAAAnRHQAAA6\nI6ABAHRGQAMA6IyABgDQGQENAKAzAhoAQGcENACAzghoAACdEdAAADojoAEAdEZAAwDojIAGANAZ\nAQ0AoDMCGgBAZwQ0AIDOCGgAAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBnBDQAgM4I\naAAAnRHQAAA6I6ABAHRGQAMA6IyABgDQGQENAKAzAhoAQGcENACAzghoAACdEdAAADojoAEAdEZA\nAwDojIAGANAZAQ0AoDMCGgBAZwQ0AIDOCGgAAJ0R0AAAOiOgAQB0RkADAOjM3ANaVd2qqp5UVW+o\nqo9X1RVVdVlVnV1VJ1dVzbtMAIDtpFpr891h1S8meUmSzyc5K8mnk9wuySOSHJzkda21R05t0+Zd\nj30xZMdF10MZylhMGYt+D23U+6OHzwKAPamqtNYWesBpEQHt+CQ3a629ZWr+bZOcm+T7kpzUWnv9\nxDIBTRnK2M8yBDSAjbERAW3upzhba2dNh7Nx/peSvHT89bh5lwsAsF1s9CCB70xNAQCYsmEBraoO\nSPK48de3bVS5AABbzQEbWNbvJ7l7kre01t65ng3PO++8fOELX1hMrQAAOrMhAa2qTknyzCQfSfJz\ns9bZtWvXdT8vLS1laWnput+f97w/zFlnfSIHHniHhdXx298+f2H7BtiKNuqqSAaG0Lvl5eUsLy9v\naJlzH8W5WwFVv5zkRUk+nOSBrbUvz1hnj6M4H/zgR+Xtbz8pyaMWVs+DDnp8vvnNV2W7jOhTxs4r\nwyhO5s1rDrNtyVGck6rqGRnC2QeTHD8rnAEAcEMLC2hV9ZwkL0hyfoZwdvGiygIA2E4WEtCq6reS\nPD/JeRlOa16yiHIAALajuQ8SqKrHJzk1yTVJ3pvkGTM6ml7UWnvlvMsGANgOFjGK84hx+l1JnrHK\nOstJBDQAgBkWcaunU1tr39Vau9E4nfX4sXmXCwCwXWz0rZ4AANgLAQ0AoDMCGgBAZwQ0AIDOCGgA\nAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBnBDQAgM4IaAAAnRHQAAA6I6ABAHRGQAMA\n6IyABgDQGQENAKAzAhoAQGcENACAzghoAACdEdAAADojoAEAdEZAAwDojIAGANAZAQ0AoDMCGgBA\nZwQ0AIDOCGgAAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBnBDQAgM4IaAAAnRHQAAA6\nI6ABAHRGQAMA6IyABgDQGQENAKAzAhoAQGcENACAzghoAACdEdAAADojoAEAdEZAAwDozAGbXQGA\nSVW12VVgg23Ea95aW3gZ28VGvQe9JnsmoAGdWfSHdiljXWVshO3yPLYTr8lmc4oTAKAzAhoAQGcE\nNACAzghoAACdEdAAADojoAEAdEZAAwDojIAGANAZAQ0AoDMCGgBAZwQ0AIDOCGgAAJ0R0AAAOiOg\nAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBnBDQAgM4IaAAAnRHQAAA6I6ABAHRGQAMA6IyABgDQGQEN\nAKAzAhoAQGcENACAzghoAACdEdAAADojoAEAdEZAAwDojIAGANAZAQ0AoDMCGgBAZwQ0AIDOCGgA\nAJ0R0AAAOiOgAQB0RkADAOiMgAYA0BkBDQCgMwIaAEBnBDQAgM4IaAAAnRHQAAA6I6ABAHRmIQGt\nqg6rqtOq6vNVdVVVXVRVL6yqWy6iPACA7eSAee+wqu6a5Jwkt0nyxiQfTXLfJE9P8uCqOra1dsm8\nywUA2C4WcQTtJRnC2dNaa49orf1Ga+2BSV6Y5Kgkv7uAMgEAto25BrTx6NkJSS5qrb14avHzklyR\n5LFVdbN5lgsAsJ3M+wja8eP0HdMLWmvfSPK+JAclud+cy2VHWd7sCrBlLG92BdhClpeXN7sKcJ15\nB7SjxumFqyz/2Dg9cs7lsqMsb3YF2DKWN7sCbCECGj2Zd0A7eJx+bZXlK/ON5gQAWMXcR3Euwo1u\nlNz0pn+QAw88Y2FlfPvb5y9s3wAA6zHvgLZyhOzgVZavzL9sekFV7XXnV175/n2r1brsvR7K6KGM\nUzegjHlafBlreQ/NoZQtWMastrIVn4cy9quENb4/Tj11np8tW1k/r8lONe+A9tFxetQqy1f6nt2g\nj1przasEADCq1tr8dlZ1lyQfT3JRku9vEzuvqlsk+UKSluTQ1tqVcysYAGAbmesggdbaJzNcYuPO\nSX5pavGpSW6W5AzhDABgdXM9gpZcdxTtnCSHJnlTrr/V01KSf0tyTGvt0rkWCgCwjcw9oCXDzdKT\n/HaSBye5dZLPJ3lDklNba6tdggMAgCzmXpxprX22tXZya+0OrbUbt9bu3Fp75mQ4q6rDquq0qvp8\nVV1VVRdV1QuryjXSdqCq+lRVXbvK4wurbHNMVf1jVV1SVVdU1Qeq6ulVtZB2zcaqqpOq6s+q6uyq\n+vrYFvZ4rZ19aRNV9fiqOreqLq+qy6rqrKp66PyfEYuynrZSVUfs4bPm2qp6zR7K0Va2gaq6VVU9\nqareUFUfHz8rLhvbz8m1yvDSjf582ZTroI337Dwnw03V35jrT4M+PcmDq+rY1tolm1E3NtVlSf5k\nxvxvTM+oqocleV2G+7v+bZJLkpyY5IVJjk3yqMVVkw3ym0mOTnJ5ks8muVuGQUYz7UubqKo/SvLM\nJJ9J8pdJbpzkMUneXFVPm3FPYfq0rrYyuiDD98+0D81aWVvZVh6V5CUZzu6dleTTSW6X5BFJXp7k\nIUkeObnBpny+tNY2/JHk7UmuTfJLU/P/eJz/F5tRL4/NeyT5VJJPrnHd70ny5SRXJrn3xPwbZ7jf\n67VJHr3Zz8ljv9vEUpK7jj8fN76ur5pXm0hyzDj/wiQHT8w/PMnF474O3+y/g8fc28oR4/LT1rF/\nbWUbPTLcN/yhM+bfNsm/j6/1Iybmb8rny4afChqPnp2Q5KK2e3p8XoZ0+tiqutlG140t46QkhyR5\nbWvtuqsXt9a+leE/6SR56mZUjPlprS231j4x/rq3ayXuS5v4xXH6u22i+0Vr7d+TvDjDh+8T97H6\nbKB1tpV9oa1sI621s1prb5kx/0tJXjr+etzEok35fNmMvjrHj9N3TC9orX0jQxo9KMn9NrJSdOEm\nVfXYqvqN8bz+0irn9n9snL5txrL3ZPjP5P5VdeDCakpv1tomvntqm7bKNm8dp8fPWMb2cMeq+oXx\n8+YXquoee1hXW9k5vjM1TTbp82Uz+qCt3GXgwlWWfyzDEbYjk7xrQ2pED1qGPgCvmpp/UVU9sbX2\nnol5q7ah1to1VXVRkh9McpcMl3Zh+1tPm/hoVR2U5A5JLh//a5728XH6A4uoLF04YXxcp6qWkzy+\ntfaZiXnayg5RVQckedz462Sw2pTPl804grZyP87VLrexMt9ozp3l9Az/cdw2wwWN75Hkv2foL/LW\nqjp6Yt2DMwS6PbWhija0k6y1TRw8sf7K/NXWT7Sh7eibGS4Dde8Mr+8tM5zOOitDX7Yzp7rYaCs7\nx+8nuXuSt7TW3jkxf1M+X1yOgC601n577EfyldbaVa21D7fWnprkBUlummTX5tYQ2A7Gz5hdrbUL\nWmtfHx9nJ/mJJP+S5PuTPGlza8lGq6pTMoy4/EiSn9vk6iTZnIC2khwPXmX5yvzLNqAu9G+lw+aP\nTsyb/m9lmja086y3Tfgc4gZaa9dkuMRCsvvnTaKtbFtV9csZLvH04STHt9amX8tN+XzZjID20XF6\n1CrLjxynq/VRY2e5eJweNDFvpV/Zbm1o7ENw5yRXJ/nkYqtGR9bVJlpr38xwDaSbV9XtZuzP59DO\ntNvnjbayvVXVM5K8KMkHM4SzL89YbVM+XzYjoJ01Tk+YvlpvVd0iwwXfvpnknze6YnRpZTTvZNg6\nc5w+eMb6D8hwSvSc1trVi6wYXdmXNnFmhv+KZ23zkHFqoNLOMuvzJtFWtqWqek6GbjTnZwhnF6+y\n6uZ8vmzSReLeluECbr88Nf8F4/yXbEa9PDbnkeGq3wfNmH9EhlG91yb5tYn5t8hw0cCrktxnYv5N\nMtyh4tokj9rs5+Ux1zaylD1ffHTdbSLJ/cf5H0tyy4n5RyT5aoZrMt5ps5+7x9zbyr0z3od6av4D\nx/ZzTZL7aSvb+5Hkt8bX9NzJ13SVdTfl82UhN0vfm6q6y/ikDk3yplx/q6elDIcSj2mtXbrhFWNT\nVNWuJP8lybsz3HLj8iR3TfLQDBfze0uSn26tfWdim4cl+YcMb5jXJrk0w203fiDJ37fWHr2BT4EF\nqKqHJ3n4+OvtMnTi/mSS947zvtJae9bE+utuExO3Yvlshtu4fHeSRyf53iRPa629ZP7PjHlbT1sZ\nL6Xx/Rm+gz43Lj86wzWpWpLfaq393owytJVtoqoen+HKAdck+bMkX5+x2kWttVdObLPxny+bmF4P\nS3JahvO030pyUYYjaAdvVp08Nq0tPCDJ32QYPXNpkm8n+VKGW4I9dg/bHZMhvF2S4b+RD2S4n+tu\n/x17bL1HhjuLXDt+iE4+rh0fu90abF/aRJLHZ/gv+hsZOveeleT/2ezn77GYtpLk5CRvHr9zLh+/\ncD+V5DVJjt1LOdrKNnhMtZdrV3m8a8Z2G/r5silH0AAAWJ3roAEAdEZAAwDojIAGANAZAQ0AoDMC\nGgBAZwQ0AIDOCGgAAJ0R0AAAOiOgwRZQVbuq6tqqesA6tlmuqmsXWa95qqojxud4+tT8V4zz77RZ\ndduqxr/bWZtdD2D9BDR2hPGL6tqquma8F+xq6501se7j97PMJ8xjP/uhjY+tZrrO+/U8tlpQXYCt\n2AZgxxPQ2Em+k6SS/PyshVV1ZJLjxvWS+X2xbdYX5OOS/OAmlT1Pv57heXx+P/axU0PK3TK0A2CL\nEdDYSb6U5LwkT6yqG81Y/qRx+uY5l1tz3t+atNY+01q7cDPKnqfW2hdbaxe21r6z97WZNP7dPrvZ\n9QDWT0BjJ2lJXpbkdkl+cnJBVR2Y5AlJ3pfkX2dtXFX3qao/raoPVNVXq+rKqrqwqv6oqm45te5y\nktPGX0+fOG16g75UVXWjqvrFqnpfVX2tqq6oqo9V1cuq6vtnV6NOqqpzq+qbYz1eU1V3mLHibqf2\nqmpprMPzqupeVfWWqrps3NdyVd1/led++6o6vaq+PNbx/Kp63OT+Zm23yr5uUVUvqKrPjn/Dj1TV\nr2SVz6PV+qBV1YlVdWZVfaGqrqqqz43P4anj8iPG5/+A8e82+RqcNbGf46vqL6vqXydegw9W1XOr\n6sYz6rPSH/C4tb4W43a3qqrfraoPjetfVlUXVNXzq+pmM9Z9/vi3uWJc95+q6oS1/p3H/ezWB21f\n67/K/q87jV9VJ1TV2VV1eVV9papOq6qDx/XuPba1S8flb6qqw2fsb83vsXH9766qU6rq/VV1yfg8\nLqqqN1bVA9fzt4LeHLDZFYAN9pokL8hwtOxNE/NPTHKbJM9K8gOrbPvkJA9PspzkHRkCxQ8neWaS\nh1TVfVtr3xjXPT3JpUkeluSNSS6Y2M/XkuHLJcn/n+THk3w6yauTfD3Jncdyzk7y8ak6/Oexrm9K\nclaS+yV5dJJ7VtW9Wmvfnlp/tVN7P5zk2UnOSfKXSQ5P8jNJzhz3c92Rt6o6NMn/SHKnJO8et7l9\nkpckeefRQ6IMAAAI9ElEQVReyrmBMfCcOZZ/QZIzknxvkt9KsrSHTW+w/6p6SpKXJvlChr/FxUkO\nTXLPDEH7LzL8/U8dfz88ya6JXXxq4udnJzlqfF5vTnKTJD8yrr9UVT/eWpvVh23Nr0VV3Xlc504Z\njuK+JEP7OSrJM8b6fnpc9/AMbezwJO9J8o9Jbp7hn4q3VdUvtNZevoe/1bTVXpv1tqU9OXGs35vH\n53Jshr/7Xavq15P80/icXpbk6CQ/leQuVXV0a22yfut5jyXJK5I8JskHk7wyyZVJ7jiW/6AMbQ22\nptaah8e2fyS5Nsmnx59fluTqJHecWP62DF/oN0nyO+P6j5vax52S1Ix9nzyu/+yp+U+YtZ+J5b83\nLn9jkgOnlh2Y5JCJ33eN616W5O5T6/71uOyRU/OXk1wzNW9pXHfW83vKOP/FU/P/apz//Kn5Rye5\nalz23DW+Dr8xrv/3U/OPSPLVcdlpU8teMc6/08S8/5Xhy/iQGWXcam9/h6nld15l/m+P5T5qav6+\nvBbnjPOfM6u+SW48Vd/vzCj34CTnJ7kiyaHraPfv2t/672H/K23820l+dGJ+ZQhYK+X87NR2Lx+X\nnbiv77Hx73FtknNX2eZWa3kOHh69PpziZCd6WZIbZfjQXzlicUKSv26tXbXaRq21T7fWZh2NOD3J\n5Ul+Yq0VqKEP3H/O8GX7i621q6fKurq1dvGMTV/UWvvwjOeTJP9xreUneW9r7VVT805Lcs3kfsaj\nfD+b4Uv2d6bq+L+TTO9jb544lvHsqX19KsmL1rmva3L9gI7JfV2ynp201i5aZdGfjNPVXtc1vRZV\ndZ8MR6fOb639wYzyL2mtfWtc954ZTsm+rrX2d1PrfS1DuLpJhqOd+2tebSlJXtNaO3vll/F9csb4\n6wWttddMrb/Sbu45OXOd77GV9b41a5v1tgPojVOc7DittXOr6oNJTq6q38lwurNy/ZfTTDX0U/uF\nDKdU/kOS78kN+03dcR3VuNu4/T+31r64ju3OmzFvpRP49+7Pflpr36mqL03t56gMgeDc1to3Z+zn\nfbl+cMUeVdUtktw1w5HMWaHo3WvZz+jVSf44yb9W1WsznAp8X2vtK+vYx0q9Dkry9CQ/neH09s1z\nw4Edq72ua30t7jdO376G6qz0AbxlVe2asfw243Qeo3Pn1ZZW29cXxun/mrFsZUTuYZMz1/Mea619\nvarenOSnquqCJK/L0C3g3NbaFeusP3RHQGOnelmGIzYPyXBU57zW2gf2ss3fZugf84kkb0jyxSTf\nyvBl/owku3Uo34OVDs+fW8c2yXAka9rKUaRZI1PXs5+VfU3u5+Bx+qVV1l9t/ix729eag2pr7YVV\ndXGGo5CnZPj7t6p6d5JntdZmhYLdjIHgXRmOGH0wQx/Fr2Q4BV5JnpfVX9e1vhbrea1vPU5PGB+z\ntCQHrWFfezOvtpSM/SpX2deelh04NX+977FHJ3lOkv+Uob9hklxVVf+Q5Fdba19ex3OArgho7FRn\nJPmDJP89yR1yww7ku6mqH87wxfHOJA9pE53Gq6oyfEmsx8qX43qOum2Gr4/T266yfLX5s6x8Ua+2\nze3Wsa+01s5IcsY4UvCYDEfATk7y9qq62yqniKc9LEM4O721doPr41XV7TMEtP218loftse1Bit/\no1Naa38+h7K3jH15j41dEk5NcmpVHZbh9PATkjw2Q7/GNd95A3qjDxo70tif5x8yBKRvZDhysicr\nl7z4/9ruI/rum+E04LRrxumsoxEfyfBlfM8xCPTqIxkGAhxdVTefsfxH1rqj1trlGUalHlaz7+aw\ntC8VbK19rbX21tbaUzIMKLhVbvjFfE1y3Zf8tJXX9fUzlh23L/WZ4X+M0wetUodZ6+7EYLEv77Hr\ntNY+21r7mwyjNz+R5Eeqar2naqEbAho72W9m+I/9Qav0r5q00mfq+MmZ4yUoXrzKNl8dp7td72n8\nAnpxkpsmeenYGX9yv99dVYfspU4LNw5eeG2G03S/Obls7NC+3qvUn57hc+cPJsPKeBmKU9a6k6o6\nfpVFK0fnJl/Pr2Y4Rbbb65DVX9e7ZDjCut9aa+/PMIrzXplxFKiqbj1efiTjqdmzkzyiqp44a39V\ndY+qus2sZVvcut5jVXVIVd1jxn5uPj6uzjC6FLYkpzjZsVprn0nymTWu/j8zdIh/RFW9b/z5tkke\nnOSjGTo9Tx8dOSfDKM1nVNWtc33fqxe11r6e4dTMfTNcE+rCqnpLhpFq35eh/9GvZv2jJKfN4y4G\nv5bkx5I8u6rum+Eoz+2TPDLJWzKE3LXe6/KPx/V/Jsn7q+odGcLfIzN09D9xjft5Q1VdnuSfk/x7\nhuf5oxmumXVehuturfinJCcleX1VvTXD5Tk+1Vp7dYbrdn08yTPHL/sLMlzq4aEZrlH3mDXWZ28e\nm+HyGb9XVT+TYUBEJTkyw2t9VMbroGXoT/WuJH9VVadkuIzEZRlOkR6d5O4ZBh6se0BE59b7Hjss\nQxv6YIb+g5/JMKjgJ8ft/nQN/3hBtxxBg93tdnPu8YjXiRkuwnmHJE/L0O/pZRm+QK6esc1lGYLI\nv2boF/PbGULZLcflV4/bPi1DeHtckl/OEDJen+S9e6rTvjyPNW5zwxlDR+tjMoTFu2forH3PJE/N\ncN2s5Pq+anve+XDx0x9P8sIMIxJPyRCs/t8kv7KHOk3X6zkZvtDvPdbjCRlOJT87yfGttWsm1n15\nkudnGKTwrAyvw8ljfa7IED7/ZnxuT0vyQ+M6P7eO+uzReBmReyf5wyS3SPJLGQanHJbkjzIRtlpr\nn0tynyT/NcPp2f801ut+GS6w+5QkH1pP+ftb/43Y1z68xy7K0Efw4gynx38l1w8w+NnW2mrtCbaE\nmn3JGYC9q6rfzXAz8we11t65t/UBWBsBDdirqrpDa+3zU/PukeE07lUZ7sqgvw/AnOiDBqzFeVX1\nsSQfztAB/8gM/bSS5MnCGcB8OYIG7FVVPTdD/54jMvShujRDB/0/aq29ZxOrBrAtCWgAAJ0xihMA\noDMCGgBAZwQ0AIDOCGgAAJ0R0AAAOvN/ANK/Lvog4zyzAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f70531aff10>"
       ]
      }
     ],
     "prompt_number": 41
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}