Fitting/hepstat_tutorial_fit.py.nbconvert.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "middle-depth",
   "metadata": {},
   "source": [
    "# Hepstat_Tutorial_Fit\n",
    "Fit example with ROOT, with following objectives:\n",
    " * Fit a histogram with a linear chi-squre fit, and compare results with by-hand calculations\n",
    " * Different fiting options\n",
    " * Compare chi-squre fit and likelihood fit\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "**Author:** Lailin XU  \n",
    "<i><small>This notebook tutorial was automatically generated with <a href= \"https://github.com/root-project/root/blob/master/documentation/doxygen/converttonotebook.py\">ROOTBOOK-izer</a> from the macro found in the ROOT repository  on Friday, April 16, 2021 at 06:48 AM.</small></i>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "variable-happiness",
   "metadata": {},
   "source": [
    "Import the ROOT libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "parental-glucose",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-15T22:48:36.664401Z",
     "iopub.status.busy": "2021-04-15T22:48:36.663650Z",
     "iopub.status.idle": "2021-04-15T22:48:38.284309Z",
     "shell.execute_reply": "2021-04-15T22:48:38.284520Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to JupyROOT 6.22/07\n"
     ]
    }
   ],
   "source": [
    "import ROOT as R\n",
    "from math import pow, sqrt\n",
    "R.gROOT.SetStyle(\"ATLAS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "presidential-faculty",
   "metadata": {},
   "source": [
    "Linear least square method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "natural-adrian",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-15T22:48:38.289495Z",
     "iopub.status.busy": "2021-04-15T22:48:38.289194Z",
     "iopub.status.idle": "2021-04-15T22:48:38.396432Z",
     "shell.execute_reply": "2021-04-15T22:48:38.396781Z"
    }
   },
   "outputs": [],
   "source": [
    "def lls(l_x=[], l_y=[], l_ye=[],):\n",
    "  \"\"\"\n",
    "  A simple linear regression, y=m*x+b\n",
    "  \"\"\"\n",
    "  sum_1, sum_x, sum_x2 = 0, 0, 0\n",
    "  sum_xy, sum_y, sum_y2 = 0, 0, 0\n",
    "  \n",
    "  for i in range(len(l_x)):\n",
    "    ye = l_ye[i]\n",
    "    x = l_x[i]\n",
    "    y = l_y[i]\n",
    "    if ye ==0: continue # ignore empty bins\n",
    "    sum_1 += pow(ye, -2)\n",
    "    sum_x += x*pow(ye, -2)\n",
    "    sum_y += y*pow(ye, -2)\n",
    "    sum_x2 += pow(x/ye, 2)\n",
    "    sum_xy += x*y*pow(ye, -2)\n",
    "    sum_y2 += pow(y/ye, 2)\n",
    "\n",
    "  m = (sum_1 * sum_xy - sum_x * sum_y) / (sum_1 * sum_x2 - pow(sum_x, 2))\n",
    "  b = (sum_y * sum_x2 - sum_x * sum_xy) / (sum_1 * sum_x2 - pow(sum_x, 2))\n",
    "\n",
    "  # Error\n",
    "  v_m = sum_1 / (sum_1 * sum_x2 - pow(sum_x, 2))\n",
    "  v_b = sum_x2 / (sum_1 * sum_x2 - pow(sum_x, 2))\n",
    "  corr_mb0 = - sum_x / sqrt(sum_1 * sum_x2)\n",
    "\n",
    "  # Error matrix\n",
    "  mat = R.TMatrixD(2,2)\n",
    "  R.TMatrixDRow(mat, 0)[0] = sum_x2\n",
    "  R.TMatrixDRow(mat, 0)[1] = sum_x\n",
    "  R.TMatrixDRow(mat, 1)[0] = sum_x\n",
    "  R.TMatrixDRow(mat, 1)[1] = sum_1\n",
    "  mat_inv = mat.Invert()\n",
    "  e_m = sqrt(R.TMatrixDRow(mat_inv, 0)(0)) \n",
    "  e_b = sqrt(R.TMatrixDRow(mat_inv, 1)(1)) \n",
    "  corr_mb = (R.TMatrixDRow(mat_inv, 0)(1)) / (e_m*e_b) \n",
    "\n",
    "  print(sqrt(v_m), e_m, sqrt(v_b), e_b, corr_mb0, corr_mb)\n",
    "  return [m, e_m, b, e_b, corr_mb]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "stretch-legislature",
   "metadata": {},
   "source": [
    "Generate some toy data, assuming the model is y = m*x + b. Then define a histogram to save the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "competitive-challenge",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-15T22:48:38.537623Z",
     "iopub.status.busy": "2021-04-15T22:48:38.537266Z",
     "iopub.status.idle": "2021-04-15T22:48:38.740309Z",
     "shell.execute_reply": "2021-04-15T22:48:38.740583Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAI8CAIAAAD0vjrdAAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dS3LbOLsGYOrUP025PU+6eiekNpOLvYJkImqU3kD/TnkzIjeSysDz2CvQGeCYh60LdTEsgPTzDFI2RVufYNl8A4DAbL1eFwAAvMz/pC4AAGAKhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACP7zSt+3aZqqql7pm0/GbDZLXQIATNB6vb78k85e6VmrqmrbtizLuq6lq31ms8Ptf8w5cU+b/Le6/DMqPvq3uvwz5vmtLv+MeX6ryz+j4qN/qyTPGN3rDv+1bTufz2ezWVVVTdO86nMBACT0WqGqaZrValWWZfi0S1d1Xb/SMwIAJPSKPVWhd2q9XvfT1XK5DB1X0hUAMCWXuPuvS1eLxSIcadu2S1eGBQGACbjokgp1Xa/X6/V6bVgQAJiYi4aqpmnqup7NZm3bbjwUOq5EKwBgpC5xz2HTNE3TLJfL/sGyLLuZVXVdd4+uVqu3swTD8DpVSW4HPV6qG1ajUHwS4628UHw6ik8l8+IPLvQ4qXWqiud+qY1OqX6W2ixlNgsnvJ1ZVpm/ZYcpPpXxFj/eygvFp6P4VBR/htdaUT0s/tl9OpClOovFYqM3CwBgLF4rVAXHZKlOuBPQtCoAYIxeq39sGnv/hdlgw+dUVbX9SsMX7nyoT+dqKopPYryVF4pPR/GpKP6c5x1vk13AMRseLxaLfu9af9J9MDBLzFs2FcUnMd7KC8Wno/hUFH+Giy6pMDrdelpH6ieq/lpcE+i0AwCGjTiHptU0zXw+L3o3bXZH+l1TXczauVTEqP8fAAB5Mvw3MmFksB+VuhseN5p03/FCqAKAV2D4b0xCkFosFv3Op5Ccuv0NO+5nBIC3QE/JybYH/oLtvquDD+mpAoDo9FSNRkhUq9Wqf7CbRDUwJ/3trBQPAG+QUHWaMJYXFjVNXAoAkBOh6jThVr7taVLDvVBheYWd58zOFekFAUCmRneJfN1taibmNbqpzKkCIDchlyS/Qp1dQKpcpafqBPu6qYrBqVTF842BRgwBYMKEqmMdMxUdAHizhKpjhQ6q7WWogi5pDUyuksYAYMKEqmOFIbyDK3luhypdXADwFghVRzlmVfR9t/hZngoA3gKh6ighGIXYtE8IXm3b9hNY0zRhevu+cUMAYBpsk3KUcHPmYrEY7rLq9k4Oyy40TRM+LfbcF2qbGgAylMmSCmdLdXl1UT9s32Z/O3W5qlOW5b4RQKEKgAwJVWc+73ibLFvNs+rZvjOFKgAyJFSd+bzjbbIJEKoAyJBQdR7b1CQ2sJT+eN/NAPCq8twDV6hKTHICgFMNXz3t/QcA/zKbzfLskICdhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAjs/ZeYDZUB4FR57l8kVCUmOQHAqWyoDAAwWUIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAE9v5LzIbKAHAqGyqzg+QEAKeyoTIAwGQJVQAAEQhVAAARCFUAABEIVQAAEQhVAAARCFUAwP95eHj49u1b+Pj6+vrr168PDw9pSxoRoQoAKIqieHh4+PTp0+/fv8Onj4+PT09PHz9+lKuOJFQBAEVRFP/973///PPPHz9+dEfu7u7++uuvf/75J2FVIzKzondCs5n2B9grrIvt7+TFXF9fPz4+bh+/urraeTxbqS6vLuopDa+j70cDvHFC1SU9PT398ccf+x79/fv3wKOXd3AjGqHqzdFTBTBAqLowPVUvZENlAKAoiuLz589PT093d3f9gzc3N+/evUtV0riYqA4AFEVR3N7e/vr168uXL92Rm5ubnz9/3t7eJqxqRISqEzRNU9d10zQRzwSATLx///7+/v7q6ip8enV19e7du/v7+w8fPqQtbCzM6TlKXdfL5bJ/pCzLnZnp+DMLc6oABplTlcrYWz7V5VVP1WFVVXU5qSzL8EHbtlVVbZzZT1TDZwIAEyNUHVDXddu2RVEsFov1et00zXq9Xq1WRVG0bdvvgmqaJiSqsiy7MxeLxfaZAMD0GH46IHSBLhaLuq77x6uqatu2P7QXjhRb/aX7jheG/wAGjX0QarzG3vKG/3LUBaaNRBUe2khaXYfWxpnbXwvAMNv6MkZC1ZDhMbu6rrcnSw0cMQIIcAzb+jJSQtWQEIO6KecDCyX0BwGHvxsAw2zry0gJVYdVVVXX9Ww2m8/ny+VyPp/PZjM39AG8kru7u36i6g5urPQNuTFResjGfo2hyyrMnQq61usWU9jZntuz2nd+/5P4wQGTNK5tfacqk4nqo7tK2vvvWP0fT9M08/m8KIqqql44qJf8LQuQlaurqz/++GPftr4S1Zty9iXyJWnsJQz/HWV7lYRuAaruyMCXh9OMGAIc4/Pnz/3t54Kbm5vPnz8nqQeOJFQNCeN93UT1Pvf0AbwS2/oyUkLVmTa6nY7JWHqqAI5hW19GSqgaEmJQf2Z6Z98aCtuh6pjVFgDoe//+/ffv38PHj4+Pf//9t0RF/oSqISctOhVGCQdCFQAwYULVkC5UbdxH0O2d3N+UJmxH07Ztf1+anWcCANNjnaoDutUTiqIoyzKsodANCO7bO/ngmYENlQEGZLJa0hs09pZPdXl1UT+sn6s62yt5Bl2uOnhmIVQBDBr7pX28xt7yQlXumqbpTzkfnm4VVM/2nSlUAQwY+6V9vMbe8kLVWyRUAQwY+6V9vMbe8qkuryaqAwBEYO+/xAb2JxrvfxEA4FWl2t1vmFCVmOQEAKcavnraUBkAYMSEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACOz9l5gNlQHgVDZUZgfJCQBOZUNlAIDJEqoAACIQqgAAIhCqAAAiEKoAACIQqgAAIhCqAAAiEKoAACIQqgAAIhCqAAAisE1NYvb+AyA3+V+A7P3HDvm/cQEgN/b+AwCYLKEKACACoQoAIAKhCgAgAqEKACACoQoAIAKhCgAgAqEKACACoQoAIAKhCgAgAtvUAJApG3kxLkJVYjZUBoBT2VCZHSQnADiVDZUBACZLqAIAiMDw32F1Xe97qKqqqqq2jzdN0zTNvkcBgOkRqg5omma5XA6csBGb6rruzg8flGXZNM1r1QcA5MHw3wHDeWggUZVlGT5o21Z/FQBMnlB1lLIs17v001LXpxVObppmvV4vFouiKNq21VkFANMmVB0QwtAxXU3d1Kt+fqrrOnRZzefz+MUBANkQqg5o27Y4LlSFM0PXVN/APHcAYDKEqsi241d3xAggAEyYUHWsWU9VVRsJqft0oE9LqAKACROqhnQxaGNGVNu28/ncuB4A0BGqhvT7llarVbjjb7VahSPL5bI7YbgXKsxV33nO7FxRXiAAZGt0l0iLfw4Jw3xVVfU7paqqWq/X4Wc2n89fuCOyDZUBYKezL5E2VM5RCFU7h/m6/qruzIHvc/wthADASAlVZ3JPHwDQJ1TFcUzG0lMFABMmVA2p63rfLX771lDYDlXHrLYAAIydUDVkuVwul8udYWg7PO27xc/4IAC8BULVkH3bIXd7J/enq4c+rbZt+51b3Znb29cAAFMyc0v/sO62zLIsQ5dV0zThbr6yLDfCVlVV3UPhzsHwabHnvtDZTPsDQGSpLq8u6od1UalvO1HtO3nfmYVQBQCvQKjKWtM0/fnmw1POm2fVs31nClUAEJ1Q9RYJVQAQXarLq21qEhtYSl/eAoCd8twDV6hKTHICgFMNXz3t/QcAMGJCFQBABEIVwAGz2SzPCRxAVoQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIhCoAgAjs/ZeYDZUB4FR5rscrVCUmOQHAqWyoDAAwWUIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAE9v5LzIbKAHAqGyqzg+QEAKeyoTIAwGQJVQAAEQhVAAARCFUAABEIVQAAEQhVAAARCFUAABEIVQAAEQhVAAARCFUAABHYpiYxe/8BwKns/ccOkhMAnMrefwAAkyVUAez18PDw7du38PH19fXXr18fHh7SlgRkS6g6R13XTdMMnNA0zcFzgMw9PDx8+vTp9+/f4dPHx8enp6ePHz/KVcBOM3N6TlVVVdu2xZ4B3bqul8tl/0hZlvvS1Wym/SFf3759+/37948fP/oHv3z5cnV19f3791RVAQelury6qJ+maZr5fB4+3m66fqIqyzJkr2J/rhKqIGfX19ePj4/bx6+urnYeBzKR6vJq+O80XaLa1jRNSFRlWa7X66Zp1uv1YrEoiqJtW0OBMC5PT0/7ktPAQ8BbpqfkBN3AX7DRdPuGBQeGC/VUQc70VMFI6anKXV3XbduWZRk6n7aF5LT9aF3Xr10b8Bo+f/785cuXjYM3NzefP39OUg+QOaHqKN3Q3sFRvKqq9h0xAgjjcnt7++vXr36uurm5+fnz5+3tbcKqgGwJVUcJU6lWq9W+E7rAtB2qts8BRuH9+/f39/dXV1fh06urq3fv3t3f33/48CFtYUCebFNzWMhJZVkOBCZgkt6/f//9+/e///67KArzqIBheqoOCFOpikP9TMOPlmW575zZuV7wmgBgBEZ3idRTdUCYSjUw8PdC7v4DgJ3OvkTaUDlH4aeyWCwODvwNnxD6uoweAsCE6anaqz9a118WoTseDlZVJS0BAELVYRt7+W0f7+eqpmn2ZSzZCwAmTKgaEiaYb+jv6FdsRaXtUHXMagsAwNgJVXtVVbXzfr1u1+SNR8MOyttfYnkqAHgLTFSPJkyxatt2YwJWSGD7NrcBAKbBhr4n63qqtpuu2zs5rBTaNE03VriznW2oDKMQbgT22wpjYUPlKWiaJky0att2uVx2AcvfYgCYPD0l8TXPqmf7ztRTBaOgpwrGJdXl1UU9JaEKRkGognEx/AcAMGKWVEhsYH8i/y0GgJ0S7po8QKhKTHICgFMNXz1tqAwAMGJCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQAS2qUnM3n8AcCp7/7GD5AQAp7L3HwDAZAlVAAARCFUAABEIVQAAEQhVAAARCFUAABEIVQAAEQhVAAARWPwTxiSsaGfN2AvT4MAx9FQBAEQgVAEARGD4LzEbKgPAqWyozA6SEwCcyobKAACTJVQBAEQgVAEARCBUAQBEIFQBAEQgVAEARCBUAQBEIFQBAEQgVAEARCBUAQBEIFQBAERg77/EbKgMAKeyoTI7SE4AcCobKgMATJaeqmM1TRP+raqqqqqDJx95JgAwDTPDTwfVdb1cLjcOLhaLuq6PObksyxDIts1m2p/ThD5tbxuAAakur4b/DuiHpLIsy7IMHy+Xy+1QtXFy+KBtW/1VADB5ekoOCB0DG/1SVVW1bVsUxWq16gJT0zTz+bz4d9dUF7P6Z/a/ufbnJHqqAA5KdXl1UR/SRaLtVtoOW13S2jh53/FCqOJ0QhXAQYb/chQ6nLqBvH0nBCE5LRaLjXN2Tr0CACbG3X9DqmcDJxx/JNwPGK04ACAnQtWQfZ1MXQdVPzBtHNn5VUIVAEyVUHWCkLG6+/sWi4WQBAAEQtWxmqbpL0C1cT/gvpWogrIs27bdec7ZS+mbqgzAtOW5wd8AE9WPVVXVYrHolqpaLpdRftjrc738qQEgZ6O7RApVJ6jrOuw/s1qtwpFu+G94HDDcGGisEAAmTKg6R+i1Kp7TElzAw8PDt2/fwsfX19dfv359eHhIWxIAfULVkNAvtfOhjW6n7dsAD34JHO/h4eHTp0+/f/8Onz4+Pj49PX38+FGuAsiHUDVkPp/P5/PhSegbtk8+ZrUFGPbf//73zz///PHjR3fk7u7ur7/++ueffxJWBUCfUHXYzlC1fTBMYB8IVXC2u7u7fqLqDt7d3SWpB4BtQtWQMHFquVxuBKNueYX+pjRhhYW2bTeWWtg+E07y9PT0+Ph46kMAXJgNfQ/o1k0oyzKM3zVN081P37d3cjh54Mzum2t/jnF9fb0zPF1dXQlVABtSXV5d1A/rolJfWZY7x/W2T953ZiFUcbSvX78+PT1tDPbd3Ny8e/fu77//TlUVQJ6Eqqz1bwMc3mK5Ozns9Dd8slDFkR4eHj5+/PjXX391uerm5ubnz5/39/cfPnxIWxtAboSqt0io4ngPDw///PNP6Je6urr6/Pnz7e2tRAWwTah6i4QqThUm+XnbAAxIdXm1oXJiAxsIunACwE557rUsVCUmOQHAqYavnqkil3WqAAAiEKoAACIQqgAAIhCqAAAiEKoAACIQqgAAIhCqAAAiEKoAACIQqgAAIhCqAAAiEKoAACKw919iNlQGgFPZUJkdJCcAOJUNlQEAJkuoAgCIQKgCAIhAqAIAiECoAgCIQKgCAIhAqAIAiECoAgCIQKgCAIhAqAIAiMA2NYnZ+w8ATmXvP3aQnADgVPb+AwCYLKEKACACw38wJsaLAbKlpwoAIAKhCgAgAqEKACACoQoAIAKhCgAgAqEKACACoQoAIAKhCgAgAot/JmZDZQA4lQ2V2UFyAoBT2VAZAGCyhCoAgAgM/x2raZrwb1VVVVUdPPnIMwGAaZiZ03NQ0zTz+XzjYFmWIWZtqOt6uVwec2ZRFLOZ9geAyFJdXg3/HdBPVGVZlmUZPm7bdnseXD9R9c/UXwUAkydUHdAlqvV6HQb11uv1arUKB/tpqWmakKjKsgwnr9frxWJRFEXbtvs6qwCAaTD8NKTrptpupa5Tqnuoqqq2bbdP3ne8MPwHAK/A8F+OBrqXuj6q7pyQnELXVF9d1/ErAwAyI1QNCYGpmx3Vt2+a1Pbx7fgFAEyPJRWGDMSgrv8pZKbuzIE56WGRhUilAQB50VN1jv6c9NS1AABZEKpOVtd1d0tg10E1PLQXstfOc2bnivqaACA7o7tEGv47wcYqoFHuLHD3HwDsdPYl0obKuauqqktUi8Vie92Ega8NNwaaUAUAEyZUHdY0zWw261ZMWK/XVkkAADYY/jugG/Ib2MKv+Pe6CcevtgAATIaeqgOOSVR926cds9oCADB2QtWQbpjvmES17xY/a34CwFtg77kh4faBsiwHJlH1B/5Ct9Ziseinse2D/e+v/QEgrlSXVxf1Icfck9lvwG7v5LIsq6pqmiZ8Wuy5L1SoAoDohKocnRqqil6u6gzMxxKqACA6oWo6mmfVs31nClUAEJ1Q9RYJVQAQXarLq3WqEhsYYZS3AGCnPPfAFaoSk5wA4FTDV097/wEAjJhQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgVAFABCBUAUAEIFQBQAQgb3/ErOhMgCcyobK7CA5AcCpbKgMADBZQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQAT2/kvMhsoAcCobKrOD5AQAp7KhMgDAZAlVAAARCFUAABEIVQAAEQhVAAARCFUAABEIVQAAEQhVAAARCFUAABEIVQAAEdimJjF7/wHAqez9xw6SEwCcyt5/AACTJVQBAEQgVJ2maZq6ro88rWmaVy8IAMjDzJyek4Rh2oFGq+t6uVz2j5RluS9dzWbaHwAiS3V51VN1goM9T/1EVZZl+KBt26qqXrMuACA9oepYdV3P5/OBE5qmCYmqLMv1et00zXq9XiwWRVG0bWsoEACmzfDTAU3TbGepnY1WVVXbttuP7jtejHn47+AwKACkYvhv9EJyCl1TfcdMbAcAxk6oOqCqqvWz1Wp1zPn7jhgBBIAJE6ri6ALTwJx0oQoAJkyoAgCIQKiKY7gXKiyvsPOc2ble53UAQC5Gd4m0oXJibqADgJ3OvkTaUHnchpf3DDcGWgIUACZMqAIAiECoiuOYdRP0VAHAhAlVkW2HqmNWWwAAxk6oimbfLX6WpwKAt0CoiiZsR9O2bX9fmm6X5e3tawCAKRnrhr5JdJsr72u0bu/ksiyrqmqaJny670tsqAwA0dlQeQqapgmDgG3bLpfLLmBNKXw8PDx8+/YtfHx9ff3169eHh4e0JQFADsbaU5Kz5ln1bN+Zo+upenh4+PTp059//vnjx49w5MuXL79+/bq/v3///n3a2gAgSHV5HdlFfWJGF6q+ffv2+/fvLlEFX758ubq6+v79e6qqAKBPqHqLRheqrq+vHx8ft49fXV3tPA4AlydUvUXDmxPl9qN5enr6448/9j36+/fvgUcBIKKDu/sJVW+OnioAiC7V5fU/l39Kxuvz589PT093d3f9gzc3N+/evUtVEgBkwpIKnOD29vbXr19fvnzpjtzc3Pz8+fP29jZhVQCQA6GKE7x///7+/v7q6ip8enV19e7du/v7+w8fPqQtDACSG9mcnokZ3ZyqjhXVAciWFdUBAEZMqAIAiECoAgCIQKgCAIhAqAIAiECoAgCIwIrqiQ3sXmTBAgDY6eDef0kIVYlJTgBwquGrZ6rIZfgPACACoQoAIAKhCgAgAqEKACACE9U5h/n1ALBBTxUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVZ8pz36UjKT6J8VZeKD4dxacy6uJTsaRCYjZUBoBT5Zn5hKrEJCcAOJUNlQEAJkuoAgCIQKgCAIhAqAIAiECoyt2Rs+0inhZxft/lq1J89GdU/EkuX9V4i487lVjx0b9VrKeLe1qeN/11hCoAgAiEKgCACIQqAIAIhCoAgAiEKgCACIQqAIAIZvaeSyjzW0MBYKSSxBuhCgAgAsN/AAARCFUAABEIVQAAEQhVAAARCFUAABH8J3UBb1fTNE3TVFVVVVXqWv7lvMLCV+17tK7rF9cVX6g5q9rOKGkULR8qzPANf3Zho2j2QstnoK7r3Bq/c1Jt+bf8QA2X+xGsubiyLDd+CqvVKnVR6/XLCtv+2vzfZhnWdkZJmbf8arXarqosy9R1vbSwzJt9vV4vFovtwhaLReq6XlpY/i3f11WbupAdTq0t85bf+Rt9+Xd++oZ4a/o/5v57NHmuemFhmf++bet+A1MX8v/OKynnlu//mSvLcqPUUReWc7Ov/x1cNl5g2lz18sIyb/m+/tssdS2bzqgt85bfGdY7F7vCpm+IN6V7U/Z/wDm8KV9e2IXfuC/U//VLXcv/ObuknFt+5yvq/pon7K96eWE5N/v6ubyNmLLz1/zCXl5Y8pdwvKwyx4Yzasu85cOf0OS94Nn9pKdt5x+UdQZv1pcXlucfjg07+4fHXlIOr2KngY63LkFevqp1pMKybfb14KvY95t+GVEKy7nl+/Lpl912Xm0ZvpC+8KKSD3C7++9yuil++ybTpZrol21hjNrAnNZuxujAOa8n28JiCcUPDNakenXZFhZdXddt25ZlOTwmlUTOtb1E27ZF71c4FaHqcgb+XoQ3d3hPXN7LCxvLn8Kqqrr/TwzParyYF5aUc8sPXEHT/uF7eWE5N3tRFFVVLRaL4TuhLlfNv5/3hYVl3vJB0zTL5bLIstqza8vwteRJqLqcyV9jyrKs63rWk/w/DZOXc8s3TbNer3f+Le4uq0nqfHlhOTd7URR1XYdb5TeOdy85VZ0vLyzzlg/m83lRFJn8t23D2bWNouWDjfIuHQcvP+L4ZoUG3zmNLu2daC8vLJPbLk4yjbv/Rt3yyaeUbji+sHE1+2Kx6BecfNJJ54zC8m/58L/T7i2UdvrghpfUlnnLZ7KewtpE9Usa+OnmEKpeUljXy9X/JjnfTrxO3eY7nVHS6Fo+w/sug5MKG1Gzb1xs8klU5xWWectvx5R8QtULaxtLyxe9hNcvz5IKEzTwt2PsoWqxWIRpj/u+PJ8/5Z1phKoRtfzGRTR1Of/vjMJG1Ozr52r74/upK/o/ZxSWectvX7/zCVUvrC3zll+tVjvLW1980aL0P+m3Y+CGz7QX+FctbKPDOR/TCFUDsmr5fBaf3BC9sKyafcMEBl4HJG/5ne+iTELVq9aWvOWHXfhPvb3/Li3beyheqbCqqtq2TXVj41uWScs3TRPmxhZFMXzn14W9UmGZNPtO4ea75XKZW3lRCkvb8v2/n/330saCNUk2AXzt2nJ+zxf/XiflAo0vVF1OeOftfCht0sq2MMauCy5lWWb1Xsq2sCjCK9p5/aiqKtxRn0S2hUW071V0xxPeLpdzbZMhVF3azvhycE28Czi7sG7r8oH/66d9aVOVf8tnG1xeUthYmn21WuV2jXxhYZm3/M6n7v6uhkdT/UReWFvmLV88F7azvEuvJHKZUUaC0ObbtyGE4wnnmrywsH1ffvx3uLxpzKnKueUzmU2y7eWF5dzswzWk/aG8vLDMW35btr8F6xNry7zlQw0753Vd+EeQ4096wrosf/DghZ1UWLhhp/8rNHBm8pe2zxhD1bhavvszt9pvFIWNq9nXvavIxgtJfqPWqYWNruW3jTRUja7lD761LKkwWcWzxWIRblLNIeafVFh3TeqO9G9KD7+K2d7t1RljqBpXyxdHGEVh42r2YKO8jQrHUtgYW37DSEPVGFt+4K11yTsTc/xJT9vGojj5vCOPL2zn23Tnl2fy0naaRqhaZ9zyO6vaMIrCxtXsnZ1zXHK47/34wkba8n2TCVXrMbR8Du/52fq4PzHE1c37Kwan/l3eCwtrnlXPolbHXlo+icybvf/rnFV5Ly8s85afsMxbPvl7XqgCAIjgf1IXAAAwBUIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFfPsgCIAAAHaSURBVABABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQBABEIVAEAEQhUAQARCFQAwGlVVVVXVNM3G8aZpwkMJanomVAEAo1FVVdu28/l84/h8Pm/bNm2omq3X64RPDwBwkpCryrLs+qvqul4ul0VRpE01eqoAgDEJWapt27quw6chUa1Wq6R16akCAMamaZowArher2ezWVEUi8UiZKyEhCoAYHzCIGD4uD8UmJBQBQCMUuijKlJPpeqYUwUAEIFQBQCMT3/1hLQrKXSEKgBgZJqmCROqwh1/3Z2AaZlTBQCMTJhNtVqtqqrKZJGqQk8VADAuYbCvLMvwQV3XZVkWGQwCClUAwGjUdR0G/vprKGwsB5qKUAUAjEYY6VssFhvHw5HwaCrmVAEARKCnCgAgAqEKACACoQoAIAKhCgAgAqEKACACoQoAIAKhCgAgAqEKACACoQoAIAKhCgAgAqEKACACoQoAIAKhCgAgAqEKACACoQoAIAKhCgAgAqEKACACoQoAIAKhCgAggv8Fis/+ZWGAkyIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = 11.\n",
    "b = 2.5\n",
    "bin0, bin1=0, 5\n",
    "nbins = bin1-bin0\n",
    "\n",
    "fl = R.TF1(\"flinear\", \"pol1\", bin0, bin1)\n",
    "fl.SetParameters(0, b)\n",
    "fl.SetParameters(1, m)\n",
    "\n",
    "hname = \"test_data\"\n",
    "hd = R.TH1F(hname, hname, nbins, bin0, bin1)\n",
    "hd.Sumw2()\n",
    "hd.GetXaxis().SetTitle(\"x\")\n",
    "hd.GetYaxis().SetTitle(\"y\")\n",
    "hd.FillRandom(\"flinear\", 200) \n",
    "\n",
    "myc = R.TCanvas(\"c\", \"c\", 800, 600)\n",
    "myc.SetFillColor(0)\n",
    "\n",
    "myc.cd()\n",
    "hd.Draw()\n",
    "myc.Draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "egyptian-spotlight",
   "metadata": {},
   "source": [
    "Do a fitting. By default ROOT uses Chi-square fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "selective-argentina",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-15T22:48:38.800431Z",
     "iopub.status.busy": "2021-04-15T22:48:38.800112Z",
     "iopub.status.idle": "2021-04-15T22:48:38.936209Z",
     "shell.execute_reply": "2021-04-15T22:48:38.936479Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "****************************************\n",
      "Minimizer is Linear / Migrad\n",
      "Chi2                      =      2.43795\n",
      "NDf                       =            3\n",
      "p0                        =      1.71809   +/-   3.1885      \n",
      "p1                        =      15.1177   +/-   1.66766     \n",
      "\n",
      "****************************************\n",
      "Minimizer is Linear / Migrad\n",
      "Chi2                      =      2.43795\n",
      "NDf                       =            3\n",
      "p0                        =      1.71809   +/-   3.1885      \n",
      "p1                        =      15.1177   +/-   1.66766     \n",
      "\n",
      "Covariance Matrix:\n",
      "\n",
      "            \t          p0          p1\n",
      "p0          \t      10.167     -3.9292\n",
      "p1          \t     -3.9292      2.7811\n",
      "\n",
      "Correlation Matrix:\n",
      "\n",
      "            \t          p0          p1\n",
      "p0          \t           1    -0.73894\n",
      "p1          \t    -0.73894           1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Info in <TCanvas::Print>: png file test_fitting_1.png has been created\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAI8CAIAAAD0vjrdAAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dQXLbRrs2UPDWP3Upnjup7ATQZhJbXoE9IThKNpDPLm9GwEZSqbqa21oB/0F/wkVIEKRIEN0NnjNwSSAkNVsU8fjtRvdqu90WAABc5n9iNwAAYAmEKgCACQhVAAATEKoAACYgVAEATECoAgCYgFAFADABoQoAYAJCFQDABIQqAIAJCFUAABMQqgAAJiBUAQBMQKgCAJiAUAUAMAGhCgBgAkIVAMAEhCoAgAkIVQAAExCqAAAmIFQBAExAqAIAmIBQBQAwAaEKAGACQhUAwASEKgCACfy/2A0AAG7CavXfD7bbqO24mmuFqqZpqqq60jdfjFX3+gKAmzHD1W8bI7itrvRTq6pq27Ysy7qupatDVqvj/X/KOdOetvhvNf9P1PjJv9X8PzHNbzX/T0zzW83/EzX+td+qKP512uAXzd/zk7vunKq2be/v71erVVVVTdNc9WcBAKnZL0otdeyvuF6oaprm8fGxLMvwaZeu6rq+0k8EAJJyU4mquGqlKlSnttttP11tNptQuJKuAGDBbi1RFfMsqdClq/V6HY60bdulK8OCALAwO4lqu11+oipmXqeqruvtdrvdbg0LAsAirVZz3NyXpllDVdM0dV2vVqu2bXceCoUr0QoA8nWDQ359c9xz2DRN0zSbzaZ/sCzLbmZVXdfdo4+Pj7ezBMP4OlVRbgc9XawbVieh8VHk2/JC4+PR+Fhe2/jBC9r1nv3RhR4XtU5V8VKX2ilK9bPUblNWq3DC7cyyuqm/t6RofBT5trzQ+Hg0PpZXNT61AlWsnr/Wiuph8c/u05Es1Vmv1zvVLAAgZTMXqBJ33b3/TslSnXAnoGlVAJCF1ApU0V2rPraMvf/CbLDxc6qq2n+m4QsHH+q7ncpwajQ+inxbXmh8PBofy9HGp5yoYvV8xr/vGZyy4fF6ve5X1/qT7oORWWLL/ntLmcZHkW/LC42PR+NjGWl8ynEqWObef7nr1tM6UT9R9dfiWkDRDgCKHBJVRBmH6Liaprm/vy96N212R/qlqS5mDS4VkfV/YgC4NbkkKsN/mQkjg/2o1N3wuNOlh44XQhUAmcglTgWG/3ISgtR6ve4Xn0Jy6vY37LifEYCs5ZWoIlIpebX9gb9gv3Z19CGVKgASl2OiUqnKRkhUj4+P/YPdJKqROem3s1I8AMuwk6i22wwSVURC1euEsbywqGnkpgDANe0nKsYJVa8TbuXbnyY1XoUKyysMnrM610RPCAB2rVZJJKrsLpHX3aZmYa5RpjKnCoCkpDOJ6uxLZKxcpVL1CofKVMXoVKri5cZAI4YAJC6dRJUjlapTnTIVHQDylcKQX9ZUqk4VClT7y1AFXdIamVwljQGQpv1JVEVh8u6rCVWnCkN4R1fy3A9VSlwApGxoApJEdQ6h6iSnrIp+6BY/y1MBkCyTqCYkVJ0kBKMQmw4Jwatt234Ca5omTG8/NG4IAFEMrpsgUV3CNiknCTdnrtfr8ZJVt3dyWHahaZrwaXHgvlDb1AAQxXiBKlz18r1Cxbq8uqgfd2izv0FdruqUZXloBFCoAmBmg0s47VyLhKozf26+XZas5kX14tCZQhUAczpxBpVQdebPzbfLFkCoAmAepxSoeicLVeew+GdkI0vp5/tqBiApy7vFL809cIWqyCQnAK5qeYmqOHb1jBW5hCoAEpX7IFR0i4xTKbNOFQAskEQ1P6EKAJZGoorC8B8ALMr+OunMQ6UKAJZDoopIpQoAlsCQX3QqVQCQPYkqBUIVAORNokqE4T8AyJhJVOkQqgAgSwpUqTH8BwD5kagSpFIVmQ2VAXgticqGygyQnAA4nTgVpLmhsuE/AMiDRJU4oQoAMiBRpc/wHwAkTZzKhUoVAKRLosqIUAUAiZKo8mL4DwBSZKn07AhVAJAWBapMGf4DgIRIVPlSqQKAJAyuWClRZUSoAoD4FKgWQKgCgJgUqBZDqIrMhsoAt0yB6jw2VGaA5ARwsySqs6W5obJQBQBzE6cWyZIKADAriWqphCoAmI9EtWCG/wBgDuLU4qlUAcDVSVS3QKgCgOuSqG6EUAUAV7STqLbbpBPV09PT58+fw8dv37799OnT09NT3CZlRKgCgGvZT1Qpe3p6+v33379//x4+/fHjx/Pz82+//SZXnUioAoDprVaZJaqiKP7zn//88ssvX79+7Y58+fLl119//euvvyK2KiMrK3pHtFrpf4CDwrrYOb5PZjqJ6u3btz9+/Ng/fnd3N3g8WbEury7qMY2vo+9XA9y4TENVponq+fn5p59+OvTo9+/fRx6d39GNaISqm6NSBTAix1CV3ZBfn0rVhSz+CQATyLRA1ff+/fvn5+cvX770Dz48PLx58yZWk/JiojoAXGR/TnqRYaIqiuLjx4///PPPhw8fuiMPDw9///33x48fI7YqI0LVKzRNU9d10zQTnglAvgbjVJFnoiqK4t27d9++fbu7uwuf3t3dvXnz5tu3bz///HPchuXCnJ6T1HW92Wz6R8qyHMxMp59ZmFMFMCrxOVVLilM7Eu/5o2JdXlWqjquqqstJZVmGD9q2rapq58x+oho/E4CsDY73ZRtCmIZQdURd123bFkWxXq+3223TNNvt9vHxsSiKtm37JaimaUKiKsuyO3O9Xu+fCUC+Blf1FKcoDP8dFUqg6/W6ruv+8aqq2rbtD+2FI8VevfTQ8cLwH8CoBAehljEh/agEe/5VDP+lqAtMO4kqPLSTtLqC1s6Z+18LwLg0t/W9kUTF2YSqMeNjdnVd70+WGjliBBDgFGlu62vIj6OEqjEhBnVTzkcWSugPAo5/NwDGJbitb9brpDMboeq4qqrqul6tVvf395vN5v7+frVauaEP4Eq+fPnST1TdwZ2VvucxOC0dBpkoPWZnv8ZQsgpzp4Ku97rFFAb7c39W++D3fxW/OGCRktrW92YnUSUyUT27q6S9/07V//U0TXN/f18URVVVFw7qRX/JAiTl7u7up59+OrStr0R1U86+RF6Sxi5h+O8k+6skdAtQdUdGvjycZsQQ4BTv37/vbz8XPDw8vH//frY2mJbOGYSqMWG8r5uo3ueePoArib6tr0lUnEeoOtNO2emUjKVSBXCKiNv6mpbOJYSqMSEG9Wemdw6tobAfqk5ZbQGAvnfv3v3xxx/h4x8/fvz555/zJKodEhWvIlSNedWiU2GUcCRUAZAsiYrLCVVjulC1cx9Bt3dyf1OasB1N27b9fWkGzwQgKaalMwnrVB3RrZ5QFEVZlmENhW5A8NDeyUfPDGyoDDBihtWSFKgGJbJO1dliXV5d1I/r56rO/kqeQZerjp5ZCFUAo659aZeoDhGqzvy5+XbZzJqm6U85H59uFVQvDp0pVAGMuOqlXaIaIVSd+XPz7bIFEKoARlzp0i5OHSVUncdEdQBuiETF9dj7L7KR/Yny/S8CQIIG32690WYq1u5+44SqyCQngBkoUC3M+NUzVuQSqgBYMgUqZiNUAbBYClTMyUR1AJZJomJmKlUALI04RRQqVQAsikRFLEIVAMshURGR4T8AFmInUYlTzEyoAiB7ClSkwPAfAHmTqEiEUAVAxiQq0mH4D4AsiVOkRqiKzIbKAGeQqG6cDZUZIDkBvJZEhQ2VAeAi4hQpM1EdgDxIVCROqAIgAxIV6TP8B0DqLJVOFlSqAEiaREUuVKoASNZugJKoSJlKFQApMomK7AhVACRHoiJHhv8ASItJVGRKqAIgFUPrYK/2Z1ZBmoSqyOz9BxAcSFREkP4FyN5/DEj/hQtwbYPXx+12+DgU9v4DgH3mpLMYQhUAcRwqUEGmhCoAIlCgYnmsUwXA3CQqFkmlCoD5iFMsmEoVADORqFg2oQqAOUhULJ7hPwCuS5ziRqhUAXBFEhW3Q6gC4FokKm6K4T8ArmInUZ0Rp2zkRV6EqshsqAws0uWJCkbYUJkBkhOwMIb8mEGaGyqbUwXAZCQqbplQBcA0JCpunOG/4+q6PvRQVVVVVe0fb5qmaZpDjwIsj0lUIFQd0TTNZrMZOWEnNtV13Z0fPijLsmmaa7UPIDYFKggM/x0xnodGElVZluGDtm3Vq4ClkqigI1SdpCzL7ZB+WupqWuHkpmm22+16vS6Kom1bxSpgeSQq6BOqjghh6JRSUzf1qp+f6roOJav7+/vpGwcQyWo1MIlKouLGCVVHtG1bnBaqwpmhNNU3Ms8dIEcKVDBIqJrYfvzqjhgBBHK3X6AqJCp4IVSdatVTVdVOQuo+HalpCVVA1gbjlEQFHaFqTBeDdmZEtW17f39vXA+4EQpUcAqhaky/tvT4+Bju+Ht8fAxHNptNd8J4FSrMVR88Z3WuSZ4gwFEKVMSS3SXS4p9jwjBfVVX9olRVVdvtNvzO7u/vL9wR2YbKQMoUqIjo7EukDZVTFELV4DBfV6/qzhz5PqffQgiQCIsmwGsJVWdyTx+wYApUcAahahqnZCyVKiALEhWcR6gaU9f1oVv8Dq2hsB+qTlltASARhvzgbELVmM1ms9lsBsPQfng6dIuf8UEgF/uJCjidUDXm0HbI3d7J/enqoabVtm2/uNWdub99DUA6BqelA6+yckv/uO62zLIsQ8mqaZpwN19Zljthq6qq7qFw52D4tDhwX+hqpf+B+EyiYmFiXV5d1I/rolLffqI6dPKhMwuhCkiARMXyCFVJa5qmP998fMp586J6cehMoQqIy5AfiyRU3SKhCohFgYoFi3V5tU1NZCNL6ctbwJVIVOQuzT1wharIJCdgZhIVCzB+9YwVuYQqgFshTsFVWacK4CZIVHBtKlUAR4ShhHwH6wdHQrJ9NpAuoQpgyRSoYDZCFcAyKVDBzIQqgAVSoIL5magOsDQSFUShUgWwHOIURKRSBbAQEhXEJVQBLIFEBdEZ/gPI3k6iEqcgCqEqMhsqA5dQoOI22VCZAZITcDaJipuV5obK5lQBZEmigtSoVAFkRpyCNKlUAeREooJkCVUA2ZCoIGWG/wAyIE5B+lSqAFInUUEWhCqApElUkAvDfwDpslQ6ZESlCiBREhXkRaUKIDmG/CBHKlUAaZGoIFMqVZHZUBnok6jgFDZUZoDkBHRMooITpbmhslAFEJ8CFSyAOVUAkUlUsAwqVQDRDI5RSFSQKaEKIA4FKlgYoQpgbgpUsEhCFcCsFKhgqUxUB5iPRAULplIFMAdxChZPpQrg6iQquAVCFcB1SVRwIwz/RWbvP1gwcQquxN5/DJCcYKkkKrieNPf+M/wHMD2JCm6QUAVw0NPT0+fPn8PHb9++/fTp09PT09Gv2klU261EBTdBqDpHXddN04yc0DTN0XOAxD09Pf3+++/fv38Pn/748eP5+fm3334bz1X7iQq4EStzel6rqqq2bYsDA7p1XW82m/6RsiwPpavVSv9Duj5//vz9+/evX7/2D3748OHu7u6PP/7YP9+QHyQi1uXVRf11mqa5v78PH+93XT9RlWUZsldxOFcJVZCyt2/f/vjxY//43d3d/nGJCtIR6/Jq+O91ukS1r2makKjKstxut03TbLfb9XpdFEXbtoYCIS/Pz8+DiWrwIYkKKFSqXqUb+At2uu7QsODIcKFKFaTsxEqVSVSQGpWq1NV13bZtWZah+LQvJKf9R+u6vnbbgGt4//79hw8fdg4+PDy8f/8+fLxaSVTA/xGqTtIN7R0dxauq6tARI4CQl48fP/7zzz/9XPXw8PD3339//PixMOQH7BGqThKmUj0+Ph46oQtM+6Fq/xwgC+/evfv27dvd3V349O7u7s2bN9++ffv5558lKmCfbWqOCzmpLMuRwAQs0rt37/74448///yzKIowj0qcAg5RqToiTKUqjtWZxh8ty/LQOatzXfCcgDNJVDCn7C6RKlVHhKlUIwN/F3L3H2RiK1HBzM6+RNpQOUXht7Jer48O/I2fEGpdRg8hW7vv7LbzA/apVB3UH63rL4vQHQ8Hq6qSlmCpBv+7K04Bg4Sq43b28ts/3s9VTdMcyliyF+TFeB/wKkLVmDDBfEd/R79iLyrth6pTVlsAUiNRAa8lVB1UVdXg/Xrdrsk7j4YdlPe/xPJUkJehIb9wSKoCxpioPpkwxapt250JWCGBHdrcBkiKAhVwNhv6vlpXqdrvum7v5LBSaNM03VjhYD/bUBmScihRhRuB/bVCLmJdXg3/TSlMqGpfhINlWRoBhPTZGhm4kErJ9JoX1YtDZ6pUQSLGE5VKFeQl1uXVRT0moQqiO2USlVAFeYl1eTVRHbhdpqUDEzKnKrKR/Yn8txiuSqKCfEXcNXmEUBWZ5ARRmJYOWRu/esaKXEIVcFsUqIArMacKuCESFXA9QhVwKyQq4KoM/wHLJ04BM1CpAhZOogLmoVIFLNbgDUASFXAlQhWwTApUwMyEKmBpFKiAKIQqYFEUqIBYTFQHlkOiAiJSqYrM3n8wCXEKboq9/xggOcHlJCq4NWnu/Wf4D8ibRAUkQqUKyNhOohKngIiEKiBLClRAagz/AfmRqIAECVVAZiQqIE2G/4BsiFNAylSqgDxIVEDihCogAxIVkD7Df5CTsKLdTa0Zm0KcuqkOB86mUgWkK4VEBXAioQpIlEQF5MXwX2Q2VIZBlkoHRthQmQGSE+yTqIBxaW6oLFQBCTHkB+TLnCogFRIVkDWhCkiCRAXkzvAfEJ9JVMACCFVATApUwGIY/gOikaiAJVGpAiIYvN9ZogKyJlQBc1OgAhZJqALmo0AFLJhQBcxEgQpYNhPVgTlIVMDiqVRFZkNlFk+cAiZnQ2UGSE4sm0QFXEOaGyob/gOuRaICbopK1amapgn/VlVVVdXRk088ExZJnAJu0Mrw01F1XW82m52D6/W6rutTTi7LMgSyfauV/ud1Qk078ZeNRAXEFevyavjviH5IKsuyLMvw8Waz2Q9VOyeHD9q2Va/idkhUwM1SKTkiFAZ26lJVVbVtWxTF4+NjF5iaprm/vy/+XZrqYlb/zP431/+8SuKVqp1ElWozgYWLdXl1UR/TRaL9XtoPW13S2jn50PFCqOL1Ug5VEhWQCMN/KQoFp24g79AJQUhO6/V655zBqVewJKuVRAXg7r9R1YuRE04/Eu4HnKxxkAaTqAACoWrMoSJTV6DqB6adI4NfJVSxMBIVQEeoeoWQsbr7+9brtZDELTPkB9AnVJ2qaZr+AlQ79wMeWokqKMuybdvBc85eSj/NqcrcCAUqYAZpbvA3wkT1U1VVtV6vu6WqNpvNJL/s7bku/9FwHokKmEd2l0ih6hXqug77zzw+PoYj3fDf+DhguDHQWCELIFEBHCJUnSNUrYqXtAQzeHp6+vz5c/j47du3nz59enp6mrMBg+smSFQAHaFqTKhLDT60U3bavw3w6JfA6Z6enn7//ffv37+HT3/8+PH8/Pzbb7/NlqsUqACOEqrG3N/f39/fj09C37F/8imrLcC4//znP7/88svXr1+7I1++fPn111//+uuva//o/QJVIVEBDBGqjhsMVfsHwwT2kVAFZ/vy5Us/UXUHv3z5ctWfOxinJCqAQULVmDBxarPZ7ASjbnmF/qY0YYWFtm13llrYPxNe5fn5+cePH6996EIKVACvZUPfI7p1E8qyDON3TdN089MP7Z0cTh45s/vm+p9TvH37djA83d3dXSNUiVNA1mJdXl3Uj+uiUl9ZloPjevsnHzqzEKo42adPn56fn3cG+x4eHt68efPnn39O+7MkKiB3QlXS+rcBjm+x3J0cdvobP1mo4kRPT0+//fbbr7/+2uWqh4eHv//++9u3bz///PNUP0WcApZBqLpFQhWne3p6+uuvv0Jd6u7u7v379x8/fpSoAPYJVbdIqOK1wiS/yV82EhWwJLEurzZUjmxkA0F5i3nsr5MOkLg091oWqiKTnIhLogJyNH71jBW5hCq4UYb8AKZl8U+4RRIVwOSEKrg5EhXANRj+g9tiEhXAlQhVcCsUqACuyvAf3ASJCuDahCpYPokKYAaG/2DJxCmA2ahUwWJJVABzEqpgmSQqgJkZ/oOlEacAohCqIrOhMtOSqIBbYENlBkhOTEiiAm6EDZWBK7JUOkBcQhVkT4EKIAXu/oO8SVQAiVCpglwNzhmQqABiEaogSwpUAKkx/AfZ2UpUAAlSqYK87KYncQogESpVkA0FKoCUqVRBBsQpgPSpVEHqJCqALKhURWbvP8YNvUBW+zOrAG6Kvf8YIDlxyGCBKs33EYCZpbn3n+E/SJEhP4DsCFWQHIkKIEeG/yAtO4lqJ04ZLwZIlkoVJGQ8UQGQMpUqSIIhP4DcqVRBfBIVwAIIVRCZRAWwDIb/ICaTqAAWQ6iCOBSoABbG8B9EIFEBLI9KFcxqcO8EiQpgAYSqyGyofFMUqAAmkeZGqEJVZJLTjVCgAphQmhsqC1VwdQpUALfARHW4LokK4EaoVJ2qaZrwb1VVVVUdPfnEM1kwcQrgpqzM6TmqaZr7+/udg2VZhpi1o67rzWZzyplFUaxW+n+xJCqAWGJdXg3/HdFPVGVZlmUZPm7bdn8eXD9R9c9Ur7o1EhXADRKqjugS1Xa7DYN62+328fExHOynpaZpQqIqyzKcvN1u1+t1URRt2x4qVrEwq9XAzjMSFcAtMPw0pitT7fdSV5TqHqqqqm3b/ZMPHS8M/y2OAhVACgz/pWikvNTVqLpzQnIKpam+uq6nbxnpkagAbpxQNSYEpm52VN+haVL7x/fjF8tjyA8ASyqMGYlBXf0pZKbuzJE56WGRhYmaRkL2ExUAN0ioOkd/TnrsthCTIT8AOob/Xq2u6+6WwK5ANT60F7LX4Dmrc036nDiHRAVwVdldIlWqXmFnFdBJ7ixw91+mJCqAazv7EhkrV6lUnaqqqi5Rrdfr/XUTRr423BhoQtVimJYOwD6VquP6Bar1em2JhFumQAXAIULVEV2iGtnCr/j3ugmnr7ZAXiQqAEYY/jvilETVt3/aKastkD6JCoBxQtWYbqTvlER16BY/a37mznZ+AJzC3nNjwu0DZVmOzKPqD/yFslZ/3tXgwf731/+JU6ACyE6sy6uL+phT7snsd2C3d3JZllVVNU0TPi0O3BcqVKVs8Jfv1wWQPqEqRa8NVUUvV3VG5mMJVclSoALIl1C1HM2L6sWhM4WqBClQAeROqLpFQlVqFKgAFiDW5dU6VZGNjDDKWzOTqABykeYeuEJVZJJTCsQpgLyMXz3t/QdxSFQATEKo4qZJVABMxfAft2t/nXQAOJtQxS1SoAJgcob/uDkSFQDXIFRxWyQqAK7E8B+3QpwC4KpUqrgJEhUA1yZUsXwSFQAzMPzHkolTAMxGpYrFkqgAmJNKVWQ2VL4SiQpgwWyozADJ6RoslQ6wbDZUhjlIVABEoVLFchjyAyAilSoWQqICIC6hiiWQqACIzvAf2TOJCoAUCFVkTIEKgHQY/iNXEhUASVGpIj+D649IVADEJVSRGQUqANIkVJENBSoAUiZUkQcFKgASJ1RFZkPlU0hUAPTZUJkBktM4cQqAfTZUhteRqADIiFBFoiQqAPJi+I/kiFMA5EilirRIVABkSqgiIRIVAPky/EcqdhKVOAVAXlSqSIJEBUDuVKqIzJAfAMugUkVMEhUAiyFUEY1EBcCSGP6L7Gb3/jOJCoCz2fuPActOToMUqAC4kL3/QKICYLGEKuYjUQGwYELV6zRNU9f1iac1TXP1BmVitRqYRCVRAbAkqxuc03OJMEw70ml1XW82m/6RsiwPpavV6ib6X4EKgDnFuryqVL3C0cpTP1GVZRk+aNu2qqprtitd+wWqQqICYKGEqlPVdX1/fz9yQtM0IVGVZbndbpum2W636/W6KIq2bW9wKHAwTklUACzVTQw/XaJpmv0sNdhpVVW1bbv/6KHjRc7Df+PDoIO3sub5RAHIj+G/7IXkFEpTfadMbF8SBSoAbpNQdURVVdsXj4+Pp5x/6MgtjACaQQXAzRKqptEFppE56csOVRZNAODGCVVMQIEKAISqaYxXocLyCoPnrM51nedxDokKgGvI7hJpQ+XIMr37r7M/5AcAkzj7EmlD5byNL+8Zbgxc3hKgEhUAdFSqOM9ugJKoALhxKlXTOGXdhMVUqkyiAoB9QtXE9kPVKastZESiAoBBQtVkDt3it6TlqaxEBQCHCFWTCdvRtG3b35em22V5f/uavOyv7VkUCS3rAADR5bqhbxTd5sqHOq3bO7ksy6qqmqYJnx76klw2VB66NXVsQ2UAiMiGykvQNE0YBGzbdrPZdAEr6/Cxn6g+ffocPnj79u2nT5+enp7mbhMApCePSklemhfVi0NnJl6p2o9T//u/T7///vsvv/zy9evXcOTDhw///PPPt2/f3r17N3f7AGBIrMtr0hf1xUs5VA3e5ff58+fv3793iSr48OHD3d3dH3/8MV/jAOAwoeoWpRmqBhf3D818+/btjx8/9h+9u7sbPA4A8xOqbtH45kRxXhCHl6F6fn7+6aefDn3h9+/fRx4FgAkd3d1PqLo5SVWqRgpUHZUqANIX6/Jq7z+K4uR10t+/f//8/Pzly5f+wYeHhzdv3lytaQCQB0sq8IqdZz5+/PjPP/98+PChO/Lw8PD3339//Pjxaq0DgDwIVTdtf5308Z1n3r179+3bt7u7u/Dp3d3dmzdvvn379vPPP1+zmQCQgYTm9NyguHOqLtkaOcwQ9OIBIEFWVGdWlyQqAGCfieq3aH/IDwC4kFB1WxSoAOBKDP/dEIkKAK5HqLoVEhUAXJXhv+UTpwBgBkJVZCO7F01yO6hEBcDyHN37LwqhKrKrLqQhUQGwSONXz1iRS6haJnEKAGZmovoCSVQAMD+hamkkKgCIwvDfolgqHQBiEaqWY85EZStlANghVC2BIT8AiM6cquxJVACQAqEqbxIVACTC8F/GTEsHgHQIVVlSoAKA1Bj+y08iiSrNfZdOpPFR5NvyQuPj0fhYsm58LF7oHU0AAAhaSURBVCpVkb1qQ+XBc9WoALg1aWY+oSqy0xd8SqRABQDR2VCZMylQAUD6hKrUKVABQBZMVE+aRAUAuVCpSpQ4BQB5UalK0b8T1ao4IVGdOCnvlNMmnN83f6s0fvKfqPGvMn+r8m38tFOJNX7ybzXVj5v2tDRv+usIVclRowKAHAlVCVmt7DwDALkSqlKhQAUAWROqkiBRAUDu3P0XnyE/AFgAlaq0SFQAkKnV6XvPMbnEbw0FgExFiTdCFQDABAz/AQBMQKgCAJiAUAUAMAGhCgBgAkIVAMAELP4ZTdM0TdNUVVVVVey2/Mt5DQtfdejRuq4vbtf0QpuTatsZTcqi50MLE3zBn92wLLq90PMJqOs6tc7vvKpt6ff8SBvm+xVsmV1Zlju/hcfHx9iN2m4va9j+16b/MkuwbWc0KfGef3x83G9VWZax23VpwxLv9u12u16v9xu2Xq9jt+vShqXf831da2M3ZMBr25Z4zw/+Rc//yo/fEbem/2vuv0aj56oLG5b439u+7i8wdkP+z3lNSrnn+29zZVnuNDXrhqXc7dt/B5edJxg3V13esMR7vq//Movdll1ntC3xnh8M653ZrrDxO+KmdC/K/i84hRfl5Q2b+YV7of6fX+y2/NfZTUq55wefUfduHrFedXnDUu727UvzdmLK4J/5zC5vWPSncLqkMseOM9qWeM+Ht9DoVfDkftPLNviGsk3gxXp5w9J849gxWB/OvUkpPItBI4W3LkHO36rtRA1Lttu3o8/i0F/6PCZpWMo935dOXXbfeW1L8In0hScVfYDb3X/z6ab4HZpMF2uiX7INI2sjc1q7GaMj51xPsg2bSmj8yGBNrGeXbMMmV9d127ZlWY6PSUWRctsu0bZt0fsTjkWoms/I+0V4cYfXxPwub1gub4VVVXX/nxif1TibC5uUcs+PXEHjvvFd3rCUu70oiqqq1uv1+J1Q87Xm3z/3woYl3vNB0zSbzaZIsrVnty3B55ImoWo+i7/GlGVZ1/WqJ/p/GhYv5Z5vmma73Q6+F3eX1SjtvLxhKXd7URR1XYdb5XeOd085Vjsvb1jiPR/c398XRZHIf9t2nN22LHo+2Gne3HFw/hHHmxU6fHAaXdw70S5vWCK3XbzKMu7+y7rno08p3XF6w/Lq9vV63W9w9EknnTMaln7Ph/+ddi+huNMHd1zStsR7PpH1FLYmqs9p5LebQqi6pGFdlav/TVK+nXgbu88HndGk7Ho+wfsug1c1LKNu37nYpJOozmtY4j2/H1PSCVUXti2Xni96Ca/fPEsqLNDIe0fuoWq9Xodpj4e+PJ238s4yQlVGPb9zEY3dnP9zRsMy6vbtS2v74/uxW/RfZzQs8Z7fv36nE6oubFviPf/4+DjYvO3sixbF/03fjpEbPuNe4K/asJ2CczqWEapGJNXz6Sw+uWPyhiXV7TsWMPA6InrPD76KEglVV21b9J4fN/Nbvb3/5pbsPRRXalhVVW3bxrqx8ZYl0vNN04S5sUVRjN/5NbMrNSyRbh8Ubr7bbDapNW+ShsXt+f77Z/+1tLNgTZRNAK/dtpRf88W/10mZofOFqvmEV97gQ3GTVrINI3ddcCnLMqnXUrINm0R4RoPXj6qqwh31USTbsAkdehbd8Yi3y6XctsUQquY2GF+Orok3g7Mb1m1dPvJ//bhPbanS7/lkg8slDcul2x8fH1O7Rl7YsMR7fvBHd++r4dFYv5EL25Z4zxcvDRts3twricwzykgQ+nz/NoRwPOJckwsbdujLT/8O81vGnKqUez6R2ST7Lm9Yyt0+3oa4v5TLG5Z4z+9L9q9g+8q2Jd7zoQ2D87pm/hWk+JtesC7LHz04s1c1LNyw0/8TGjkz+lM7JMdQlVfPd29zj4dl0bC8un3bu4rsPJHoN2q9tmHZ9fy+TENVdj1/9KVlSYXFKl6s1+twk2oKMf9VDeuuSd2R/k3p4U8x2bu9OjmGqrx6vjhBFg3Lq9uDnebttDCXhuXY8zsyDVU59vzIS2vOOxNT/E0v286iOOm8Ik9v2ODLdPDLE3lqg5YRqrYJ9/xgq3Zk0bC8ur0zOMclhfveT29Ypj3ft5hQtc2h51N4za+2p73FMK1u3l8xOvVvfhc2rHlRvZi0dRyk56NIvNv7f85JNe/yhiXe8wuWeM9Hf80LVQAAE/if2A0AAFgCoQoAYAJCFQDABIQqAIAJCFUAABMQqgAAJiBUAQBMQKgCAJiAUAUAMAGhCgBgAkIVAMAEhCoAgAkIVQAAExCqAAAmIFQBAExAqAIAmIBQBQAwAaEKAGACQhUAwASEKgCACQhVAAATEKoAACYgVAEATECoAgCYgFAFADABoQoAYAJCFQDABIQqAIAJCFUAABMQqgAAJiBUAQDZqKqqqqqmaXaON00THorQphdCFQCQjaqq2ra9v7/fOX5/f9+2bdxQtdputxF/PADAq4RcVZZlV6+q63qz2RRFETfVqFQBADkJWapt27quw6chUT0+PkZtl0oVAJCbpmnCCOB2u12tVkVRrNfrkLEiEqoAgPyEQcDwcX8oMCKhCgDIUqhRFbGnUnXMqQIAmIBQBQDkp796QtyVFDpCFQCQmaZpwoSqcMdfdydgXOZUAQCZCbOpHh8fq6pKZJGqQqUKAMhLGOwryzJ8UNd1WZZFAoOAQhUAkI26rsPAX38NhZ3lQGMRqgCAbISRvvV6vXM8HAmPxmJOFQDABFSqAAAmIFQBAExAqAIAmIBQBQAwAaEKAGACQhUAwASEKgCACQhVAAATEKoAACYgVAEATECoAgCYgFAFADABoQoAYAJCFQDABIQqAIAJCFUAABMQqgAAJiBUAQBMQKgCAJjA/wc3i6xUmzf39AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "my1 = R.TF1(\"myFunc1\", \"pol1\", bin0, bin1)\n",
    "my1.SetLineColor(4)\n",
    "\n",
    "frp1 = hd.Fit(\"myFunc1\", \"S\")\n",
    "frp1.Print(\"V\")\n",
    "myc.Update()\n",
    "myc.Draw()\n",
    "myc.SaveAs(\"test_fitting_1.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "incorporated-female",
   "metadata": {},
   "source": [
    "Do a linear regression BY HAND"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "upset-columbia",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-15T22:48:38.972976Z",
     "iopub.status.busy": "2021-04-15T22:48:38.972628Z",
     "iopub.status.idle": "2021-04-15T22:48:39.080515Z",
     "shell.execute_reply": "2021-04-15T22:48:39.080888Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.6676596448093777 1.6676596448093777 3.188503842606061 3.1885038426060603 -0.7389364642331806 -0.7389364642331806\n",
      "\n",
      "Input [x]: [0.5, 1.5, 2.5, 3.5, 4.5]\n",
      "Input [y]: [8.0, 29.0, 41.0, 60.0, 62.0]\n",
      "Input [y_err]: [2.8284271247461903, 5.385164807134504, 6.4031242374328485, 7.745966692414834, 7.874007874011811]\n",
      "Estimated m (p1): 15.117729551767999 +/- 1.6676596448093777\n",
      "Estimated b (p0): 1.7180860595587057 +/- 3.1885038426060603\n",
      "Estimated correlation (m, b): -0.7389364642331806\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "l_x, l_y, l_ye=[], [], []\n",
    "for i in range(1, nbins+1):\n",
    "  l_x.append(hd.GetBinCenter(i))\n",
    "  l_y.append(hd.GetBinContent(i))\n",
    "  l_ye.append(hd.GetBinError(i))\n",
    "\n",
    "[est_m, est_err_m, est_b, est_err_b, est_corr_mb] = lls(l_x, l_y, l_ye)\n",
    "print(\"\\nInput [x]:\", l_x)\n",
    "print(\"Input [y]:\", l_y)\n",
    "print(\"Input [y_err]:\", l_ye)\n",
    "print(\"Estimated m (p1): {0} +/- {1}\".format(est_m, est_err_m))\n",
    "print(\"Estimated b (p0): {0} +/- {1}\".format(est_b, est_err_b))\n",
    "print(\"Estimated correlation (m, b): {0}\\n\\n\".format(est_corr_mb))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lesser-world",
   "metadata": {},
   "source": [
    "Now do a likelihood fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "proper-pottery",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-15T22:48:39.115356Z",
     "iopub.status.busy": "2021-04-15T22:48:39.114946Z",
     "iopub.status.idle": "2021-04-15T22:48:39.255734Z",
     "shell.execute_reply": "2021-04-15T22:48:39.256071Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " FCN=1.21602 FROM MIGRAD    STATUS=CONVERGED     157 CALLS         158 TOTAL\n",
      "                     EDM=4.672e-09    STRATEGY= 1      ERROR MATRIX ACCURATE \n",
      "  EXT PARAMETER                                   STEP         FIRST   \n",
      "  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE \n",
      "   1  p0           1.93421e+00   3.65636e+00   2.06054e-03  -2.05085e-06\n",
      "   2  p1           1.52263e+01   1.81514e+00   1.02238e-03  -4.08006e-05\n",
      "                               ERR DEF= 0.5\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Info in <TCanvas::Print>: png file test_fitting_2.png has been created\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAI8CAIAAAD0vjrdAAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dTXLbxto2YPDUN3XJnjup7ATQZhL/rMCeEBwlG8hxypshsJFUql7Nba2A3wARDkwS4B+I7gaua+CSQEhqtiji9tON7tVut8sAALjNf0I3AABgDoQqAIARCFUAACMQqgAARiBUAQCMQKgCABiBUAUAMAKhCgBgBEIVAMAIhCoAgBEIVQAAIxCqAABGIFQBAIxAqAIAGIFQBQAwAqEKAGAEQhUAwAiEKgCAEQhVAAAjEKoAAEYgVAEAjECoAgAYgVAFADACoQoAYARCFQDACIQqAIAR/L/QDQAAFmC1+veD3a7nQPJWu/s8laqqiqK4x3eek1X7ggKAZenGj/GvhneKN8PuFaqKoqjrOs/zsiylqz6r1en+P+eccU+b/bea/idq/OjfavqfGOe3mv4nxvmtpv+JGn/Bt1qtVvsBav9Lmu8xfc+P7r5zquq6fnx8XK1WRVFUVXXXnwUAxOVgQKYvUc3DvUJVVVXb7TbP8+bTNl2VZXmnnwgARGRhiSq7a6WqqU7tdrtuutpsNk3hSroCgNnaS1S73V6i2u3mlqiyaZZUaNPVer1ujtR13aYrw4IAMCv7NardQcSasDETmnSdqrIsd7vdbrczLAgAM7Ra7SWq2Q/5dU0aqqqqKstytVrVdb33UFO4Eq0AIFUHcWpRiSq735IKXVVVVVW12Wy6B/M8b2dWlWXZPrrdbpezBEOoez5HofGhpNv4dFueaXw4Gh/KxY2PqUAVqufv+FObutReUaqbpfabslo1JyxnltWy/t5iovFBpNvyTOPD0fhQLmt8TIkqC9fz99qmpln8s/10IEu11uv1XjULAIhaZHEqrPvu/XdOlmo1dwKaVgUACVjeMlQn2fsvpOG9/yIvGi+orB2ZdBufbsszjQ9H40MZavyxi9fEierk5rlzm1PFSbP9e4uexgeRbsszjQ9H40PpbfwZiSr4k57bnCoAYG4M+Q1KOERPoFkMYvicoigOBzqbLzz6UFfS/4kBYEEOC1S7/XXSs2gSlUpVjB4fH885rZucumtuNR8sapEIAGYoqUQV0KQrqien3U7nTN1E1d2Kx5x9AFJ1bMjvcC8/iSoz/He1qqqaOlbbge2RbmmqjVlHV4o3/AdA1A7SUxIFqhmuqD5vzc2c3ajUrne616V9xzOhCoCYpZmosnCXV8N/12iC1Hq97hafmuS0Xq/3TracKQCJWa2O7I5syO8UlZKLHQ78NQ5rVycfUqkCIDrpr5ugUpWMJlFtt9vuwXYS1cCcdPcAAhC79BNVQJZUuEwzltfsaRi4KQAwrlOTqMSpYSpVl2lu5TucJjVchWqWVzh6zupaIz0hADh3EtXEkrtEqlRd4B5lKnOqAAgs1iG/qy+RoXKVStUF+spU2eBUquzlxkAjhgBE57BAFUeiSpFK1bnOmYoOAMmItUCVLpWqczUFqsNlqBpt0hqYXCWNARALBao7EKrO1QzhnVzJ8zBUKXEBEJeDRHV4Rpa5I+piQtVZzlkVve8WP8tTARCLY3f57Z2iQHU1oeosTTBqYlOfJnjVdd1NYFVVNdPb+8YNAWAip4b87DxzI9uknKW5OXO9Xg+XrNq9k5tlF6qqaj7Neu4LtU0NABO5pEDVXPXSvUKFury6qJ/Wt9nfUW2uauV53jcCKFQBMIULh/yEqit/brpdFq3qRfGi70yhCoD7umoGlVB15c9Nt8tmQKgC4I6unZMuVF3H4p8AMDsW9gxBqApsYH+idP+LAEBIpxLVDC4vAXdNHiBUBSY5ATCmBSSq7NTVM1TkEqoAiFTqM3umZsgvNIt/AkD6JKoICFUAkDiJKg6G/wAgZXuJarc7ODBlaxZNqAKANB0tUElU4Rj+A4AEGfKLj0oVAKTm1FLpmUQVglAFAOlQoIqYUAUAiVCgiptQBQApuHZ3ZCYjVAFA3Az5JUKoCsyGygAMUaA6xobKHCE5AdBLouoR54bK1qkCgChJVKlRqQKAyIhTaVKpAoCYSFTJUqkCgDi4yy9xQhUAROBUohKn4idUAUBQx25Vk6hSJFQBQDiG/GbERHUACESimheVKgCY3OGQ32537Ng0rWEcQhUATOuMRCVOpUioCszefwDLcnTIT6K6kL3/OEJyAliKM2ZQZRLVeeLc+0+oAoD7O7VOekOiSppQBQB3dkaiEqdmQKgCgHuyl99iCFUAcDedRCVOzZ5QBQB3oEC1PFZUB4CxSVSLpFIFAKMaTFTi1IypVAHAeBJPVE9PT58/f24+fvPmzadPn56ensI2KSFCFQCMYbXam5aeYqL67bffvn371nz6/fv35+fnX3/9Va46k1AFADebxSSq//73vz///PNff/3VHvny5csvv/zy559/BmxVQla2SQlotdL/AL2azUYSeJ+cRaLKsuzNmzffv38/PP7w8HD0eLRCXV5d1EMa3pzIrwZYuARC1ant/GJu+57n5+fXr1/3Pfrt27eBR6d3cnc/oWpxVKoABsQequZSoGqpVN3IkgoAcKFTBaoswUSVZdm7d++en5+/fPnSPfjhw4dXr16FalJaTFS/QFVVZVlWVTXimQAkZqaJKsuyjx8//vPPP+/fv2+PfPjw4e+///748WPAViVEqDpLWZar1erx8XGz2Tw+Pq5Wq6IobjwTgMT8uGhC1rNuQqKJKsuyt2/ffv369eHhofn04eHh1atXX79+/emnn8I2LBXm9JxWFEVd183HeZ53P96rRZVludlszjmzYU4VwIC45lTNt0B1KK6ev1yoy6tK1QllWTbZaL1e73a7qqp2u912u82yrK7rblSqqqpJVHmet2eu1+vDMwFIzKkCVTajRMXVVEpOaNL6er0uy7J7vClfdUtQbUFrr0v7jmcqVQCDoqiXLKlA1Yqi52+gUhWjNjDtJarmob2k1Ra09s48/FoAhsWyA50CFZcQqoYMj9mVZXk4CX3giBFAgHPEsgPdGctQSVR0CVVDmhiU53n7ad9CCd1BwOHvBsCwKHagm93CnkzA4p+nFUXRva2vnY0uJAHcw5cvXw7X7/7y5cvDw8Pvv/8+RQs6iUqc4nwmSg/Z21qoKVm1CyVknUl8beo62p+Hs9qPfv+L+MUBsxR4BzoFqizLopmontxVUqXqXN1fT1VVj4+PWZYVRXFjvSr4SxYgKg8PD69fv+7bgU6iWpSrL5G3pLFbmFN1lsNVEtoFqNojA1/enGZpdYBzvHv3rrtTSuPDhw/v3r27408dTFTmpHMOoWpIM97XTlTvck8fwJ0E2IHuVKKCcwhVV9orO52TsVSqAM4x6Q50P27nd3QvPziTUDWkiUHdmemtvjUUDkPVOastAND19u3b9ka/79+///HHH/dKVN3PTKLiNkLVkIsWnWpGCQdCFQARkagYm1A1pA1Ve/cRtHsndzelabajqeu6uy/N0TMBCOnHIb/MtHRGYp2qE9rVE7Isy/O8WUOhHRDs2zv55JkNGyoDDLjLakkKVGeIZJ2qq4W6vLqon9bNVa2+FdXbXHXyzEyoAhg08qX9YO0iiaqPUHXlz023yyZWVVV3yvnwdKtG8aLvTKEKYMCYl3aJ6hJC1ZU/N90umwGhCmDAOJd2cepyQtV1TFQHYL4kKiZk77/ABvYnSve/CABRODUnPZOokhVqd79hQlVgkhPA+BSo5m746hkqcglVAMyLAhWBCFUAzIhlqAjHRHUA5kKiIiiVKgDSJ04RAZUqABInUREHoQqAlElURMPwHwBpOrVugjjFxIQqABKkQEV8DP8BkBqJiiipVAGQDkulEzGhKjB7/wGcS4GKF/b+4wjJCeAsEhUdce79Z04VANGTqEiBShUAceskKnGKmAlVAERqLy9JVETO8B8AUTLkR2pUqgCIz2CiEqeIk1AFQGRMoiJNQhUA0TDkR8qEKgDiMFigyiQqoidUARCaAhWz4O4/AIKSqJgLlSoAAjkVp7KsOUGqIg1CVWA2VAYW6owCVZR75i5C/BcgGypzRPwvXIDxGfLjNnFuqCxUATChg6udhT2ZDaEKgKkoUDFr7v4DYBISFXOnUgXAnZ0a8sskKmZBqALgnhSoWAzDfwDcjUTFkghVANyHRMXCGP4D4A4Gd0cWp5gloQqAUSlQsVSG/wAYj0TFgqlUBWbvP2A+BhPVFW9p3gbpY+8/jvCWAcyESVRMyN5/AMyRIT/IskyoAuAmgwWqTKJiSYSq08qy7HuoKIqiKA6PV1VVVVXfowAzYcgPOoSqE6qq2mw2AyfsxaayLNvzmw/yPK+q6l7tAwhCgQoOWFLhhOE8NJCo8jxvPqjrWr0KmJVTBSqJimUSqs6S5/numG5aamtazclVVe12u/V6nWVZXdeKVcAcrFaG/KCPUHVCE4bOKTW1U6+6+aksy6Zk9fj4OH7jAKb0Y5w6XIZKomLhhKoT6rrOzgtVzZlNaaprYJ47QBoUqOAMQtXIDuNXe8QIIJAky1DBeYSqc606iqLYS0jtpwM1LaEKSI8hPzibUDWkjUF7M6Lqun58fDSuB8yZIT+4kFA1pFtb2m63zR1/2+22ObLZbNoThqtQzVz1o+esrjXKEwQ4zpAfEUjuEmnxzyHNMF9RFN2iVFEUu92u+Z09Pj7euCOyDZWB6AwmKm9aTObqS2SoXKVSNaQJVUeH+dp6VXvmwPc5/xZCgMAkKriWUHUl9/QBM3RqWjowQKgaxzkZS6UKiJdp6XAzoWpIWZZ9t/j1raFwGKrOWW0BICTT0mEMJqoPafbyq6pqICq18jw/usef8UEgaiZRwUhUqob0bYfc7p3cna7e1LTquu4Wt9ozD7evAQjPJCoYz8ot/cPa2zLzPG/G76qqau7my/N8L2wVRdE+1Nw52Hya9dwXulrpfyAQQ37MV6jLq4v6aW1U6jpMVH0n952ZCVVAKINz0jOJisQJVVHrTqsqimJ4ynn1onjRd6ZQBQTgLj/mTqhaIqEKmJQCFcsQ6vLq7j+AZVCggjsTqgIb2J9IEQsYhznpzE7AXZMHCFWBSU7AfSlQMUfDV89QkUuoApgpBSqYlsU/AU5YrVZxjjUMkahgcipVALNjyA9CEKoAZkSBCsIx/AcwFxIVBKVSBTALg4lKnIIJqFQBpE+iggioVAGkzJAfREOlCiBZEhXERKgCSJNEBZEx/BeYvf+Aix28b5hExdLEuR6vUBWY5ARcRoEK7P0HwK0Gl0rPJCoISqgCSMGpOJVJVBCaUAUQPQUqSIFQBRA3uyNDIoQqgIi9JCpxCuJnnSqAWElUkBSVKoD4GPKDBAlVAJHpT1TiFMTM8B9ATCQqSJZKFUA0TKKClAlVABEwiQrSJ1QFZkNlQKKCS9lQmSMkJ1i0wd2RvT1Anzg3VDZRHSAQiQrmRaUKYHKDcSqTqCBNKlUA05KoYKaEKoAJSVQwX4b/ACZxMHPWJCqYGaEK4P4UqGABhCqAOxtchiqTqGAuhCqAuzkVpzKJCmZEqAK4DwUqWBihCuAO7DwDyyNUBWbvP5ihl79rcQruxN5/HCE5wdxIVHB/ce79J1QBjMSQHyybFdUBej09PX3+/Ln5+M2bN58+fXp6ejp+an+i2u0kKlgEoeoaZVlWVTVwQlVVJ88BIvf09PTbb799+/at+fT79+/Pz8+//vrrkVw1mKiAhViZ03Opoijqus56BnTLstxsNt0jeZ73pavVSv9DvD5//vzt27e//vqre/D9+/cPDw+///77/w6ZRAWRCXV5dVG/TFVVj4+PzceHXddNVHmeN9kr689VQhXE7M2bN9+/fz88/vDw8O9xk6ggSkJVGro3FOx1XZu3uhGqjVnb7bYoisPvpv8hTs/Pz69fv+57dO/vVqKCqAhVCWgH/hp7Xdc3LDgwXChUQcyOVqoO/2JNooLYhLq8mqh+rrIs67rO83y9Xh89oUlOh4+WZXnvtgH38O7du/fv33ePHBaoJCqgJVSdpaqqZhTv5A19h2N87RE3A0JaPn78+M8//7S5ypAfMEyoOkszWWq73fad0Aamw1B1eA6QhLdv3379+vXh4WEnUQFnsKL6aU1OyvN8IDABs/T27dvf//ije0ScAvqoVJ3QTKXKTtWZhh/N87zvnNW1bnhOwNl+/FuTqGBKyV0iVapOaBdEuNP3d/cfxGtwGapMooI7u/oSGSpXqVQNaX4r6/X65MDf8AlNrcvoIaTk1MKeEhWwR6WqV3e0rrssQndhzyzLiqKQlmBWFKiAqwhVp+3t5Xd4vJurqqrqy1iyFyTAzjPAtYSqIc0E8z3dHf2yg6h0GKrOWW0BiIJEBdxAqOpVFMXR+/Xa7fz2Hm12UD78EstTQRpeEtWxIb/mIakKGGKi+miaKVZ1Xe9NwGoSWN/mNkAU+hOVAhVwJhv6XqytVB12Xbt3crNSaFVV7Vjh0X62oTKEd8aQX3MjsL9WSEWoy6vhvzE1E6rqF83BPM+NAEKk+hOVBAVcSqVkfNWL4kXfmSpVENLZiUqlCtIS6vLqoh6SUAXBXDKJSqiCtBj+W6iBpfS9g8NdWDcB0hfnHrhCVWCSE0xKooJZGL56hopcQhWwGKalA/ckVAELoEAF3J9QBczd0FLpEhUwGqEKmK/BOJVJVMCohCpgphSogGkJVcDsmEEFhCBUAfNia2QgkP+EbgDAeCQqIByVKmAWDPkBoQlVQPqs6glEwPAfkDiJCoiDSlVgNlSGm5hEBYtkQ2WOkJzgSiZRwYLFuaGy4T8gQRIVEB+VKiA1JlEBURKqgHQoUAERM/wHJEKiAuImVAEpkKiA6Bn+A+ImTgGJUKkCIiZRAelQqQJi1b+wZyZRAfERqiAlzYp2i1gzNqal0hfR4cDNhCogMgpUQJqEqsDs/Qc/iKlABUTL3n8cITnBv8xJB84W595/QhUQAQUqIH2WVABCk6iAWVCpAsIx5AfMiFAFBNKfqMQpIEWG/4AQJCpgdlSqgMmZRAXMkVAFTMgkKmC+DP8BU5GogFlTqQImYRIVMHdCFXBnClTAMhj+A+5JogIWQ6UqMBsqM2cSFXAfNlTmCMmJeRKngHuKc0Nlw3/A2CQqYJFUqoBR9S/smUlUwKwJVeeqqqr5tyiKoihOnnzmmTArlkoHFmxlTs9JZVluNpu9g+v1uizLc07O87wJZIdWK/3PZZqJApG+bFarTIEKiECoy6s5VSd0Q1Ke53meNx9vNpvDULV3cvNBXdfqVczcatWXqHY7iQpYCpWSE5rCwF5dqiiKuq6zLNtut21gqqrq8fEx+7E01cas7pndb67/uUiMlSpDfkBkQl1eXdSHtJHosJcOw1abtPZO7jueCVVcLrpQ1V+gAgjF8F+MmoJTO5DXd0KjSU7r9XrvnKNTr2AOJCqADnf/DSleDJxw/pHmfsDRGgcBGfIDOCBUDekrMrUFqm5g2jty9KuEKpLXv7CnOAUsnFB1gSZjtff3rddrIYkFsU46wCCh6lxVVXUXoNq7H7BvJapGnud1XR895+r9iSKaqswSSFTA5OLcNXmAiernKopivV63S1VtNptRftm7a93+o+FcEhUQQnKXSKHqAmVZNvvPbLfb5kg7/Dc8DtjcGGiskFs8PT19/vy5+fjNmzefPn16enq6+099WdUzy7JVtjucRCVRAbSEqms0VavsJS3BvT09Pf3222/fvn1rPv3+/fvz8/Ovv/5631ylQAVwCaFqSFOXOvrQXtnp8DbAk18C5/vvf//7888///XXX+2RL1++/PLLL3/++ee9fqREBXAhK3oPaWZNHd1hpt2Upu3AoxvaHD2z+/31P+d48+bN9+/fD48/PDwcPX4ry1ABKbOieryOFp8ODzYT2A+PD98YCCc9Pz/3JaeBh64nUQFcRaga0kyc2mw2e8GoXV6huylNU6Cq63pvqYXDM+EiDw8Pr1+/vvSha/w4LX3vQYkKYJjhpxPadRPyPG8GAauqauen9+2d3Jw8cGb7zfU/5/j06dPz8/OXL1+6Bz98+PDq1as//vhjnJ9hqXRgLkJdXl3UT2ujUlee50fH9Q5P7jszE6o429PT06+//vrLL7+0uerDhw9///33169ff/rppxF+gAIVMCNCVdS6twEOb7Hcntzs9Dd8slDF+Z6env7888+mLvXw8PDu3buPHz9KVACHhKolEqq4VDMePebLZrXKjsWpTKICkhXq8mrvv8AG9rqRt7gvBSogWXFuCyhUBSY5EYZEBaRs+OoZKnIJVbA8PUN+4hTALaxTBQsjUQHch0oVLIYhP4B7EqpgASzsCXB/QhXMWn+cyiQqgFGZUwXzJVEBTEiogpmSqACmZfgPZkecAghBpQrmRaICCESoghmRqADCMfwHc2EZKoCghKrAbKjMOCQqYElsqMwRkhO3MuQHLI8NlYGxWSodIBpCFSTLkB9ATNz9B2mSqAAio1IFidll/yaqwziVSVQA4QhVkJI2MilQAcRGqIJ0KFABREyoghSYQQUQPRPVIXoSFUAKVKogbj1DfuIUQGyEKoiVAhVAUoSqwOz9x3H9iSrLVtmxueoAy2HvP46QnDiif8gvzvcRgInFufefieoQGZOoANKkUgXRMIkKIGVCFcShJ1HtxSnjxQDREqogtM7Y/3CiAiBmQhUEZcgPYC5MVIdwJCqAGVGpghD6h/wyiQogTUIVTM4kKoA5EqpgWob8AGbKnCqYkEQFMF8qVTCJob38JCqAORCqArOh8iIoUAGMKs6NUIWqwCSn+evZyy+TqACuFeeGykIV3I0CFcCSCFXnqqqq+bcoiqIoTp585pnMlkQFsDArw08nVVX1+Pi4dzDP8yZm7SnLcrPZnHNmlmWrlf6fqZ4hP79tgAmEurxaUuGEbqLK8zzP8+bjuq4Ph2y7iap7pnrVgqxWEhXAMglVJ7SJarfbNYN6u91uu902B7tpqaqqJlHled6cvNvt1ut1lmV1XfcVq5gVQ34AC2b4aUhbpjrspbYo1T5UFEVd14cn9x3PDP/NjAIVQBwM/8VooLzU1qjac5rk1JSmusqyHL9lxEaiAlg8oWpIE5ja2VFdfdOkDo8fxi9mxSQqALIss6TCsIEY1NafmszUnjkwJ71ZZGGkphGHnklU4hTAAglV1+jOSQ/dFgLp3PspUQGQGf67QlmW7S2BbYFqeGivyV5Hz1lda9TnxIU6BSqJCuBOkrtEqlRdYG8V0FHuLHD3X2L6C1SZRAUwqqsvkaFylUrVuYqiaBPVer0+XDdh4GubGwNNqEqeRAVAP5Wq07oFqvV6bYmEhTKJCoBBQtUJbaIa2MIv+3HdhPNXWyANClQAnMHw3wnnJKquw9POWW2BePXPSc8kKgA6hKoh7UjfOYmq7xY/a34m7GVVz8M4lUlUAPzI3nNDmtsH8jwfmEfVHfhrylrdeVdHD3a/v/6PVP/WyJk4BRC3UJdXF/Uh59yT2e3Adu/kPM+Loqiqqvk067kvVKiKVH+i8usCiJ9QFaNLQ1XWyVWtgflYQlWMbOQHkDihaj6qF8WLvjOFquhIVADpE6qWSKiKi0QFMAuhLq/WqYLeSVTiFADnE6oCG5i2pYg1EYkKIDUBd00eIFQFJjkFZsgPIEHDV08bKsPkJCoAxqNSxSJZiQqAsQlVLI9JVADcgVDFwhjyA+A+hCoWw5AfAPdkojrLIFEBcGcqVcydOAXAJFSqmDWJCoCpqFQxU/1xKpOoALgDoYo5UqACYHJCVWD2/htfz6IJmUQFMBf2/uMIyWlMClQAy2DvP7gniQqAoFSqmAXrpAMQmkoV6ZOoAIiAUEXiJCoA4mD4j2T1TKISpwAIQqWKNElUAERGpYoEGfIDID4qVaRGogIgSipVpMNKVABETKgiESZRARA3oYoUGPIDIHpCVWA2VD7BkB8AB2yozBGS0xCJCoBj4txQWagiVob8AEiKUEV8FKgASJBQRWR6ClSZRAVA3IQqYmLID4BkCVVEY7VSoAIgXUIVEVCgAiB99v4jNIkKgFlQqSIccQqAGVGpIhCJCoB5UalicpahAmCOhCqm1ZOoxCkAUidUBbasDZUN+QEwBhsqc8QMk9NRhvwAGE+cGyqbqM79SVQALIBQdZmqqsqyPPO0qqru3qD4vQz5HU6ikqgAmJPVUoafRtJUFAc6rSzLzWbTPZLneV+6Wq1m3f8KVACEEOryqlJ1gZOVp26iyvO8+aCu66Io7tmuKElUACyMUHWusiwfHx8HTqiqqklUeZ7vdruqqna73Xq9zrKsrutlDQUa8gNgeWY9/DSGqqoOs9TRTiuKoq7rw0f7jmcpD/8NDYNaNwGAoAz/Ja9JTk1pquucie3zIVEBsFRC1QlFUexebLfbc87vOzLzEcDV6miiMuQHwEIIVeNoA9PAnPQ5h6r+SVQAsBBCFTcz5AcAQtVYhqtQzfIKR89ZXes+z+Nyq9VhgSqTqAC4WXKXSHv/BZbo3X9Zlu2yfxPV/vFUnxAAcbn6Emnvv7QNL+/Z3Bg4pyVAm5e5RAUALZUqLmQGFQAco1I1jnPWTZhDpUqiAoAeQtXIDkPVOastJKBnGapMogKALMuEqhH13eI3h+WpenZHtrAnALSEqtE029HUdd3dl6bdZflw+5pk9CzsmWXRLOsAABFIdUPfINrNlfs6rd07Oc/zoiiqqmo+7fuS2DdU7ilQ/Xss5fUgAJgxGyrPQVVVzSBgXdebzaYNWEmGj/5E9enT5+aDN2/efPr06enpadKGAUCU4q6UpKl6UbzoOzPeSlXPnPT/+7+n33777eeff/7rr7+aI+/fv//nn3++fv369u3bqRsJAMeEurzGelFfhhhDVX+BarfLPn/+/O3btzZRNd6/f//w8PD7779P1EIAGCRULVF0oWowUWVZ9ubNm+/fvx9+3cPDw9HjADA9oWqJhjcnmvpXc2oZqufn59evX/d99bdv3wYeBYARndzdT6hanIgqVect7KlSBUD8Ql1e7f3HBZvPvJ54O6MAAAnCSURBVHv37vn5+cuXL92DHz58ePXq1T3bBwAJsKTCsvVsPtO3VPrHjx//+eef9+/ft0c+fPjw999/f/z48c4NBYDYCVUL1rNU+kDF9O3bt1+/fn14eGg+fXh4ePXq1devX3/66ad7NhQAEhDNnJ5FCjmn6rbdkZsZgl48AETInComtFod23nG7sgAcD2hamFuK1ABAH2EqiVRoAKAuxGqFuNYohKnAGAsQtUCGPIDgPsTquZOgQoAJiFUBTawe9Gtt4MqUAEwUyf3/gtCqArsXgtpSFQAzNfw1TNU5BKqZkecAoAQbFMzLxIVAAQiVM2IRAUA4Rj+m4WXwePzt0YGAMYlVKUvRKKylTIA7BGqEmfIDwDiYE5VyiQqAIiGUJUsiQoAYmL4L0GmpQNAfISq1ChQAUCUDP8lJaZEFee+S2fS+CDSbXmm8eFofChJNz4UlarALthQebU6jFOZGhUAyxNn5hOqAjtrwaeYClQAEJwNlbmKAhUApECoituxRCVOAUCEhKpYGfIDgKQIVVFSoAKA1FhSITKr1XWJ6sxJeeecNuL8vulbpfGj/0SNv8j0rUq38eNOJdb40b/VWD9u3NPivOmvJVTFxJAfACTL8F8c+uPUarXKjt39BwBERaUqAgpUAJA+oSoKEhUApM7wX3h7iUqcAoAUCVWBHUyZ+t+dDWftYAMAyxPnbYArV+6AXl4T7a8gxpcIACQnSLwRqgAARmCiOgDACIQqAIARCFUAACMQqgAARiBUAQCMwDpVwVRVVVVVURRFUYRuyw+ua1jzVX2PlmV5c7vG17Q5qrZd0aQker5pYYQv+KsblkS3Z3o+AmVZxtb5rYvaFn/PD7Rhul/Bjsnleb73W9hut6Ebtdvd1rDDr43/ZRZh265oUuQ9v91uD1uV53nodt3asMi7fbfbrdfrw4at1+vQ7bq1YfH3fFfb2tANOeLStkXe80f/oqd/5YfviKXp/pq7r9HguerGhkX+93ao/QsM3ZD/ua5JMfd8920uz/O9pibdsJi7ffdjcNl7gmFz1e0Ni7znu7ovs9Bt2XdF2yLv+aNhvTXZFTZ8RyxK+6Ls/oJjeFHe3rCJX7g36v75hW7Lv65uUsw9f/QZte/mAetVtzcs5m7fvTRvL6Yc/TOf2O0NC/4UzhdV5thzRdsi7/nmLTR4FTy63/S8HX1D2UXwYr29YXG+cew5Wh9OvUkxPIujBgpvbYKcvlW7kRoWbbfvBp9F31/6NEZpWMw93xVPXfbQdW2L8Il0NU8q+AC3u/+m007x65tMF2qiX7QNI2kDc1rbGaMD59xPtA0bS9P4gcGaUM8u2oaNrizLuq7zPB8ekwoi5rbdoq7rrPMnHIpQNZ2B94vmxd28JqZ3e8NSeSssiqL9/8TwrMbJ3NikmHt+4Aoa9o3v9obF3O1ZlhVFsV6vh++Emq41P/7cGxsWec83qqrabDZZlK29um0RPpc4CVXTmf01Js/zsixXHcH/0zB7Mfd8VVW73e7oe3F7WQ3SztsbFnO3Z1lWlmVzq/ze8fYph2rn7Q2LvOcbj4+PWZZF8t+2PVe3LYmeb+w1b+o4OP2I42I1HX50Gl3YO9Fub1gkt11cZB53/yXd88GnlO45v2Fpdft6ve42OPikk9YVDYu/55v/nbYvobDTB/fc0rbIez6S9RR2JqpPaeC3G0OouqVhbZWr+01ivp14F7rPj7qiScn1fIT3XTYualhC3b53sYknUV3XsMh7/jCmxBOqbmxbKj2fdRJet3mWVJihgfeO1EPVer1upj32fXk8b+WteYSqhHp+7yIaujn/c0XDEur23Utru+P7oVv0rysaFnnPH16/4wlVN7Yt8p7fbrdHm7ebfNGi8L/p5Ri44TPsBf6uDdsrOMdjHqFqQFQ9H8/ik3tGb1hU3b5nBgOvA4L3/NFXUSSh6q5tC97zwyZ+q7f339SivYfiTg0riqKu61A3Ni5ZJD1fVVUzNzbLsuE7vyZ2p4ZF0u1HNTffbTab2Jo3SsPC9nz3/bP7WtpbsCbIJoD3blvMr/nsx3VSJuh8oWo6zSvv6ENhk1a0DSN1bXDJ8zyq11K0DRtF84yOXj+KomjuqA8i2oaNqO9ZtMcD3i4Xc9tmQ6ia2tH4cnJNvAlc3bB26/KB/+uHfWpzFX/PRxtcbmlYKt2+3W5ju0be2LDIe/7oj27fV5tHQ/1Gbmxb5D2fvTTsaPOmXklkmlFGGk2fH96G0BwPONfkxob1ffn532F685hTFXPPRzKb5NDtDYu524fbEPaXcnvDIu/5Q9H+FewubFvkPd+04ei8rol/BTH+pmeszfInD07sooY1N+x0/4QGzgz+1PqkGKrS6vn2bW7bL4mGpdXtu85VZO+JBL9R69KGJdfzhxINVcn1/MmXliUVZit7sV6vm5tUY4j5FzWsvSa1R7o3pTd/itHe7dVKMVSl1fPZGZJoWFrd3thr3l4LU2lYij2/J9FQlWLPD7y0prwzMcbf9LztLYoTzyvy/IYdfZke/fJIntpR8whVu4h7/mir9iTRsLS6vXV0jksM972f37BEe75rNqFql0LPx/CaX+3Oe4thXO28v2xw6t/0bmxY9aJ4MWrr6KXng4i827t/zlE17/aGRd7zMxZ5zwd/zQtVAAAj+E/oBgAAzIFQBQAwAqEKAGAEQhUAwAiEKgCAEQhVAAAjEKoAAEYgVAEAjECoAgAYgVAFADACoQoAYARCFQDACIQqAIARCFUAACMQqgAARiBUAQCMQKgCABiBUAUAMAKhCgBgBEIVAMAIhCoAgBEIVQAAIxCqAABGIFQBAIxAqAIAGIFQBQAwAqEKAGAEQhUAwAiEKgCAEQhVAAAjEKoAgGQURVEURVVVe8erqmoeCtCmF0IVAJCMoijqun58fNw7/vj4WNd12FC12u12AX88AMBFmlyV53lbryrLcrPZZFkWNtWoVAEAKWmyVF3XZVk2nzaJarvdBm2XShUAkJqqqpoRwN1ut1qtsixbr9dNxgpIqAIA0tMMAjYfd4cCAxKqAIAkNTWqLPRUqpY5VQAAIxCqAID0dFdPCLuSQkuoAgASU1VVM6GqueOvvRMwLHOqAIDENLOpttttURSRLFKVqVQBAGlpBvvyPG8+KMsyz/MsgkFAoQoASEZZls3AX3cNhb3lQEMRqgCAZDQjfev1eu94c6R5NBRzqgAARqBSBQAwAqEKAGAEQhUAwAiEKgCAEQhVAAAjEKoAAEYgVAEAjECoAgAYgVAFADACoQoAYARCFQDACIQqAIARCFUAACMQqgAARiBUAQCMQKgCABiBUAUAMAKhCgBgBEIVAMAI/j8Cy5TQXXmd4gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "my2 = R.TF1(\"myFunc2\", \"pol1\", bin0, bin1)\n",
    "my2.SetLineColor(2)\n",
    "frp2 = hd.Fit(\"myFunc2\", \"LS\")\n",
    "my1.Draw(\"same\")\n",
    "myc.Update()\n",
    "myc.Draw()\n",
    "myc.SaveAs(\"test_fitting_2.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "quality-zimbabwe",
   "metadata": {},
   "source": [
    "Plot the NLL scan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "czech-unemployment",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-15T22:48:39.306405Z",
     "iopub.status.busy": "2021-04-15T22:48:39.306081Z",
     "iopub.status.idle": "2021-04-15T22:48:39.450004Z",
     "shell.execute_reply": "2021-04-15T22:48:39.450248Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "****************************************\n",
      "Minimizer is Minuit / Migrad\n",
      "MinFCN                    =      1.21602\n",
      "Chi2                      =      2.47176\n",
      "NDf                       =            3\n",
      "Edm                       =    4.672e-09\n",
      "NCalls                    =          158\n",
      "p0                        =      1.93421   +/-   3.65636     \n",
      "p1                        =      15.2263   +/-   1.81514     \n",
      "\n",
      "Covariance Matrix:\n",
      "\n",
      "            \t          p0          p1\n",
      "p0          \t      13.369     -5.1922\n",
      "p1          \t     -5.1922      3.2947\n",
      "\n",
      "Correlation Matrix:\n",
      "\n",
      "            \t          p0          p1\n",
      "p0          \t           1    -0.78233\n",
      "p1          \t    -0.78233           1\n",
      "\n",
      "****************************************\n",
      "Minimizer is Minuit / Migrad\n",
      "MinFCN                    =      1.21602\n",
      "Chi2                      =      2.47176\n",
      "NDf                       =            3\n",
      "Edm                       =    4.672e-09\n",
      "NCalls                    =          158\n",
      "p0                        =      1.93421   +/-   3.65636     \n",
      "p1                        =      15.2263   +/-   1.81514     \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Info in <TCanvas::Print>: png file test_fitting_3.png has been created\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAI8CAIAAAD0vjrdAAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dS3LiyP72cemNnhIyc3CwEyWzU70Qm8sCuuyJxMg+CwAcLOTfZyaxEYejrDloBbyD7MpQC0nchDKV+n4GHdXi4hQC9JCZ+qV7OBwcAAAA3Ob/6W4AAACADQhVAAAANSBUAQAA1IBQBQAAUANCFQAAQA0IVQAAADUgVAEAANSAUAUAAFADQhUAAEANCFUAAAA1IFQBAADUgFAFAABQA0IVAABADQhVAAAANSBUAQAA1IBQBQAAUANCFQAAQA0IVQAAADX4Q3cDmhPHcRzHYRje8nAhhBCi8A4Vz1zxKAAAYAf3cDjobkNDXNd1HOeK/RVCbLfb7JYgCI4jlHz+QoX3BwAANunK8F8cx9c90HVdlah835f/WCwWuZB09fMDAAA7dCJUhWE4Ho+veKAaszscDofDIY7jw+Ego9ViscjeU4Yq3/cPReimAgDAejaHqjiOXdd1XTcXgM5/uOyjiqIotz33DwAAgA5NVL+UykzHc8yjKMolKvm/zEYHAKCzujJRPY5jOQJ4/v7KiednzjGXd46iiFwFAEA32Tz8VwsVksIwDMOweshP1lxwMxgiBACgI+ipKqU6n8IwPFlSoaKegu/7RCsAAKzXXKiSxTOdf08/aqwq5tWhSpEX/WXLK2RnrKurC7MjgKrAVdkYYkUUAwAAV9PTZ1RYAqBGURSp8k4VgiC4dzMu2t/sFX+5QglBEMjtURSpO/u+X1hPQe174V85vz33uCcNaFFTaUCLmqq9AS1qqvYGtKipNOBOTa3XHXuqwjDM1TKQCUN15MienuzImu/7YRjeo+/qlp6q44fIm84Z11N/t3AOu+ue+/rf4540oEVNpQEtaqr2BrSoqdob0KKm0oA7NbVedympkB0OC4LgnDE+OQd8u92Ox2OjJiGpfqncxsVikZtoVSibILkwEAAAi9UfquREokv7nNSUI9m/5bquIeUJTGgDAAAw311KKsjamNfFkTAMD4dDEASGdFad04xzqi0QzgAAsBslFUqpOWEn51SpzrnjXFXxJE6rRp2tbECLmkoDWtRU7Q1oUVO1N6BFTaUBrZhTpbn4pyHdUc7v4g7ZwgfZEcnsPY+3y39st9vcPeM4lomqcGIWAACwieaeKtWRc+8ZVCd7qgov6MtewCiDUfZ/sxFKlaRS91TrMVf/Uft6Cq3cKcfS/bJypxxL98vKnXIs3S8rd8pp1X7paqoRL5AMNHdtyXWhyikqDOGUFPPM5iql+krGFr1Bz2flTjmW7peVO+VYul9W7pRj6X5ZuVNOq/aLUGX0asQqQp0sD3H+PZ1TFdVNODRXaNGn7iJW7peVO+VYul9W7pRj6X5ZuVOOYftVfQJ1NJ1D9b9AqoNHe0uaZ9QbtC5W7pRj6X5ZuVOOpftl5U45lu6XlTvltGq/rCr+eVL8mxovyy4LAwAA0DoNRbnjFOU4ju/74rcG2mCgFqX+81m5U46l+2XlTjmW7peVO+VYul9W7pTTqv2ysKdKpqjcLO8zV60BAABol7tEuVwlguMOKnPKU+nVotQPAEBb2FP8UwixWCx83z8cDofDQa7fIleekXfYbreu6wohiFYAAMAaNUe5k+WgcmWfLl132TL0VAEAUDt76lSdU8nzuKKmyUWq7odQBQBA7ewZ/pPDfCcrZB4OhyiKfN+XW2T/Vge55XQ3DQAAQ1WcPTWeQO8V5eRUqnNmTcVxHIbhdrvtYJ8NPVUAANTOnuE/nI9QBQBA7ewZ/gMAAOggQhUAAEANCFUAAAA1qDlUyaqetzxDHMeu64ZhWE+DAAAAGlFzqAqC4OqC6WEYuq47Ho9lRdB6GwYAAHBXd5keL4RQi/35vi9+O75n/Ju6f6eqgHL1HwAAtbOwpEI2Wp0jCIKudVARqgAAqJ2FoUpStT0Lb63ux7JeddVX8hYAAIVOlk23M1TlqIlW3UxROfRUAQBQO2t7qlCBUAUAQO2oqA4AANBif9z7D5wsrMA4IAAAsMB9QxUXABZKkmS1Wq3Xa8dx+v3+ZDKZz+eDwUB3uwAAwPXMGv5bLBau615aNbRdkiR5fn7e7Xb7/d5xnP1+n6bp09NTkiS6mwYAAK5331AVx3EQBPLfvu9HUXQ4HA6HQxRFURT5vi+3yy3qnuPx+K6t0mu1Wj0+Pn58fKgt6/V6NBotl0uNrQIAADe67/T4OI5lQir7K2EYLhYL3/dV75SsPGFxXfV+vy/7qHI8zyvcDgAALmLn1X8yUUVRVHaHMAx9399utypUyf4qW0cA0zQtS05pmrr/1nDbAABoC7eSrlY1Maequs9J3porCqo9VMlC8Dc+/HgvPM97eHgofIjneYd/u/qvAwBgt0MlXa26e0mFk7Tnp0Kyj+2KXJW94HGxWDhHlzROJpM0TeWlf8psNuv1eje0FwAAaNZET1V1NJERRPVmyYyld0LV1TnPdV2VqOQ0fMdxFotF9hWYz+dfX1/T6VRtmc1mn5+f8/n8uj8KAACMUN2BdiN1QV8QBMe3Zuda5bao6wSbp9p86YujUtTJjd/f3z9//vQ8z3Ec+d9fv37V0HQAAHA43DvelLn79PjscJhMGEKIOI6zRUHltX7ySkB5t+bHBNWFilnnvzjq4cfXLVZczygvT2AFQAAAaqTrxNrEMjUqLckglauxnksbWhLV7XIT7bOiKGrjHgEAgIs0F+XkvCLZR+X7vvitmb9+kZPltY7J7qhLl9lRUZrOKgAA6mJtT5XShUX9VEaUO2tsagQAALVreu2/OI7tHgsTQriuu1gsFovFeDx2XbcLaRIAYL3c7B0cayhUyajhuu54PJZRQ7IscIzHY/me830/W1Khor9K1X6tLg5rTrlYAECnJEny+vra7/eFEP1+/+XlJUmSZv70pWdG7afIu4eqOI7dTPWmnOrA0RbZvje5PrTskDv8LtCQXYcnR12H6Vx4/WcTOwYA6LYkSZ6fn3e7nVxjbb/fp2n69PTUTK66urRBA20rdPdQpeoUBEGQrT4VRZHsy9lut23PVdn258KT6oo72ScnayvU2i4AAG6yWq0eHx8/Pj7UlvV6PRqNlsulxlYZ677T41WRqrK/oqotGNX1Uu/VfxX7mLs8gWsAAQBG6ff7so8qx/O8wu2G0HU+vW9PlUxU2crpOWEYyv4qOyZXtb3LDQAAJU3TsuRUcVOXNTFRvTpq2BREbrywkRFAAIA5PM97eHi49KYua7qkgq3khHQ5zJej1t5puk0AANxmMplMp9PcxtlsNplMtLTHcE2EquqhvcIgYjhZ1TO7X2UT0s+fqA4AgGnm8/nX11c2V81ms8/Pz/l8rrFV5rr6esVzyP4bx3Gy1/1lqf6buzbjUmoSWNkd5K2yeoKidtZxnCAIcv9b9jxnbgQAQIvv72/HcTzPk//966+/fv36pbtRJ+g6k959enx2klAQBHIGlZx7pPqoLl0y795OXv0nd+p47Wd1oV9Wxd4VXp7ANYAAAKPIE1Mcx22ZBq3rTNrEX62efG1aorqd2p2Ta/9VvzJEKwCAdmb+1D95XZe1ocpxnDiOwzDM1lX3fT83LamDyt6pZr6DAQAd1MZTks09VShDqAIAmKyl5yM7i3/iOgcKVgEA0DaEKgAAgBr8Ue/TXd2/0sbeRQAALNbSsT+N6KkyFCOAAAC0S809VRVrJwMAgLagm+oKvGQ6nXzL8p4GAGjR6hOQrsbX3FOFS1WM8bX33QwAaDXzE5WZM2TuNadKFrN3Xdd1XSFEbjmX7N3kfe7UDPNVLCGku2kAABiqeg0+Xa26S6gSQozHY1U/fbvdjsfjjhdPvw7T1QEADTO/m8pY9Yeq7HI0QRAEQSD/vVgsyFUAAMBW9adR2bPi+352yC8Mw8Vi4ThOFEXZNYbjOB6Px05X5w+d82uAXwwAgCZZcN6xZJkaFaRyk6jCMPR933EcGaFwPkYAAQCNsSBRaXSXUCXzU+FNjuMwCHgdNagKAAAM1GhFdTm/So4D4kxJkry8vMiLKPv9/svLS5IkuhsFALAQ3VQ3qjlUyflSZX0qqo+KzqozJUny/Py82+3k/+73+zRNn56eyFUAAJjmLqEq+48cuY7NYrEoq1yFrNVq9fj4+PHxobas1+vRaLRcLjW2CgBgH7qpblf/K6gu9HMcJwgCIUQuYKk7ZEcDu3kgq2egHw6Hfr+/3++Pb/I8r3A7AADXaVeoOnkJl5Z9ucsrKITIjgAe/4ls8Cq7TxdUv4PTNH14eCi7dbfbVdwKAMBF2hWqqllSUkGK4ziKosJrAKUwDKvvAMdxPM8ri00VNwEAcCmbEpVG91pQuWK9v9wd4jhmflWZyWSSpul6vc5unM1mvV5PV5MAAEAhkqlOJ38ZJEny9PQ0Go1UrprNZp+fn5vNZjgcNtJGAIC1ttut7/v2dVNZNfxXIQxD+qXONxgMNpuN53me58ktvV6PRAUAuEWSJK+vr/1+X1ZAlFt0N8oGTUc513VzywJ22UVROo5jIYR9vycAAE2SFRCz9Xqm0+nX19dmsxkMBnrbVpeu9FThamWlvwAAOB8VEO+HniqdrovSdFYBAK7WhQqI9FR1lFtOd9MAALZJ07QsOVXcZKCKs6fGE+i9SioYSNZuuGjZwepyD7mnqnjm47LyCn1OAIDGyDKHZT1VLaqAWH321JWrOjT8J1/ii/Y3Vxo+J/dUFYcwCILCyHV1/yQjgACA67y8vJRVQHx/f9fVqnp1ZfjvcDhoSVT3/qPMEgMAtMJ8Pv/6+ppOp2qLrIA4n881tsoOnRj+O15q8EyymyqKopNX3slQxRx8AIDhZAXE4XDoeV6app7nUQGxLjaHqjiOx+Px7c9jYC2Dw+HACCAA4DqyHtV+v5cVEHU3xx5c/VcP2UHFWxMAYD71s5zTVr1sDlVCiMNvURRd+vCLBvLkQCHvTgAAOsvmUHUjNU0qDMNs9YuK5KRWklHuN8VKjgDe6ckBALZi9sj9EKpO2G63uUnu2+22LC0tFotcCYbxeEz3FQAAXUCoKqViUxAEx8OI2Snw2YAVRZG6s+/7juNst9uKuqDVNWENLBcLAGivdnVTte4USagqJYTwfT9Xt1MIoXJVdrvv+77vHw6HbL9UHMcyV1UUdDhcy2EEEABgtVtOkVq0KbHeQpVXqGV/ZaX1c6pSqb9bWOzq9l8M7frNAQDQqDunDF17Sk/VNWQ8qljBJndP524l1+msAgDAEIQqAADs151uKo1srqh+iziOZd9SxRxzOV9K3UcIUXGhH9cAAgBgN0JVMSGEnAtVGJXkxPPs0N52u1U5LEtlsvuFqgNL1gAAKnGaaAbDf/+Q4Sl3QZ/z79IJ6p7yH+rO8h/HpRPiOJbxKwiCe7QZAACYoyvR9eTVf3K6d/aCvux6zL7vCyFkj5TcclxqIXuT87v7Sm6p+KO1vP78BAEAVOjaaULX/nblVb4iVDn/zlVZuUQlZXOVUl12ofrCvYsOTdc+MACAM1l5gjh55TuhylDxb+K3ijtnJ1GdnEdV4xvdys8MAOB2HTxB0FPVRYQqAMBddfPsQPFP3IQqoAAA6EWoAgDATt3sptKIUAUAAFADQpU9GAEEACh0UzWPUGWhc1Z6BgAA9SJU2SNJkpeXF9d1hRD9fv/l5SVJEt2NAgA0Sv6upptKC0KVJZIkeX5+3u128n/3+32apk9PT+QqAOiCJEleX1/7/b78XS236G5U5xCqNHPLXfQ8q9Xq8fHx4+NDbVmv16PRaLlc1t1kAIBZ1O/q/X7vOM5+v59Op3b/rq44e2qcXkz3oE41ds/2+335WcrxPK9wOwDAGq+vr7vdLvu72nGc6XTqed7b25uuVmlERfUuquuop2n68PBQdutut6u4FQDQdvyuzqGiOq7neV5ZbKq4CQBggTRNy5JTxU24hz90NwD1mEwmaZqu1+vsxtls1uv1dDUJANAA+eO5rKeK39VNoqfKEvP5/Ovrazqdqi2z2ezz83M+n2tsFQCgAZPJJPv9L81ms8lkoqU9ncWcKp3qHfRNkmS5XK7X6zRNPc9L0/TXr1/D4bCu5wcAmClJkqenp9FopMYr5O/qzWbTzbMAE9W76E5HPY5jIQSV3wCgO5IkGQ6H8he153mTyWQ+n3czUTmEqm6691EnVwFAd8jvfPm7WndbNCNUdVF1gbLbDw2hCgA6omtf+CcrfGp5Nbj6T7NOfQYAAKhF9dlTV1F1rv6z2eFw0FitHwDQjK51UxmLUAUAAFADQpXl6KwCALvRTWUOQhUAAEANCFX2o7MKAIAGEKoAAGgrxv6MQqgCAACoAaGqExgBBAD70E1lGop/FovjOI7jMAwvfUjZrWVPVZF1+KgAAFDIzJ4CQm4xebQuenGEENvttuzWwqdq+EcGv2kAwBp8pVfQ9eIw/FegosMJAACgEKEqLwzD8Xh8xQNlN1UURYcidTfzGsysAgA70E1lJo7KP+I4Ps5SF704V4wYNv+p4HMIABbgy7waw38AAAAtRqj6hxBCDdVFUXTpw9syDYsRQABoO7qpjEWoqocMVb7vh2HoZgghNLcMAAA0glBVp+12u1gscltc1zWqH4vOKgBoL7qpTEaoqoeKTUEQHA8jVlxO6F6rgZ0CAECj1p0iCVX1EEL4vh8EQbZyuhBC5aqyiuqF9RfOcUtr6awCgHaRVXu61k2l5RR5C5apqUdZZpJha7vdGjUCCABohSRJVqvVer3e7/cPDw9yy2Aw0N0uFKOn6u7kXPWKFWwAADiWJMnz8/Nut9vv947j7Pf76XT69PSUJInupqEYoaqj1AggaQ8AzLRarR4fHz8+PtSW9Xo9Go2Wy6XGVqECw381iONYju6VDQI6juP7fmPtOYf8odPv92WX8mQymc/ndCkDgDnkqN/xRs/z3t7etDQJ1eipqoEQYrFYLBaLwolTssiCUQWrZJfyZDJRXcppmtKlDADmSNP0OFGdvAl6EaquIYQQQmT7pWRH1HHpBJWlKjqxmkeXMgAYzvM8OTP9opugF8N/1ziehxSGoUxUruv6vi+EiONY3S0IgqabWIkuZQAw32QySdN0vV5nN85ms16vp6tJqEZPVT2yJalkXfVsojKqm4ouZQBohfl8/vX1ld0ym80+Pz/n87mmFuGEbpURa0D8m/it4s7VFTjvd2jk/PTj7Z7nEaoAwBxJkgyHQ8/z0jT1PE9eVDQcDnW3S7+TJay1xBtClU66auO+vLyUdSm/v7833x4AQCF1mpC/1XU3pzV0nV4JVTrpOupJkjw9PY1GI5WrZJfyZrPhBxAAmKNr69LUhVDVRRo/LUmSLJfL9Xotu5TTNP316xeJCgDMQaK6GqGqi0z4wMguZRNaAgDI4pv5aoSqLjLqA2NUYwCg4/hOvoWuV4+SCgAAADUgVAEAYBa6qVqKUIV/HA6Hk2U/AABAGZap0awix/AzBQA6iG6qc5jZC0Co0syoT47srDKqSQAAHKs+VemKXAz/AQBgCn7ZthqhCv/CzCoAAK5DqAIAwAh0U7UdoQoAAKAGhCrkMQIIAM2jm8oChCoAAIAaEKpQgM4qAGgS3VR2IFQBAADUgOKfmlFRHQA6jm6qK5g5nEJPlWaHctobZuZbFgDssN1udTehxSrOnhpPoIQqAACakyTJ6+trv98XQvT7fdd1v7+/dTcK9SBUoRSdVQBQryRJnp+fd7vdfr93HEf+9+npKUkS3U1DDQhVAAA0ZLVaPT4+fnx8ZDeORqPlcqmrSagRk+N0asXkxFY0EgBaod/vy96pHM/zCrfjOrrOXPRUAQDQhDRNy5JTxU1oEUIVAABN8Dzv4eHh0pvQItSpwglyuvrhcNhut77v624OALTYZDJJ03S9Xmc3zmazXq+nq0moET1VOEFek6Ku/n15eeEqFQC4znw+//r6ym6ZzWafn5/z+VxTi1CnDs1BjuM4juMwDG9/EiGEECJ3U8UzF97fOVUQ1oRDI6/+zV6rMp1Ov76+NpvNYDDQ2zYAaKMkSYbDoed5aZp6njeZTObz+XA41N2uljlZ8UfLObRDoUoegBv3Vz5JEATHEariABfe32nDhXWvr6+73S539e90OvU87+3tTVerAKC91De//Imuuzl20nV67cqcqjiOb3+Sind/Lc9voPV6fXxBynq9JlQBwI1IVPbpxJyqMAzH4/GNTxLHccU6TTJU+b5fuAjRjWOOunD1LwDUy/wBCtzC5p6qOI5vz1JKjU/VFvIS37I6dVz9CwBAVid6qm6nplKV3UH2VNnXlzuZTKbTaW7jbDabTCZa2gMA7UU3lfVsDlVCCDUAF0XR1c8jB+98368YxZMjg/aFKnn1bzZXcfUvAACFbA5VtYjjeLFYOOdNRZeXcrgZbZ/APhgMNpuN53me5zmO43nearXabDZc/QsAF6GbqgsIVSfIqVRndnQtFovcZPbxeNz27qvBYPD29rbf76MokvOrSFQAABwjVFWRecj3/epglO2OiqJIjTnKRV22223FuKF7rXr28BLyRZCr1jT/1wGgveimuk6LTpESoapUGIay2+mcITzf92U9hWz8iuNY5io5gFiosATDOW7fQQBAA0hUV2vdKZJQVUxNpTpn4E8IIZevOb5J9VG1fXJVFp1VAAAcs7lO1S2yYeg4D6k1BE9W9VQdVyxHAAAdRDdVpxCqTigcudtut6qGQjejkuys4psCAACFUFUsDMPCXih1cZ+cLCUTlbxndcDqZvYCgC7jx2fXEKqKyWlSx9vd36XVs5FLLgtYOFCo7mZfqKKzCgCALCaq/0P2M1238rF81HHpBDXbvWJ9GwCAlfjZ2UH0VP0jV7TzIkII3/e32+1isVgsFjJCye4reYfrspr56KwCAEAhVNVDXtwnU1R2brvv+zYVUwAAnIMfnN3EUa9ZdhLVyXlU1dWe2nJo+O4AAGm73crLmPhivLeT5RK1vP4cdZ3s+NTZsRcAcLUkSVar1Xq93u/3Dw8P+/3++/t7MBjobld36ToxMVEdt6LAOoAuS5Lk+fl5t9vJJeflf5+enpIk0d00NI1QBQDA9Var1ePj48fHR3bjaDRaLpe6mgRdGLjRyaaBM5v2BQDO1+/3Ze9Ujud5hdvRAIb/AABomTRNy5JTxU2wFaEK9WBmFYAO8jzv4eHh0ptgK+pUAQBwvclkkqbper3ObpzNZr1eT1eToAs9VagNnVUAOmg+n399fWW3zGazz8/P+XyuqUXQhlClmVtOd9MAAKcNBoPNZuM4jud58r+9Xm+z2QyHQ91Ns1nF2VPjCZQrtnSy8oo5uVOqrDAAWE99mcsly3Q3B1z9ByvIYnf9fl8I0e/3X15eKH8HoDtIVB1HqEJtZFnhyWSiygqnaUpZYQB2s3LMAdchVKE2x2WF1+s1ZYUBAB1BvtbJst83lBUG0DWWfY1bgzlVaDfKCgMAOo5QhXpQVhhA19BNhRwqqqM2lBUGAHQZPVWojSwrPJ1OsxspKwzASnRT4RihSjMDC8JeTZYV9jxPlRV2HIeywgDsQ6LSjorqyLP4Y6nKClu8jwA6i282w+k6QLwtdOrCx7IL+wigU/haMx8lFWCnw+HQxqFMAAAuRagCAOBcdFOhAqEKd0dnFQCgCwhVAACchW4qVCNUoQl0VgFor+12q7sJaAdCFQAABZIkeX197ff7Qoh+v++67vf3t+5GwWiEKs0MrF12J3RWAWiRJEmen593u51cD17+9+npKUkS3U2D45ha/JNQpdmhnO6mAUB3rVarx8fHj4+P7MbRaLRcLnU1CVkVZ0+NJ1Dm3BWL4ziO4zAMb38SIYSsLX6sg3MeO7jLANqo3+/L3qkcz/MKt8MoVFQ3i+w8vPHFkU8SBEFZOOtgwujgLgNonTRNHx4eym7d7XYVt8IEVFQ3SBzHtz9JWe9UxzGzCoD5PM8ri00VNwF/6G6AccIwXCwWNz5JHMdcggsA7TWZTNI0Xa/X2Y2z2azX6+lqEsxHT9U/4jiWlwzcnqgcxxmPx7c/ia3orAJgvvl8/vX1ld0ym80+Pz/n87mmFqEFCFX1U1OpdDcEAHClwWDw999///z50/M8x3E8z+v1epvNZjgc6m4azEWo+ocQQl2KGUXR1c8j56T7vn/jlYN2o7MKQCu8vb3t9/soivb7/fv7O4kK1QhVdYrjWI4e1jLVHQCgS/byMS48wpmYqF4nOZXqoo6uqzts2l6YQHZWtX0vAAD307oxDXqqaiN/yvi+f9FvmuqasAaWi70HrpQEYBR+8hmidadIQlU9wjCUyYCBv/N9f3+7rqsWK315eWFRLQDakahwNUJVDdRUqltmuHeNXKzU+b1M6X6/T9OUxUoBAO3FnKoaqAv95GJ/uVvVGoJcD5glFyvNblmv19PpdLlcvr296WoVgI6jmwq34N1TII5jOeX8zBdHCHHOrKAoinLTrbr86WWxUgCm6fJ3smV0HUp6qmoQhmFhL5RKWr7vO1yUm5GmaVlykjextBYAoHVI5QUu7akqo0qrlw38dflXET1VAIzS5S9k+9BT1Sayz0kIwTSpq7FYKQDAMlz9d43tdrvdbmupnuCWu/3JTSYXK51Op9mNLFYKQAu6qVqn4uyp8QRKqNLMwNplzRgMBpvNxvM8tVip4zgsVgqgeSSqNjKz+CfvJJ34JEtxHMsRVV4QAM3jm8c+uo4pPVXQj+siATRMXZ1NokKNCFUwyOFwsH4yGQCNkiR5fX1Vq2O5rvv9/a27UbAHoQoA0AlydazdbqdWx3Ich9WxUCNCFcxCZxWAO5GrY318fGQ3jkaj5XKpq0mwDGPJOjGWX4iXBcA9UHO4O5ioDvyDzioAtTu5OlbD7YGVCFUAAPt5nle2qGjFTcBFWKZGs4oumS4PgcnOqi6/AgBqx+pYNjFzQIOeKs0MLAhrCAYBAdRLro6V3TKbzVgdq6XMrKhOqAIAdMJgMPj7779//vypVsfq9Xqsjj6ZutEAAB+vSURBVIUaMcKiEyNcJ/ESAaiR+kpRq2PBSrrOHZyxdCIxnMRLBKAufJ90ByUVgALMrAJQCxIVGkCoAgAAqAGhCqajswrAjeimQjMIVWiT7XaruwkAWoZEhcYQqjRzy+lumkG+v79d1+33+0KIfr//8vLCqvIA0GUVZ0+NJ1BClWYG1i4zTZIkz8/PjuPIxbn2+32apk9PT+QqACfRTWUrin8C11itVo+Pj9kt6/V6NBotl0tdTQIA4BgRXid+Qp2j3+8XLiDveR4LywOowHdsZ1GnCiiQpmlZcqq4CQBIVGgeoQpG8zzv4eHh0psAAGjeH7obAJwwmUzSNF2v19mNs9ms1+vpahIAw9FNBS3oqYLp5vP519fXdDrNbvz8/JzP57qaBMBAqo4diQq6EKpgusFgsNlsPM/zPM9xHM/z/vrrr//973/D4VB30wDolyTJ6+trto6d7hahu4jzOvFz6lJxHAshHF46AI7j/K5j9/j4+PHxoTb++PFjs9kMBgONDYNeXP3XUQYWhDWZTFQOCwICcBzndx27bKJyHIc6dl1gZkX1Dv3cj+M4juMwDK9+rDyjq/N6TsUzCyEKH0V3y414AYGOo44dCuk6O3TonCSj66X7K4Q4XsQ3iqLjkFQRjYMgKIxcZIIb8QICXZamaUVdld1uR9WVzmL4777iOL7iUa7rqkTl+77v+/Lf4/E494TXPT9uxCAg0GXUsYNpOlGnKgzDxWJxxaPkP3L9UvIsPh6PsylYhirf90lXDZO5iv4qoJuoYwej2NxTFcexnLB2RaJyHEc+yvf93EhfFEW1NA8AcCNZxy67ZTabUccOutgcqmpxPBdKZaxsp5T8d9kcdtwVg4BAZw0Gg7///vvnz5+qjl2v19tsNtSxgxZdGTeJ43g8HjuXTFSXceo4VKmnyg4LypN64QT2Coxb1YVXEugs9fFX12gDXP13X1eEqjLqesDsU8lQFQRBHMfZqwWrYxZRoEa8mEAH8cFHIV1vjE5MVL9d/JsKTIUzq44nb43HY2avN4MZ60DX8JGHaZhTdRZ5/WA2UWX7n7KZKYqiw2+yBMN2u62oC1pdE9bAcrEAADSjdadIQtVZwjAMgiBbpyqXk2QVq8PhkAtb8iEVlx8ernWP3Ww7ZqwD3UE3VRe07hTZlTdljXOqVNWrc6alF85qV/hSuAf1qm63W5WDAdiEL09U0/UOoafqYqqP6pxlBAvrL+DeXl9f+/2+EKLf77+8vCRJortFAAD7EaqKqZnphbeqyVKNtglnSJLkzz//fHt7k2up7vf7NE2fnp7IVYA16KaCsQhVpcbjsRy5K6OGlsIwDMOwui+K6inNWK1Wj4+P2S3r9Xo0Gi2XS11NAlAjEhVM1pV35xVzqtzyep6qKpUcAZSVqwpLJ6gJWIV/l2+H2vX7fdlHleN5XuF2AO3C1ybOwZwqzYQQQojjaVLHnVXqPipsyS3HpRPiOJaJKgiCmpuLImmaliWnipsAtAWJCobryhv0ZE+V7HzK9japhzi/l1XOFv9U3VSSKrPu/I5Q2TtX/NGOvP6NoacKsIy6jJcvTJxP17uFiuqlhBBRFMlctd1us9PSc4nK+b3mlLxPtirVyXLqFXWV+Pq4wmQySdN0vV5nN85ms16vp6tJAK6QJMlqtVqv1/v9/uHhYTKZ6G4RzGJmVUKC/2nZywDlKGHFnbODgycnp/PDq3ZJkjw9PY1Go2yu+s9//sOq9UCLJEny/Pz8+Pj48fGhNv748WOz2QwGA40NQ1voOr1yUteJUHUPSZIsl8v1ep2mqed5k8nkv//9L68z0CKvr6+73S6bqBzHmU6nnue9vb3pahVahFDVRYSqu5Jjsg6vM9A2TI7EjQhVXcTJvjG81EBbpGn68PBQdutut6u4FZAoqQDcEWstA23heV5ZbKq4CTABV/8BAMzCZbxoKXqq0BV0VgFtMZ/Pv76+sltms9nn5+d8PtfUIuAshCp0CLkKaIXBYPD333///PnT8zzHcTzP6/V6FEaB+Zi9qxOzp7XgZQcMl/2Qqst4gfNx9V8XVfeacGjuhFAFmIxPKM5xctiBZWq6iO+O5slBQF55wEB8NnGm6veJrpkezKlCFzG5CgBQO0IVAMAIdFOh7QhV6Cg6qwCjkKhgAUIVuiubq7bbrd7GAF1GooIdCFXoutfX136/L4To9/svLy9JkuhuEQCglQhV6K4kSf7888+3tze57v1+v0/T9OnpiVwFNIluKliDUIXuWq1Wj4+P2S3r9Xo0Gi2XS11NArogO9pOooJNeDfrRPFPvfr9vuyjyvE8r3A7gFskSbJardbr9X6/f3h4mEwm7+/vfNHhOmYW/6SnSrNDOd1Ns1yapmXJqeImANdJkuT5+Xm326nR9vf39x8/fjDajutUnD01nkAJVegoz/MeHh4uvQnAdeRo+8fHR3Yjo+2wDMN/OjGZQK+Xl5c0TdfrdXbjbDbr9Xrv7++6WgVYidF2NEnX6ZWeKnTXfD7/+vqaTqfZjZ+fn/P5XFeTACsx2o6OIFShuwaDwWaz8TzP8zzHcTzP++uvv/73v/8Nh0PdTQOswmg7OoLhJ50Y/jNHHMdCCPlvjgtQO0bb0SRdX+OcPHTi5G0mjgtQuyRJhsPhdDpVuWo2m31+fm42G/qGUTvmVAGmYK1loHaDwcBxnOxoe6/XI1HBMvwi14keEZNxdIAaZT9Q2dF24B4Y/jNLHMdxHIdhePVj5VdG9RcHFdUNR64CasFHCbUzs6I6b/Ri8mhd+uIIIbJrWklRFJVFK75oDMcBAm7H5wjNY06VQeI4vuJRruuqROX7vu/78t/j8fi6J4R2TK4CbkSiQqf8obsBxgnDcLFYXPEo+Y9cv5Q8JY/HY75WWkrmKg4fcAU+O+ga3vH/iON4PB7nNp7/4sjw5Pt+rlNKPW3hU/GN0xYcKeAKfHCgC8N/Njie2K56rRgBtMbxtDkAx0hU6CBC1T+EEIffoii69OFBEARBcDwhnSxlBzkI+Pr62u/3hRD9fv/l5SVJEt3tAgxFokI38b4vUD1mdxF1PSDDf62WJMnz8/P//d//qS3T6fTr62uz2ciShkDHbbdbdXUO32zQjuE/e8gCV0IIdT3gFV1fMMpqtXp8fMxuWa/Xo9FouVzqahJggiRJcj24JCp0Ge/+Ajf2VOWqVVXXqbqqgY5DadBm9fv9/X5/vN3zvMLtQBfIHtzHx8ePjw+18cePH/Tgoi6tO0vSU1W/MAyDIMjWqaqozH64VkM7A8dJ07QsOVXcBFhP9uBmE5XjOPTgokatO0XSU1WgxjlVqupVYX8V/eRtQU8VcIzPBYzFnCo7qT6q65YRhCEmk8l0Os1tnM1mk8lES3sA7ejBBY4RqmoQ/1Z4qxwHpLhRq83n86+vr2yums1mq9VqPp9rbBWgked5Dw8Pl94E2I1laupxcrhQTbFCGw0Gg81ms1wuPc9L09TzvF6v9+vXr8fHRwZw0VmTySRN0/V6nd04m816vZ6uJgF6MaenwBVzquQVCmUTpxzHCYLgeASQOVVtFMexOsocQXRZkiTD4XA6napcNZvNPj8/N5vNcDjU2zZ0HHOq2kQIIYQ4DknHqweq+5RVVUDrZA/l4XC45YpfoNWGw+H397fneZ7nOY4je3BJVOgyhv+ucTxBKooimahc1/V9XwgRx7G6W+EKNrCDzFX0V6Fr1Nv+7e3t7e0t24MLdBahqh5CCJWrttttNnUVDvzBJuQqdM3xG55EBTjMqapd9jJAOUpYcefqkSMOTbuQq9ARvNVhgpNTL7S8S/ls6MR3k2U4oOgC3ucwHxPVgdZj3jqsR6ICKhCqgDrlchVFX9FquTcwiQqoRqgC6pckyevra7/fF0L0+/2Xl5ckSXQ3CjhX4RuYRAWcRKgCavb9/T0cDne7nVz+bL/fp2n69PRErkIrJEny/PycewPLklS6mwaYjl8eOvHLz0qvr6+73e7j4yO7cTqdep739vamq1XAmXgDwwK6Tq+c1HUiVFmp3+/Ln/g5nucVbgeMwhsYFuDqP8AGaZqWnXgqbgIMwRsYuAU9JTpR/NNK/NBHq/EGRiuYWfyTZWo0IznZZzKZpGm6Xq+zG2ezWa/X09Uk4Hy8gdEK1WdPXSUD6anSiTlVVkqS5OnpaTQaqdPSbDZbrVa/fv0aDod62waclCTJcDicTqfZN/Dn5+dms+ENjLZgThVgicFgsNlsPM/zPM9xHM/zer3er1+/Hh8fdTcNOE1WT8i9gUlUwDnoKdGJnirrxXGcXVSbIw7D5d6iuTcw0BaUVOgiTrEdxEGHsXhzwhoM/wGdwKLLMJDruiQq4HaEKqBp5CoYRcYpEhVwO0IVoAG5Crpst9vs/9JBBdSIUAXocZyrcmc7oEZJkry+vvb7fSFEv99/eXlJkoREBdSLUKWZW05303B3MlcVnu10Nw1WSZLk+fl5t9vJkuj7/T5NU1k6QXfTgCtVnD01nkAJVZodyuluGprw/f09HA5zZ7unpydyFWq0Wq0eHx8/Pj7UlvV6PZ1Ol8ulxlYBt6g4e2o8gdL3qxN973h9fd3tdtmzneM40+nU87y3tzddrYJlWM4PXUOdqi4iVIGzHe4tTdOHh4eyW3e7XcWtQEtRpwronDRNy5JTxU3ARTzPK4tNFTcBuMIfuhsAdJc8pZX1VHG2Q10mk0mapmqBZGk2m/V6PV1NAqxEqAJ04myHBry/vzuOM51O1TttNpt9fn5uNhut7QJsw/AfoNN8Pv/6+ppOp2rLbDZbrVbz+Vxjq2ATObnk+/vb8zzP8xzH8Tyv1+ttNpvhcKi7dYBVCFWAToPBYLPZ5M52v379enx81N002EBN1x0MBm9vb/v9Poqi/X7//v5OogJqx9VnOlUXKOPQdE0cx0II9b9cHIob8RaCxU5W+KSkwn3FcRzHcRiGVz/WcRwhRPa0l1XxzGWP4isP1XiH4Gq8edBl1Km6O5lqL93fMAwXi0VuYxAExxGqIjUX3t/hWw9nuO59i07Zbre+72e38N2CjqNO1X3JfqZLCSFUovJ9X31tLRaLXM/Tdc8PnCSXXGDpZRwrWzWSRAXo0olQFYbheDy+9FFxHMtTl+/7h8NBjgAeDgcZrbbbbTZIyX/Lex67bswRUFh6GTllaySTqACNbA5VcRzL1aqPx+/OoZJQrhdK/S9RCU1i6WVkla2R/PPnT42tAjrO5lB1I9lNFQTB8U1yY3YIRk1jb6Zt6KDVajWZTHIn0dFotFwuNbYKuqzX69w63HJjrpAsgCZ1paM4jmM5Anj+/spZLFEUHUclNXtdPVvFnav/REdef9yOpZehsEYyUI2J6saR06EKQ5KaQXW8XQjhZjCBHbVg6WVksUYyYCZC1cXCMJQDf8dzqhaLRe6yrPF4zJggbsdJFDmTyeR442w2K9wOoBmEqsuoIgu+76u0lO2OiqJIXfSnrhOsmNLuXuueewkTTSaT7BKB2e3NNwbavb+///jxI7dq5OfnJ6tGwiatO0USqs4VhqHrumr2em5cT1axyg0XxnEsc1XF5YeFJRjOcafdhLEKl17+z3/+89///ldjq9A8ec44HA7Hq0ayRjIs07pTZFcmSl8xUf34sdJFs9HVYwsfxUR1XCRJkuVyuV6v0zT1PG8ymcznc1mayKHqunWO66Q7JV8auVUjAeg6vdJTdUK2cGgQBIeSqetlCocIgesMBoO3t7f9fh9F0X6/f39/l90S8pcZg8J2qKiTXnaeIFEBhiBUVVGlE2ScotonDFF4Ei3LVaxp0yKFddKfnp5knKIzEjDcH7obYK44jmWiOjneJ8OWEKLibvyURANUrjocDkmSrFar9Xq93+8fHh7kWOFgMNDdRlQprJPuOA510oFWoKeqlIxKQRCczEMyfhX2Y6mNhCo0Qw0FFnZ4sKaN4QrrpDu/oxUAw3VlovTJieoy9AghVAySv/irQ5W8ST15EATZaFW2XWGiOu7k9fX17e0tt3E6nXqed7wdhqBOOlAXXafXrpzUT4YqGaF835czynNX/JVRzyaEUDNX5MqAcRyrLRV/tCOvPxrGmjYtxYEDaqHr9MqcqmKXXqwnL2mWKSpblUqltDIVV2yRt3Cdk2va0OFhrMIDN5vNer1e840BTGbm9c70lNQsO4nq5DwqeqpwJ3R4tJHrut/f309PT6PRSE2iknXSqeoJXIThvy4iVOFOXl5e0jQtnN3MW06jwnqezu/f3PLQlJV4bbqtQJsRqrqIUIU7SZKkrMPj8fGRd13DKspbVFTDp046cDUqqgOozWAwKFsYThZcKJuOQKXQ2pXV80ySpLqkJ4kKaB16SnSipwoNKOvwyA05USn0Tl5fX3e7XWH1KT7+wJ0w/NdFhCpoJydHPz8/Zwt5T6fTr6+vzWZDrrodFw0AzSNUdRGhCiYoHAqkUmgtqOcJaEGo6iJCFUxAV8pd8fICzWOiOgANTlYKbbg9rVM9td913f1+P51Oc9tns9lkMrlnuwBoQKjSzC2nu2noBM/zykagKm5yOn+dYJIkr6+v/X5fCNHv919eXnKLVctP8eFw+P7+/vr6yuYqWd5iPp833mrAHhVnT40nUEKVZodyupuGrphMJoVdKWmaHn83nQwTXVBRJcHJxCn5Ka4ob6F3L4BWqzh7ajyBEqqArpvP54VdKb9+/VJFrWS6qg4T3bFarbIXSzqOs16vR6PRcDgsLD01GAze3t72+30URfv9/v39nUQFWImJ0joxUR2GOGdpFJmrJpNJruRSB68TZO45YDiu/usiQhVMU700SnfCRNkifQ5VEoA24Oo/APpVJKpbrhNsy6z2c2aMyalRhaqn9gOwHqEKwFmqE0PhTe2a1X5yxpiaXvbz50+qJAA4RqgCcK6y6wT/+uuv4yuZTZvVfrK3rGz6+XK5zF7Qdzgcyqb2UyUB6DhCFYBzVYQJFThUuqrIKNV/5fyxwnPueX5v2Xq9Pl72eL1ev7+/5y7oo0oCgGLVlR5wV7z+aJ3v7++fP3+qMPHXX3/JygvHyr5zPM8re+aXlxc5jPjw8PDz58/v7+/b7/nnn39mR+Wm0+mPHz+O7189J2y325W9IFEUld0EQBddp1d6qjQzsCAsUOHMkktpmpY9Q+Gs9vPHCi8aVTynt0x+3Cqmi1VPJquY2g/gfsysqM4l/TpRUgEWK6u/kKXe/6+vr7vd7pwKWOff89I2vLy8pGm6Xq+zt85ms16v9/7+Xv0kAIxCnaouIlTBYudklJM/KI8rYJ1TK+uc36nHBaWSJHl6ehqNRqrNcsYYk6WA1qFOFQCrnHOJnJyFUNGZJNcfzKqolaXuo+Y3XLRWNNPPAdyIUAXgLs7PKBWTlo5ntZ9zT7WxrAZEWUEpFukDcAuGn3Ri+A8dUb36jXPJfKaLZj4xogd0E8N/AKx18hK588tpXlR4kxE9AE0iVAHQ7/z0c2lOYkQPQGMYftKJ4T/g2MmxwivuCaBTKKnQRdUXfnNoAAAodLJyipZz6B/N/0lkkZwAALhU9dlTV1F1QlWxOI7jOA7D8OrHOo4jhGBsAgCAjmD4r5gMuZe+OGEYLhaL3MYgCMrCGXOqAACoHXOqDBLH8Xg8di4MVUKI7XYr/+37vuM42f+VfVc5hCoAAGpHnSpThGEoE9VF4jiWEcr3/cPhIEcAD4eDSleFoQoAAFiDnpJ/qN6prPNfHNVNdfwQOZJY2FlFTxUAALWjp6rdZKIKguD4JrlRDQUCAAArEar+IYRQq7FGUXT1k9TaKAAA0BqUVKhHRTejHPWTk6sAAICtCFX3FYahHPirqKpw3TMzGQsAYDddNTyvRqi6IzV73ff9spFBshEAAIWuPkXqSmPMqbqLMAxd11Wz16mnAACA9eipqlmuNEMURcxeBwCgC+ipqlO2cGgQBIfDgUQFAEBH0FNVG7XwX8VifwAAwFaEqnrEcSwTFeN9AAB0E8N/9ZBdU0EQkKgAAOgmeqquIZOTEEIN86lVaCou9OtI3rJ1QUMr98vKnXIs3S8rd8qxdL+s3CnH3v2qEaHqGrmF/FSQWiwWchCwUOF7saKWBu9dAAAKmVkXlFBVg1vKUJGcAAC4VPXZU1fkoitPJyu7Uq3cKcfS/bJypxxL98vKnXIs3S8rd8pp1X7paioT1QEAAGpAqAIAAKgBoQoAAKAGhKp2OH/O3UWz8+70tHobYOVO3e9p9TbAyp2639Pe4zk7frDu9Frd6Tm1vwLna1FT60WoAgAAqAGhCgAAoAaEKgAAgBpQ/FMzKqoDAHApM6dYEao0IzkBAHApMyuqM/wHAABQg9aUnLeSmb2XAAC0nZZ4Q6gCAACoAcN/AAAANSBUAQAA1IBQBQAAUANCFQAAQA0IVQAAADWg+Kc2cRzHcSyEEELobstZrmuwfFTZrWEY3tyuRsndaUWzr2iqBQdLtr9Fn6yrG8zBal6XD5YShmFbjpceBzTO9/3cUYiiSHejqtzS4OPHtvrt16JmX9HUVh+sKIqO2+z7vu52lbqxwa0+WEEQHLc5CALd7Sp1Y4NbfbCy1I7oboi5eGmalvsCVf82Nlfd2GBrvk0OmbOg7oacdl1T23uwsgHF9/3cjuhuXYHbG9zeg5UNKLl9NzNX3d7g9h6srOybVndbzMVL0yj10comEpPfprc3+PjhLZX9YtXdlhOubmp7D1bh/qpzgIH9Vbc3uO0HKxdHCr9qDHF7g43dtYu0MQg2j5emUYUfzoPBH7nbG9z2T2DhGI3uRhW7vakm712Fim45lS+bb1WFWhps4H6do2IHy75t9KqlwS09WFnm9/4agqv/mqMmKpZNSzRtumLrGoxuqpgCrKbTVtynea1rcI3kflUMh5m2461r8D2EYbjdbn3fL5xbhixCVXMqPnvynbrdbptrzRlub7AFXzdCCPUTpLAryBw3NrW9B6vitGfmNUq3N7i9B0sIEQRBxe8x0w7Z7Q1u78GS4jheLBZO+3ekGYSq5nT2q9/3/TAM3Qwz97fj2nuw4jg+HA6F3/jqXGjUXtze4PYerDAM5TX5ue3q1TBtF25vcHsPljQejx3HMfwnpTmoU9Uc0zqiTqqrwdvtNvdU2+3Wdd0oitrytdIdNh0s9Qu7+torc1zaYAsOlgyRcq8dxwmCwPCWX93glh4s2Tbf901upFman8bVWfIFL5zVaOa1+rc3WJ0bsk/S3utyzTxMha5oqmUHq0VXa0oXNdiOg5Xr/DBtivqx6xrc3oN1PEnfzCs/jMLwH+5ICCHnNmZnJAgh1BcKU93NYc3BiuPYdV3Vl3D496XgBrqiwXYcLDldSVV+WiwWruvqblSV6xrc3oMl35MM/F1Gd6rrEPk5bFFP1V0bLJ/cwAJCFcw8TIXqbWqLDpb5lSRzam9wiw5WlslFxQrV0mCTD1bhe5KeqpPoqWpa6y6guFOD5Qh96+aZdVMrDpbs75GNlGcCYzsApDs1uBUH65jsBHLa0/JaGmzswcp+7YcZ2To72f+FwkT15gghyj48Zr41W9dgdFYcx/IaJd/3W/HmbF2D6yJ3tnDWsxBCjYGao3UNrlfZDqrtTGDPIVQ1rTCmnKwvp9HVDVYLs1f8/jZzlzuo7QerdQHllga3+mCp6/PbcjK+scHtPViFrVKnA3lrWw5io/SOPnaNfM2PV3eR2w2cAnJjg8sefv4zGMXuOVXtPVitm+dxe4Pbe7Aqmmfmcby9we09WMfMPEZG4aVplMr+Jzca4qIGy4tisl8QFfc0dpcr2BSqbDpYsnm+70fldLfxXy5tsE0HS52Vc/uo3rGmJYxLG2zTwTpGqDqJl6Zpzm9BEMirc838KlHOb7A6Vagt2Wtx5RdN6y7OyrIpVNl0sJwz6G7jv1zaYJsO1iGz+7Llucbrbl2Bixps2cHKIVSdxEvTtMKaHyZ/rs5v8PG3SdnDDd/lMnaHqkNrD1Zhm3N0t/FfLm2wTQdLKpyvY2ZlAen8Btt3sLIIVSe5B+Mr41lJzV50DK78lnVjg+PfxG+1tg514mC1SHsPVvYrpRUtv73B7T1YOB+hCgAAoAYU/wQAAKgBoQoAAKAGhCoAAIAaEKoAAABqQKgCAACoAaEKAACgBoQqAACAGhCqAAAAakCoAgAAqAGhCgAAoAaEKgAAgBoQqgAAAGpAqAIAAKgBoQoAAKAGhCoAAIAaEKoAAABqQKgCAACoAaEKAACgBoQqAACAGhCqAAAAakCoAgAAqAGhCgAAoAaEKgAAgBoQqgAAAGpAqAIAAKgBoQoAAKAGhCoAAIAaEKoAAABqQKgCAADtI4RwXdd13TAMdbflH+7hcNDdBgAAgHPFcTwejx3HiaLIcRz5bxPyDKEKAAC0hkpUKsAcb9GF4T8AANAacrAvCAK1RQiRvUkjQhUAAGiN7XbrZIKUJDPWYrHQ0iSFUAUAANpB9UXlQlXuf3UhVAEAgHaI47j6DnpHAAlVAACgHeTY3zFDeqr+0N0AAACAy5SlqJNdWXdFqAIAAC1T1mWlF6EKAAC0zHFJKtd1Hd3jgMypAgAA7eD7vu4mVCFUAQCAdlNTqeipAgAAOK2sYgKhCgAA4AIqM+Wu8pP/q31wkFAFAABaQ65Ik+uykhcDsvYfAADAucIw9H1/u93mlqwJgkB7CVD3+KJEAAAAkwkhsqWqgiDQ3k3lEKoAAEBLyalU2juoFEIVAABADZhTBQAAUANCFQAAQA0IVQAAADUgVAEAANSAUAUAAFADQhUAAEANCFUAAAA1IFQBAADUgFAFAABQA0IVAABADQhVAAAANSBUAQAA1IBQBQAAUANCFQAAQA0IVQAAADUgVAEAANTg/wMEw/6Zohj9RgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "myc.Clear()\n",
    "frp2.Print(\"V\")\n",
    "gr=R.TGraph()\n",
    "smin, smax=0, 4\n",
    "fr2=R.TFitResult(frp2.Get())\n",
    "fr2.Print()\n",
    "fr2.Scan(0, gr, smin, smax)\n",
    "gr.Draw()\n",
    "gr.GetXaxis().SetTitle(\"#theta\")\n",
    "gr.GetYaxis().SetTitle(\"-2log(#lambda)\")\n",
    "myc.Draw()\n",
    "\n",
    "myc.SaveAs(\"test_fitting_3.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "colored-portugal",
   "metadata": {},
   "source": [
    "Draw all canvases "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "sharp-transcription",
   "metadata": {
    "collapsed": false,
    "execution": {
     "iopub.execute_input": "2021-04-15T22:48:39.455029Z",
     "iopub.status.busy": "2021-04-15T22:48:39.454661Z",
     "iopub.status.idle": "2021-04-15T22:48:39.602507Z",
     "shell.execute_reply": "2021-04-15T22:48:39.602793Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAI8CAIAAAD0vjrdAAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dS3LiyP72cemNnhIyc3CwEyWzU70Qm8sCuuyJxMg+CwAcLOTfZyaxEYejrDloBbyD7MpQC0nchDKV+n4GHdXi4hQC9JCZ+qV7OBwcAAAA3Ob/6W4AAACADQhVAAAANSBUAQAA1IBQBQAAUANCFQAAQA0IVQAAADUgVAEAANSAUAUAAFADQhUAAEANCFUAAAA1IFQBAADUgFAFAABQA0IVAABADQhVAAAANSBUAQAA1IBQBQAAUANCFQAAQA0IVQAAADX4Q3cDmhPHcRzHYRje8nAhhBCi8A4Vz1zxKAAAYAf3cDjobkNDXNd1HOeK/RVCbLfb7JYgCI4jlHz+QoX3BwAANunK8F8cx9c90HVdlah835f/WCwWuZB09fMDAAA7dCJUhWE4Ho+veKAaszscDofDIY7jw+Ego9ViscjeU4Yq3/cPReimAgDAejaHqjiOXdd1XTcXgM5/uOyjiqIotz33DwAAgA5NVL+UykzHc8yjKMolKvm/zEYHAKCzujJRPY5jOQJ4/v7KiednzjGXd46iiFwFAEA32Tz8VwsVksIwDMOweshP1lxwMxgiBACgI+ipKqU6n8IwPFlSoaKegu/7RCsAAKzXXKiSxTOdf08/aqwq5tWhSpEX/WXLK2RnrKurC7MjgKrAVdkYYkUUAwAAV9PTZ1RYAqBGURSp8k4VgiC4dzMu2t/sFX+5QglBEMjtURSpO/u+X1hPQe174V85vz33uCcNaFFTaUCLmqq9AS1qqvYGtKipNOBOTa3XHXuqwjDM1TKQCUN15MienuzImu/7YRjeo+/qlp6q44fIm84Z11N/t3AOu+ue+/rf4540oEVNpQEtaqr2BrSoqdob0KKm0oA7NbVedympkB0OC4LgnDE+OQd8u92Ox2OjJiGpfqncxsVikZtoVSibILkwEAAAi9UfquREokv7nNSUI9m/5bquIeUJTGgDAAAw311KKsjamNfFkTAMD4dDEASGdFad04xzqi0QzgAAsBslFUqpOWEn51SpzrnjXFXxJE6rRp2tbECLmkoDWtRU7Q1oUVO1N6BFTaUBrZhTpbn4pyHdUc7v4g7ZwgfZEcnsPY+3y39st9vcPeM4lomqcGIWAACwieaeKtWRc+8ZVCd7qgov6MtewCiDUfZ/sxFKlaRS91TrMVf/Uft6Cq3cKcfS/bJypxxL98vKnXIs3S8rd8pp1X7paqoRL5AMNHdtyXWhyikqDOGUFPPM5iql+krGFr1Bz2flTjmW7peVO+VYul9W7pRj6X5ZuVNOq/aLUGX0asQqQp0sD3H+PZ1TFdVNODRXaNGn7iJW7peVO+VYul9W7pRj6X5ZuVOOYftVfQJ1NJ1D9b9AqoNHe0uaZ9QbtC5W7pRj6X5ZuVOOpftl5U45lu6XlTvltGq/rCr+eVL8mxovyy4LAwAA0DoNRbnjFOU4ju/74rcG2mCgFqX+81m5U46l+2XlTjmW7peVO+VYul9W7pTTqv2ysKdKpqjcLO8zV60BAABol7tEuVwlguMOKnPKU+nVotQPAEBb2FP8UwixWCx83z8cDofDQa7fIleekXfYbreu6wohiFYAAMAaNUe5k+WgcmWfLl132TL0VAEAUDt76lSdU8nzuKKmyUWq7odQBQBA7ewZ/pPDfCcrZB4OhyiKfN+XW2T/Vge55XQ3DQAAQ1WcPTWeQO8V5eRUqnNmTcVxHIbhdrvtYJ8NPVUAANTOnuE/nI9QBQBA7ewZ/gMAAOggQhUAAEANCFUAAAA1qDlUyaqetzxDHMeu64ZhWE+DAAAAGlFzqAqC4OqC6WEYuq47Ho9lRdB6GwYAAHBXd5keL4RQi/35vi9+O75n/Ju6f6eqgHL1HwAAtbOwpEI2Wp0jCIKudVARqgAAqJ2FoUpStT0Lb63ux7JeddVX8hYAAIVOlk23M1TlqIlW3UxROfRUAQBQO2t7qlCBUAUAQO2oqA4AANBif9z7D5wsrMA4IAAAsMB9QxUXABZKkmS1Wq3Xa8dx+v3+ZDKZz+eDwUB3uwAAwPXMGv5bLBau615aNbRdkiR5fn7e7Xb7/d5xnP1+n6bp09NTkiS6mwYAAK5331AVx3EQBPLfvu9HUXQ4HA6HQxRFURT5vi+3yy3qnuPx+K6t0mu1Wj0+Pn58fKgt6/V6NBotl0uNrQIAADe67/T4OI5lQir7K2EYLhYL3/dV75SsPGFxXfV+vy/7qHI8zyvcDgAALmLn1X8yUUVRVHaHMAx9399utypUyf4qW0cA0zQtS05pmrr/1nDbAABoC7eSrlY1Maequs9J3porCqo9VMlC8Dc+/HgvPM97eHgofIjneYd/u/qvAwBgt0MlXa26e0mFk7Tnp0Kyj+2KXJW94HGxWDhHlzROJpM0TeWlf8psNuv1eje0FwAAaNZET1V1NJERRPVmyYyld0LV1TnPdV2VqOQ0fMdxFotF9hWYz+dfX1/T6VRtmc1mn5+f8/n8uj8KAACMUN2BdiN1QV8QBMe3Zuda5bao6wSbp9p86YujUtTJjd/f3z9//vQ8z3Ec+d9fv37V0HQAAHA43DvelLn79PjscJhMGEKIOI6zRUHltX7ySkB5t+bHBNWFilnnvzjq4cfXLVZczygvT2AFQAAAaqTrxNrEMjUqLckglauxnksbWhLV7XIT7bOiKGrjHgEAgIs0F+XkvCLZR+X7vvitmb9+kZPltY7J7qhLl9lRUZrOKgAA6mJtT5XShUX9VEaUO2tsagQAALVreu2/OI7tHgsTQriuu1gsFovFeDx2XbcLaRIAYL3c7B0cayhUyajhuu54PJZRQ7IscIzHY/me830/W1Khor9K1X6tLg5rTrlYAECnJEny+vra7/eFEP1+/+XlJUmSZv70pWdG7afIu4eqOI7dTPWmnOrA0RbZvje5PrTskDv8LtCQXYcnR12H6Vx4/WcTOwYA6LYkSZ6fn3e7nVxjbb/fp2n69PTUTK66urRBA20rdPdQpeoUBEGQrT4VRZHsy9lut23PVdn258KT6oo72ScnayvU2i4AAG6yWq0eHx8/Pj7UlvV6PRqNlsulxlYZ677T41WRqrK/oqotGNX1Uu/VfxX7mLs8gWsAAQBG6ff7so8qx/O8wu2G0HU+vW9PlUxU2crpOWEYyv4qOyZXtb3LDQAAJU3TsuRUcVOXNTFRvTpq2BREbrywkRFAAIA5PM97eHi49KYua7qkgq3khHQ5zJej1t5puk0AANxmMplMp9PcxtlsNplMtLTHcE2EquqhvcIgYjhZ1TO7X2UT0s+fqA4AgGnm8/nX11c2V81ms8/Pz/l8rrFV5rr6esVzyP4bx3Gy1/1lqf6buzbjUmoSWNkd5K2yeoKidtZxnCAIcv9b9jxnbgQAQIvv72/HcTzPk//966+/fv36pbtRJ+g6k959enx2klAQBHIGlZx7pPqoLl0y795OXv0nd+p47Wd1oV9Wxd4VXp7ANYAAAKPIE1Mcx22ZBq3rTNrEX62efG1aorqd2p2Ta/9VvzJEKwCAdmb+1D95XZe1ocpxnDiOwzDM1lX3fT83LamDyt6pZr6DAQAd1MZTks09VShDqAIAmKyl5yM7i3/iOgcKVgEA0DaEKgAAgBr8Ue/TXd2/0sbeRQAALNbSsT+N6KkyFCOAAAC0S809VRVrJwMAgLagm+oKvGQ6nXzL8p4GAGjR6hOQrsbX3FOFS1WM8bX33QwAaDXzE5WZM2TuNadKFrN3Xdd1XSFEbjmX7N3kfe7UDPNVLCGku2kAABiqeg0+Xa26S6gSQozHY1U/fbvdjsfjjhdPvw7T1QEADTO/m8pY9Yeq7HI0QRAEQSD/vVgsyFUAAMBW9adR2bPi+352yC8Mw8Vi4ThOFEXZNYbjOB6Px05X5w+d82uAXwwAgCZZcN6xZJkaFaRyk6jCMPR933EcGaFwPkYAAQCNsSBRaXSXUCXzU+FNjuMwCHgdNagKAAAM1GhFdTm/So4D4kxJkry8vMiLKPv9/svLS5IkuhsFALAQ3VQ3qjlUyflSZX0qqo+KzqozJUny/Py82+3k/+73+zRNn56eyFUAAJjmLqEq+48cuY7NYrEoq1yFrNVq9fj4+PHxobas1+vRaLRcLjW2CgBgH7qpblf/K6gu9HMcJwgCIUQuYKk7ZEcDu3kgq2egHw6Hfr+/3++Pb/I8r3A7AADXaVeoOnkJl5Z9ucsrKITIjgAe/4ls8Cq7TxdUv4PTNH14eCi7dbfbVdwKAMBF2hWqqllSUkGK4ziKosJrAKUwDKvvAMdxPM8ri00VNwEAcCmbEpVG91pQuWK9v9wd4jhmflWZyWSSpul6vc5unM1mvV5PV5MAAEAhkqlOJ38ZJEny9PQ0Go1UrprNZp+fn5vNZjgcNtJGAIC1ttut7/v2dVNZNfxXIQxD+qXONxgMNpuN53me58ktvV6PRAUAuEWSJK+vr/1+X1ZAlFt0N8oGTUc513VzywJ22UVROo5jIYR9vycAAE2SFRCz9Xqm0+nX19dmsxkMBnrbVpeu9FThamWlvwAAOB8VEO+HniqdrovSdFYBAK7WhQqI9FR1lFtOd9MAALZJ07QsOVXcZKCKs6fGE+i9SioYSNZuuGjZwepyD7mnqnjm47LyCn1OAIDGyDKHZT1VLaqAWH321JWrOjT8J1/ii/Y3Vxo+J/dUFYcwCILCyHV1/yQjgACA67y8vJRVQHx/f9fVqnp1ZfjvcDhoSVT3/qPMEgMAtMJ8Pv/6+ppOp2qLrIA4n881tsoOnRj+O15q8EyymyqKopNX3slQxRx8AIDhZAXE4XDoeV6app7nUQGxLjaHqjiOx+Px7c9jYC2Dw+HACCAA4DqyHtV+v5cVEHU3xx5c/VcP2UHFWxMAYD71s5zTVr1sDlVCiMNvURRd+vCLBvLkQCHvTgAAOsvmUHUjNU0qDMNs9YuK5KRWklHuN8VKjgDe6ckBALZi9sj9EKpO2G63uUnu2+22LC0tFotcCYbxeEz3FQAAXUCoKqViUxAEx8OI2Snw2YAVRZG6s+/7juNst9uKuqDVNWENLBcLAGivdnVTte4USagqJYTwfT9Xt1MIoXJVdrvv+77vHw6HbL9UHMcyV1UUdDhcy2EEEABgtVtOkVq0KbHeQpVXqGV/ZaX1c6pSqb9bWOzq9l8M7frNAQDQqDunDF17Sk/VNWQ8qljBJndP524l1+msAgDAEIQqAADs151uKo1srqh+iziOZd9SxRxzOV9K3UcIUXGhH9cAAgBgN0JVMSGEnAtVGJXkxPPs0N52u1U5LEtlsvuFqgNL1gAAKnGaaAbDf/+Q4Sl3QZ/z79IJ6p7yH+rO8h/HpRPiOJbxKwiCe7QZAACYoyvR9eTVf3K6d/aCvux6zL7vCyFkj5TcclxqIXuT87v7Sm6p+KO1vP78BAEAVOjaaULX/nblVb4iVDn/zlVZuUQlZXOVUl12ofrCvYsOTdc+MACAM1l5gjh55TuhylDxb+K3ijtnJ1GdnEdV4xvdys8MAOB2HTxB0FPVRYQqAMBddfPsQPFP3IQqoAAA6EWoAgDATt3sptKIUAUAAFADQpU9GAEEACh0UzWPUGWhc1Z6BgAA9SJU2SNJkpeXF9d1hRD9fv/l5SVJEt2NAgA0Sv6upptKC0KVJZIkeX5+3u128n/3+32apk9PT+QqAOiCJEleX1/7/b78XS236G5U5xCqNHPLXfQ8q9Xq8fHx4+NDbVmv16PRaLlc1t1kAIBZ1O/q/X7vOM5+v59Op3b/rq44e2qcXkz3oE41ds/2+335WcrxPK9wOwDAGq+vr7vdLvu72nGc6XTqed7b25uuVmlERfUuquuop2n68PBQdutut6u4FQDQdvyuzqGiOq7neV5ZbKq4CQBggTRNy5JTxU24hz90NwD1mEwmaZqu1+vsxtls1uv1dDUJANAA+eO5rKeK39VNoqfKEvP5/Ovrazqdqi2z2ezz83M+n2tsFQCgAZPJJPv9L81ms8lkoqU9ncWcKp3qHfRNkmS5XK7X6zRNPc9L0/TXr1/D4bCu5wcAmClJkqenp9FopMYr5O/qzWbTzbMAE9W76E5HPY5jIQSV3wCgO5IkGQ6H8he153mTyWQ+n3czUTmEqm6691EnVwFAd8jvfPm7WndbNCNUdVF1gbLbDw2hCgA6omtf+CcrfGp5Nbj6T7NOfQYAAKhF9dlTV1F1rv6z2eFw0FitHwDQjK51UxmLUAUAAFADQpXl6KwCALvRTWUOQhUAAEANCFX2o7MKAIAGEKoAAGgrxv6MQqgCAACoAaGqExgBBAD70E1lGop/FovjOI7jMAwvfUjZrWVPVZF1+KgAAFDIzJ4CQm4xebQuenGEENvttuzWwqdq+EcGv2kAwBp8pVfQ9eIw/FegosMJAACgEKEqLwzD8Xh8xQNlN1UURYcidTfzGsysAgA70E1lJo7KP+I4Ps5SF704V4wYNv+p4HMIABbgy7waw38AAAAtRqj6hxBCDdVFUXTpw9syDYsRQABoO7qpjEWoqocMVb7vh2HoZgghNLcMAAA0glBVp+12u1gscltc1zWqH4vOKgBoL7qpTEaoqoeKTUEQHA8jVlxO6F6rgZ0CAECj1p0iCVX1EEL4vh8EQbZyuhBC5aqyiuqF9RfOcUtr6awCgHaRVXu61k2l5RR5C5apqUdZZpJha7vdGjUCCABohSRJVqvVer3e7/cPDw9yy2Aw0N0uFKOn6u7kXPWKFWwAADiWJMnz8/Nut9vv947j7Pf76XT69PSUJInupqEYoaqj1AggaQ8AzLRarR4fHz8+PtSW9Xo9Go2Wy6XGVqECw381iONYju6VDQI6juP7fmPtOYf8odPv92WX8mQymc/ndCkDgDnkqN/xRs/z3t7etDQJ1eipqoEQYrFYLBaLwolTssiCUQWrZJfyZDJRXcppmtKlDADmSNP0OFGdvAl6EaquIYQQQmT7pWRH1HHpBJWlKjqxmkeXMgAYzvM8OTP9opugF8N/1ziehxSGoUxUruv6vi+EiONY3S0IgqabWIkuZQAw32QySdN0vV5nN85ms16vp6tJqEZPVT2yJalkXfVsojKqm4ouZQBohfl8/vX1ld0ym80+Pz/n87mmFuGEbpURa0D8m/it4s7VFTjvd2jk/PTj7Z7nEaoAwBxJkgyHQ8/z0jT1PE9eVDQcDnW3S7+TJay1xBtClU66auO+vLyUdSm/v7833x4AQCF1mpC/1XU3pzV0nV4JVTrpOupJkjw9PY1GI5WrZJfyZrPhBxAAmKNr69LUhVDVRRo/LUmSLJfL9Xotu5TTNP316xeJCgDMQaK6GqGqi0z4wMguZRNaAgDI4pv5aoSqLjLqA2NUYwCg4/hOvoWuV4+SCgAAADUgVAEAYBa6qVqKUIV/HA6Hk2U/AABAGZap0awix/AzBQA6iG6qc5jZC0Co0syoT47srDKqSQAAHKs+VemKXAz/AQBgCn7ZthqhCv/CzCoAAK5DqAIAwAh0U7UdoQoAAKAGhCrkMQIIAM2jm8oChCoAAIAaEKpQgM4qAGgS3VR2IFQBAADUgOKfmlFRHQA6jm6qK5g5nEJPlWaHctobZuZbFgDssN1udTehxSrOnhpPoIQqAACakyTJ6+trv98XQvT7fdd1v7+/dTcK9SBUoRSdVQBQryRJnp+fd7vdfr93HEf+9+npKUkS3U1DDQhVAAA0ZLVaPT4+fnx8ZDeORqPlcqmrSagRk+N0asXkxFY0EgBaod/vy96pHM/zCrfjOrrOXPRUAQDQhDRNy5JTxU1oEUIVAABN8Dzv4eHh0pvQItSpwglyuvrhcNhut77v624OALTYZDJJ03S9Xmc3zmazXq+nq0moET1VOEFek6Ku/n15eeEqFQC4znw+//r6ym6ZzWafn5/z+VxTi1CnDs1BjuM4juMwDG9/EiGEECJ3U8UzF97fOVUQ1oRDI6/+zV6rMp1Ov76+NpvNYDDQ2zYAaKMkSYbDoed5aZp6njeZTObz+XA41N2uljlZ8UfLObRDoUoegBv3Vz5JEATHEariABfe32nDhXWvr6+73S539e90OvU87+3tTVerAKC91De//Imuuzl20nV67cqcqjiOb3+Sind/Lc9voPV6fXxBynq9JlQBwI1IVPbpxJyqMAzH4/GNTxLHccU6TTJU+b5fuAjRjWOOunD1LwDUy/wBCtzC5p6qOI5vz1JKjU/VFvIS37I6dVz9CwBAVid6qm6nplKV3UH2VNnXlzuZTKbTaW7jbDabTCZa2gMA7UU3lfVsDlVCCDUAF0XR1c8jB+98368YxZMjg/aFKnn1bzZXcfUvAACFbA5VtYjjeLFYOOdNRZeXcrgZbZ/APhgMNpuN53me5zmO43nearXabDZc/QsAF6GbqgsIVSfIqVRndnQtFovcZPbxeNz27qvBYPD29rbf76MokvOrSFQAABwjVFWRecj3/epglO2OiqJIjTnKRV22223FuKF7rXr28BLyRZCr1jT/1wGgveimuk6LTpESoapUGIay2+mcITzf92U9hWz8iuNY5io5gFiosATDOW7fQQBAA0hUV2vdKZJQVUxNpTpn4E8IIZevOb5J9VG1fXJVFp1VAAAcs7lO1S2yYeg4D6k1BE9W9VQdVyxHAAAdRDdVpxCqTigcudtut6qGQjejkuys4psCAACFUFUsDMPCXih1cZ+cLCUTlbxndcDqZvYCgC7jx2fXEKqKyWlSx9vd36XVs5FLLgtYOFCo7mZfqKKzCgCALCaq/0P2M1238rF81HHpBDXbvWJ9GwCAlfjZ2UH0VP0jV7TzIkII3/e32+1isVgsFjJCye4reYfrspr56KwCAEAhVNVDXtwnU1R2brvv+zYVUwAAnIMfnN3EUa9ZdhLVyXlU1dWe2nJo+O4AAGm73crLmPhivLeT5RK1vP4cdZ3s+NTZsRcAcLUkSVar1Xq93u/3Dw8P+/3++/t7MBjobld36ToxMVEdt6LAOoAuS5Lk+fl5t9vJJeflf5+enpIk0d00NI1QBQDA9Var1ePj48fHR3bjaDRaLpe6mgRdGLjRyaaBM5v2BQDO1+/3Ze9Ujud5hdvRAIb/AABomTRNy5JTxU2wFaEK9WBmFYAO8jzv4eHh0ptgK+pUAQBwvclkkqbper3ObpzNZr1eT1eToAs9VagNnVUAOmg+n399fWW3zGazz8/P+XyuqUXQhlClmVtOd9MAAKcNBoPNZuM4jud58r+9Xm+z2QyHQ91Ns1nF2VPjCZQrtnSy8oo5uVOqrDAAWE99mcsly3Q3B1z9ByvIYnf9fl8I0e/3X15eKH8HoDtIVB1HqEJtZFnhyWSiygqnaUpZYQB2s3LMAdchVKE2x2WF1+s1ZYUBAB1BvtbJst83lBUG0DWWfY1bgzlVaDfKCgMAOo5QhXpQVhhA19BNhRwqqqM2lBUGAHQZPVWojSwrPJ1OsxspKwzASnRT4RihSjMDC8JeTZYV9jxPlRV2HIeywgDsQ6LSjorqyLP4Y6nKClu8jwA6i282w+k6QLwtdOrCx7IL+wigU/haMx8lFWCnw+HQxqFMAAAuRagCAOBcdFOhAqEKd0dnFQCgCwhVAACchW4qVCNUoQl0VgFor+12q7sJaAdCFQAABZIkeX197ff7Qoh+v++67vf3t+5GwWiEKs0MrF12J3RWAWiRJEmen593u51cD17+9+npKUkS3U2D45ha/JNQpdmhnO6mAUB3rVarx8fHj4+P7MbRaLRcLnU1CVkVZ0+NJ1Dm3BWL4ziO4zAMb38SIYSsLX6sg3MeO7jLANqo3+/L3qkcz/MKt8MoVFQ3i+w8vPHFkU8SBEFZOOtgwujgLgNonTRNHx4eym7d7XYVt8IEVFQ3SBzHtz9JWe9UxzGzCoD5PM8ri00VNwF/6G6AccIwXCwWNz5JHMdcggsA7TWZTNI0Xa/X2Y2z2azX6+lqEsxHT9U/4jiWlwzcnqgcxxmPx7c/ia3orAJgvvl8/vX1ld0ym80+Pz/n87mmFqEFCFX1U1OpdDcEAHClwWDw999///z50/M8x3E8z+v1epvNZjgc6m4azEWo+ocQQl2KGUXR1c8j56T7vn/jlYN2o7MKQCu8vb3t9/soivb7/fv7O4kK1QhVdYrjWI4e1jLVHQCgS/byMS48wpmYqF4nOZXqoo6uqzts2l6YQHZWtX0vAAD307oxDXqqaiN/yvi+f9FvmuqasAaWi70HrpQEYBR+8hmidadIQlU9wjCUyYCBv/N9f3+7rqsWK315eWFRLQDakahwNUJVDdRUqltmuHeNXKzU+b1M6X6/T9OUxUoBAO3FnKoaqAv95GJ/uVvVGoJcD5glFyvNblmv19PpdLlcvr296WoVgI6jmwq34N1TII5jOeX8zBdHCHHOrKAoinLTrbr86WWxUgCm6fJ3smV0HUp6qmoQhmFhL5RKWr7vO1yUm5GmaVlykjextBYAoHVI5QUu7akqo0qrlw38dflXET1VAIzS5S9k+9BT1Sayz0kIwTSpq7FYKQDAMlz9d43tdrvdbmupnuCWu/3JTSYXK51Op9mNLFYKQAu6qVqn4uyp8QRKqNLMwNplzRgMBpvNxvM8tVip4zgsVgqgeSSqNjKz+CfvJJ34JEtxHMsRVV4QAM3jm8c+uo4pPVXQj+siATRMXZ1NokKNCFUwyOFwsH4yGQCNkiR5fX1Vq2O5rvv9/a27UbAHoQoA0AlydazdbqdWx3Ich9WxUCNCFcxCZxWAO5GrY318fGQ3jkaj5XKpq0mwDGPJOjGWX4iXBcA9UHO4O5ioDvyDzioAtTu5OlbD7YGVCFUAAPt5nle2qGjFTcBFWKZGs4oumS4PgcnOqi6/AgBqx+pYNjFzQIOeKs0MLAhrCAYBAdRLro6V3TKbzVgdq6XMrKhOqAIAdMJgMPj7779//vypVsfq9Xqsjj6ZutEAAB+vSURBVIUaMcKiEyNcJ/ESAaiR+kpRq2PBSrrOHZyxdCIxnMRLBKAufJ90ByUVgALMrAJQCxIVGkCoAgAAqAGhCqajswrAjeimQjMIVWiT7XaruwkAWoZEhcYQqjRzy+lumkG+v79d1+33+0KIfr//8vLCqvIA0GUVZ0+NJ1BClWYG1i4zTZIkz8/PjuPIxbn2+32apk9PT+QqACfRTWUrin8C11itVo+Pj9kt6/V6NBotl0tdTQIA4BgRXid+Qp2j3+8XLiDveR4LywOowHdsZ1GnCiiQpmlZcqq4CQBIVGgeoQpG8zzv4eHh0psAAGjeH7obAJwwmUzSNF2v19mNs9ms1+vpahIAw9FNBS3oqYLp5vP519fXdDrNbvz8/JzP57qaBMBAqo4diQq6EKpgusFgsNlsPM/zPM9xHM/z/vrrr//973/D4VB30wDolyTJ6+trto6d7hahu4jzOvFz6lJxHAshHF46AI7j/K5j9/j4+PHxoTb++PFjs9kMBgONDYNeXP3XUQYWhDWZTFQOCwICcBzndx27bKJyHIc6dl1gZkX1Dv3cj+M4juMwDK9+rDyjq/N6TsUzCyEKH0V3y414AYGOo44dCuk6O3TonCSj66X7K4Q4XsQ3iqLjkFQRjYMgKIxcZIIb8QICXZamaUVdld1uR9WVzmL4777iOL7iUa7rqkTl+77v+/Lf4/E494TXPT9uxCAg0GXUsYNpOlGnKgzDxWJxxaPkP3L9UvIsPh6PsylYhirf90lXDZO5iv4qoJuoYwej2NxTFcexnLB2RaJyHEc+yvf93EhfFEW1NA8AcCNZxy67ZTabUccOutgcqmpxPBdKZaxsp5T8d9kcdtwVg4BAZw0Gg7///vvnz5+qjl2v19tsNtSxgxZdGTeJ43g8HjuXTFSXceo4VKmnyg4LypN64QT2Coxb1YVXEugs9fFX12gDXP13X1eEqjLqesDsU8lQFQRBHMfZqwWrYxZRoEa8mEAH8cFHIV1vjE5MVL9d/JsKTIUzq44nb43HY2avN4MZ60DX8JGHaZhTdRZ5/WA2UWX7n7KZKYqiw2+yBMN2u62oC1pdE9bAcrEAADSjdadIQtVZwjAMgiBbpyqXk2QVq8PhkAtb8iEVlx8ernWP3Ww7ZqwD3UE3VRe07hTZlTdljXOqVNWrc6alF85qV/hSuAf1qm63W5WDAdiEL09U0/UOoafqYqqP6pxlBAvrL+DeXl9f+/2+EKLf77+8vCRJortFAAD7EaqKqZnphbeqyVKNtglnSJLkzz//fHt7k2up7vf7NE2fnp7IVYA16KaCsQhVpcbjsRy5K6OGlsIwDMOwui+K6inNWK1Wj4+P2S3r9Xo0Gi2XS11NAlAjEhVM1pV35xVzqtzyep6qKpUcAZSVqwpLJ6gJWIV/l2+H2vX7fdlHleN5XuF2AO3C1ybOwZwqzYQQQojjaVLHnVXqPipsyS3HpRPiOJaJKgiCmpuLImmaliWnipsAtAWJCobryhv0ZE+V7HzK9japhzi/l1XOFv9U3VSSKrPu/I5Q2TtX/NGOvP6NoacKsIy6jJcvTJxP17uFiuqlhBBRFMlctd1us9PSc4nK+b3mlLxPtirVyXLqFXWV+Pq4wmQySdN0vV5nN85ms16vp6tJAK6QJMlqtVqv1/v9/uHhYTKZ6G4RzGJmVUKC/2nZywDlKGHFnbODgycnp/PDq3ZJkjw9PY1Go2yu+s9//sOq9UCLJEny/Pz8+Pj48fGhNv748WOz2QwGA40NQ1voOr1yUteJUHUPSZIsl8v1ep2mqed5k8nkv//9L68z0CKvr6+73S6bqBzHmU6nnue9vb3pahVahFDVRYSqu5Jjsg6vM9A2TI7EjQhVXcTJvjG81EBbpGn68PBQdutut6u4FZAoqQDcEWstA23heV5ZbKq4CTABV/8BAMzCZbxoKXqq0BV0VgFtMZ/Pv76+sltms9nn5+d8PtfUIuAshCp0CLkKaIXBYPD333///PnT8zzHcTzP6/V6FEaB+Zi9qxOzp7XgZQcMl/2Qqst4gfNx9V8XVfeacGjuhFAFmIxPKM5xctiBZWq6iO+O5slBQF55wEB8NnGm6veJrpkezKlCFzG5CgBQO0IVAMAIdFOh7QhV6Cg6qwCjkKhgAUIVuiubq7bbrd7GAF1GooIdCFXoutfX136/L4To9/svLy9JkuhuEQCglQhV6K4kSf7888+3tze57v1+v0/T9OnpiVwFNIluKliDUIXuWq1Wj4+P2S3r9Xo0Gi2XS11NArogO9pOooJNeDfrRPFPvfr9vuyjyvE8r3A7gFskSbJardbr9X6/f3h4mEwm7+/vfNHhOmYW/6SnSrNDOd1Ns1yapmXJqeImANdJkuT5+Xm326nR9vf39x8/fjDajutUnD01nkAJVegoz/MeHh4uvQnAdeRo+8fHR3Yjo+2wDMN/OjGZQK+Xl5c0TdfrdXbjbDbr9Xrv7++6WgVYidF2NEnX6ZWeKnTXfD7/+vqaTqfZjZ+fn/P5XFeTACsx2o6OIFShuwaDwWaz8TzP8zzHcTzP++uvv/73v/8Nh0PdTQOswmg7OoLhJ50Y/jNHHMdCCPlvjgtQO0bb0SRdX+OcPHTi5G0mjgtQuyRJhsPhdDpVuWo2m31+fm42G/qGUTvmVAGmYK1loHaDwcBxnOxoe6/XI1HBMvwi14keEZNxdIAaZT9Q2dF24B4Y/jNLHMdxHIdhePVj5VdG9RcHFdUNR64CasFHCbUzs6I6b/Ri8mhd+uIIIbJrWklRFJVFK75oDMcBAm7H5wjNY06VQeI4vuJRruuqROX7vu/78t/j8fi6J4R2TK4CbkSiQqf8obsBxgnDcLFYXPEo+Y9cv5Q8JY/HY75WWkrmKg4fcAU+O+ga3vH/iON4PB7nNp7/4sjw5Pt+rlNKPW3hU/GN0xYcKeAKfHCgC8N/Njie2K56rRgBtMbxtDkAx0hU6CBC1T+EEIffoii69OFBEARBcDwhnSxlBzkI+Pr62u/3hRD9fv/l5SVJEt3tAgxFokI38b4vUD1mdxF1PSDDf62WJMnz8/P//d//qS3T6fTr62uz2ciShkDHbbdbdXUO32zQjuE/e8gCV0IIdT3gFV1fMMpqtXp8fMxuWa/Xo9FouVzqahJggiRJcj24JCp0Ge/+Ajf2VOWqVVXXqbqqgY5DadBm9fv9/X5/vN3zvMLtQBfIHtzHx8ePjw+18cePH/Tgoi6tO0vSU1W/MAyDIMjWqaqozH64VkM7A8dJ07QsOVXcBFhP9uBmE5XjOPTgokatO0XSU1WgxjlVqupVYX8V/eRtQU8VcIzPBYzFnCo7qT6q65YRhCEmk8l0Os1tnM1mk8lES3sA7ejBBY4RqmoQ/1Z4qxwHpLhRq83n86+vr2yums1mq9VqPp9rbBWgked5Dw8Pl94E2I1laupxcrhQTbFCGw0Gg81ms1wuPc9L09TzvF6v9+vXr8fHRwZw0VmTySRN0/V6nd04m816vZ6uJgF6MaenwBVzquQVCmUTpxzHCYLgeASQOVVtFMexOsocQXRZkiTD4XA6napcNZvNPj8/N5vNcDjU2zZ0HHOq2kQIIYQ4DknHqweq+5RVVUDrZA/l4XC45YpfoNWGw+H397fneZ7nOY4je3BJVOgyhv+ucTxBKooimahc1/V9XwgRx7G6W+EKNrCDzFX0V6Fr1Nv+7e3t7e0t24MLdBahqh5CCJWrttttNnUVDvzBJuQqdM3xG55EBTjMqapd9jJAOUpYcefqkSMOTbuQq9ARvNVhgpNTL7S8S/ls6MR3k2U4oOgC3ucwHxPVgdZj3jqsR6ICKhCqgDrlchVFX9FquTcwiQqoRqgC6pckyevra7/fF0L0+/2Xl5ckSXQ3CjhX4RuYRAWcRKgCavb9/T0cDne7nVz+bL/fp2n69PRErkIrJEny/PycewPLklS6mwaYjl8eOvHLz0qvr6+73e7j4yO7cTqdep739vamq1XAmXgDwwK6Tq+c1HUiVFmp3+/Ln/g5nucVbgeMwhsYFuDqP8AGaZqWnXgqbgIMwRsYuAU9JTpR/NNK/NBHq/EGRiuYWfyTZWo0IznZZzKZpGm6Xq+zG2ezWa/X09Uk4Hy8gdEK1WdPXSUD6anSiTlVVkqS5OnpaTQaqdPSbDZbrVa/fv0aDod62waclCTJcDicTqfZN/Dn5+dms+ENjLZgThVgicFgsNlsPM/zPM9xHM/zer3er1+/Hh8fdTcNOE1WT8i9gUlUwDnoKdGJnirrxXGcXVSbIw7D5d6iuTcw0BaUVOgiTrEdxEGHsXhzwhoM/wGdwKLLMJDruiQq4HaEKqBp5CoYRcYpEhVwO0IVoAG5Crpst9vs/9JBBdSIUAXocZyrcmc7oEZJkry+vvb7fSFEv99/eXlJkoREBdSLUKWZW05303B3MlcVnu10Nw1WSZLk+fl5t9vJkuj7/T5NU1k6QXfTgCtVnD01nkAJVZodyuluGprw/f09HA5zZ7unpydyFWq0Wq0eHx8/Pj7UlvV6PZ1Ol8ulxlYBt6g4e2o8gdL3qxN973h9fd3tdtmzneM40+nU87y3tzddrYJlWM4PXUOdqi4iVIGzHe4tTdOHh4eyW3e7XcWtQEtRpwronDRNy5JTxU3ARTzPK4tNFTcBuMIfuhsAdJc8pZX1VHG2Q10mk0mapmqBZGk2m/V6PV1NAqxEqAJ04myHBry/vzuOM51O1TttNpt9fn5uNhut7QJsw/AfoNN8Pv/6+ppOp2rLbDZbrVbz+Vxjq2ATObnk+/vb8zzP8xzH8Tyv1+ttNpvhcKi7dYBVCFWAToPBYLPZ5M52v379enx81N002EBN1x0MBm9vb/v9Poqi/X7//v5OogJqx9VnOlUXKOPQdE0cx0II9b9cHIob8RaCxU5W+KSkwn3FcRzHcRiGVz/WcRwhRPa0l1XxzGWP4isP1XiH4Gq8edBl1Km6O5lqL93fMAwXi0VuYxAExxGqIjUX3t/hWw9nuO59i07Zbre+72e38N2CjqNO1X3JfqZLCSFUovJ9X31tLRaLXM/Tdc8PnCSXXGDpZRwrWzWSRAXo0olQFYbheDy+9FFxHMtTl+/7h8NBjgAeDgcZrbbbbTZIyX/Lex67bswRUFh6GTllaySTqACNbA5VcRzL1aqPx+/OoZJQrhdK/S9RCU1i6WVkla2R/PPnT42tAjrO5lB1I9lNFQTB8U1yY3YIRk1jb6Zt6KDVajWZTHIn0dFotFwuNbYKuqzX69w63HJjrpAsgCZ1paM4jmM5Anj+/spZLFEUHUclNXtdPVvFnav/REdef9yOpZehsEYyUI2J6saR06EKQ5KaQXW8XQjhZjCBHbVg6WVksUYyYCZC1cXCMJQDf8dzqhaLRe6yrPF4zJggbsdJFDmTyeR442w2K9wOoBmEqsuoIgu+76u0lO2OiqJIXfSnrhOsmNLuXuueewkTTSaT7BKB2e3NNwbavb+///jxI7dq5OfnJ6tGwiatO0USqs4VhqHrumr2em5cT1axyg0XxnEsc1XF5YeFJRjOcafdhLEKl17+z3/+89///ldjq9A8ec44HA7Hq0ayRjIs07pTZFcmSl8xUf34sdJFs9HVYwsfxUR1XCRJkuVyuV6v0zT1PG8ymcznc1mayKHqunWO66Q7JV8auVUjAeg6vdJTdUK2cGgQBIeSqetlCocIgesMBoO3t7f9fh9F0X6/f39/l90S8pcZg8J2qKiTXnaeIFEBhiBUVVGlE2ScotonDFF4Ei3LVaxp0yKFddKfnp5knKIzEjDcH7obYK44jmWiOjneJ8OWEKLibvyURANUrjocDkmSrFar9Xq93+8fHh7kWOFgMNDdRlQprJPuOA510oFWoKeqlIxKQRCczEMyfhX2Y6mNhCo0Qw0FFnZ4sKaN4QrrpDu/oxUAw3VlovTJieoy9AghVAySv/irQ5W8ST15EATZaFW2XWGiOu7k9fX17e0tt3E6nXqed7wdhqBOOlAXXafXrpzUT4YqGaF835czynNX/JVRzyaEUDNX5MqAcRyrLRV/tCOvPxrGmjYtxYEDaqHr9MqcqmKXXqwnL2mWKSpblUqltDIVV2yRt3Cdk2va0OFhrMIDN5vNer1e840BTGbm9c70lNQsO4nq5DwqeqpwJ3R4tJHrut/f309PT6PRSE2iknXSqeoJXIThvy4iVOFOXl5e0jQtnN3MW06jwnqezu/f3PLQlJV4bbqtQJsRqrqIUIU7SZKkrMPj8fGRd13DKspbVFTDp046cDUqqgOozWAwKFsYThZcKJuOQKXQ2pXV80ySpLqkJ4kKaB16SnSipwoNKOvwyA05USn0Tl5fX3e7XWH1KT7+wJ0w/NdFhCpoJydHPz8/Zwt5T6fTr6+vzWZDrrodFw0AzSNUdRGhCiYoHAqkUmgtqOcJaEGo6iJCFUxAV8pd8fICzWOiOgANTlYKbbg9rVM9td913f1+P51Oc9tns9lkMrlnuwBoQKjSzC2nu2noBM/zykagKm5yOn+dYJIkr6+v/X5fCNHv919eXnKLVctP8eFw+P7+/vr6yuYqWd5iPp833mrAHhVnT40nUEKVZodyupuGrphMJoVdKWmaHn83nQwTXVBRJcHJxCn5Ka4ob6F3L4BWqzh7ajyBEqqArpvP54VdKb9+/VJFrWS6qg4T3bFarbIXSzqOs16vR6PRcDgsLD01GAze3t72+30URfv9/v39nUQFWImJ0joxUR2GOGdpFJmrJpNJruRSB68TZO45YDiu/usiQhVMU700SnfCRNkifQ5VEoA24Oo/APpVJKpbrhNsy6z2c2aMyalRhaqn9gOwHqEKwFmqE0PhTe2a1X5yxpiaXvbz50+qJAA4RqgCcK6y6wT/+uuv4yuZTZvVfrK3rGz6+XK5zF7Qdzgcyqb2UyUB6DhCFYBzVYQJFThUuqrIKNV/5fyxwnPueX5v2Xq9Pl72eL1ev7+/5y7oo0oCgGLVlR5wV7z+aJ3v7++fP3+qMPHXX3/JygvHyr5zPM8re+aXlxc5jPjw8PDz58/v7+/b7/nnn39mR+Wm0+mPHz+O7189J2y325W9IFEUld0EQBddp1d6qjQzsCAsUOHMkktpmpY9Q+Gs9vPHCi8aVTynt0x+3Cqmi1VPJquY2g/gfsysqM4l/TpRUgEWK6u/kKXe/6+vr7vd7pwKWOff89I2vLy8pGm6Xq+zt85ms16v9/7+Xv0kAIxCnaouIlTBYudklJM/KI8rYJ1TK+uc36nHBaWSJHl6ehqNRqrNcsYYk6WA1qFOFQCrnHOJnJyFUNGZJNcfzKqolaXuo+Y3XLRWNNPPAdyIUAXgLs7PKBWTlo5ntZ9zT7WxrAZEWUEpFukDcAuGn3Ri+A8dUb36jXPJfKaLZj4xogd0E8N/AKx18hK588tpXlR4kxE9AE0iVAHQ7/z0c2lOYkQPQGMYftKJ4T/g2MmxwivuCaBTKKnQRdUXfnNoAAAodLJyipZz6B/N/0lkkZwAALhU9dlTV1F1QlWxOI7jOA7D8OrHOo4jhGBsAgCAjmD4r5gMuZe+OGEYLhaL3MYgCMrCGXOqAACoHXOqDBLH8Xg8di4MVUKI7XYr/+37vuM42f+VfVc5hCoAAGpHnSpThGEoE9VF4jiWEcr3/cPhIEcAD4eDSleFoQoAAFiDnpJ/qN6prPNfHNVNdfwQOZJY2FlFTxUAALWjp6rdZKIKguD4JrlRDQUCAAArEar+IYRQq7FGUXT1k9TaKAAA0BqUVKhHRTejHPWTk6sAAICtCFX3FYahHPirqKpw3TMzGQsAYDddNTyvRqi6IzV73ff9spFBshEAAIWuPkXqSmPMqbqLMAxd11Wz16mnAACA9eipqlmuNEMURcxeBwCgC+ipqlO2cGgQBIfDgUQFAEBH0FNVG7XwX8VifwAAwFaEqnrEcSwTFeN9AAB0E8N/9ZBdU0EQkKgAAOgmeqquIZOTEEIN86lVaCou9OtI3rJ1QUMr98vKnXIs3S8rd8qxdL+s3CnH3v2qEaHqGrmF/FSQWiwWchCwUOF7saKWBu9dAAAKmVkXlFBVg1vKUJGcAAC4VPXZU1fkoitPJyu7Uq3cKcfS/bJypxxL98vKnXIs3S8rd8pp1X7paioT1QEAAGpAqAIAAKgBoQoAAKAGhKp2OH/O3UWz8+70tHobYOVO3e9p9TbAyp2639Pe4zk7frDu9Frd6Tm1vwLna1FT60WoAgAAqAGhCgAAoAaEKgAAgBpQ/FMzKqoDAHApM6dYEao0IzkBAHApMyuqM/wHAABQg9aUnLeSmb2XAAC0nZZ4Q6gCAACoAcN/AAAANSBUAQAA1IBQBQAAUANCFQAAQA0IVQAAADWg+Kc2cRzHcSyEEELobstZrmuwfFTZrWEY3tyuRsndaUWzr2iqBQdLtr9Fn6yrG8zBal6XD5YShmFbjpceBzTO9/3cUYiiSHejqtzS4OPHtvrt16JmX9HUVh+sKIqO2+z7vu52lbqxwa0+WEEQHLc5CALd7Sp1Y4NbfbCy1I7oboi5eGmalvsCVf82Nlfd2GBrvk0OmbOg7oacdl1T23uwsgHF9/3cjuhuXYHbG9zeg5UNKLl9NzNX3d7g9h6srOybVndbzMVL0yj10comEpPfprc3+PjhLZX9YtXdlhOubmp7D1bh/qpzgIH9Vbc3uO0HKxdHCr9qDHF7g43dtYu0MQg2j5emUYUfzoPBH7nbG9z2T2DhGI3uRhW7vakm712Fim45lS+bb1WFWhps4H6do2IHy75t9KqlwS09WFnm9/4agqv/mqMmKpZNSzRtumLrGoxuqpgCrKbTVtynea1rcI3kflUMh5m2461r8D2EYbjdbn3fL5xbhixCVXMqPnvynbrdbptrzRlub7AFXzdCCPUTpLAryBw3NrW9B6vitGfmNUq3N7i9B0sIEQRBxe8x0w7Z7Q1u78GS4jheLBZO+3ekGYSq5nT2q9/3/TAM3Qwz97fj2nuw4jg+HA6F3/jqXGjUXtze4PYerDAM5TX5ue3q1TBtF25vcHsPljQejx3HMfwnpTmoU9Uc0zqiTqqrwdvtNvdU2+3Wdd0oitrytdIdNh0s9Qu7+torc1zaYAsOlgyRcq8dxwmCwPCWX93glh4s2Tbf901upFman8bVWfIFL5zVaOa1+rc3WJ0bsk/S3utyzTxMha5oqmUHq0VXa0oXNdiOg5Xr/DBtivqx6xrc3oN1PEnfzCs/jMLwH+5ICCHnNmZnJAgh1BcKU93NYc3BiuPYdV3Vl3D496XgBrqiwXYcLDldSVV+WiwWruvqblSV6xrc3oMl35MM/F1Gd6rrEPk5bFFP1V0bLJ/cwAJCFcw8TIXqbWqLDpb5lSRzam9wiw5WlslFxQrV0mCTD1bhe5KeqpPoqWpa6y6guFOD5Qh96+aZdVMrDpbs75GNlGcCYzsApDs1uBUH65jsBHLa0/JaGmzswcp+7YcZ2To72f+FwkT15gghyj48Zr41W9dgdFYcx/IaJd/3W/HmbF2D6yJ3tnDWsxBCjYGao3UNrlfZDqrtTGDPIVQ1rTCmnKwvp9HVDVYLs1f8/jZzlzuo7QerdQHllga3+mCp6/PbcjK+scHtPViFrVKnA3lrWw5io/SOPnaNfM2PV3eR2w2cAnJjg8sefv4zGMXuOVXtPVitm+dxe4Pbe7Aqmmfmcby9we09WMfMPEZG4aVplMr+Jzca4qIGy4tisl8QFfc0dpcr2BSqbDpYsnm+70fldLfxXy5tsE0HS52Vc/uo3rGmJYxLG2zTwTpGqDqJl6Zpzm9BEMirc838KlHOb7A6Vagt2Wtx5RdN6y7OyrIpVNl0sJwz6G7jv1zaYJsO1iGz+7Llucbrbl2Bixps2cHKIVSdxEvTtMKaHyZ/rs5v8PG3SdnDDd/lMnaHqkNrD1Zhm3N0t/FfLm2wTQdLKpyvY2ZlAen8Btt3sLIIVSe5B+Mr41lJzV50DK78lnVjg+PfxG+1tg514mC1SHsPVvYrpRUtv73B7T1YOB+hCgAAoAYU/wQAAKgBoQoAAKAGhCoAAIAaEKoAAABqQKgCAACoAaEKAACgBoQqAACAGhCqAAAAakCoAgAAqAGhCgAAoAaEKgAAgBoQqgAAAGpAqAIAAKgBoQoAAKAGhCoAAIAaEKoAAABqQKgCAACoAaEKAACgBoQqAACAGhCqAAAAakCoAgAAqAGhCgAAoAaEKgAAgBoQqgAAAGpAqAIAAKgBoQoAAKAGhCoAAIAaEKoAAABqQKgCAADtI4RwXdd13TAMdbflH+7hcNDdBgAAgHPFcTwejx3HiaLIcRz5bxPyDKEKAAC0hkpUKsAcb9GF4T8AANAacrAvCAK1RQiRvUkjQhUAAGiN7XbrZIKUJDPWYrHQ0iSFUAUAANpB9UXlQlXuf3UhVAEAgHaI47j6DnpHAAlVAACgHeTY3zFDeqr+0N0AAACAy5SlqJNdWXdFqAIAAC1T1mWlF6EKAAC0zHFJKtd1Hd3jgMypAgAA7eD7vu4mVCFUAQCAdlNTqeipAgAAOK2sYgKhCgAA4AIqM+Wu8pP/q31wkFAFAABaQ65Ik+uykhcDsvYfAADAucIw9H1/u93mlqwJgkB7CVD3+KJEAAAAkwkhsqWqgiDQ3k3lEKoAAEBLyalU2juoFEIVAABADZhTBQAAUANCFQAAQA0IVQAAADUgVAEAANSAUAUAAFADQhUAAEANCFUAAAA1IFQBAADUgFAFAABQA0IVAABADQhVAAAANSBUAQAA1IBQBQAAUANCFQAAQA0IVQAAADUgVAEAANTg/wMEw/6Zohj9RgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ROOT import gROOT \n",
    "gROOT.GetListOfCanvases().Draw()"
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}