3_systematics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Joint spectrum fit with systematics\n",
    "\n",
    "In this notebook, we re-execute the joint spectrum fit including systematics\n",
    "on the per-instrument energy scales, as described in the paper and scripted\n",
    "in ``joint_crab/fit_systematics.py``. Please see the paper for a description\n",
    "of the likelihood function and parameters.\n",
    "\n",
    "First we illustrate the effect of the energy scale on the spectral model.\n",
    "After the fit, we show the parameter correlations, and we show a plot\n",
    "illustrating the energy scale shifts for each instrument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import astropy.units as u\n",
    "import matplotlib.pyplot as plt\n",
    "from multinorm import MultiNorm\n",
    "\n",
    "from joint_crab.fit_models import Log10ParabolaEnergyScale\n",
    "from joint_crab.fit_systematics import SystematicsSpectrumFit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Energy scale shift\n",
    "\n",
    "Let us see how the log-parabolic spectrum is modified by an energy scale shift of $\\pm20\\%$.\n",
    "\n",
    "For negative energy scale factor `z` the spectrum is shifted to lower energies,\n",
    "for positive `z` to higher energies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEOCAYAAAC0BAELAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzs3XdcldUfwPHPuewtWxRRERWcKOTeI8X9c1vulSOttDTLLFempaW5E2eO1DJ37j1w7z1QcKGoIMjm/P64ZCZwGXFlnffrdV917z3Pc7/3UfnynPE9QkqJoiiKomQ1TXYHoCiKouRNKsEoiqIoeqESjKIoiqIXKsEoiqIoeqESjKIoiqIXKsEoiqIoeqESjKIoiqIXKsEoiqIoeqESjKIoiqIXKsEoiqIoemGY3QFkByFES6CllZVVv1KlSmXqHJGRkVhYWGRtYPmUupZZQ13HrKGuY9pOnjz5RErpmFY7kZ9rkfn6+soTJ05k6ti9e/dSr169rA0on1LXMmuo65g11HVMmxDipJTSN612qotMURRF0QuVYBRFURS9yJcJRgjRUggxPywsLLtDURRFybPy5SC/lHIjsNHX17dfdseiKErOERcXh6WlJZcvX87uUHIEU1NTXF1dMTIyytTx+TLBKIqipCQ4OBhnZ2dcXV0RQmR3ONlKSkloaCjBwcEUL148U+fIl11kiqIoKYmOjsbGxibfJxcAIQT29vZER0dn+hzqDkZRdJBSEhkXSVhsGBGxEUTFRxGdEE1sQiyJMhEpJUIIjDRGGBsYY2xgjJWRFZbGllgbW2NqaJrdX0HJIJVc/vFfr4VKMEq+lpCYwL2IewSGBxL0IoigF0E8iHhAyMsQQl6G8DT6KfEyPtPntzCywMHMAXtTewpZFsLFwgVXK1eKWRejuE1xbE1ts/DbKAqcPHmSnj17EhUVRbNmzZg+fXqyRLF8+XImT54MgKWlJXPmzKFixYpZHku+TDB/r+T38PDI7lCUtygmIYarT69yMfQiF59c5Nqza9wKu0VMQsyrNmaGZhS2LIyTuRMeth7Ym9pTwKQANiY2WBlbYWZohqmhKcYaYzRCAwKQEJcYR2xCLNEJ0UTERhARF0F4bDihUaGERoUSEhXCqUenePTyEQky4dXn2ZrYUsquFJ62nnjae1LBoQJFrIqo36KVTBs4cCDz58+nWrVqNGvWjL/++gs/P79/tSlevDj79u3D1taWrVu30r9/fwICArI8lnyZYNQssvwhMi6Sk49OcuLRCU4/Os3F0IvEJcYBYGdqh6edJ50KdsLD3IXiRja4Gphhn5CIiI2A+GiIi4L4eIiPgcgQEE/A0AQMTcHIDExttA8zW7BwAqO0u8PiE+N5EPmAwLBAbofd5mbYTa4+vcrKKyuJTYx9FVsFhwr4FvTlnYLvUNq2NAYaA71eKyVnmDt3LnPnzgUgLCyMYsWKsWfPnnQf/+DBA8LDw6levToA3bt3588//0yWYGrUqPHq/6tVq0ZwcHAWRJ9cvkwwSt4kpeTqs6vsC9rH4fuHOff4HPEyHkONIWWti9PV4R0qSkPKRkfjHHYf8eAIvFgHr93B/CemBcDKBQq4gW1RsC0GDqXBsRRYu4JGg6HGkCJWRShiVYTarrVfHRqXGMet57c49+QcZ0POcubxGfYG7wXA2tia6oWqU6twLWoVroWDmUPWxKvoNHbjRS7dD8/Sc5YpZM3XLcum+v6AAQMYMGAAcXFxNGjQgGHDhtGpUyeuXr2arO2wYcPo3r37v167d+8erq6ur567urpy7949nTH5+/snS0BZRSUYJVdLSEzg5KOTbL+znX3B+3gY+RAALxMHuhs6Uj3sCd6Pb2OaeOufg6wKaX/4F6kCVgW1ScHCUZsgzAqAidU/dykaI/i7uyoxQZuM4mMgNhKiw7SPl6EQEQIRjyD8Pjy/C3cOQ+yLfz7T2BKcy4FLRSjkDa5VwL7Eq3MbaYwobVea0nal6VCqAwCPIh9x7OExAh4EcPj+YbYFbgOggmMFGrk1opFbI4pYF9H7NVbevo8++ogGDRrQsmVLWrZsme7jUqotqau7dc+ePfj7+3Pw4MFMxZkWlWCUXEdKybkn59h4cyM77+wgNPopphhQI9GAQU+fUftlJA4Jd7V3DS4VoXQ7cPIER0+wK5GurqwsCFKbeB5fhSdXIeQKPDwHp3+FY/O0bczsoEhVKF4bitcBp7Kg+WflgLOFMy1LtKRliZb/ujvbHbSbaSenMe3kNLzsvGju3pymxZrq/zvlM7ruNPRp8eLF3Llzh5kzZwLovIN5//338fHxAaBVq1YMHDjwX91dwcHBFCpUKMXPOXfuHH379mXr1q3Y29vr4ZuoBKPkIo8iH/HnjT/ZeO137rx8gIkU1Hn5kiYREdSOTcS8kA+yQkue2lfitJkXQbGWPHgexZOIGJ5eiuPZy5eER50iMjaBl7HxxMQlkiAlCYna3/oMNQIjAw0mhhosTAyxMDHA2tQIe0sT7C2McbQyoVABM1xsTClia46NuY7VzUKAhYP2UazmP68nJsCT6xB8DIIC4M4RuLZV+565PZRoCCXfBY+GYG732ukEnnaeeNp58kHFD7gXcY+dd3ay9fZWfjjxA1NPTKW0aWmibkfRwK0BJgYm+vgjUPTs5MmT/PDDDxw4cABN0i8bv/32m85jzpw586/nVlZWHD16lKpVq7J06VKGDBmS7Ji7d+/Stm1bli1bRma3LEkPlWCUHC1RJnIw+ABrzi1g/5MzJAK+UdH0iYikrmkxntu/yxknb0ZFuXPlSTx3bkcSHZcIXH91DhNDDfYWxthZGmNtakQBc2MsTAwwNtBgaKDBIOmmIT5BEpcgiY5PIDImnpcxCdwJfcmpu895GhlD4hu9D3YWxhR3sMDdwQJPF2u8XKwo42JNAXPj1L+QxkB7N+XkCZWT+s/D7kHgAbi5B27shPOrQWigaE0o0xo8W4C1y79OU9iyMD3K9qBH2R4EhgWy+fZmfrvwGyP2j8DGxIZWJVrRqXQniloX/e9/CMpbM3PmTJ4+fUr9+vUB8PX1ZcGCBRk6x5w5c15NU/bz83s1vvL35IEBAwYwbtw4QkNDGTRoEACGhoZkdusSXXLtfjBCCHfgS8BGStk+6bU2QHPACZglpdyu6xxqP5icIaVr+SL2BX+ensfK62sJSojEPj6BNpFRVBEeXJU1WPXMk6tR1oD2ZqGonTkeTpYUs7egmIMFRezMcbExpaCNKdammauj9LrEREloZCwPwqK4/zyKu09fcvvJS24/ieD6owhCI2NftS1mb453kQJ4FylAleL2eBa0QqNJ57TjxAS4fxquboXLG7XdawgoVgvKd9AmHLMCKR66e89uTEub8sf1P9h1ZxfxMp6ahWryntd71CpcSzutWtHp8uXLuLq6YmVlld2h5BiXL1/Gy8vrX6+ldz+YbLmDEUIsBFoAIVLKcq+93hSYDhgAC6SU36V2DinlLaCPEGLta6/9CfwphLAFfgB0Jhgl53kYeo1lhyewNvQ0LwV4R8fQ4qUTN0JrsjDehxnCglLOVniXLUD3IjaULWRDKWdLzI31+1dZoxE4WpngaGVCBdfkP+BDXkRz5cELLtwP42zQc47cCuXPM/cBKGBuRJVidtQu5Ujdko642Zvr+CADcPXVPhp+pR3Dufin9q5m41DY8qn2jqZydyhe919jNhqhoUahGtQoVIPHLx+z9vpa1l5dy+Bdg3G3cad7me60KNFCdZ8pb012dZEtBmYCS/9+QQhhAMwCGgPBwHEhxAa0yWbSG8f3llKG6Dj/6KRzKblE4LVNzD82ja3xIUigeqQG49DKHIp6l3DXotSoZc/04nb4FLXNkjuSrOZkZYqTlSl1Sv2zi2zws5cE3HpKwO1QDt8MZfulRwC4O1jQqIwz75ZxppKbLQa67m4cS0O9kVB3hPbO5uwqOPcbXPxDOx3apxdU7gEW/x6kdTR3ZGDFgfQt15e/Av9i6aWlfHPkG2admUWPsj3oUKoD5kY6Ep2iZIFs6yITQhQDNv19ByOEqA58I6VskvR8FICU8s3k8uZ51r7WRSaA74AdUsqdacWgusiyWXwsV47NY87ZBew1S8BIQrkwWyIiW1G05LvU93SkloeD7jGNXEJKye0nkey79pjdV0I4eiuUuASJg6UxfuVcaOVdCB832/R1pcVFw5VNcHKxduzGwATKteOE4Tv4tuyd6ucHPAxgwbkFBDwMwNbElu5lu/Oe53sq0bxGdZEl91+6yHJSgmkPNJVS9k163g2oKqX8MJXj7YGJaO94FkgpJwkhhgI9gOPAGSnl3BSO6w/0B3B2dvZZtWpVpuKPiIjA0tIyU8fmdyIukue31nEg/iA7LAwxkuAeXhRboy5ULuhCSVsNmjxeKuVlnOTckwROPornTEgCcYlgZyqoUciQWoUNKWiRvvES88i7FL63hYIP92CQGM1T24oEFWnDM9tK/6zfecOt6FtsC9vGpehLWGmsaGzTmFpWtTASOe/O8G2zsbGhePHiGBioygl/u3HjBm9uzli/fv1cl2A6AE3eSDBVpJTJ59hlEXUH83Y9eHCPC5u/5ejLbfxhbYqQgoqaSnhb/I8hbdukfyA8j4mIiWfX5Uf8efoe+649JlGCb1FbOldxo0UFF0yN0vHDLuo5N9d8RYmQ7RDxEApWgDqfacdrNCknqzMhZ5h5ZiYBDwIoaFGQIZWG0MK9Rb6eDKDuYJL7L3cwOelvUjDw+rJkV+C+Pj5IbZn89kgpCbhwjQ0/9mP9ylqMMdjLWmszqllVZ13bbSzqsZSKDunsGsqjLE0Mae1dmEW9qnB0VEM+9/Pk6ctYPl1zlmqTdjFh0yXuhr7UfRKzAgS5tYOPz0PrWdpKA6u7wZwacOEPSExMdoi3kzcL3l2A/7v+2Jva8+XBL+m4sSNH7h/R0zdV8puclGCOAyWFEMWFEMZAZ2CDPj5ISrlRStnfxsZGH6dXgJj4BP44fJllkwfyeHMjZlsdZJa9JWXtK7GuzXrmtPuFojYprzDOz5ysTRlQtwS7htVlRd+q1Chhz+LDgdT7YQ8fLDvB8cCnKZYDecXQGCp1hQ+PQ9sFIBNhbS+YXweubdNWGHhDFZcqrGi+gil1phARF0H/Hf0ZunsoQeFBevymSlaJiYmhU6dOeHh4ULVqVQIDA7M7pFeya5rySqAe4CCECAa+llL6CyE+BLahnTm2UEp5UU+fr8r160lETDyrjtzg+f55NNKsZY69CYfNbShhWYR51UdTo1CNtE+iIISghocDNTwceBgWzdIjgSwPuMu2i4+o7FaAwfU9aODplHqdKY0BVOgA5drC+bWw91tY0RHcqsO7E7TToF9vLjT4FfejoVtDll5ayvxz82m9vjU9y/akX4V+mBma6f9LK5ni7++Pra0tN27cYNWqVYwcOTLN1f9vS65daJkV1BhM1gmPjmPhgVsEHlrNALmM7XbRLLKxwdTQjA99PqJT6U4YalL+fUZdy/R5GRvP2pPBzNt3i3vPo/AsaMWQBiXxK1cQjUbovo4JcXB6GeyZpN16oOz/oNE32qKfKQh5GcKPJ39k061NFLYszKgqo6hbpK6evlnOkd1jMJkp19+kSRO++eYbqlevTnx8PAULFuTx48dZtqdQrltomd3UHUzWeREdx8KDgew7uJdPExZSxeIWnzg5EWxgQkv3FgzzHa7Ky2cRc2NDulcvRpcqbmw4c5/Ze28weMUpPAta8XGjUpjo+mXRwAh8e2urARz+Wfu4sgVqDoVan4Cxxb+aO5k7Man2JNqWbMuEoxP4cPeHNHRryBdVv8DJ3EnP3zSH2Po5PDyftecsWB78Ul0/nqly/ffu3aNIEe3wtaGhITY2NoSGhuLgkP3/7vJlglEbjv13MfEJ/Hr0Lkt2n6VX7Ar8jXfyvbMjG82dKGrlxoLqY6jqUjW7w8yTjAw0tPNxpU2lwmw6d5/pO68z4NeTFLfRYOL2hBoldPxgMbGC+l+AT0/Y8TXs/x7OrIAmE6FMm2RTm98p+A5rW65lyaUlzD07lzZ/tmG473Dalmyrdt3Uo4yU689oif63KV8mGHUHk3lSSjacvc+UrVeo9GIP602Xc9o6mrYFi/OMBPqV680HFT9Q5UjeAgONoLV3YZqXd2Hd6XtM2nSe934JoF5pR0b5eVG6oI5uHutC0O4X7V3N1s9gTU/waATNvgc79381NTIwom/5vjQu2phvDn/DN0e+YWvgVsbVGEchyzw8UUPHnYY+ZaRcf/fu3XF1dSUoKAhXV1fi4+MJCwvDzs4uWfvskC8TjLqDyZxTd58xftMlHgddY7rFUkqZnuY71xJsNIildIHizK45Hi97r7RPpGQpQwMNHXyLYBN2g0AjN2btuYnf9P28V9WNYY1LY2ehoxJC0erQby8cXwC7J8Ds6tr1MzU/0narvd7Uuij+TfxZe20tU09Mpe2Gtox8ZyRtPNrkmN+Yc7vMlOtv1aoVS5YsoXr16qxdu5YGDRrkmD+PnDRNWcmhHr+IYdjqM7SffZBaT9aw13wU8aa3aFfCiy2GCQyoOICVzVeq5JLNjA0E/euUYN9n9ehevRgrjwVR9/s9LDp0m/iE5OtgXjEwhGoD4MNj2r1odo+HX+rD/TPJmmqEho6lO/JH6z8oY1+GMYfHMGT3EEKjQvX4zfKP18v1e3t707dv3zSP6dOnD6GhoXh4eDBt2jS++y577rxSki9nkb3WRdbv+vXrabZPSX6Y+RSfkMivR+8wdfs1XOKDWGi7BOeIc0x3r8gS+Yxi1sX4tta3lHcs/58+Jz9cy7fhzet4I+QFYzde4sD1J3i5WDOhTTl8itqmfaLLG2HzcIh8or2Tqfc5GCbv8kyUiay4vIIfT/6IpbElE2pOoLZr7Sz8Rm9fds8iy4nyykr+t0YttEzbhXthtJl9iLEbL/CZ7V7+MvuSeBnM+2WqsEQ+o1PpTqxuufo/JxdFfzycrFjauwpz3q/M85extJtzmM9/P0fYyzjdB3q1hMEBULEzHJwGvzRIcTaVRmjoWqYrq1qswt7MnkG7BjEpYBIxCTF6+kZKbpMvE4ySupex8UzcfIlWMw+S+Pw+x9xm0f35bDa7VaCjiwP3E6P4qf5PjK42Wi2+ywWEEPiVd2HnsLr0r+POmpPBNJy2jy3nH+iuCGBmC21mQ5dVEBEC8+vDganaDdHeUNK2JCubr6SrV1dWXFlBty3duBN+R4/fSsktVIJRXjl6K5SmPx3glwO3GVfqNpuMRmD17BzfVGrGF/IhZRzKsrblWhq6NczuUJUMsjAx5ItmXqwfXJOCNiYMWn6KD5ad5PGLNO42SvvBoKPg2Qx2jYMlrSAsOFkzEwMTRlYZyc8NfuZ+5H06buzIlltb9PRtlNwiXyYYVezy317GxvP1+gt0nn8UUxlNQPkNdL3zJUH2brzv5cPvzy/Qr3w/Fry7gIIWBbM7XOU/KFfYhj8H1WSUnyd7rz3m3R/3seHsfd13Mxb20GEJtJ6t3fRsTk24tD7FpvWK1GNty7WUsi3FyAMj+TbgW+IS0uiSU/KsfJlg1BjMP07ffUaz6QdYcuQOn1VK5C/zMThf/409Pp3obBHHo9hnzG44m6GVh6Za6kXJXQwNNHxQtwRbhtbGzd6CoStPM3jFKZ5FxqZ+kBBQ6X0YcEC7TmZ1d9j8qXbzszcUtCjIwqYL6VGmByuvrKTntp48jHyox2+k5FT5MsEoEJeQyI87rtF+7hHiEiQ7G95j8I3+yOhn/FynL0OfHsHNuiirW6zO9TODlJR5OFny+4DqfNakNDsuPaLp9P0cvP5E90H2JaD3Nqj+IRz/BfwbQ+jNZM2MNEZ8+s6nTK07lRvPbtBpUyeOPzyup2+i5FQqweRDQU9f0nHeEabvuk678vbsKbUWj0Of8aJwJYaUr8v8oG208WjDUr+leXultoKhgYbB9T1YN6gmliaGdPUPYMKmS8TG61g3Y2isLS3TZRU8vwvz68GVzSk2fbfYu6xssRJrY2v6b+/PqiurdHfHKRmWnnL9QUFB1K9fHy8vL8qWLcv06dPfSmz5MsHk5zGYTefu02z6AW48imBBayemhI/A+Nxy7lb/gK4FjDgScorRVUczrsY4Ve4lHylX2IZNQ2rTrVpRFhy8Tfu5h7kTGqn7oNJ+2i4z+xKw6j3YOTbFWWbuNu6saL6CGoVrMDFgImOPjFXjMlno9XL9n3zyCSNHjkzWxtDQkKlTp3L58mWOHj3KrFmzuHTpkt5jy5cJJj+OwUTHJfDFuvN8uOI0JZws2fU/SaP9HeHpLY76jafL00M8jX7KvMbz6OTZKceUmlDeHjNjA8a3Kcfcrj4EPomk+YyDbDibxqayBdyg119QuYd2zcyvbeHl02TNrIytmFF/Bv3K9+P367/Tf0d/nkc/19M3yb3mzp2Lt7c33t7eFC9enPr166d5zPr16+nRowcA7du3Z9euXcnuEl1cXKhcuTIAVlZWeHl5ce/evaz/Am9Qo7b5wN3QlwxcfpKL98P5oE5xRhTYg8H60eBQmt9r9GTC+XkUsynGjAYzKGJVJO0TKnla03IFKVfYmo9WnWHoytOcCHzKl829MDE0SPkAI1NoNUO7idnm4doyM51XgnOZfzUz0BgwtPJQ3Au48/Whr3lvy3vMbDgTdxv3lM+bzSYfm8yVp1ey9Jyedp6MrJL8DuNvb6Ncf2BgIKdPn6ZqVf1XO1cJJo/bcekRw1afQQALu1agwY1JsH05iaWbM72ENwvPzaJmoZr8UPcHLI0tsztcJYdwtTVnVf9qTN56hQUHb3M2OIzZ71emcAEdi2srdwdHT/itKyxoBG3naasCvKGFewtcLV35aM9HdN3clZ/q/0QVlyp6/Da5j77K9UdERNCuXTt++uknrK2tsyRWXVSCyaMSEyU/7bzGjN03KF/Yhrn/c6Pwtj4QFEBMnc/4QvOM7ZeX0bFUR0ZVHaWmICvJGBloGN2iDD5Fbfls7Tla/nyQme9V0r3fTJEq0H8vrHpfm2gajoFaw5LtM+Pt5M2K5isYvHMwH+z8gPE1x9PCvYVev09G6brT0Cd9leuPi4ujXbt2vP/++7Rt21bv3wNUgsmTwqLiGPbbGXZdCaGDjysTahpisro5RIQQ9r85DL2/jVMhpxjuM5weZXuo8RZFJ7/yLpQuaEX/ZSfp5n+ML5p50btmsdT/3lgXgl5bYP1g7er/0JvQ4ift7LPXFLYszBK/JXyy9xNGHRjFg4gH9C3fN1//fdRXuX4pJX369MHLy4thw4bpLf435ctB/rzs1uMI/jfrEPuuPWZ867JMqfwUkyVNIS6aB52X0eP2as4/Oc+UOlPoWa5nvv7HrKSfu6Ml6wbVoKGnE+M3XWL4mrNExyWfMfaKkRm084e6n8OZ5bCsTYqD/zYmNsxtNJdmxZsx4/QMJgZMJCGFmWj5RVaW679//z7NmjUD4NChQyxbtozdu3e/mkSwZYv+S/nkyzuYvLqj5cHrTxi0/CSGBhqW961K1fBtsHwIOJTiRovv+SDga6LiopjXeB7vFHwnu8NVchkrUyPmdvXh5903+HHnNQKfRDKvmy+OVqlMZxcC6o8Cew9YPwgWNoH314Jt0X81MzYw5rva3+Fs7syii4t4Gv2U72p/h7GBjo3S8qhFixZl+BhTU1PWrFmT7PVChQq9SiK1atXKlvVH+fIOJi9OU156JJAei47hYmPG+kE1qBq0EP4cCEVrcrbVD/Q4NBIpJYv9FqvkomSaRiP4qFFJ5rxfmUsPwmkz6xCX7ofrPqhCB+i2DiIeaVf+p7CRmRCCYb7D+NT3U3bc2cHAnQOJiI3Q07dQ3pZ8mWDykoREydiNFxmz/iL1Szvy+4AqFDn8BeyZABU6c7jBcPrt+wQbExuW+i2llG2p7A5ZyQP8yruwdkANEhIl7eceZs+VEN0HFKsFvbeDgTEsagY3d6fYrEfZHnxb61tOPTpF3+191VqZXC5TCUYIYSGESGVSvPK2vIyN54NlJ1l0KJDeNYszr0s5LNf3gZOLofZwdlZuz+C9H+Nm5cZSv6W4Wrlmd8hKHlKusA3rP6yJu6MFfZYc59ejaewB4+QJfXaAXXFY3hHOr02xWcsSLfmp/k9cf3adnn/1JORlGslLybHSlWCEEBohxHtCiM1CiBDgCvBACHFRCPG9EKKkfsNU3vT4RQyd5x9l95VHjG1VljGNXTFY0QGubIKmk9lYrBKf7v+MsvZlWdh0IQ5mOqaWKkomOVub8lv/6tQr7cToPy8wactlEhN19PVbu0DPzeD6DvzeFwLmpdisbpG6zGk0hweRD+i+tTvBL5LvQaMvqlbaP/7rtUjvHcweoAQwCigopSwipXQCagNHge+EEF3/UyRKugU+iaT93MNce/SCX7r70qOiJSxpAXePQNtfWGNnz5cHv8TH2Yf5jedjbaz/BVVK/mVhYsj8bj50rebGvP23GLb6jO5imWYFoNsfULoZbB0B+6ZACj/IqrhUYcG7C3gR+4Je23pxN/yuHr+FlqmpKWFhYSrJoE0uoaGhmJqaZvoc6Z1F1khKmaw6nZTyKfA78LsQwijTUSjpdi74Ob0WHSdRSlb2q0alAtGwqCU8vwOdV7IyMZRvj3xL7cK1mVZvGqaGmf/LoSjpZWigYXzrcrjYmPH9tquERsYyp6sPliap/IgxMoOOS2HDENgzEWLCofH4ZAsyyzuWx7+JP/2396fnXz1Z0GSBXkvLuLq6cvbsWSIi1AQD0CZcV9fMd62nK8GklFwy0yYrCSHcgS8BGyll+6TXvICPAAdgl5RyztuMSd8O3XhCv6UnsLMwZmnvKrgbPoFFrSHyCXT9nV9fBjL5+GTqF6nPD3V/yJfTPJXsI4RgcH0PnKxM+PyP83Sef4Qlvapgb5nKNGYDQ2g9C0ws4fDPEPMCmk8Dzb+Hdz3tPPFv4k+/7f3o9Vcv/N/1x8NWP0sMjIyMiIiIwNfXVy/nz28yPMgvhPjP9ROEEAuFECFCiAtvvN5UCHFVCHFDCPG5rnNIKW9JKfu88dplKeUAoCOQp/6G/HXhAb0WHcfNzpw/BtbAXfMIFvpB1HPdcQ2rAAAgAElEQVTovoGlkbeYfHwyjdwaMbXuVJVclGzTwbcIv3T34UZIBB3mHeH+86jUG2s04DcFag/XTk75c1CKJf9L2pZkYdOFGAgD+mzvw83nyTc5U3KeNBOMEGL1a481QNpLS9O2GGj6xucYALMAP6AM0EUIUUYIUV4IsemNh5OOeFsBB4FdWRBnjrD6eBCDlp+iXGFrfutfHaeYu9qpngkx0HMzv4Zf5vsT39O4aGOm1J2CkYHqrVSyVwNPZ5b2rsrj8BjazznMrcc6upyE0NYsq/8lnFsFf/SDhPhkzdxt3PFv4o+BMKD3tt4qyeQC6bmDCZdSdkx6dAB2/tcPlVLuB96sG1EFuJF0ZxILrAJaSynPSylbvPFIdd6ilHKDlLIG8P5/jTMnWHToNiN+P0etko782rcqNhE3YHEzkInQczMrn59/decyuc5kjDQquSg5Q5XidqzsX42Y+EQ6zD2S9oLMuiOg0Tdw4XdY2wtS2JSsuE1xFjRZgEZo6LOtD7ee39JL7ErWEGnNlhBCFJdS3hZCOEgpnwgh7JIG9//bBwtRDNgkpSyX9Lw90FRK2TfpeTegqpTyw1SOtwcmAo2BBVLKSUKIekBbwAQ4J6WclcJx/YH+AM7Ozj6rVq3KVPwRERFYWuq3vP2mW7GsvRaHj7MBAyqaYBMVhPeZ0Uih4WzFCexIvMOqp6sob1ae3o69MRS5s/LP27iW+UFOvY4PIhL5/kQ0MQmST31NKW6jewmda9B6PG4u5LFDdS6V+RSZQqXvh3EPmfFwBkIIPnb+GEcjxyyLN6dex5ykfv36J6WUaQ9DSCnT9QA2pLdtOs9XDLjw2vMOaBPF38+7AT9n5We+du6WwHwPDw+ZWXv27Mn0sWlJTEyUU7ddkUVHbpJDV56ScfEJUoZckXJKCSm/LyXl4+tyw40NstzicnLgjoEyJj5Gb7G8Dfq8lvlJTr6Od0MjZc3vdslyY/6SJwJD0z7g8Cwpv7aW8rduUsbHptjk2tNrstbKWrLxmsby3ot7WRZrTr6OOQVwQqbjZ21GBvn1XXY3GHh9O0VXII39WjNH5uBaZFJKvt92lRm7b9DR15VpHb0xfHYTlrQEoYGem9gReZvRh0ZTtWBVfqz/oxrQV3K8InbmrBlQHQcrE7r5HyPgVqjuA6oPgibfwqX1qY7JlLQtybzG84iIjaDv9r5qxX8OlJEEo++VR8eBkkKI4kIIY6AzsEEfHySEaCmEmB8WFqaP02ealJLJf11l9t6bvFfVje/aVsDg2S1Y3EI75tJjIwdjHjFi/wgqOFRgRoMZmBikMgVUUXIYFxszfutfDRcbU3otPp6OJDNYuzbm4jptNebE5Is3y9iXYW7juYRGhfLBjg9U7bIcJlvuYIQQK4EjQGkhRLAQoo+UMh74ENgGXAZWSykvZtVnvi4n3sFIKfnuryvM3XeTrtXcmNC6HJrwYFjaGhJiofsGTiVG8vGejylZoCSzGs3C3Mg8u8NWlAxxsjZlZf9qFCpglr4kU3MoNBgN536DzZ+kuOK/gmMFfm7wM3fD7zJo1yAi4yL1FL2SURlJMKOy6kOllF2klC5SSiMppauU0j/p9S1SylJSyhJSyolZ9Xlvyml3MH93i83bd4tu1YoyvnU5NJGPYGkriA6Hbuu4amTAh7s+xMXChbmN56ryL0qu5WRlyop+VSlUwIyei45zPDCNOUN1Pvtnncy2L1ItK/ND3R+4FHqJobuHEpMQo5/glQxJd4KRUl5Iu1XukNPuYKbvuv6qW2xc67KIqGewtA28eARd13LX0o4PdnyAuZE58xvPx840+X7bipKbOFmZsrJfNVwKmNJr0XFO332m+4AGX0G1QXB0Nuz5NsUm9d3qM77meI49PMbI/SPz9c6YOUWGVvILIToIIayS/n+0EOIPIURl/YSWP8zac4Ofdl6ng48rE1qXQ8RGwvL28PQWdFnJEwd3+u/oT6JMZP6783GxdMnukBUlSzhambCibzXsLY3pvvAY54N19CgIoR30r9QN9k+BI7NTbNayREs+r/I5u+7uYkLABFW0MptltFTMV1LKF0KIWkATYAmQ6+p95ZQusoUHb/P9tqv8r1JhvmtXAU1iLKx6D+6fhvYLiXD1YeDOgTyNfsrsRrP1WuRPUbJDQRtTVvSrhrWpEd0WBnDloY7FmEJAy+ng1Qq2jYLTy1Ns9r7X+/Qt35e119Yy52yu+/GUp2Q0wfx9z9kcmCOlXA/kujmyOaGLbPXxIMZtukTTsgX5vn0FDEjU7o9xex+0mklcqSZ8vPdjbjy7wbR60yjnUC7bYlUUfSpcwIyV/aphYqihm/8xAp/oGKTXGEC7BeBeDzZ8CFc2p9hsaKWh/M/jf8w5O4fVV1frJW4lbRlNMPeEEPPQFpPcIoQwycQ58r1N5+7z+R/nqFPKkeldvDHUCNjyKVzeAE2+JdG7C18e+pKABwGMrTmWWoVrZXfIiqJXbvbm/NqnKvEJiby/IIAHYToKZBqaQKflUKgSrO0Nd48mayKEYEz1MdRxrcPEgInsDdqrv+CVVGU0OXREO424qZTyOWAHfJblUelZdnaR7b0awserzuBT1JZ5XX0wMTSA/d/DiYVQ8yOoPpjpp6az9fZWPqr8Ea1KtHrrMSpKdijpbMXS3lUJi4qj64IAQiN0zAQzsYT31oCNK6zoCCGXkzUx1BjyfZ3v8bLzYsT+EVx4kmfmKeUaGUowUsqXUso/pJTXk54/kFJu109o+pNdXWSn7j5j4K+nKOVshX/PdzAzNoCTS7QbLlXsAo3GsvrqahZeWEiHUh3oU65P2idVlDykvKsN/j18CX4WRa/Fx4mISb6C/xULe+j6Bxiawa/tICz5tsrmRubMbDgTO1M7Bu8aTFB4kB6jV96kurfekuuPXtB78XGcrE1Y0rsK1qZGcG07bPoYPBpBq585cO8g3wZod6P8ouoXCKHv6jyKkvNUdbdn1nuVuXg/nAHLThITr2O6sW1R6Pq7drOy5R20+yO9wcHMgbmN5pIgExi0axBhMTlj/Vt+kJ79YBoLIX4RQngnPe+v/7D06213kd1/HkU3/2MYGWhY1rsqjlYm2plia3pCwfLQYQlXw27x6b5PKWVbih/q/oBhChVkFSW/aFTGme/alufgjScMX32WhEQd040LloNOy+DJdfitK8THJmtSzKYY0+tP517EPT7Z+wlxKWwFoGS99NzBDEI7ztJVCNEA8NZvSPr3NrvIwl7G0WPhMSJj4lnauwpu9ubw7A6s6ATm9vDeGh4nRDF412AsjS2Z2XCmKgGjKGh3xhzl58mmcw+YsPmS7jUt7vWg9UwIPKCdXZZCWx9nH8bWGMvxh8cZe2SsWiPzFqTn1+THSQP6nwohvgPe0XNMeUZ0XAL9lp7gTuhLFvd+By8Xa4h6pr2Vj4+G7huIMrNm6F+9CY8NZ6nfUpzMU92sU1Hynf513HkYHs2iQ4EUsjGjXx0da8EqdoawINg9AQoUhQZfJmvSskRLgl4EMefsHIpaF6VfhX56jF5JT4J5NdFcSvm5EGKIHuPJMxISJcNWn+FY4FNmdKlEjRIO2h36VnfXrtLvto5Ex1KM3vcZF0MvMqPBDDztPLM7bEXJUYQQfNW8DCHhMUzcchlnG1NaVSyU+gG1P9X2EOyfAnbu4N0lWZOBFQdyJ/wOM07PwN3GnYZFG+rxG+RvaXaRJS2mRAjhkPT8Z30HlRd8u+UyW84/5MtmXtp/EFLC5mFwez+0mgHFazPn7By239nOcN/h1CtSL7tDVpQcSaMRTO1YkSrF7fh09VmO6qrALAQ0nwbF68CGIRB4KIUmgnE1x1HBoQKjDo7icmjyKc5K1sjILLKFeoviLdP3IP+Sw4H4H7xNzxrF/rmlP/wznFqqrQrr/R7bArcx9+xc2ni0oXuZ7nqJQ1HyClMjA37p5oubvTkfLDvJzccRqTc2NIaOS8G2GPz2PoTeTNbExMCE6Q2mY21szZDdQ3gS9UR/wedjOWlHy7dGn4P8uy4/YuzGizTycuKrFmW0L17ZAjvGQJk2UH80l0IvMfrgaLwdvfmq2ldqOrKipIONuRGLer6DoUbQa9Fx3QsxzWzh/dWA0E6oSWX68s8NfiY8NpyP93xMbELy2WfKf5OTdrTM9S7cC+PDFacpW8iGGV0qYaAR8OiidsvXQt7QZg5PYp4ydPdQCpgWUNsdK0oGFbEzZ0EPXx6FR9N36Qmi43SskbFzh06/wrNAbUmZFLZd9rL3YkLNCZx9fJaJARPVzLIsli/vYPThUXg0fZecwNbcCP8evpgbG0LkE1jZGYwtofMK4gyMGL53OGExYcyoPwMHM4fsDltRcp1Kbrb81Mmb03ef89nac7qTQrGa0GIa3NwF20en2OTdYu/Sv0J//rj+B6uurtJT1PlTtuxomddExSbQd8kJwqPj8O/5Dk7WptrFXqu7Q0QIdF4B1oWYfHwyp0JOMbbGWLzsvbI7bEXJtfzKuzCiaWk2nr3PjF03dDeu3F27WVnAHG1pphQM9h5MPdd6TD42mevR1/UQcf6UL3e0zEqJiZLha85w4X4YMzpX0q51Afjrc7hzCFrNBFcf/rj+B79d/Y2eZXvSzL1Z9gatKHnAwLolaFu5MD/uvMamc/d1N248Hko0hM3D4W5Asrc1QsOk2pNws3Zj4eOFPIx8qKeo85eM7mjpK4RYJ4Q4JYQ4J4Q4L4Q4p6/g9CUrZ5H9tOs6W84/ZJSfJ43KOGtfPLkYTvhDzY+hQgfOPz7PhKMTqO5SnY8qf/SfP1NRFO1040lty+Nb1Jbhq89yNij5QP4rBobQ3l9bfXl1NwhPnpAsjS35qf5PxMt4PtnzCTEJOiYRKOmS0WKXy4FFQDugJdAi6b+5SlbNItty/gEzdl2nvY8r/WonTUe+GwCbP9X+ttRwDE+jn/LJ3k9wMndiSp0pqsaYomQhE0MD5nXzwdHKhP7LThASHp16YzNb6LISYiO1Ncvikrd1t3Gnq0NXLoReYFLAJD1Gnj9kNME8llJukFLellLe+fuhl8hyuDvhCQxffZZKbgWY+L9y2qnG4Q+0vx3ZuEJ7f+KRjNg3gucxz5lWbxoFTAtkd9iKkufYW5owv5sv4VHxDPg1jerLTl7wv7lw7yRsGZ5izbKK5hXpV74fv1//nTXX1ugx8rwvownmayHEAiFEFyFE278feoksBwuNiGHGqRhszIz+2TTs70H9mAjtb0lmtsw8PZOAhwGMrjaaMvZlsjtsRcmzyhSyZlrHipy6+5zR6y7onlnm1RLqjIDTv2q7s1Mw2HswNQvVZFLAJC4+uaifoPOBjCaYXmirKTdF2zX2dzdZvrLrcggvYiXzu/toZ4wBbPsCgo9Bm1ng5MXuu7vxv+BPh1IdaOPRJnsDVpR8wK+8C0MbeLDmZDBLDgfqblzvc+0+TFs+g6Djyd420BgwqfYk7M3sGbZ3mNpDJpMymmAqSil9pZQ9pJS9kh699RJZDtbxnSJMqm1GBdekLq+zq+D4L1D9Qyj7P4JeBDH64GjK2pfl8yqfZ2+wipKPfNyoFI28nBm/+TIBumqWaQyg7S9gXeif5QRvsDW1ZVrdaYREhTDqwCgSZaIeI8+bMppgjgohVF8PYG+WdOkenoeNH0Gx2tBoLDEJMQzfOxwE/FD3B7VSX1HeIo1GMK1TRYramTN4xSkehEWl3tjcTrvSP+ppqiv9yzuWZ8Q7Izhw7wALzi/QY+R5U0YTTC3gjBDianZPUxZCuAsh/IUQa9943UIIcVIIof+uu+gw+K2bdnZK+4VgYMjkY5O5/PQy39b6FlcrV72HoCjKv1mbGjG/uw9RsQkM+PWU7kF/lwrQ4iftRmV7JqbYpHPpzvgV92PWmVkcf5i8O01JXUYTTFOgJPAu/2GashBioRAiRAhx4Y3XmyYlrxtCCJ19S1LKW1LKPim8NRJYndGYMkxK+HMQPL8LHRaDpRObb21mzbU19CrXS5XfV5Rs5OFkxdSO3pwNes43G9IYpPfuApV7wMFpcPWvZG8LIfi6+te4WbkxYv8IVXk5AzKaYMYBYa9NTw4Hvs7E5y5Gm6xeEUIYALMAP6AM0EUIUUYIUV4IsemNR4rbPgohGgGXgEeZiClDXIM3wJVN0HgcuFUjMCyQcUfGUcmpEkMrDdX3xyuKkoam5QoyqF4JVh4LYvWJIN2N/aZAwQqwrj+mUcl/fFgYWfBD3R94EfuCUQdGkZCo465IeSWjCaZC0vbJAEgpnwGVMvqhUsr9wNM3Xq4C3Ei6M4kFVgGtpZTnpZQt3ngkH5HTqg9UA94D+gkhMvr90ufuUdxvLQHPFlB9MDEJMXy671OMDIzUYkpFyUGGNS5FjRL2fPXnBS7e1zETzMhUu4cMUPbiZIhPvoq/tF1pRlUZxdEHR/nl/C/6CjlPyehPQo0QwjYpsSCEsMvEOVJTGHj914xgoGpqjYUQ9sBEoJIQYpSUcpKU8suk93oCT6RMPu1DCNEf6A/g7OzM3r17MxyozfOLuFl4cMnhPRL27WN16GquRlzlA8cPuHL8Cle4kuFz5mcRERGZ+nNQ/k1dx5R1dpNcDpb0/OUQ39Qww8Io9cLwDiUGUe7iJIIX9eZGyX7J3reTdvha+DL7zGwMHxjiYeqhz9BzvYwmh6nA4aSBdQl0RPtDPiuk9Kee6mopKWUoMCCV9xbrOG4+MB/A19dX1qtXL0NBatVj754y1Ktfnx13dnDgzgF6lOnBh+98mIlzKXv37iVzfw7K69R1TF1hz2d0nn+Edfet+KW7j45N/uoR/PwCrvc24lqrs3ZR5huqxFWh48aOrHqxirX11qoKHTpkqAtJSrkUbR2yR8BjoK2UclkWxRIMFHntuSuQRonUzMmSYpdC8CDiAV8f/ppy9uVUEUtFycF8itoyys+LnZcf4X/wts62N0v0gEKV4c/B2s3K3mBhZMGUulMIjQ7lq8NfqU3KdMjwGIWU8pKUcqaU8mcp5aUsjOU4UFIIUVwIYQx0BjZk4fmzVIJM4PMDn5MoE5lSZwpGBkbZHZKiKDr0qlmMpmUL8t3WK5y88yzVdlJjBB0WaZ+s7a0tA/WGsvZlGeYzjL1Be1lxZYW+Qs719DMIngYhxErgCFBaCBEshOgjpYwHPgS2AZeB1VJKvRQByopqytvCtnEq5BSjq42miHWRtA9QFCVbCSGY3L4CLgVMGbLiFM8ikyeOV2yLQeuftUUx90xIsUlXr67Uca3DtBPTuPr0qn6CzuUyuh/MsBQefYQQ3hk5j5Syi5TSRUppJKV0lVL6J72+RUpZSkpZQkqZVWM7We7EwxP8FfYXrUq0ooV7vivFpii5lo2ZEbPf8+FJRCzD15zV3b1VpjX49IJD0+HGrmRvCyEYX3M81ibWjNw/kuh4HVsF5FMZvYPxRTuwXjjp0R+oB/wihBiRtaHpz38dgzHUGFLStCRfVP0iiyNTFEXfyrva8EUzT3ZfCUlzPIamk8DRC9Z9AC+Sr4+xM7VjYs2J3Ay7ydQTU/UUce6V0QRjD1SWUg6XUg5Hm3AcgTpAzyyOTW/+axeZt5M3Q5yHYGFkkcWRKYryNvSoUYx3yzgz+a8rnAvWsROmkZl2PCbmBfw5ABKTF7ysUbgG3cp0Y9XVVewL2qfHqHOfjCYYN+D1jss4oKiUMgrINfuLZuWWyYqi5D5CCKa0r4CjpQkfrjhNeHRc6o2dvKDJRLi5GwLmpNjk48ofU9q2NGMOj1GlZF6T0QSzAm1F5a+FEF8Dh4CVQggLtCVacoWs2jJZUZTcq4C5MTO6VOLe8yi+TGuTMt8+ULoZ7PxGW0H9DcYGxkyuM5nIuEi+Pvy1mrqcJN0JRmhXJi0G+gHPgTBggJRynJQyUkr5vn5CVBRF0Q/fYnZ80qgkG8/e5/dT91JvKAS0+llbOf33vhCXfBuAEgVK8InPJ+wP3q+2Wk6S7gQjtSn5TynlSSnldCnlT1LKE3qMTW9UF5miKH8bWM+DqsXtGLP+AoFPIlNvaOEAbebA4yuw/asUm3Tx7EJ1l+r8cOIHAsMC9RNwLpKZDcfe0Uskb5HqIlMU5W8GGsGPnbwxMtAwdNVp4hN1dG95NIRqg7U72F7fkextjdAwvuZ4jDRGjDowivjE5JuY5ScZTTD1gSNCiJvZveGYoihKVilUwIzJ7cpzLjiMddd1DPgDNByjnbq8fjC8fLMoPDhbOPNV9a+4EHoh3++CmdEE4weUABrwHzYcy26qi0xRlDc1LedClypF2HI7joBboak3NDKFdr9ok8vGj7SbD755rmJN8Svux7yz87gUmmvmP2W5jBa7vJPSQ1/B6YvqIlMUJSWjm5fByVwwbPVZ3VOXC5aHBqPh8gY4uzLFJl9W/RJbU1u+PPglMQm5ZhVHlspoqRghhOgqhBiT9NxNCFFFP6EpiqK8XRYmhvSvYMLD8Gi+Xp9GKcQaQ6BoTdg6Ep4n3zHTxsSGsTXGcuP5DWadnqWniHO2jHaRzQaqA12Snr9Au82xoihKnlCigAFDG5Rk3el7bDirY8cQjQG0mQ0yEdYPSnGVf23X2rQv1Z7FFxdzJuSMHqPOmTKaYKpKKQcD0fBqy2TjLI9KURQlGw2uX4LKbgUYve48D8N0FLG0LaZd5X97v3ZmWQo+9f0UFwsXvjr0Vb4riJnRBBMnhDAgaadJIYQjkDxt53BqkF9RFF0MDTRM7ehNXILks7VpVF2u3AM8GsOOr+HJjWRvWxhZMK7mOALDA/n59M96jDrnyWiCmQGsA5yEEBOBg8C3WR6VnqlBfkVR0lLcwYIvmntx4PoTfg24m3rDv1f5G5okFcRMSNakqktVOpXuxLJLyzj16JQeo85Z0pVgksrEIKVcDowAJgEPgDZSyjWvt1EURckrulZ1o04pR77dfJnbulb5W7tAs+8h+DgcmZlik2E+wyhkWYivDn1FVHzyUjN5UXrvYPYIIYYIIdyklFeklLOklDOBm0KIBkKIJUAPPcapKIry1gkhmNKuAkYGgmGrz5Cga5V/+Q5QujnsngiPryV729zInHE1xnH3xV1mn5mtx6hzjvQmmKZAAtrKyfeFEJeEELeA62hnlP0opVyspxgVRVGyTUEbU8a1Lsfpu89ZcOBW6g2FgBY/grF50qyy5F1lVVyq0KFUB5ZeWsr5x8mrMuc16UowUspoKeVsKWVNoCjQEO3GY0WllP2klPlv/p2iKPlGa+9CNCnrzNQd17j+6EXqDa2cwS/trjJHM0fGHB5DbEJsim3yiowO8iOljJNSPpBS6tgGLmdTs8gURckIIQQT2pTHwtiAT9ecJT5Bx+TZ8u3Bs4W2q+zJ9WRvWxpbMqb6GG48v8H8c/P1GHX2S+8gf/W8NIivZpEpipJRjlYmjG9TjrPBYczbn0ZXWfOp2pplG4akuACzjmsdWrq3xP+8P1efXtVj1NkrvXcwPYCTQohVQoieQoiC+gxKURQlJ2pRoRDNy7swfed1runsKisITb+Du0fgeMoVlUe8MwJrE2u+OfwNCSmM1+QF6R2DGSClrAx8A9gCi4UQR4QQ3woh6iQtvlQURcnzxrYui6WpIZ+tPad7VlnFLuDRSLvN8rPkNYELmBbg8yqfcyH0AssvL9dfwNkoo9WUr0gpf5RSNkVbsv8g0AEI0EdwiqIoOY2DpQnftCrL2aDn+B9Ma1bZT9r/6ijrX8e1DjPPzCT4RbAeo84eGR7k/5uUMkpKuUVKOURK6ZuVQSmKouRkLSu40LiMM1O3X+PW44jUGxYoAo2+gVt74OyqZG8LIfiq2lcIBOOPjtddkiYXSu8gf1l9B6IoipJbCCGY2KYcJoYaRv5+jkRdXWW+faBINdg2CiIeJ3u7oEVBPvb5mMP3D7Pp1iY9Rv32pfcOZtnf/yOE6Pv6G0II8yyNKJ2EEO5CCH8hxNrXXqsnhDgghJgrhKiXHXEpipI/OFmb8lWLMhwPfMbyAB37Lmo00HI6xEZqk0wKOpbqSAWHCnx//HueR+faFSDJpDfBvD5FedAb7x3I6IcKIRYKIUKEEBfeeL2pEOKqEOKGEOJzXeeQUt6SUvZ582UgAjAF8l6HpqIoOUp7H1dql3Tgu61XuP9cR30xJ0+oPRzOr4HrO5K9baAxYEz1MbyIfcG0k9P0GPHbld4E8/r935vrYTIzjrMYbfmZf06qnYk2C/ADygBdhBBlhBDlhRCb3ng4pXLeA1JKP2AkMDYTcSmKoqSbEIJv/1eeRAlf/XlB9xhKrU/A0RM2fQIxycdtStuVpnvZ7qy7sY7jD4/rMeq3J73JoWDS+pdKJE8wGR6VklLuB56+8XIV4EbSnUkssApoLaU8L6Vs8cYjJJXz/r2i6RlgktG4FEVRMqqInTnD3y3FrishbDz3IPWGhibQcgaEBcHeSSk2GVBxAIUtCzPuyLg8UUbGMJ3tvgF8gV6AqxDiInAl6eGQRbEUBl7f2DoYqJpaYyGEPTARqCSEGCWlnCSEaAs0AQoAKRYCEkL0B/oDODs7s3fv3kwFGxERkeljlX9T1zJrqOuYNTJzHd2lxN1Gw5e/n0bz6CqWxqkXPinl0gSXI7M5GVuCCCv3ZO+3Nm/N7JDZjNk4hqYFmqZwhtwjXQlGSvmvgjlCCFegAlAe2J9FsaT0J5Lq3ZGUMhQY8MZrfwB/6PqQpO8yH8DX11fWq1cvw4EC7N27l8weq/ybupZZQ13HrJHZ6+jiGU7Lnw+yL8yO7ztUTL1h1Yowswq+D5ZC812g+fc69XrU4+a+m+y4u4NBDQbhZu2W4VhyivROU/YTQhT6+7mUMjhpDcxkKWXXLIolGCjy2nNX4H4WnftfVLFLRVGympeLNf3quLPmZPD/27vz+Cjra4/jn5OdVWTHJCg7hkVSQJGyowherXrFitZeWxGXCqItWpWL1VpcKih137BYe8Va6oIWRVAjCKfULzMAABXxSURBVCKbICICYZOI7DuoEPjdPybUmMlMZsI8s2S+79drXmae5zfzHM4rePg9y/kxd832wAOrnQgD74NNn8L8Z8sdcmvXW0lPTefeT+5N6GdjQr0GczEw1cw2mtkMM3vQzK4ws/YRbBOzAGhlZs3MLAMYAkyN0Hf/iJpdiogXRvZvxcn1qjP6tc/57nCQ/mLtL4YW/eH9e2DP1367G1ZvyIj8EczZNId3N7zrYcTeCrUX2dUlT+uPB1YB64C+wHwgyA3g5TOzycDHQBszKzKzoc65YmA4MB1YAbzinFse7neHeHzNYEQk4rLSUxl7YQfWbT/AY+8XBh54rOPy0eKAz8Zc2uZSTq17Kg/Mf4D9h4J0C4hj4d5i/Gvn3A0li48NBXri60cWFufcZc65Js65dOdcjnNuYsn2ac651s65Fs65seF+bxjH1wxGRDzRo1V9LsrP5ulZa4IvTla3GfQcBV+8Aatn+u1OS0ljTLcxbP92O08sTcwllsMtMHvNrPOxN865RUDryIbkPc1gRMRLo//rVKpnpDG6omdjfnoj1GsF034Hh/0f1OzQoAMXt76Yl1a8xKpdqzyM2BvhFpihwBNm9lczG25mTwOHPYjLU5rBiIiX6tfM5LZBbZm/bif/Wux/jeU/0jJ9p8p2rYfZ5T/BPzJ/JLUyajF23tiEu+Af1oqWzrlVQHdgGtAI37WScz2MzxOawYiI1y7tkkvnk0/k3mkr2HUgyEOTzXtDh0tgzgTY7n/dpk5WHW76yU0s3ro44Zphhr2iJfBLfC1ZxjjnJpQ8j5JQNIMREa+lpBhjL2rPnm8Pc//bXwYfPGAspGXB27eUu27MRa0uomP9joxbOI69h/Z6FHHkaUVLERGPtG1cm6t7NOMfCzeyaEPZ7lil1GoEfUfDmvdhhf/TGSmWwuhuo9n13S6eXPKkhxFHVlKuaKlTZCISLTf2b0WTE7IY/drnFB85Gnhg16uhUXt453Zfa/8y8urlMbj1YCZ/OZnVu1Z7GHHkhN0J2cwyzKxaIq9oqVNkIhItNTLTuPO8PL7cvI+/fRzkscHUNDh3HOz9GmY9WO6QG/NvpGZGTe6bf19CXPAPq8CY2UjgG6DQzFaY2XBvwhIRqToGtm9Mr9YNeGjGKrbu/S7wwJPPhNMuh7mPwXb/WUqdrDqM6DSCBZsXMH39dA8jjoxQ7yKbYGb/A4wETnXOZQO9gDwzu8fLAEVEEp2ZcffP2nGo+Chjp60IPvjsuyG9Grz9+3Iv+A9uPZi2ddsybuE4Dh4+6FHEkRHqDOZDoCW+1vxzzWwx8CCwBhhiZnU8is8TugYjItHWrH4NruvdnDeWbGLe2iA339ZsCH1uhzXvwcppfrtTU1K544w72HJwC88te87DiI9fqHeRveacuxOYB1wAnAW8ABQDdYECMwvSeCe+6BqMiMTC9X1akl2nGn94Y3nwC/6nD/OtfvnO7eU+4Z/fMJ9zm53LC8tfoGhf/K4OH+5F/uHA3/E1vcwH2gPLnHOd8C1zLCIiAVTLSGXMeXms3FLRBf90GPRn2L0B5j5a7pCbO99Makoq4xeO9yja4xfWejAlT/KfAUwBqgGfARcBlCxzLCIiQZzTrhG9Wjfg4Rmr2Lbv+8ADm/eGvAt8LWR2b/Tb3bhGY4a2H8rMr2byyTfx+aRI2OvBAP8G+uBr0/8BkDiPlZbQNRgRiRUz4w/n5/Fd8REeeCeEJ/xxMOPOcndf2e5Ksmtmc//8+yk+Whz5YI9TTNaDiTVdgxGRWGrRoCZDezRnyqIiFn+1K/DAOrnw05tg+auwYa7f7qy0LEZ1GUXh7kKmrJriYcSVE5P1YEREkt3wfi1pWCuTu6cu5+jRYC39R0LtbN9ty0f9V8ns37Q/XRt35fElj8ddn7KkXA9GRCTWamamcdugtiwt2sOUxUHuBMuoDmf/ETZ/Bkv+z2+3mXFr11vZ8/0enl76tIcRhy/cAnMVVWA9GBGReHBhp2zym9bhz++sZN93Qf5X2v5iyO0G7/0RvvO/dty2blsuanURL614ifV71nsXcJjCbXa5miqwHoyISDxISTHuOr8d2/d/z6PvB3mU0AwG3Q8HtsOsceUOGZE/gozUDMYvip/blsNudumcO+Kc+2cirwcjIhIvTsutw8+75PDXOetYu21/4IEn5UOny+GTp2DnWr/d9avVZ1jHYRRsLGDeN/M8jDh0YReYqkC3KYtIPLnlnLZkpqVyb0V9yvqNgZT0gLct/zLvl2TXzGbcgnEcKeeGgGirsMCY2dlm9qyZdSp5f433YXlLtymLSDxpUCuTG/q2ZOaKrcxevS3wwNpNoMfNsOJNWO9/A29maiYjfzKSlbtWMnWN/8Jl0RbKDOY3wC3AFWbWD+jkbUgiIsnnqh6n0LRude5564vgfcq6D4faOb4+ZeXMUgaeMpCO9Tvy6KePxrzbcigFZptzbrdzbhQwAOjqcUwiIkknMy2VO85ty6ot+5k8/6vAA9Or+Vr6b/4Mlr7st9vMuKXrLWz7dhuTlk/yLuAQhFJg/n3sB+fcbcDfvAtHRCR5ndOuMd2a1+WhGavYc7CC25azu8D795S7vHKnhp0YcPIAJi2fxJYDWzyMOLgKC4xz7g0za1fqffmtPUVE5LiYGWPOy2P3t4d57AP/FS1LDYRz7oV93/hWvyzHTZ1vovhoMU8sfcKjaCsW6l1kLx77wcyuLr3DzKpHNCIRkSTW7qQTuKRzDpPmrmfDDv/ZyX80PcPXbXnOX2DfZr/dubVyGdJ2CK8Xvs7qXUGKlYdCLTBW6ufflNk3O0KxhMXMmpvZRDObUmpbipmNNbNHzezKWMQlInK8fjegDempKdz/dgXdls+6C44cgg/Glrv72o7XUiO9Bg8vejjiMYYi1AJTuhObldkX9rM0Zva8mW01s8/LbB9oZivNrNDMbgsakHNrSxpulnYBkI2vfU38LvMmIhJEo9pZXNurBW9/vpkF63cGHli3OZx+DXz6d9iy3G/3CZknMKzDMGZ/PTsma8aEWhwam9mvzCwf/wITpA1oQJOAgaU3mFkq8DgwCN/qmJeZWZ6ZdTCzt8q8Ggb43jbAx8653wLXVyIuEZG4MKxXMxrXzuJPb30RvNtyr1GQWTvgw5eXn3o5TWo0YfzC8Rx1QW5/9kCoBeYuoAswAcgxs+Vm9i8zGwvUD/egzrlZQNmyfDpQWDIzOQS8DFzgnFvmnDuvzGtrgK8uAo4trhD7x1hFRCqpekYao85pw9KiPbz52aYgA+v6ikzhTFjzgd/uzNRMRuSPYMXOFUxbN83DiP2lhTLIOfdM6fdmlgN0BDoAsyIUSzZQel3QInzLM5fLzOoBY4F8M7vdOXcf8CrwqJn1DBRXSSeCawAaNWpEQUFBpYLdv39/pT8rP6ZcRobyGBnxlMe6ztG0Vgp/fH0p1XasIiO17Akkn5QjbTg9syGHX/stizqPB/vx3KGGq0FOeg4Pzn2QzA2ZpFt6NMIPrcCU5ZwrwlcAIlkOy8tcwHlhSZPN68psOwiUvS5T9nPPmNk3wPm1atXq3KdPn0qECgUFBVT2s/JjymVkKI+REW95zMjZzhUTP2F9elOu6dUi8MAGY8l6dRh96m6D0y712531dRbXzryWLY22cEXeFR5G/IN4anZZBOSWep8DBJkXVp56kYlIoujRqj69WzfgsfcL2X3wUOCB7QdDk9N8D18e/s5vd/fs7nRr0o2nP3uafYf2eRjxD+KpwCwAWplZMzPLAIYAnnRrUzdlEUkkd5x7Kvu/Lw6+ZkxKCpx9D+zZCPOfKXfITZ1vYvf3u6PWQiYmBcbMJgMfA23MrMjMhjrnioHhwHR8C5m94pzzv+8uAjSDEZFE0qZxLS7pnMvfPl7PVzuCNLBs3htangWzx8O3u/12t6vXjkGnDOLFL15k28EgXZsjJCYFxjl3mXOuiXMu3TmX45ybWLJ9mnOutXOuhXOu/CeHIkAzGBFJNL8d0JrUFGPcuyuDDzzrLt+yynMmlLt7RP4IDh85HJVZTDydIosazWBEJNE0qp3F0B7NmLp0E8uKgvzjuHEH6HAJzHsS9vpfxs6tncsj/R7hhk43eBitT1IWGBGRRHRt7xacWD2dB96poIVMv9G+tWIK7it3d8+cnlRP976NZFIWGJ0iE5FEVDsrnRH9WvFR4XZmrQpyDeXEU6Dr1b4WMtsqOKXmoaQsMDpFJiKJ6hfdmpJbtxr3v/1lxS1k0mv4bluOkaQsMJrBiEiiykxLZdSANnzxzV6mLg3yqGCN+r7llVe8CV8vil6ApSRlgdEMRkQS2fkdTyKvSW3Gz1jJoeIgDSzPvAGq14OZd0cvuFKSssCIiCSylBTj1oFt2LjzW15e8FXggZm1oOcoWPchrC2IWnzHqMCIiCSg3q0bcEazujzy3moOfF8ceGCXq6B2jm8W4yqzukrlJWWB0TUYEUl0ZsbvB7Vl+/5DPP/RusAD07Og7+2wabHvekwUJWWB0TUYEakKftL0RAbkNeLpWWvZeSBII8yOQ6B+a9/Sykejt1RWUhYYEZGq4pZz2nDwUDFPfbgm8KDUNOh7B2z7EpZNiVpsKjAiIgmsVaNaXJifzQtz17N5j3+b/v849QJfG5mCe+HI4ajEpgIjIpLgbj6rNUed49H3VwcelJIC/cbArvXw6YtRiSspC4wu8otIVZJbtzpDujblHws2smHHgcADWw2AnNPhwwfLXZQs0pKywOgiv4hUNSP6tSQt1ZgwM8gsxgz6j4F9m2DhRM9jSsoCIyJS1TSsncWV3U/h9SVfs3JzkCWRm/WCM4dDdmfPY1KBERGpIq7v3YKaGWk8NKOCDsrnjIWm3TyPRwVGRKSKqFM9g6E9mzF9+Zbgi5JFiQqMiEgVclWPZtSpnl7xLCYKkrLA6C4yEamqamelc22vFnywchuLNuyKaSxJWWB0F5mIVGVXdj+Z+jUzYj6LScoCIyJSlVXPSOP6Pi2ZU7iDuWu2xywOFRgRkSroF2c0pVHtTCbMWI2Lcpv+Y1RgRESqoKz0VG7o25L563cyp3BHTGJQgRERqaIu7ZpLkxOyeHjmqpjMYlRgRESqqMw03yxm0YZdzFod/WsxCVtgzKy5mU00symltvU0s6fM7DkzmxvL+ERE4sHPu+SSXacaD8+I/iwmJgXGzJ43s61m9nmZ7QPNbKWZFZrZbcG+wzm31jk3tMy22c6564C3gBciH7mISGLJSEtheL+WLNm4m4KV26J67FjNYCYBA0tvMLNU4HFgEJAHXGZmeWbWwczeKvNqWMH3Xw5M9iJwEZFEM7hzDjknVmPCe9G9oywmBcY5NwvYWWbz6UBhyczkEPAycIFzbplz7rwyr62BvtvMmgJ7nHN7vfsTiIgkjvTUFG7o25KlG3dTsCp6s5i0qB2pYtnAxlLvi4AzAg02s3rAWCDfzG53zt1Xsmso8Ncgn7sGuAagUaNGFBQUVCrY/fv3V/qz8mPKZWQoj5FRVfPY4KijXpZxz6uLoFsWZub5MeOpwJT3pw04l3PO7QCuK2f7H4IdxDn3jJl9A5xfq1atzn369Ak3TgAKCgqo7Gflx5TLyFAeI6Mq5/F3Nb7ijteWYSe1o0+biq40HL94uousCMgt9T4H2OTFgdSLTESS0eDOOWTXqcZfonQtJp4KzAKglZk1M7MMYAgw1YsDqZuyiCSjjLQUftO3BZ9+tZvZUXguJla3KU8GPgbamFmRmQ11zhUDw4HpwArgFefcci+OrxmMiCSrSzrn0qx+DdZs2+/5sWJyDcY5d1mA7dOAaV4f38zOB85v2bKl14cSEYkrGWkpvHtzL9JTvZ9fxNMpsqjRDEZEklk0igskaYERERHvJWWB0UV+ERHvJWWB0SkyERHvJWWBERER7yVlgdEpMhER7yVlgdEpMhER7yVlgREREe/FU7PLqDn2oCWw18xWl2w+ASh9zqzs+7Lb6gNe91ooL4ZIf7aicYH2h7O9otx6nct4zmOgffGYx0BxRfJzymPkPufl3+2TQ4rAOaeXr+nbM8Hel90GLIx2TF58tqJxgfaHs72i3Hqdy3jOY6g5i4c8Hk8ulcfo5vF4chnu9mAvnSL7wZsVvA+0zUvHc7xQP1vRuED7w9keSm69FM95DLQvHvN4PMdUHiNzzHA+F42/20FZSWWSMJnZQudcl1jHURUol5GhPEaG8hg5msFU3jOxDqAKUS4jQ3mMDOUxQjSDERERT2gGIyIinlCBERERT6jAiIiIJ1RgPGBmF5rZs2b2hpkNiHU8icrMmpvZRDObEutYEpGZ1TCzF0p+F38R63gSlX4PK08Fpgwze97MtprZ52W2DzSzlWZWaGa3BfsO59zrzrlhwK+ASz0MN25FKI9rnXNDvY00sYSZ1/8GppT8Lv4s6sHGsXDyqN/DylOB8TcJGFh6g5mlAo8Dg4A84DIzyzOzDmb2VplXw1If/d+SzyWjSUQuj/KDSYSYVyAH2Fgy7EgUY0wEkwg9j1JJSdmLLBjn3CwzO6XM5tOBQufcWgAzexm4wDl3H3Be2e8wMwPuB952zi32NuL4FIk8ir9w8goU4SsyS9A/Jn8kzDx+Ed3oqg790oUmmx/+JQi+v7jZQcaPAM4CBpvZdV4GlmDCyqOZ1TOzp4B8M7vd6+ASWKC8vgpcbGZPEpt2KImm3Dzq97DyNIMJjZWzLeATqs65R4BHvAsnYYWbxx2ACnTFys2rc+4A8OtoB5PAAuVRv4eVpBlMaIqA3FLvc4BNMYolkSmP3lBeI0N5jDAVmNAsAFqZWTMzywCGAFNjHFMiUh69obxGhvIYYSowZZjZZOBjoI2ZFZnZUOdcMTAcmA6sAF5xzi2PZZzxTnn0hvIaGcpjdKjZpYiIeEIzGBER8YQKjIiIeEIFRkREPKECIyIinlCBERERT6jAiIiIJ1RgRDxgZkfMbEmpV7lLE5jZmWb2bKn3HUp9ZqeZrSv5eWaQY31kZv3LbBtlZmpXJDGlXmQi3vjWOdcphHEDgXeOvXHOLQM6AZjZJOAt51xFC11NxvfU+Xultg3B13RVJGY0gxGJrf5AwNlJWWZ2m5nNN7PPzOzOks3/BH5mZuklY1oC9YB5EY9WJAwqMCLeqFbmFJnfyqZmVh847JzbE8oXmtm5QFPgDHyznO5m1t05txXfmi/HluceArzs1KZDYkynyES8EcopsgHAu2F85wB8qy1+WvK+JtAamMsPp8n+XfLfy8OKVsQDKjAisTMIeCiM8Qb8yTk3sZx9rwJ/NrMuQIpz7rNIBChyPHSKTCQGSpbV7ojv1FaopgNDzaxGyXfklJxmwzm3F/gIeA54KcLhilSKZjAi3qhmZqWLxzvOudK3KncGPg3nOolzbpqZtQXm+eoT+/CdCtteMmQy8Aow+LgiF4kQtesXiQEz+1+g0Dn3cqxjEfGKCoyIiHhC12BERMQTKjAiIuIJFRgREfGECoyIiHhCBUZERDyhAiMiIp5QgREREU/8Pwps/xLai11LAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1175fb128>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for z in [-0.2, 0, +0.2]:\n",
    "    model = Log10ParabolaEnergyScale(z=z)\n",
    "    model.plot([0.01, 50] * u.TeV, energy_power=2, label=f\"z={z}\")\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel(\"E / TeV\")\n",
    "plt.ylabel(\n",
    "    r\"$E^2\\,d\\phi/dE\\,/\\,(\\mathrm{erg}\\,\\mathrm{cm}^{-2}\\,\\mathrm{s}^{-1})$\"\n",
    ")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Execute fit\n",
    "\n",
    "Let's run the fit.\n",
    "The ``delta`` values used are 15% for IACTs and 5% for Fermi-LAT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'fermi': 0.05, 'magic': 0.15, 'hess': 0.15, 'veritas': 0.15, 'fact': 0.15}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "SystematicsSpectrumFit.delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 9.13 s, sys: 39.5 ms, total: 9.17 s\n",
      "Wall time: 9.19 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "fit = SystematicsSpectrumFit()\n",
    "result = fit.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit results\n",
    "\n",
    "We use `MultiNorm` to have nice table representations of the values,\n",
    "errors and parameter correlation matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "FitResult\n",
       "\n",
       "\tbackend    : minuit\n",
       "\tmethod     : minuit\n",
       "\tsuccess    : True\n",
       "\tnfev       : 594\n",
       "\ttotal stat : 441.23\n",
       "\tmessage    : Optimization terminated successfully."
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = result.model.parameters.to_table()\n",
    "multi_norm = MultiNorm(\n",
    "    mean=list(t[\"value\"]),\n",
    "    names=list(t[\"name\"]),\n",
    "    cov=result.model.parameters.covariance,\n",
    ").drop([\"reference\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean</th>\n",
       "      <th>err</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>name</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>amplitude</th>\n",
       "      <td>3.880091e-11</td>\n",
       "      <td>3.536017e-12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>2.510771e+00</td>\n",
       "      <td>3.763238e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta</th>\n",
       "      <td>2.423330e-01</td>\n",
       "      <td>2.759517e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>z_fermi</th>\n",
       "      <td>4.365007e-03</td>\n",
       "      <td>4.511147e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>z_magic</th>\n",
       "      <td>-2.684315e-03</td>\n",
       "      <td>6.089172e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>z_veritas</th>\n",
       "      <td>-5.560264e-02</td>\n",
       "      <td>6.427045e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>z_fact</th>\n",
       "      <td>-1.029063e-01</td>\n",
       "      <td>6.579716e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>z_hess</th>\n",
       "      <td>9.929812e-02</td>\n",
       "      <td>7.293912e-02</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   mean           err\n",
       "name                                 \n",
       "amplitude  3.880091e-11  3.536017e-12\n",
       "alpha      2.510771e+00  3.763238e-02\n",
       "beta       2.423330e-01  2.759517e-02\n",
       "z_fermi    4.365007e-03  4.511147e-02\n",
       "z_magic   -2.684315e-03  6.089172e-02\n",
       "z_veritas -5.560264e-02  6.427045e-02\n",
       "z_fact    -1.029063e-01  6.579716e-02\n",
       "z_hess     9.929812e-02  7.293912e-02"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "multi_norm.parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>name</th>\n",
       "      <th>amplitude</th>\n",
       "      <th>alpha</th>\n",
       "      <th>beta</th>\n",
       "      <th>z_fermi</th>\n",
       "      <th>z_magic</th>\n",
       "      <th>z_veritas</th>\n",
       "      <th>z_fact</th>\n",
       "      <th>z_hess</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>name</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>amplitude</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>-0.489970</td>\n",
       "      <td>0.192608</td>\n",
       "      <td>-0.400052</td>\n",
       "      <td>-0.814372</td>\n",
       "      <td>-0.759988</td>\n",
       "      <td>-0.794817</td>\n",
       "      <td>-0.867305</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>-0.489970</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.626201</td>\n",
       "      <td>0.025503</td>\n",
       "      <td>0.252517</td>\n",
       "      <td>0.285266</td>\n",
       "      <td>0.620496</td>\n",
       "      <td>0.589589</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>beta</th>\n",
       "      <td>0.192608</td>\n",
       "      <td>0.626201</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.042951</td>\n",
       "      <td>-0.214565</td>\n",
       "      <td>-0.190374</td>\n",
       "      <td>0.093339</td>\n",
       "      <td>0.010345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>z_fermi</th>\n",
       "      <td>-0.400052</td>\n",
       "      <td>0.025503</td>\n",
       "      <td>0.042951</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.440457</td>\n",
       "      <td>0.362461</td>\n",
       "      <td>0.293658</td>\n",
       "      <td>0.331468</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>z_magic</th>\n",
       "      <td>-0.814372</td>\n",
       "      <td>0.252517</td>\n",
       "      <td>-0.214565</td>\n",
       "      <td>0.440457</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.646915</td>\n",
       "      <td>0.610670</td>\n",
       "      <td>0.680751</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>z_veritas</th>\n",
       "      <td>-0.759988</td>\n",
       "      <td>0.285266</td>\n",
       "      <td>-0.190374</td>\n",
       "      <td>0.362461</td>\n",
       "      <td>0.646915</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.581340</td>\n",
       "      <td>0.643209</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>z_fact</th>\n",
       "      <td>-0.794817</td>\n",
       "      <td>0.620496</td>\n",
       "      <td>0.093339</td>\n",
       "      <td>0.293658</td>\n",
       "      <td>0.610670</td>\n",
       "      <td>0.581340</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.741487</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>z_hess</th>\n",
       "      <td>-0.867305</td>\n",
       "      <td>0.589589</td>\n",
       "      <td>0.010345</td>\n",
       "      <td>0.331468</td>\n",
       "      <td>0.680751</td>\n",
       "      <td>0.643209</td>\n",
       "      <td>0.741487</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "name       amplitude     alpha      beta   z_fermi   z_magic  z_veritas  \\\n",
       "name                                                                      \n",
       "amplitude   1.000000 -0.489970  0.192608 -0.400052 -0.814372  -0.759988   \n",
       "alpha      -0.489970  1.000000  0.626201  0.025503  0.252517   0.285266   \n",
       "beta        0.192608  0.626201  1.000000  0.042951 -0.214565  -0.190374   \n",
       "z_fermi    -0.400052  0.025503  0.042951  1.000000  0.440457   0.362461   \n",
       "z_magic    -0.814372  0.252517 -0.214565  0.440457  1.000000   0.646915   \n",
       "z_veritas  -0.759988  0.285266 -0.190374  0.362461  0.646915   1.000000   \n",
       "z_fact     -0.794817  0.620496  0.093339  0.293658  0.610670   0.581340   \n",
       "z_hess     -0.867305  0.589589  0.010345  0.331468  0.680751   0.643209   \n",
       "\n",
       "name         z_fact    z_hess  \n",
       "name                           \n",
       "amplitude -0.794817 -0.867305  \n",
       "alpha      0.620496  0.589589  \n",
       "beta       0.093339  0.010345  \n",
       "z_fermi    0.293658  0.331468  \n",
       "z_magic    0.610670  0.680751  \n",
       "z_veritas  0.581340  0.643209  \n",
       "z_fact     1.000000  0.741487  \n",
       "z_hess     0.741487  1.000000  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "multi_norm.correlation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot\n",
    "\n",
    "The spectral results with the modified likelihood are given in Figure 4 of the paper.\n",
    "\n",
    "Here we plot the results for the energy-scale factors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.004365007495221839"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "multi_norm.mean[\"z_fermi\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfkAAAFpCAYAAACf0FGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAHx1JREFUeJzt3XFwnHd95/HP1846aWIJiA2uFaAie8EMJIJEKplyPcZJO4AgTtpjSzmwS9peNqlmoLimU3oZEgWutIA9HqiOcxTN1cFOJxCVu7NCBaEkopS0AclJ1smBIUtMg9d3xQbcFWrw4nzvDz22d5/Ikizt7k/70/s1o9Gz+/yeZ7/PVyt99Dz77LPm7gIAAPFZEboAAADQGIQ8AACRIuQBAIgUIQ8AQKQIeQAAIkXIAwAQKUIeAIBIEfIAAESKkAcAIFKEPAAAkTovdAGLtXbtWu/s7Kzb+kqlkk6cOKFVq1apo6OjbuttVT/96U910UUXhS5jyaAftejHGfSiFv04oxG5MjExcdTdXzzXOGv1a9f39PT4+Ph43dZnZqenW7039TA2NqaNGzeGLmPJoB+16McZ9KIW/TijEbliZhPu3jPXOA7XAwAQKUIeAIBIEfIAAESKkAcAIFKEPAAAkSLkAQCIVMu/T77ebr/9dh06dEj1fO89AGD5CpkrhHxKf38/7+8EANRNyFzhcD0AAJEi5AEAiBQhDwBApHhNPmXTpk06duyY1qxZo5GRkdDlAABaXMhcIeRT7r///tAlAAAiEjJXOFwPAECkCHkAACJFyAMAEClCHgCASBHyAABEipAHACBShDwAAJEi5AEAiBQhDwBApLjiXcqdd96pgwcPasOGDaFLAQBEIGSuEPIp+Xyez5MHANRNyFzhcD0AAJEi5AEAiBQhDwBApHhNPqW7u1vlclltbW2amJgIXQ4AoMWFzBVCPmX//v2hSwAARCRkrnC4HgCASBHyAABEipAHACBShDwAAJEi5AEAiBQhDwBApAh5AAAiRcgDABApQh4AgEhxxbuUffv26cCBA7riiitClwIAiEDIXCHkUzZt2qS2tjY+Tx4AUBchc4XD9QAARIqQBwAgUoQ8AACR4jX5lI6ODp04cUKrVq1SqVQKXQ4AoMWFzBVCPuXIkSOhSwAARCRkrnC4HgCASBHyAABEipAHACxasVhUX1+f2tvbtWLFCr3tbW9TX1+fisVi6NKWtbqGvJl1mtkT9VwnAGBpGx0dVVdXl4aGhlQul+Xumpqa0tDQkLq6ujQ6Ohq6xGWLPXkAwIIVi0XlcjlNTU2pUqnUzKtUKpqamlIul2OPPpBGhPxKM7vLzJ40swfM7BfMLGtmXzSzCTP7mpm9SpLM7LfM7Akze9zM/j657zVm9g0ze8zMCmZ2WQNqBADUwY4dO54X7mmVSkU7d+5sUkWo1oiQv0zSf3P310j6iaS3SxqU9F5375b0AUmfTsbeJunN7v5aSdcn990i6ZPu/jpJPZJ+0IAaAQB1sHfv3nmF/J49e5pUEao14n3yT7v7Y8n0hKROSW+QdJ+ZnRpzfvL965J2m9nnJH0+ue8fJd1qZi+V9Hl3/276AcwsLykvSevWrdPY2FgDNkMNW28rmZycpA9V6Ect+nHGcu3F5OTkvMaVy+Vl2Z+0ZvfA3L1+KzPrlHS/u1+e3P6ApA5J/8nd159lmaslvU3S70p6nbsfM7Nsct/7Jf1nd3/wbI/Z09Pj4+Pj9dyG09P17E2rGhsb4xP5qtCPWvTjjOXai/b2dpXL5XmNO378eBMqWnoakStmNuHuPXONa8aJd/8q6Wkz+y1JsmmvTaaz7v6Iu98m6aikl5nZpZK+5+6fkrRPUlcTajxtfHxcu3btUj3/cQCAWG3evFmZTGbWMZlMRlu2bGlSRUtPyFxp1tn175b0+2b2uKQnJd2Q3P8JMzuQvO3u7yU9Lum3JT1hZo9JepWkzzSpRklSd3e3NmzYoO7u7mY+LAC0pG3bts0r5Ldu3dqkipaekLlS19fk3f2QpMurbm+vmv2WGcb/xxlW8+fJFwBgictmsxoeHlYul1OlUqk5CS+TySiTyWh4eFjZbDZglcsX75MHACxKb2+vCoWC8vn86SveXXjhhcrn8yoUCurt7Q1d4rJFyAMAFi2bzWpgYEDHjx/XyZMn9YUvfEEDAwPswQfGR82mcHY9AKCeQuYKe/IAAESKkAcAIFKEPAAAkSLkAQCIFCEPAECkCHkAACJFyAMAEClCHgCASBHyAABEipAHACBShHzK4cOHdd999+nw4cOhSwEARCBkrhDyKR0dHVq7dq06OjpClwIAiEDIXCHkAQCIFCEPAECkCPmUUqmko0ePqlQqhS4FABCBkLnC58mnXHLJJaen+Tx5AMBihcwV9uQBAIgUIQ8AQKQIeQAAIkXIAwAQKUIeAIBIEfIAAESKkAcAIFKEPAAAkSLkAQCIFCEPAECkCPkUd9dDDz3EJW0BAHURMlcIeQAAIkXIAwAQKUIeAIBIEfIpExMTOnjwoCYmJkKXAgCIQMhc4fPkU3p6ek5Pc/IdAGCxQuYKe/IAAESKkAcAIFKEPAAAkSLkAQCIFCEPAECkCHkAACJFyAMAEClCHgCASBHyAABEipBPWb9+vdasWaP169eHLqVpisWi+vr61N7erhUrVqi9vV19fX0qFouhSwOAlhcyVxoe8mb2PjP7lpndc47Lvd/MLmxUXWdTKpU0PDysUqnU7IcOYnR0VF1dXRoaGlK5XJa7q1wua2hoSF1dXXrkkUdClwgALS1krjRjT75P0lvd/d3nuNz7JTU95JeTYrGoXC6nqakpVSqVmnmVSkVTU1Pq7+9njx4AWlRDQ97Mdkm6VNI+M/sTM3vYzB5Nvm9Ixqw0s+1mdsDMCmb2XjN7n6QOSQ+Z2UONrHE527Fjx/PCPa1SqWjnzp1NqggAUE8NDXl3v0VSSdI1kv67pDe6+5WSbpP00WRYXtIrJF3p7l2S7nH3T51azt2vaWSNy9nevXvnDPmTJ09qz549TaoIAFBP1uiPvTOzQ5J6JP2CpE9JukySS8q4+6vM7G8k7XL3L8+0nLsfnWGdeU3/c6B169Z133vvvXWr9+GHH9azzz6rCy64QG94wxvqtt6l6Nprr53Xxx6amR588MEmVLT0TU5OavXq1aHLWDLoxxn0ohb9OKMRuXLNNddMuHvPXOOaGfLbJe1390+ZWaekMXfvNLPPS/q0u//dTMvNFPLVenp6fHx8vJ71np6O/fPk29vbVS6X5zXu+PHjTaho6RsbG9PGjRtDl7Fk0I8z6EUt+nFGI3LFzOYV8s18C90LJB1Opm+suv8BSbeY2XmSZGYXJ/eXJbU1rbplaPPmzcpkMrOOWblypbZs2dKkigAA9dTMkP+4pD83s69LWll1/5Ckf5ZUMLPHJb0ruX9Q0ign3jXOtm3b5gz5TCajrVu3NqkiAEA9ndfoB3D3zmTyqKRXVs36UDL/55L+KPmqXu4vJf1lo+tbzrLZrIaHh5XL5VSpVGpOwstkMspkMrrtttuUzWYDVgkAWCiueLfM9fb2qlAoKJ/P11zxLp/Pq1Ao6Oqrrw5dIgBggRq+J4+lL5vNamBgQAMDA8+b98wzzwSoCABQD+zJAwAQKUIeAIBIEfIAAESKkAcAIFKceJdy1VVXqVwuq62N6/AAABYvZK4Q8ikTExNcjhEAUDchc4XD9QAARIqQBwAgUoQ8AACRIuRTBgcHNTIyosHBwdClAAAiEDJXOPEu5eabbz49nc/nA1YCAIhByFxhTx4AgEgR8gAARIqQBwAgUoQ8AACRIuQBAIgUIQ8AQKQIeQAAIkXIAwAQKUIeAIBIccW7lOuuu07Hjh3TmjVrQpcCAIhAyFwh5FNGRkb4PHkAQN2EzBUO1wMAEClCHgCASBHyAABEipBP6e/v1+7du9Xf3x+6FABABELmCifepdxxxx2npwl6AMBihcwV9uQBAIgUIQ8AQKQIeQAAIkXIAwAQKUIeAIBIEfIAAESKkAcAIFKEPAAAkSLkAQCIFFe8S7nppptUKpXU0dERuhQAQARC5gohnzI4OMjnyQMA6iZkrnC4HgCASBHyAABEipAHACBShHxKPp/X9u3blc/nQ5cCAIhAyFzhxLuUu+666/T04OBgwEoAADEImSvsyWPZKxaL6uvrU3t7u1asWKH29nb19fWpWCyGLg0AFqWhIW9mf2tmL0y++hr5WMBCjI6OqqurS0NDQyqXy3J3lctlDQ0NqaurS6Ojo6FLBIAFa0jI27QV7v5Wd/+JpBdKIuSxpBSLReVyOU1NTalSqdTMq1QqmpqaUi6XY48eQMuaNeTN7GPVe+Bm1m9m28zsj83sm2ZWMLM7knmdZvYtM/u0pP2SXmZmh8xsraS/kJQ1s8fM7BNmttrMvmJm+83sgJndkKzjIjP7gpk9bmZPmNlvN27Tsdzt2LHjeeGeVqlUtHPnziZVBAD1Ndee/L2SqoP2HZJ+KOkySa+X9DpJ3Wb2xmT+Bkmfcfcr3f37Vct9UFLR3V/n7n8s6VlJv+nuV0m6RtIOMzNJb5FUcvfXuvvlkr64yO0Dzmrv3r3zCvk9e/Y0qSIAqK9Zz65390fN7CVm1iHpxZJ+LKlL0pskPZoMW63p0P9nSd9393+ax+OapI8m/xw8J+kSSeskHZC03cw+Jul+d//ajAub5SXlJWndunUaGxubx0Oeu0att5VMTk5G24fJycl5jSuXy6d7EHM/FoJ+nEEvatGPmTW7J+busw8w+4im995/UdIRSZ2SvuPud6bGdWo6mC+vuu+QpB5N/yNwep6Z3SipV9Jmd68k4za6+yEzu1jSWyXdIukBd//wbPX19PT4+Pj4/LZ2HqYPKEybqzfLQczX8W9vb1e5XJ7XuOPHj0uKux8LQT/OoBe16McZjcgVM5tw9565xs3nxLt7Jb1TUk7SsKQvSfo9M1udPNAlZvaSOdZRltRWdfsFkv4lCfhrJP1Ssq4OSVPuvlfSdklXzaM+YEE2b96sTCYz65hMJqMtW7Y0qSIAqK85Q97dn9R0QB929yPu/oCkv5b0j2Z2QNPB3zbHOo5J+npyMt0nJN0jqcfMxiW9W9K3k6FXSPqGmT0m6VZJ/3WB2wXMadu2bfMK+a1btzapIgCor3ld8c7dr0jd/qSkT84w9PLUuM6q6Xelxv7KDMsf0vSRgmBuv/12HTp0SJ2dnSHLQBNks1kNDw8rl8upUqnUnISXyWSUyWQ0PDysbDYbsEoArS5krnBZ25T+/n5eS1pGent7VSgUtHPnTu3Zs0eTk5NavXq1tmzZoq1btxLwABYtZK4Q8lj2stmsBgYGNDAwELoUAKgrrl0PAECkCHkAACLF4fqUTZs26dixY1qzZo1GRkZClwMAaHEhc4WQT7n//vtDlwAAiEjIXOFwPQAAkSLkAQCIFCEPAECkCHkAACJFyAMAEClCHgCASBHyAABEipAHACBShDwAAJHiincpd955pw4ePKgNGzaELgUAEIGQuULIp+TzeT5PHgBQNyFzhcP1AABEipAHACBShDwAAJHiNfmU7u5ulctltbW1aWJiInQ5AIAWFzJXCPmU/fv3hy4BABCRkLnC4XoAACJFyAMAEClCHgCASBHyAABEipAHACBShDwAAJEi5AEAiBQhDwBApAh5AAAixRXvUvbt26cDBw7oiiuuCF0KACACIXOFkE/ZtGmT2tra+Dx5AEBdhMwVDtcDABApQh4AgEgR8gAARIrX5FM6Ojp04sQJrVq1SqVSKXQ5AIAWFzJXCPmUI0eOhC4BABCRkLnC4XoAACJFyAMAEClCHgCASBHyAABEipAHACBShDwAAJEi5AHMqlgsqq+vT+3t7VqxYoXa29vV19enYrEYujQAc1gyIW9mf2tmLwxdB4AzRkdH1dXVpaGhIZXLZbm7yuWyhoaG1NXVpUceeSR0iQBmsWRC3t3f6u4/CV0HgGnFYlG5XE5TU1OqVCo18yqViqamptTf388ePbCEnXPIm1mnmX3bzIbM7Akzu8fMft3Mvm5m3zWz1ydfD5vZo8n3DcmyF5rZ58ysYGafNbNHzKwnmXfIzNYm07+TjHnczPbUd5NnNz4+rl27dml8fLyZDwssOTt27HheuKdVKhXt3LmzSRUBrSlkrix0T/7fSfqkpC5Jr5L0Lkm/KukDkv6LpG9LeqO7XynpNkkfTZbrk/Rjd++S9BFJ3ekVm9lrJN0q6Vp3f62kP1xgjQvS3d2tDRs2qLv7eaUBy8revXvnDPmTJ09qz56m/h8OtJyQubLQa9c/7e4HJMnMnpT0FXd3MzsgqVPSCyTdbWaXSXJJmWS5X9X0Pwdy9yfMrDDDuq+VNOzuR5NxP0oPMLO8pLwkrVu3TmNjYwvcjJlNTk7WfZ2til7UWk79mJycnNe4crm8bHoym+X03JgP+lErVD8WGvI/q5p+rur2c8k6PyLpIXf/TTPrlDSWzLd5rNs0/Y/BWbn7oKRBSerp6fGNGzfOs+z5GRsbU73X2aroRa3l1I/Vq1erXC7POa6trW3Z9GQ2y+m5MR/0o1aofjTqxLsXSDqcTN9Ydf8/SHqHJJnZqyVdMcOyX5H0DjNbk4y7uEE1ApjF5s2blclkZh2zcuVKbdmypUkVAThXjfqo2Y9r+nD9H0l6sOr+Tyf3FyQ9Kqkg6Xj1gu7+pJn9maSvmtnJZNyNDarzeczOHGxwn/WAAhC1bdu26e677571dflMJqOtW7c2sSqg9YTMlXMOeXc/JOnyqts3nmXeK6sW+1Dy/VlJm939WTPLanqv/fvJsp1V67lb0t3nWhuA+slmsxoeHlYul1OlUqkJ+0wmo0wmo9tuu03ZbDZglQBm0+z3yV8o6R/M7HFJ/1PSH7j7iSbXAGCeent7VSgUlM/na654l8/nVSgUdPXVV4cuEcAsGnW4fkbuXpbU08zHBLA42WxWAwMDGhgYeN68Z555JkBFAOZryVzxDgAA1BchDwBApAh5AAAiRcgDABApQh4AgEgR8gAARIqQBwAgUoR8yuHDh3Xffffp8OHDcw8GAGAOIXOFkE/p6OjQ2rVr1dHREboUAEAEQuYKIQ8AQKQIeQAAIkXIp5RKJR09elSlUil0KQCACITMlaZ+QE0ruOSSS05P83nyAIDFCpkr7MkDABApQh4AgEgR8gAARIqQBwAgUoQ8AACRIuQBAIgUIQ8AQKQIeQAAIkXIAwAQKUIeAIBIEfIp7q6HHnqIS9oCAOoiZK4Q8gAARIqQBwAgUoQ8AACRIuRTJiYmdPDgQU1MTIQuBQAQgZC5wufJp/T09Jye5uQ7AMBihcwV9uQBAIgUIQ8AQKQIeQAAIkXIAwAQKUIeAIBIEfIAAESKkAcAIFKEPAAAkSLkAQCIFFe8S1m/fr1OnDihVatWhS4FABCBkLnCnnxKqVTS8PCwSqVS6FIAABEImSuEPAAAkSLkAQCIFCEPoC6KxaL6+vrU3t6uFStWqL29XX19fSoWi6FLA5atRYe8mb3PzL5lZvfUo6Cq9d5iZr9Tz3XOx8jIiB5++GGNjIw0+6GBljU6Oqquri4NDQ2pXC7L3VUulzU0NKSuri6Njo6GLhEIJmSu1OPs+j5Jve7+9FwDzew8d//5fFbq7rsWXdkCXH/99dU1hCgBaCnFYlG5XE5TU1PPm1epVFSpVJTL5VQoFJTNZgNUCIQVMlcWtSdvZrskXSppn5ndamb/w8y+aWaPmtkNyZgbzew+MxuR9ICZbTSzr5rZ58zsO2b2F2b2bjP7hpkdMLNssly/mX1g0VsIoKF27NihSqUy65hKpaKdO3c2qSIApywq5N39FkklSddIukjSg+7+y8ntT5jZRcnQX5H0Hne/Nrn9Wkl/KOkKSVskvdLdXy9pSNJ7F1MTgObau3fvvEJ+z549TaoIwCn1vBjOmyRdX7X3fYGklyfTX3b3H1WN/aa7H5EkMytKeiC5/4Cm/0GYlZnlJeUlad26dRobG1t89TNo1HpbyeTkJH2oQj9qTU5OanJycl5jy+Vy1L3juVGLfsys2T2pZ8ibpLe7+8GaO82ulvTT1NifVU0/V3X7ufnU5O6DkgYlqaenxzdu3LjAkmfXqPW2krGxMfpQhX7UGhsb0+rVq1Uul+cc29bWFnXveG7Uoh8za3ZP6vkWui9Jeq+ZmSSZ2ZV1XDeAJWrz5s3KZDKzjslkMtqyZUuTKgJwSj1D/iOSMpIKZvZEchtA5LZt2zavkN+6dWuTKgJwyqIP17t7Z9XNm2eYv1vS7qrbY5LGqm5vnGmeu/cvtjYAjZfNZjU8PKxcLnf6LXOnZDIZZTIZDQ8P8/Y5IACueAdg0Xp7e1UoFJTP52uueJfP51UoFNTb2xu6RGBZ4qNmAdRFNpvVwMCABgYGQpcCIEHIp1x11VUql8tqa2sLXQoAIAIhc4WQT5mYmOCtHwCAugmZK7wmDwBApAh5AAAiRcgDABApQj5lcHBQIyMjGhwcDF0KACACIXOFE+9Sbr75zPV88vl8wEoAADEImSvsyQMAEClCHgCASBHyAABEipAHACBShDwAAJEi5AEAiBQhDwBApAh5AAAiRcgDABAprniXct111+nYsWNas2ZN6FIAABEImSuEfMrIyAifJw8AqJuQucLhegAAIkXIAwAQKUIeAIBIEfIp/f392r17t/r7+0OXAgCIQMhc4cS7lDvuuOP0NEEPAFiskLnCnjwAAJEi5AEAiBQhDwBApAh5AAAiRcgDABApQh4AgEgR8gAARIqQBwAgUoQ8AACR4op3KTfddJNKpZI6OjpClwIAiEDIXCHkUwYHB/k8eQBA3YTMFQ7XAwAQKUIeAIBIEfIAAESKkE/J5/Pavn278vl86FIAABEImSuceJdy1113nZ4eHBwMWAkAIAYhc4U9eQAAIkXIAwAQKUIeAIBIEfIAAESKkAcAIFKEPAAAkSLkAQCIFCEPAECkCHkAACJl7h66hkUxsx9K+n6dV7tW0tE6r7NV0Yta9KMW/TiDXtSiH7Xq3Y9fcvcXzzWo5UO+Ecxs3N17QtexFNCLWvSjFv04g17Uoh+1QvWDw/UAAESKkAcAIFKE/Mz4+Lkz6EUt+lGLfpxBL2rRj1pB+sFr8gAARIo9eQAAIrVsQ97MLjazL5vZd5PvLzrLuC+a2U/M7P7U/bvN7Gkzeyz5el1zKq+/OvTiFWb2SLL8Z81sVXMqb4xz6Md7kjHfNbP3VN0/ZmYHq54bL2le9fVhZm9JtuEpM/vgDPPPT37WTyU/+86qeX+a3H/QzN7czLobZaH9MLNOM/u3qufCrmbXXm/z6MUbzWy/mf3czHKpeTP+zrSyRfbjZNVzY19DCnT3Zfkl6eOSPphMf1DSx84y7tckbZJ0f+r+3ZJyobdjifTic5LemUzvkvQHobep0f2QdLGk7yXfX5RMvyiZNyapJ/R2LGL7V0oqSrpU0ipJj0t6dWpMn6RdyfQ7JX02mX51Mv58Sa9I1rMy9DYF7EenpCdCb0OTe9EpqUvSZ6r/Rs72O9OqX4vpRzJvstE1Lts9eUk3SLo7mb5b0m/MNMjdvyKp3KyiAllwL8zMJF0raXiu5VvIfPrxZklfdvcfufuPJX1Z0luaVF+jvV7SU+7+PXc/IeleTfekWnWPhiX9WvJcuEHSve7+M3d/WtJTyfpa2WL6EZs5e+Huh9y9IOm51LIx/s4sph9NsZxDfp27H5Gk5PtCDqn+mZkVzGynmZ1f3/KaajG9WCPpJ+7+8+T2DyRdUuf6mm0+/bhE0jNVt9Pb/VfJIbgPteAf+7m2rWZM8rM/runnwnyWbTWL6YckvcLMHjWzr5rZf2h0sQ22mJ/vcn1uzOYCMxs3s38ys4bsHJ3XiJUuFWb2d5J+cYZZt9Zh9X8q6f9q+hDNoKQ/kfThOqy3IRrYi5kCbMm/ZaMO/Zhtu9/t7ofNrE3S30jaoulDda1iPj/Ts41pyefDHBbTjyOSXu7ux8ysW9L/MrPXuPu/1rvIJlnMz3e5Pjdm83J3L5nZpZIeNLMD7l6sU22SIg95d//1s80zs/9nZuvd/YiZrZf0L+e47iPJ5M/M7K8kfWARpTZcA3txVNILzey8ZA/mpZJKiyy34erQjx9I2lh1+6Wafi1e7n44+V42s7/W9CG9Vgr5H0h6WdXtmX6mp8b8wMzOk/QCST+a57KtZsH98OkXXn8mSe4+YWZFSa+UNN7wqhtjMT/fs/7OtLBFPd/dvZR8/56ZjUm6UtOv8dfNcj5cv0/SqbM73yPpf5/Lwskf/1OvSf+GpCfqWl1zLbgXyR+xhySdOmv0nHu5BM2nH1+S9CYze1Fy9v2bJH3JzM4zs7WSZGYZSdep9Z4b35R0WfKuiVWaPpEsfeZvdY9ykh5Mngv7JL0zOdv8FZIuk/SNJtXdKAvuh5m92MxWSlKyt3aZpk84a1Xz6cXZzPg706A6m2XB/Uj6cH4yvVbSv5f0f+peYeizE0N9afr1sq9I+m7y/eLk/h5JQ1Xjvibph5L+TdP/tb05uf9BSQc0/Qd8r6TVobcpYC8u1fQf8qck3Sfp/NDb1KR+/F6yzU9J+t3kvoskTUgqSHpS0ifVgmeXS3qrpO9oeq/i1uS+D0u6Ppm+IPlZP5X87C+tWvbWZLmDknpDb0vIfkh6e/I8eFzSfkmbQm9LE3rxy8nfh59KOibpyapln/c70+pfC+2HpDckGfJ48v33G1EfV7wDACBSy/lwPQAAUSPkAQCIFCEPAECkCHkAACJFyAMAEClCHgCASBHyAABEipAHACBS/x/RKSwuU9WdgwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1176f96d8>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "names = [\"fermi\", \"magic\", \"veritas\", \"fact\", \"hess\"]\n",
    "\n",
    "for y, name in enumerate(names):\n",
    "    x = multi_norm.mean[f\"z_{name}\"]\n",
    "    ax.plot(x, y, ls=\"-\", marker=\"o\", ms=10, color=\"k\")\n",
    "\n",
    "ax.set_ylim([-1, 5])\n",
    "ax.set_yticks(list(range(5)))\n",
    "ax.set_yticklabels(names)\n",
    "ax.grid()\n",
    "ax.axvline(0.15, color=\"k\", lw=3, ls=\"--\")\n",
    "ax.axvline(-0.15, color=\"k\", lw=3, ls=\"--\");"
   ]
  }
 ],
 "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}