fit.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
       "    return false;\n",
       "}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "IPython.OutputArea.prototype._should_scroll = function(lines) {\n",
    "    return false;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import glob\n",
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from bisect import *\n",
    "from decimal import Decimal\n",
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "27 27\n",
      "[  55.73411261 -131.19114854  583.84576629 -440.51423473   99.91044255]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<function matplotlib.pyplot.show(*args, **kw)>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU1fnH8c8TCJsQWaIh7NaCGjcUxLUa1FpEK9Za1BZFa8WqtbXaitYWXOta21p/hWpVcAHElcW9SEoViYCi7IoLm5EdQlgCSZ7fH3fQEBIySSZzZybf9+s1r8yce+fmmzF5PNx77jnm7oiISGpJCzuAiIjEnoq7iEgKUnEXEUlBKu4iIilIxV1EJAU1DjsAQGZmpnfr1i2mx9yyZQv77LNPTI8ZL8oeDmWPv2TNDYmRffbs2Wvdfb/KtiVEce/WrRuzZs2K6THz8vLIzc2N6THjRdnDoezxl6y5ITGym9nSqrbptIyISApScRcRSUEq7iIiKUjFXUQkBam4i4ikIBV3EZEUVG1xN7POZjbVzBaY2Xwz+02kva2ZvWVmn0a+tom0m5k9ZGZLzOxjMzu6vn8IERHZXTTj3EuAG9z9AzNrBcw2s7eAS4Ep7n6Pmd0E3AQMBc4EukcexwIjIl9FRBqkMfnLmDBnZaXbcjpkMPyHh8b8e1bbc3f3Anf/IPJ8M7AQ6AgMAEZHdhsNnBt5PgB40gMzgNZmlh3z5CIiSWLCnJUsKCiM6/es0R2qZtYNOArIB7LcvSCy6WsgK/K8I7C83NtWRNoKEBFpoHKyM3j2yuPj9v2iLu5m1hJ4AbjO3QvN7Jtt7u5mVqMlncxsCDAEICsri7y8vJq8vVpFRUUxP2a8KHs4lD3+kjU31Cz7xo3bAOL7s7p7tQ8gHXgDuL5c22IgO/I8G1gcef4v4KLK9qvq0atXL4+1qVOnxvyY8aLs4VD2+EvW3O41yz5w5HQfOHJ6zDMAs7yKuhrNaBkDHgMWuvuD5TZNBAZHng8GJpRrvyQyauY4YJN/e/pGRETiIJrTMicCFwNzzWxOpO0PwD3AeDO7HFgKDIxsexXoDywBtgKXxTSxiIhUq9ri7u7vAFbF5tMq2d+Ba+qYS0RE6iAh5nMXEUkmY/KXMTp/GyMWvxfV/gsKCsnJzqjnVLvT9AMiIjU0Yc5Klm0ui3r/nOwMBvTsWI+J9qSeu4hILXRplRbXces1pZ67iEgKUs9dRKSGGvlO9vVCWP85FG8OHtsLg68l26CsFLwseHzzvHT39l3buhwH391jbEqdqbiLSEra22RdFZmXkVG2kXala2lXtoZ2pWtpW7qWdqVraFO2nhZlW2juW4NH2VbGsDN440MxCHrSb1XcRUSitWuyrl2jVFqVbaLLzi/pVLKUzNLVtCtdExTwsrW0LV1HY0p2e/9O0lnXKJONaW1Y3yiTbdaCrWkt2GbBI6NZY/qfdAw0bRU8mmVA0wxo3AzSGoE1AkuDtLTgqzWKtJd/XtUo87pTcReR1LJjC6xZRO7WN7h4ny85u/kGWLUAtqz+dp9GTSGjA+zbCTIOC55ndIw8gufp+2TS3oz2VXybvLw86Jkbhx+odlTcRSR5bVkHy96DgjmweiGsmg8bvgScq4BimsL2HOj+fdg/B7JyYL9DoFX7eu01JwIVdxFJHpu/hqXvwtLp8OW7sGZh0G6NoN2BkH0kHHkRZOXw67d3sLpRFuOGnBRu5pCouItI4tq4LFLI3wm+rv8saG/SEjofC0f8BLqeCNk9Ib3Zbm9dNS26u0dTlYq7iCSO7YXw6Zvw6VtBD31TZN2fZq2h6wnQ+7KgmLc/AhqpfO2NPh0RCVfRGlj8KiycBF/8F0p3QItM6HYinHBtUMz3z2HMzBXB0Ma5xcDMag8bxnwuiUTFXUTib+NyWDSZnh8+Bf9dGNzQ07or9BkCh5wDnY4JhhCWU3FoY3XCmM8lkai4i0h8rFkc9M4XTgpGtwCN9+kKJ/8eDvkhZB1W7QiWeK9DmsxU3EWk/hStgY/GwIfPwNrFQVunY+D7t8PBZzNr7nJyc3NDjZiqVNxFJLbKyuCLPJg9Gha9AmU7ocvx0P8BOPis4CahbywPK2XKq7a4m9njwNnAanc/LNLWExgJNANKgKvd/f3Ieqt/J1hmbytwqbt/UF/hRSSBbP4a5jwTFPWNS6F5Wzj2Sjj6EtjvoD12T4YFL5JZND33UcDDwJPl2u4DbnP318ysf+R1LnAm0D3yOBYYEfkqIqmorBQ+extmj4LFrwUzH3b7Hpw2LDiP3rhplW/dteBF69bRfauGfoG0pqJZQ3WamXWr2Azs+l/ovsBXkecDgCcj66jOMLPWZpbt7gUxyisiiaDwK/jwafjgyWAseotMOOFXcPTg4E7RKCX6ghfJzII6XM1OQXGfXO60zCHAGwQLZ6cBJ7j7UjObDNwTWVQbM5sCDHX3WZUccwgwBCArK6vXuHHjYvID7VJUVETLli1jesx4UfZwKHv1mm9dSZdlL5K1Ko80L2F9m54UZJ/B2sw+eFp6jY51d/42SktL+eMJ+sxrq2/fvrPdvXdl22p7QfUq4Lfu/oKZDQQeA06vyQHc/RHgEYDevXt7rK+Y5+XlJe1VeGUPh7LvRcHH8M6DMP/l4FTLMZfDcVfRtu0BtK3lIUcsfo+NGzfqM68ntS3ug4HfRJ4/B/w78nwl0Lncfp0ibSKSjJbNgP/9JZgSoGlGsLDEcVdDy/0q3b0mC2QsKCikQ/NYhpXyalvcvwJOAfKAU4FPI+0TgV+Z2TiCC6mbdL5dJMm4w2dT4H8PBvO7tGgHp/4JjvkFNN/71c+a3EWak53BIS2KYpVaKohmKORYgpEwmWa2AhgOXAH83cwaA9uJnDsHXiUYBrmEYCjkZfWQWUTqQ1kZLJoU9NQLPgoWruh3bzCUsUmLqA9Tk7tI8/LyahlWqhPNaJmLqtjUq5J9HbimrqFEJI7KymD+i/Dfe2HtJ9Duu3DOw3DEBYyZ/TUTnvgo6kNpLHri0B2qIg3Zl+/Cm3+Erz4I5nY5/wnIGRCs74km60pmKu4iDdHaJfCf4bBocnD65Uf/gsMH7jETI2iyrmSl4i7SkGxZF5x+mfUYNG4WXCg97uoanVOX5KDiLtIQ7NwO+SPZkXc/jUq2MqXFmTzfchCbFraBhVWfU9c59OSl4i6Sytxh3gvwn9tg0zLmNj2W20suoNm+OVG9XefQk5eKu0iqWjo9uFi6cja0PxwGTOC+t5rSDHQOvQFQcRdJNVvWwht/gI+fhVYd4NwRcMSFkYul0U2vK8lPxV0kVbjDR+OCwl68OVi+7qTrdbG0gVJxF0lyY/KX8dZ7X9Bu+nAO3zGHxek5PNL2z6z4tBt8uvvFUl0gbThU3EWSWelOeOdBRmx/mrK0dB7NuJYpLc7Ebc/x6qALpA2JirtIsloxGyb9mp9unsc7jfpw0nVPckVGNleEnUsSgoq7SLIp3gxv3wn5/4JW2dzfZhj/2Z7DGxnZYSeTBFL5v91EJDEtfg3+79igsPe5Aq7JZ1azE8JOJQlIPXeRZLB1PbxyPcx/CfY7BC4fBZ37hJ1KEpiKu0iCqbia0eHFH3D1xgfIKCvk+ZaDmdjofEpfLWXXmHWtaCSVUXEXSTC7ptk9sn0zLtz8BGdveYkVjbtwT9s7WJp+4B77a0UjqYyKu0gC+kHmeh4o+wdsmQ/HXEGn79/OfXu5GUkrGklF1V5QNbPHzWy1mc2r0H6tmS0ys/lmdl+59pvNbImZLTazH9RHaJGU5U6/LS/z57XXwpbV8NPn4KwHdJep1Fg0PfdRwMPAk7sazKwvMAA40t2LzWz/SHsOcCFwKNAB+I+Z9XD30lgHF0k5m7+Gl6/mssIpzG7ah15XjYGW+4WdSpJUNGuoTjOzbhWarwLucffiyD6rI+0DgHGR9i/MbAnQB81WJA1YxQuklem9/T2u3PQ3mpZt5z5+wbw2P+ZZFXapAwvWtK5mp6C4T3b3wyKv5wATgH7AduB37j7TzB4GZrj705H9HgNec/fnKznmEGAIQFZWVq9x48bF5AfapaioiJYtW8b0mPGi7OGor+x3529j2eYyurTa8yxoU9/OlTuf4qzSKXxq3bi3ybUsT+vI8R0ak9s5Pervkayfe7LmhsTI3rdv39nu3rvSje5e7QPoBswr93oe8A/ACHrmX0SePwwMKrffY8D51R2/V69eHmtTp06N+THjRdnDUV/ZB46c7gNHTt9zw8oP3B862n34vu5v/sl9Z3Gtv0eyfu7Jmts9MbIDs7yKulrb0TIrgBcjB3/fzMqATGAl0Lncfp0ibSKyizvMHgWv3Qj77AeDJ8IBJ4edSlJMbacfeBnoC2BmPYAmwFpgInChmTU1swOA7sD7sQgqkhJ2boMJ18Dk66Db9+CX76iwS72otuduZmOBXCDTzFYAw4HHgccjwyN3AIMjvfj5ZjYeWACUANe4RsqIBNZ/AeMvhq/nwilDg0dao7BTSYqKZrTMRVVsGlTF/ncBd9UllEgii2b0S3kLCgoZ1GYhPPJg0PDT56DHGfWUTiSgO1RFamjX9ADRrGhkXsqtLV/ixxvGQvsj4IKnoE23+g8pDZ6Ku0gt5GRn8OyVx+99py3r4MVfwNdvw1GDoP8DkK4ZviQ+VNxF6sPK2TB+MBSthh8+BL0Gh51IGhgVd5FYKj/MsWV7+Pnr0PHosFNJA6TiLhIrO7fBKzfAnGfgwNPgx/+GFm3DTiUNlIq7CMEImNH52xixuPppkCq9mLr5axh7EXz1gYY5SkJQcRchGAGzbHMZrVtXv29OdgYDenb8tqHgo6Cwb9sIF46Bg8+qv6AiUVJxF4no0iqt+hEwFS2cDC9eAc3bBufXs4+on3AiNVTb6QdEGjZ3eOdv8Owg2P8QuOJtFXZJKOq5i9RUyY5gbpg5z8Ch58G5/9T4dUk4Ku4iNbFlXdBbXzYdTrkJcm8Cs7BTiexBxV0kWqsXwdgLoLAAfvwYHH5+2IlEqqTiLhKNJf+B5y6Dxs3gslehU+WL34gkCl1QFalO/iPwzEBo3TW4cKrCLklAPXeRqpSWwOs3wcxH4aD+cN6j0DQ51/uUhkfFXaQyO7bCC5fD4lfhhGvh9Nt0x6kkFRV3kYq2rocxF8CKmcE0vX2uCDuRSI1Ve87dzB43s9WRJfUqbrvBzNzMMiOvzcweMrMlZvaxmWk6PEkuG5bCY2cEUwoMfFKFXZJWND33UcDDwJPlG82sM3AGsKxc85kEi2J3B44FRkS+isRVbZbCO6npUnjsOijZBpe8DF1PqMeEIvWr2p67u08D1ley6a/AjYCXaxsAPOmBGUBrM8uOSVKRGti1FF60ftL2c/6249bgvPrP31Bhl6RXq3PuZjYAWOnuH9nud+d1BJaXe70i0lZQyTGGAEMAsrKyyMvLq02UKhUVFcX8mPGi7HW3ceM2OjSHqw4qrnbf/VdN4+BFf6eoWXs+yLmN4gWrYMGqOKSMnUT53GsqWXND4mevcXE3sxbAHwhOydSauz8CPALQu3dvz83Nrcvh9pCXl0esjxkvyl53u+Zlz82tZpbH6Q/Dwr9A1xP5uNPVnPT9s+OQLvYS5XOvqWTNDYmfvTY3MR0IHAB8ZGZfAp2AD8ysPbAS6Fxu306RNpHEUlYGb9wCb94Ch5wDg16kJF1j2CV11Li4u/tcd9/f3bu5ezeCUy9Hu/vXwETgksiomeOATe6+xykZkVCV7AjmYH/vYegzBH4yCtKbhZ1KJKaiGQo5FngPOMjMVpjZ5XvZ/VXgc2AJ8ChwdUxSisTK9kJ45nyY9zycNhzOvE83J0lKqvacu7tfVM32buWeO3BN3WOJ1IOiNfD0ebB6AZw7Enru9VdbJKnpDlVpGDatgCcHwKaVcNE46P79sBOJ1CsVd0l9a5fAU+fC9k1w8UvQtYbrpIokIRV3SW1fzwsKu5fB4EnQoWfYiUTiQvO5S8rqvmMBjOoPaelw2esq7NKgqLhLSjq8+AP+uP5maNEOfv467Ncj7EgicaXiLqln4SSGrh/OqkbZQY+9TdewE4nEnYq7pJY5Y2H8YL5MP5Db2t0PrbLCTiQSCl1QldSR/wi89ns44GTu2PpbitOah51IJDQq7pIU9jo/uzs/KhrHhUWjmdn0eP6+7Xd8tGo7Odkq7tJw6bSMJIUq52d352eb/82FRaOZ1vxU/trmFnZaE3KyMxjQs2P8g4okCPXcJWnkZGfw7JXlbkAqK4NXrofZL8Axv+DkM+/n5DT1V0RAxV2SVVkpTLwW5jwDJ14Hp98Kuy8cI9KgqbhL8iktgZd/CXOfg9yb4ZShKuwiFai4S3Ip3QkvXA4LJsBpw+B7N4SdSCQhqbhL0mjsO2D8JbD4VTjjLjjhV2FHEklYKu6SFNK9mBs23AFfz4L+D0CfK8KOJJLQVNwl8e3YwtD1wzl0x0fww79Dr0vDTiSS8KJZZu9xM1ttZvPKtd1vZovM7GMze8nMWpfbdrOZLTGzxWb2g/oKLg1E8WZ45iccuuNjRux7gwq7SJSiGRQ8CuhXoe0t4DB3PwL4BLgZwMxygAuBQyPv+aeZaYFKqZ3tm+Cp82DZDP7R+kamtTg97EQiSaPa4u7u04D1FdredPeSyMsZQKfI8wHAOHcvdvcvCBbK7hPDvNJQbF0fLIv31Yfwk1FMb54bdiKRpBKLc+4/B56NPO9IUOx3WRFp24OZDQGGAGRlZZGXlxeDKN8qKiqK+THjpaFnT9+xiSM/Gk6LrcuZf+hNrFudwcaNGwHq9XNp6J97GJI1NyR+9joVdzO7BSgBnqnpe939EeARgN69e3tubm5douwhLy+PWB8zXhp09qLVMPocKC6An43n8O+eBsCIxe8BkJtbf+ufNujPPSTJmhsSP3uti7uZXQqcDZzm7h5pXgl0Lrdbp0ibyB4qzvS4b+l6hq2/iczSVdzX5jbmT2kBU4KivqCgkJzsjLCiiiSdWs2yZGb9gBuBc9x9a7lNE4ELzaypmR0AdAfer3tMSUXlZ3psXbqO4etuJLN0Nfe0uZP5TXdf71SzPIrUTLU9dzMbC+QCmWa2AhhOMDqmKfCWBXN6zHD3X7r7fDMbDywgOF1zjbuX1ld4SX452Rk8e1FXGP0raLQRLnmZW7vW36kXkYai2uLu7hdV0vzYXva/C7irLqGk4WhbugZGXR2cax/0AnQ5LuxIIilBd6hKaNqVrmHYuhuhUREMehG6HBt2JJGUoeIu4di4nGHrbiSjbBNcOhE6HxN2IpGUomVrJP42LoNRZ9GqrJA7296twi5SD9Rzl/jasBRGnw3bN3Fn27v5vEmPsBOJpCT13CV+NnwJo86C7YVwyQQVdpF6pOIu8bH+c3jiLNhRBIMnQoejwk4kktJ0Wkbq37rPYPQPYedWuGQiZB8RdiKRlKfiLjEzJn8Zo/O3fTMPDED7kpUMWzeUdN/Bne3uYenELYCmFBCpbzotIzEzYc5Klm0u++Z1+5KVDF93I+ns5PZ297I0/Tu77a8pBUTqj3ruElNdWqXx7JXHw9olMOpSaGYw+HUeyMoJO5pIg6Keu8Te2k+DUTFlJTB4Eqiwi8SdirvEVKeyr2DU2eClcOlkFXaRkKi4S8x0KFnOfcW3B4V98GTY/5CwI4k0WCruEhtrFjNs3Y0YHinsB4edSKRBU3GXulu9CEadjQE3Nh2mwi6SADRaRupm9aJgrhhL47Z297G8qFXYiUQE9dylLlYv/KawM3gyXzXuXP17RCQuqi3uZva4ma02s3nl2tqa2Vtm9mnka5tIu5nZQ2a2xMw+NrOj6zO8hGj1wmBUjDWCS1+B/TQJmEgiiabnPgroV6HtJmCKu3cHpkReA5xJsCh2d2AIMCI2MSWhrFoQFPa0xsFwx8zuYScSkQqiWUN1mpl1q9A8gGDRbIDRQB4wNNL+pLs7MMPMWptZtrsXxCqwxM+Y/GVMmLNyt7auOz/nj+tvpoTG3N7uTgpeWAusBYK5Yjo0DyGoiOzBgjpczU5BcZ/s7odFXm9099aR5wZscPfWZjYZuMfd34lsmwIMdfdZlRxzCEHvnqysrF7jxo2LzU8UUVRURMuWLWN6zHhJlOx3529j2eYyurQK/oF3YNkX3FN8F8XWhBubDOOrtPZ7vOeodqX0+2742WsjUT732kjW7MmaGxIje9++fWe7e+9KN7p7tQ+gGzCv3OuNFbZviHydDJxUrn0K0Lu64/fq1ctjberUqTE/ZrwkSvaBI6f7wJHTgxcrZrvf3dn9wUPd131e5XsSJXttKHv8JWtu98TIDszyKupqbUfLrDKzbIDI19WR9pVA+SETnSJtksxWzIInz4VmrYOLp20PCDuRiFSjtsV9IjA48nwwMKFc+yWRUTPHAZtc59uTWo8dC4LC3qJtUNjbdA07kohEodoLqmY2luDiaaaZrQCGA/cA483scmApMDCy+6tAf2AJsBW4rB4yS5wcvGMeN63/E7TpEIyKyegQdiQRiVI0o2UuqmLTaZXs68A1dQ0lCeCL/3Hz+ltYl7Y/HS97FVrtefFURBKX7lCVPX02FZ75CWsaZXFbu/tU2EWSkOaWkd0t+Q+M+xm0PZDb+BObG7UOO5GI1IJ67vKtT96EsRcFd5wOnqTCLpLE1HNvYCq76xSg1/b3uH7DXSxNP4C7fBhbnv6EBQWF5GRnhJBSROpKPfcGZsKclSwoKNytrc+2d7h+w518kX4gd7a9my1pwbS9OdkZDOjZMYyYIlJH6rk3QDnZGTx75fHBizljYcKfofMxdP/ZczzRbN9ww4lITKjn3pDN/De8/Evo9j24+CVQYRdJGSruDdW7D8ErN0CPfvDT8dBkn7ATiUgM6bRMQ+PO+UVPw1vPwKE/gvMehUbpYacSkRhTcW9I3Ll486OcveVF6DkIznkI0hqFnUpE6oFOyzQUZWUw+becveVFXmtxDpzzDxV2kRSm4t4QlJYEF05nP8FL+1zAqIyrIE3/6UVSmU7LpLqSYnj+57BoMpz6J8YtPDHsRCISB+q+pbIdW4PpBBZNhn73wsm/CzuRiMSJeu6panshjL0Qlk6Hcx6Goy8OO5GIxJGKe5KrbK6YVmWbGLp+GAfsXML/tR7K9JnfhZnvAWi+GJEGQqdlklzFuWIyS1Zx29ob6Lrzcx5s80emN8/dbX/NFyPSMNSp525mvwV+ATgwl2BZvWxgHNAOmA1c7O476phT9uKbuWK+ngdPXwrp2+CSSdzY9fiwo4lISGrdczezjsCvgd7ufhjQCLgQuBf4q7t/F9gAXB6LoFKNL9+BJ/qDpcFlr4MKu0iDVtfTMo2B5mbWGGgBFACnAs9Hto8Gzq3j95Bq9Nn2Djx1XrAc3uVvQlZO2JFEJGQWrGldyzeb/Qa4C9gGvAn8BpgR6bVjZp2B1yI9+4rvHQIMAcjKyuo1bty4WueoTFFRES1btozpMeOlJtnnvjuBa3Y+weaMHsw9/E+UpLeq53R711A+90STrNmTNTckRva+ffvOdvfelW5091o9gDbA28B+QDrwMjAIWFJun87AvOqO1atXL4+1qVOnxvyY8RJV9rIy9yl3uA/P8Jl//r578ZZ6zxWNlP/cE1SyZk/W3O6JkR2Y5VXU1bqcljkd+MLd17j7TuBF4ESgdeQ0DUAnYM813aRuSktg0q9h2v283fwH/KXNMGjSIuxUIpJA6lLclwHHmVkLMzPgNGABMBU4P7LPYGBC3SLKbnZshfEXwwdPwsm/51/7XkeZaQIwEdldrYu7u+cTXDj9gGAYZBrwCDAUuN7MlhAMh3wsBjkFYOt6eOpcWPwa9H8ATv0jmIWdSkQSUJ3Gubv7cGB4hebPgT51OW5DNiZ/GaPztzFi8Xu7tbcrXcMf1t9CVkkBD7W+hfc/PAw+fE93nIpIpTT9QIKZMGclyzaX0br1t20H7ljM7zbcTjPfxl1t72Jh0yO+2aY7TkWkMiruCahLq7TgjlOAuc/DhKGQsT9cNIlbsw4NN5yIJAUV90RVVgZ5f4Zp90OXE+CCp2CfzLBTiUiSUHFPQE19Ozx3CSycBEddDGc9CI2bhB1LRJKIinuCaVe6mhuKh8Oi5dDvHjj2lxoRIyI1puKeSJbl8+e1v6ax74BBz8F3Tw87kYgkKc3nnijmjIHRZ7PNWnBd0ztU2EWkTtRzD1tZKfznVpj+EBxwMrdsvZaVhaVhpxKRJKeee5i2F8K4nwaF/ZhfwKAX2ZIW7qyOIpIa1HMPy/rPYexPYe0nwVQCfa4IO5GIpBAV9zDMfR4mXQdpjeDil+A7p4SdSERSjIp7PRuTv4wJc4JZj5uWbefSwn9y6rY3WZyew0OthrL2rSbAt/PILCgopEPzkMKKSMpQca9nE+asZEFBIf0y13LdhrvJLl3Biy0v5LmWF1c6VW9OdgaHtCgKIamIpBIV9/rmznUZ/+XyDf+G5q3hvJc57zu5nLeXt+Tl5cUpnIikKhX3+rR1PTdsuIM+xdODcevnjoSW+4WdSkQaAA2FrC/LZsDI73F0cT5PtfoF/PQ5FXYRiRsV91grKw1mcnyiPzRqzLB2DzK55fmQpo9aROKnThXHzFqb2fNmtsjMFprZ8WbW1szeMrNPI1/bxCpswissCJbBe/tOyBkAV07jsyYHhZ1KRBqgunYn/w687u4HA0cCC4GbgCnu3h2YEnmd2tzh4/Ew8kRYPhPO+Qec/zg02zfsZCLSQNX6gqqZ7QucDFwK4O47gB1mNgDIjew2GsgjWDQ7Na35BF65Hr78H3Q4Krhouv/BYacSkQbO3L12bzTrCTwCLCDotc8GfgOsdPfWkX0M2LDrdYX3DwGGAGRlZfUaN25crXJUpaioiJYtW8b0mOWllRbTdel4Oi9/mbK0pnz+nYv5qsMZUGHs+t352wC4+djo70yq7+z1SdnDkazZkzU3JEb2vn37znb33pVudPdaPYDeQAlwbOT134E7gI0V9ttQ3bF69erlsTZ16tSYH/Mbi193/+th7sMz3F+4wn3zqip3HThyug8cOb1Gh6/X7PVM2cORrNmTNbd7YmQHZnkVdbUu49xXACvcPT/y+nmC8+urzCzb3QvMLBtYXYfvkVg2rYDXhu7uQSAAAAa4SURBVMKiyZDZAwZPggNODjuViMgeal3c3f1rM1tuZge5+2LgNIJTNAuAwcA9ka8TYpI0TKU7YcYIyLsHvIw5PX7DfYWnU/pmOuXnhanMgoJCcrIz4pNTRCSirneoXgs8Y2ZNgM+BywhG4Iw3s8uBpcDAOn6PcC2bAZOvh9XzoUc/OPNe7h5fwIKvC8nJTq/27TnZGQzo2TEOQUVEvlWn4u7ucwjOvVd0Wl2OmxDWfgr/+wt8NBYyOsEFz8DBZ0UWqy4gJzuDZ688PuyUIiKV0twyFa2YBe/8FRa9Ao2bwom/gVOGQpN9wk4mIhI1FXcIbkJaMgXe/VswXr1Zazj5d9DnSs0HIyJJqWEX99ISWPAyvPM3WDUXWnWAM+6CXoOhqdYyFZHk1TCL+46tMOcZmP4P2LgUMg+CAf+Ew38CjZuEnU5EpM4aVnEvWg2zR0P+CNi6Djr1gX53Q48zNWujiKSUlCzuY/KXMTp/GyMXvct3dn7KUcXvc3Tx+xy481MAPmh6DBPa3cyikkPhvwb/za/miLvT2HURSXRJXdxvmzSfBV8V7ta2T9lmWiyfxpBGH3J68Vz2LdtEGWl8kn4wY1tdysymx7MyvWudvq/GrotIokvq4g6AO11KvuDo4vfpuX0mB+1cSFqTMramtaLFYf2h+xmkHXgqB7doy8HARWHnFRGJg6Qu7sMPWgGfXAebvwoaso+E7r+D7mfw/pJCcvsm/71UIiK1kdTFnYwO0Kk39PhBsAB1q/bfbvssL7RYIiJhS+7i3v5wuOCpsFOIiCQcjf8TEUlBKu4iIilIxV1EJAWpuIuIpCAVdxGRFKTiLiKSglTcRURSkIq7iEgKMncPOwNmtoZgMe1YygTWxviY8aLs4VD2+EvW3JAY2bu6e6XLxSVEca8PZjbL3StbvDvhKXs4lD3+kjU3JH52nZYREUlBKu4iIikolYv7I2EHqANlD4eyx1+y5oYEz56y59xFRBqyVO65i4g0WCruIiIpKOmLu5n1M7PFZrbEzG6qZPv1ZrbAzD42sylmVrfVsWOouuzl9vuxmbmZJcywq2iym9nAyGc/38zGxDtjZaL4feliZlPN7MPI70z/MHJWxsweN7PVZjaviu1mZg9FfraPzezoeGesShTZfxbJPNfMppvZkfHOWJnqcpfb7xgzKzGz8+OVrVrunrQPoBHwGfAdoAnwEZBTYZ++QIvI86uAZ8POHW32yH6tgGnADKB32Llr8Ll3Bz4E2kRe758kuR8Broo8zwG+DDt3uWwnA0cD86rY3h94DTDgOCA/7Mw1yH5Cud+VMxMle3W5y/1evQ28CpwfduZdj2TvufcBlrj75+6+AxgHDCi/g7tPdfetkZczgE5xzliVarNH3AHcC2yPZ7hqRJP9CuD/3H0DgLuvjnPGykST24GMyPN9ga/imG+v3H0asH4vuwwAnvTADKC1mWXHJ93eVZfd3afv+l0hgf5Oo/jMAa4FXgAS4Xf8G8le3DsCy8u9XhFpq8rlBD2bRFBt9sg/qzu7+yvxDBaFaD73HkAPM3vXzGaYWb+4pataNLlvBQaZ2QqCnti18YkWEzX9e0hUifR3uldm1hH4ETAi7CwVJfcC2TVgZoOA3sApYWeJhpmlAQ8Cl4YcpbYaE5yaySXohU0zs8PdfWOoqap3ETDK3f9iZscDT5nZYe5eFnawhsDM+hIU95PCzhKlvwFD3b3MzMLOsptkL+4rgc7lXneKtO3GzE4HbgFOcffiOGWrTnXZWwGHAXmRX5r2wEQzO8fdZ8UtZeWi+dxXEJw33Ql8YWafEBT7mfGJWKlocl8O9ANw9/fMrBnBBFEJ9U/uKkT195CozOwI4N/Ame6+Luw8UeoNjIv8jWYC/c2sxN1fDjdW8p+WmQl0N7MDzKwJcCEwsfwOZnYU8C/gnAQ577vLXrO7+yZ3z3T3bu7ejeA8ZCIUdojicwdeJui1Y2aZBKdpPo9nyEpEk3sZcBqAmR0CNAPWxDVl7U0ELomMmjkO2OTuBWGHioaZdQFeBC5290/CzhMtdz+g3N/o88DViVDYIcl77u5eYma/At4guGL9uLvPN7PbgVnuPhG4H2gJPBf5v+sydz8ntNARUWZPSFFmfwM4w8wWAKXA78PujUWZ+wbgUTP7LcHF1Us9MiQibGY2luB/mJmRawLDgXQAdx9JcI2gP7AE2ApcFk7SPUWRfRjQDvhn5O+0xBNgxsUocicsTT8gIpKCkv20jIiIVELFXUQkBam4i4ikIBV3EZEUpOIuIpKCVNxFRFKQiruISAr6f9GOQSWbh2qyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ey = np.array([50., 51., 51., 57., 66., 75., 83., 92., 100., 109., 118., 129., 141., 151., 156., 160., 164., 169., 173., 176., 178., 182., 186., 191., 191., 189., 186., 178., 171., 167., 166., 165., 176., 191., 203., 201., 199., 193., 184., 175., 173., 171., 174., 188., 202., 206])\n",
    "ey = np.array([50., 51., 51., 57., 66., 75., 83., 92., 100., 109., 118., 129., 141., 151., 156., 160., 164., 169., 173., 176., 178., 182., 186., 191., 191., 195.,201.])\n",
    "err =  np.empty(len(ey))\n",
    "er.fill(10.)\n",
    "bin_centres = np.arange(start=0.175,stop=0.175+0.05*len(ey),step=0.05)\n",
    "print (len(ey), len(bin_centres))\n",
    "\n",
    "def p3(x, *p):\n",
    "    a0,a1,a2,a3,a4 = p\n",
    "    return a0+x*(a1+x*(a2+x*(a3+x*a4)))\n",
    "\n",
    "# p0 is the initial guess for the fitting coefficients\n",
    "p0 = [1., 0.5, 0.,0.,0.]\n",
    "\n",
    "coeff, var_matrix = curve_fit(p3, bin_centres, ey, p0=p0, sigma=err)\n",
    "print(coeff)\n",
    "# Get the fitted curve\n",
    "hist_fit = p3(bin_centres, *coeff)\n",
    "# print(hist_fit)\n",
    "\n",
    "\n",
    "plt.step(bin_centres, ey, label='', where='mid')\n",
    "plt.plot(bin_centres, hist_fit, label='')\n",
    "plt.grid()\n",
    "plt.show"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13 13 14\n"
     ]
    }
   ],
   "source": [
    "chi2 = np.array([1.1,1.6832680195,2.1419517475,2.4151718822,2.8056427254,2.883004237,2.9962983052,3.07511163,3.2028987131,3.3467331513,3.2110779156,3.4331270667,3.6743670329])\n",
    "err =np.array([0.014508,0.00843219,0.0225587,0.0419361,0.0702273,0.0938654,0.114591,0.136352,0.166318,0.195104,0.188938,0.244216,0.30525])\n",
    "bin_edges = np.array([0.1,1.1,2.1,3.1,4.1,5.1,6.1,7.1,8.1,9.1,10.1,11.1,12.1,13.1])\n",
    "bin_centres = (bin_edges[:-1] + bin_edges[1:])/2\n",
    "\n",
    "print len(chi2),len(err),len(bin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.68177776  0.74609577 -0.08035491  0.00315399]\n",
      "[1.10118872 1.68274116 2.13386206 2.47347532 2.72050485 2.89387457\n",
      " 3.01250838 3.09533021 3.16126396 3.22923354 3.31816288 3.44697587\n",
      " 3.63459643]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOW9x/HPjyTshDXVsCioGImIgAt6RQ244YqtexfXCre1rbW9XbRebW2R2tpae2lFq1VpqRVBxVoXsDAKCiggZTUgFBGI7CEMhJDld/84gyIEMkkmOTOT7/v1ysskczLzfRr49vDMc55j7o6IiKSXZmEHEBGRxFO5i4ikIZW7iEgaUrmLiKQhlbuISBpSuYuIpCGVu4hIGlK5i4ikIZW7iEgaygzrhTt06ODHHHNMWC+fUDt37qRNmzZhx0iIdBlLuowDNJZkFdZY5s2bt9ndc2o6LrRyP+yww5g7d25YL59QkUiEgoKCsGMkRLqMJV3GARpLsgprLGb2UTzHaVpGRCQNqdxFRNKQyl1EJA2p3EVE0pDKXUQkDancRUTSkMpdRCQNhbbOXUQk1f1tzhomL1hX65/L75rNvZce3wCJPqMzdxGROpq8YB1Li0rCjlEtnbmLiNRDfm42z448PewYB9CZu4hIY6msgPceh/ULGvylVO4iIo3hwzdg7GD45/dh4YQGfzlNy4iINKRNy2HKT2DFFOjYE675Kxx3SYO/rMpdRKQh7NoKkV8G0zDN28B5P4dBIyGzRaO8vMpdRCSRKvYEhf7mA1BWAifdCAV3Qdsat2BPKJW7iEgiuMPy12DK3bDlQzhqCFxwPxyWH0oclbuISH1tWAKv3wWrItC5N3x5AvQ+H8xCi6RyFxGpo+zKYq6OjoOxr0GLbBj2AJxyC2RkhR2t5nI3s5bAW0CL2PET3f3e/Y65Efg1sPc63DHu/nhio4qIJImKMnqseZ6HNz1Hcy+DQSPg7B9B605hJ/tUPGfuZcBQd4+aWRYw08xedffZ+x33rLt/K/ERRUSShDssewmm3sPR21Yzr8Wp/DX7Vh668Oqwkx2gxnJ3dweisS+zYh/ekKFERJLO+gXBvPpHb8MX8vl3v5/xq02Dwk51UHFdoWpmGWa2ANgITHX3OdUcdoWZLTSziWbWI6EpRUTCUlIEL34THiuATR/Axb+FkTPY1ql/2MkOyYIT8zgPNusAvAB8290X7/P9zkDU3cvMbCRwjbsPrebnRwAjAHJyck6aMKHhL8FtDNFolLZt24YdIyHSZSzpMg7QWMLSrLKMHh9P5og1kzCvYG33S/noyKuozGwDBGP5vyUZANw5qFWj5RoyZMg8dz+5puNqVe4AZnYPsMvdHzzI4xnAVndvf6jnycvL88LCwlq9drKKRCIUFBSEHSMh0mUs6TIO0FganTssngRT74WStcFWAefdB52P/txhkUiERwqDq00bc1dIM4ur3ONZLZMDlLt7sZm1As4DHtjvmFx3L4p9eRmwrA6ZRUTC9fF78PqdsPY9OLwffHEs9Doz7FR1Es9qmVzg6dgZeTNggru/bGb3AXPd/SXgO2Z2GVABbAVubKjAIiIJVV4Ky/4B88fB6hnQ9jC4bAz0/zI0ywg7XZ3Fs1pmITCgmu/fs8/ndwJ3JjaaiEjd1XQLvCPLVzJ012sMLp1OW4/ySUYu09vdyGutL2P3e63hvXcP+fzFxaWsLy0jPzc70dETQleoikha2nsLvH3Lt1XVTs4ojTC09DWOLl/BHrJ4t+VgprUextLmJ+BWu1tc5OdmM7x/t0RHTwiVu4ikrfzcbJ4dcRqsmRVMuyx5ESpK4bC+MPBXND/hKga37sTgOjx38OZw8t1eby+Vu4ikpfaV2zir9A0Y861gl8bm7eDEa2Hg9dB1QKibejUGlbuIpI+qSlg5DeY/zR83vkImldDjNBj8PTj+8uCmGU2Eyl1EUt+2j+D9v8KC8VCyDlp34ZU2lzO99TAeuiX59n1pDCp3EflUTStM6qu4uJRHCmcxvH83vjzoiPo9WUUZfPDPYC59VST43jHnwLDRcOyFjH9iXr3zpjKVu4h8qroVJom2tKgEoO7lvnEZzP8L/PsZKN0K7XtAwZ3BuvQO2tZqL5W7iHxOfm52g11Ov+8l+7WyuwSWTg7O0te+C82y4LiLgzdHjypI6YuNGorKXUSST2UFrJsHq6bDyunBdgBeCV3y4PxRwaqXNl3CTpnUVO4iEj73YLniyulBoa+eCWUlgAXLFs+4HfIuhO6npP0SxkRRuYtIOHZuDt4IXTUdVkaCHRgBOhwJfb8ERw2BXmcl1a3rUonKXUQaR3kpHbcu4Msli+lXNh9+vTL4fsv2QYmf+T04egh0OircnGlC5S4iDaOqCjYs+myqZc1sTqzYzfFkUti8Dwy9G44aCl376w3RBqByF5HE2b72szJf9Sbs2hx8P6cPnHwzC3d25sGNAylr1opnz0refVnSgcpdROqupAjWzQ2KfNX04E1RCPZEP+acYN78qALIzgVgayRC2eY6LIWUWlO5i0h8dm2F9fNh3fvBf9e/DztiN2DLag1HngEn3xwU+hf6aFVLyFTuIimiobcGAD67OrVsB6xf8FmJr5sPxR99dmDn3sGboF0HQNeBwbx5Zvxn5EuLSrjm0VkNMILPv0ay3kijMajcRVJEQ20NkOVl9CxfydHlK7it1XJOLFkFo1cDHhzQ/gjoNiA4K+82EHJPDFa41FFj3dwimW+k0RhU7iIppN5bA1SWw8alwZn43rPyjcugqiJ4vO1hsTPx64Ii7zog4VeCfnnQEfXfNExqpHIXSVdVlbB5RVDg6+cHhf7JIqgsCx5v2eGzqz/3Tq9kd9VceZpQuYukg71FXrQgmCsvWgBFC6F8Z/B4VptgXvzUW4Mi7zYQOvZSkacxlbtIqvlckb8flPkniz4r8sxWcPgJMOArkNs/KPIux+pCoSZG5S6SzKoqab1zDSx4hhu2v8pR5Stg9Ef7nJG3jhX5V4Mz89z+QZFn6K92U6c/ASLJoqoSNi+PLUF8Pzgz/2QRp5bvgvdgqLXgo8yjPyvyrgN0Ri4HpXIXCUNVJWwq/Pwc+SeLoHxX8HhWazi8Hwy8nmXbW9Bn6HXc+PwW3DJ49iJdti81U7mLNIbKCvhkYbBP+eqZsGZWbL9yPlfkdB0Qm1rp/ekZ+YZIhD5f6INbw170I+lF5S5ST9VdOWpeSc+KVRxftpD8Pf+mz57FtPbgrHxdRneWthhMYft8VmX1Zn1md7wiA9YQfLANePfT59p7U+mmfsWl1I7KXaSeJi9YxwdFxVzYZVOszBfSZ8+iT8t8fUY33ml1Nkuan8jS5idQnNG5Tq/T1K+4lNpRuYvURVVlMEe+eiY/2PoSfWwxbTbHVrB0PgZ6Xg09z4Qjz6Brdi5dgXPr+FKRSISCAs2zS+3UWO5m1hJ4C2gRO36iu9+73zEtgHHAScAW4Bp3X53wtCJhqaqEDYs/mzP/6G3YvR2ArhndmNXyLM698AroOfjT7W1FwhTPmXsZMNTdo2aWBcw0s1fdffY+x9wCbHP3Y8zsWuAB4JoGyCvSOKqq9inzGZ8rczodBfnDoedZ0PMM7ngm2C3x3H46u5bkUWO5u7sD0diXWbEP3++w4cBPY59PBMaYmcV+ViQ1lJdC4Suw+Pmg1HcXB9//tMyDaRba7z/v/dEBTyUStrjm3M0sA5gHHAP8wd3n7HdIN+BjAHevMLPtQGdgcwKzitRaTXugm1dy/J6FnFk6jVN3v01r38WWZl34d4tBLGnfj6Ut+rE1IweKCD5mfbqk5VNaxSLJyGpzcm1mHYAXgG+7++J9vr8YGObua2NfrwQGufvm/X5+BDACICcn56QJEybUfwRJIBqN0rZt27BjJES6jGXvOEbPKWXNjiqOaNfsswfdOcrXMLRyBkMq3qYL29hJK2ZkDGJaxpksataHKmt28CevxuldMynokZXgUQTS5XcCGksiDBkyZJ67n1zTcbUqdwAzuwfY5e4P7vO914GfuvssM8sEPgFyDjUtk5eX54WFhbV67WQVrGYoCDtGQqTLWPaOY+/dfp4deXpw8+ZFz8HCCcGe5s2yoPd50O9qOHYYZLUKOXX10uV3AhpLIphZXOUez2qZHKDc3YvNrBVwHsEbpvt6CbgBmAVcCUzTfLskg9ZVUQbtnglPjQrm0XHoMQgu/g0c/yVo3SnsiCINIp4591zg6di8ezNggru/bGb3AXPd/SXgCeAvZvYhsBW4tsESi9Skoowum2bDhD/z6IZXaE45ZB4DQ34CJ1wJnXqFnVCkwcWzWmYhMKCa79+zz+e7gasSG02kFqqq4OM5sPBZWPICfXcXQ5sc3mh9MTNaDWX0t67XjSmkSdEVqpLaNi0PCn3RBCheE2zCddwlLKQP/S7/Dk8//l5wnIpdmhiVu6SeHRtg8aSg1IsWgDWDo4bAkLvhuIuhRVu2RiK6YYU0afrTL6nBHVZMgTmPwqrp4FXB1rgXjIa+V0C7w8JOKJJUVO4SmpouMNqrT9kirt3xJMeVL2VzsxzeanM1M1sOZR1HwEJg4Spg1ed+RtvkSlOncpfQTF6w7pDl26t8BdfueIr+ZfPY2qwTf8r+NtNbn0+lxX+xkLbJlaZK5S6hys/NDi4w2temQpg+CpZOhlYd4bz76HTKrdzavDW3xvm82iZXmjqVuySP4jUQ+SX8+5lg1cvZP4LTb4OW7cNOJpJyVO4SvuhGeOtBmPvnYOXLoG/Amd+DNl3CTiaSslTuEprWVVEujU6Eh1+CijIY8FU4+4fQvnvY0URSnspdDhDvKpa6alG1m2G7JvPwjglkszNYylhwF3Q5psFeU6SpUbnLAWpaxVJXGV7OObte5YroM3So2sa8Fqey4eQfcNF55yf0dURE5S4HUe0qlrqqqgyuJo2MhpI1wd2MzrmHk444LTHPLyIHULlLw3GHZf+Aab+AzYWQeyJc8hAcfY72ehFpYCp3STz3YIuAf90H69+HLsfCVU8H9yFVqYs0CpW7JNbH7walvnoGtO8Bw/8I/a7RJl4ijUx/4yQxNiwNSn35q9AmBy78FZx0I2S2CDuZSJOkcpf6cQ92apxyd3BV6dD/hdO+Ac3bhJ1MpElTuUvdlRbD5Nvgg5fh2Ath+B+gTeewU4kIKnepq7XzYOKNULIezh8V7AGjN0tFkobKXWrHHWY/AlPvgXa5cPPr0P3ksFOJyH5U7hK/0m3w4m1Q+E/Iuxgu/0OwJa+IJB2Vu8Rn7Vx47ibYURTc2u60b2gaRiSJqdzl0Nxh1h/gjXshuyvc8jp0OynsVCJSA5W7HNyurfDiN4O168ddEqyGadUh7FQiEgeVu1Sr956lMPbrEN0Awx6AQSM1DSOSQlTu8nlexaXR57h2x1PQobumYURSlMpdPrNrKycsGkXBjrnMbjmY00aO1zSMSIpSuUtgzWyYeDMdd2zgiexvMqX1pTyrYhdJWc3CDiAhq6qCmQ/BkxdBRnPmD3yAKW0u0/y6SIrTmXtTtnMLvDASPpwK+ZfDZb8nOvv9sFOJSALUeOZuZj3MbLqZLTWzJWZ2ezXHFJjZdjNbEPu4p2HiSsJ89A6MHQz/eRMu/g1c9RS0bB92KhFJkHjO3CuA77v7fDNrB8wzs6nuvnS/42a4+yWJjygJVVUFbz8E00ZBxyPh628Et78TkbRSY7m7exFQFPt8h5ktA7oB+5e7JLvopmAaZuW/4PgvwaUPQ8vssFOJSAMwd4//YLOewFtAX3cv2ef7BcAkYC2wHvgfd19Szc+PAEYA5OTknDRhwoR6RE8e0WiUtm3bhh3jkNoXLyZ/6YNklUdZ0fvrFOVeUO2bptFolP9bkgHAnYNaNXbMhEmF30m8NJbkFNZYhgwZMs/da9yKNe5yN7O2wJvAKHd/fr/HsoEqd4+a2UXAw+7e+1DPl5eX54WFhXG9drKLRCIUFBSEHaN6VVUw8zcw/X7o2CuYW8/td9DDI5EIjxQGt8Z7duTpjRQy8ZL6d1JLGktyCmssZhZXuce1FNLMsgjOzMfvX+wA7l7i7tHY568AWWbWpZaZJdHc4ZX/gWm/CKZhRr55yGIXkfQRz2oZA54Alrn7bw9yzOGx4zCzU2PPuyWRQaWW3OH1u2DuE3DG7XDF49CiXdipRKSRxLNa5gzga8AiM1sQ+95dwBEA7j4WuBL4hplVAKXAtV6byXxJLHd446cw+48w6Btw7s90UZJIExPPapmZwCGbwd3HAGMSFUrqKTIa3v4dnHwLDButYhdpgnSFaor525w1TF6w7qCPXx79O9fteIpprS7gsbVX4I/NrtXzFxeXsr60jPxcLZEUSWXaWybFTF6wjqVFJdU+dkl0EtfteIq3Wp3DY+2/g1vdfr35udkM79+tPjFFJGQ6c09B+bnZBy5TnPMovPonOP6LnPWlxzkro26/2mB5V+ougRSRgM7c08HcJ+HVHwa3wvvSn6COxS4i6UPlnureHw8vfxd6XwBXPgkZWWEnEpEkoHJPZQufg8m3wdFD4epxkNk87EQikiRU7qlqyYvBJmA9B8M14yGrZdiJRCSJqNxT0Mm7Z8GkW6D7KXDd36F567AjiUiSUbmnmP673+OObaMgtz985TlokR477IlIYqncU8nK6Xx/232syewJX52kvdhF5KBU7qli9Ux45jqKMrszqvP90KpD2IlEJIlpQXQqWDMHxl8NHY/kF/ZTos10xi4ih6Yz92S3bh6MvxLaHQ7XT6YkQ2fsIlIzlXsyK/o3/OWL0LoT3PCPoOBFROKgck9WG5bAuMuhRXZQ7O21kZeIxE/lnow2LYdxwyGzJdzwEnQ4IuxEIpJiVO7JZstKePpSwIJi73RU2IlEJAVptUwy2bY6KPaqcrjxn9Cld9iJRCRFqdyTRfHHQbHv2Qk3vgxf6BN2IhFJYZqWSQYlRUGxl26H61+Ew08IO5GIpDiduYctuhHGXQY7N8HXXoSuA8JOJCJpQOWeIDXduLo67aq2c8+WH/GFyiLu7zSKwlcqgFmH/JmlRSW6ebWI1EjTMglyqBtXV8e8ku9t+wWHV6znVx1/RmHzvnH9nG5eLSLx0Jl7AlV74+qDifwSIovgi49yz4nXNmwwEWlydOYehtUz4c0HoN+1oGIXkQagcm9su7bCpFuhYy+4+MGw04hImtK0TGNyD25ovXMTfP0NaNEu7EQikqZU7o3p3ceg8BUY9kvo2j/sNCKSxjQt01iKFsKUu+HYYTDov8NOIyJprsZyN7MeZjbdzJaa2RIzu72aY8zMfm9mH5rZQjMb2DBxU1RZFCbeBK07w/A/glnYiUQkzcUzLVMBfN/d55tZO2CemU1196X7HHMh0Dv2MQh4JPZfAXj1h8Fujzf8A9p0DjuNiDQBNZ65u3uRu8+Pfb4DWAbsfxXNcGCcB2YDHcwsN+FpU9HCCbBgPJz1A+h1ZthpRKSJqNWcu5n1BAYAc/Z7qBvw8T5fr+XA/wNoerashJfvgCNOh7N/FHYaEWlC4l4tY2ZtgUnAd909/uvsP/8cI4ARADk5OUQikbo8TdKJRqMUF5cCfDomqypn4Pwf0bIK5na9hbIZM0NMGL9oNJoWv5d0GQdoLMkq2ccSV7mbWRZBsY939+erOWQd0GOfr7vHvvc57v4Y8BhAXl6eFxQU1DZvUopEInTo0AKAgoLY9gOv/wSiK+Ga8Zze55IQ09VOJBIhHX4v6TIO0FiSVbKPJZ7VMgY8ASxz998e5LCXgOtjq2ZOA7a7e1ECc6aW5VNg1hg45VZIoWIXkfQRz5n7GcDXgEVmtiD2vbuAIwDcfSzwCnAR8CGwC7gp8VFTREkRvPjfcFhfOP8XYacRkSaqxnJ395nAIRdmu7sDtyUqVKoyr4Tnb4XyUrjySchqGXYkEWmitP1AAl0enQCfzIDLxkDOsWHHEZEmTNsPJEjeniVcFf0L9L0SBnw17Dgi0sSp3BMgszzKt4sfYFPGYXDJQ9peQERCp3KvL3fyCsfQsXILv+/wY2ip+5uKSPhU7vU19wlyNs/imXY3sbJ5XthpREQAvaFaP58shtfuYkungfyz+ZfCTiMi8imdudfVnp0w8WZo1YEPjrsdN/1PKSLJQ2fudfXaj2HzcvjaC5R/rDdQRSS56HSzLhZPgvnjYPAdcPSQsNOIiBxA5V5bW/8D//gudD8FhtwVdhoRkWqp3Gujshwm3QIYXPEEZGSFnUhEpFqac6+NaT+HdfPgqqeh45FhpxEROSiducfrwzfg7YfhpJvg+MvDTiMickgq93js2AAv/Dfk9IFho8NOIyJSI03L1KSqCl4YCWU74IZ/QFarsBOJiNRI5V6Tdx6GVdPhkt/BF/qEnUZEJC6aljmUj9+Dab+A/MvhpBvDTiMiEjeV+8GUFsOkm6FdV7j0YW3jKyIppUlMy/xtzhomL1hXq5+5fdtoBu1ey72dH2TFuGWHPLa4uJT1pWXk52q7XxFJDk3izH3ygnUsLSqJ+/h+ZfP4r91vMrHtV1jRPD+un8nPzWZ4/251jSgiklBN4swdgvJ9duTpNR9YWQFj74COPbnmtt9wTWaLGn8kEolQUBDHc4uINJImU+5xm/80bFoGV/8F4ih2EZFk1CSmZeJWWgzTR8GRg6HPpWGnERGpM5X7vmY8CLu2wgWjtDpGRFKayn2vLSth9ljo/xXo2j/sNCIi9aJy32vqPZDRHM7537CTiIjUm8od4D8z4IOX4cw7oN3hYacREak3lXtVJbx+J7TvAad/K+w0IiIJoaWQC/4GnywK7qykHR9FJE007TP3sh3wr/ug+6nQ94qw04iIJEyN5W5mfzazjWa2+CCPF5jZdjNbEPu4J/ExG8jMh2DnxuAGHFr6KCJpJJ5pmaeAMcC4Qxwzw90vSUiixlK8Bt4ZAydcDd1PDjuNiEhC1Xjm7u5vAVsbIUvjmnovWDM4996wk4iIJJy5e80HmfUEXnb3vtU8VgBMAtYC64H/cfclB3meEcAIgJycnJMmTJhQ19y1MnpOKQB3DgreMM3evoyB7/+Y1Udew+peX67380ejUdq2bVvv50kG6TKWdBkHaCzJKqyxDBkyZJ671zzd4O41fgA9gcUHeSwbaBv7/CJgRTzPeeyxx3pjuXrsO3712HeCLyor3R8tcH8wz70smpDnnz59ekKeJxmky1jSZRzuGkuyCmsswFyPo2PrvVrG3UvcPRr7/BUgy8y61Pd5G8yi52D9fDjnXmjeJuw0IiINot7lbmaHmwVLTczs1Nhzbqnv8zaIPTvhjZ9C1wHQ75qw04iINJgaV8uY2TNAAdDFzNYC9wJZAO4+FrgS+IaZVQClwLWxfzokn3f+D3ashyufgGZNe4m/iKS3Gsvd3a+r4fExBEslk1qnyk3w9sOQfzkc+V9hxxERaVBN5vT1uh1PQVUFnPezsKOIiDS4JlHuR+8p5KzSf8Hpt0HHnmHHERFpcOlf7u5cX/Ioxc06wODvhZ1GRKRRpH+5L3mB48qX8my7G6BldthpREQaRXqXe/lumHovH2X2Ynqr88NOIyLSaNK73Gf/Abav4enskbhlhJ1GRKTRpG+579gAM34LeRezpIVueC0iTUv6lvu0n0NFGZz/87CTiIg0uvQs96KF8P5fYdBI6Hx02GlERBpd+pW7O7x+F7TqCGf9IOw0IiKhSL9y/+CfsHoGDLkLWnUIO42ISCjSq9wrymDK3ZBzHJx0U9hpRERCE889VFPHu4/Btv/AVyZBRnoNTUSkNtLnzH3nZnjz13DMedD73LDTiIiEKn3Kffr9sCcKF4wKO4mISOjSo9w3LoN5T8LJN0NOXthpRERCl/rlvnfpY4t2UHBn2GlERJJC6pf7iqmwchqc/SNo0znsNCIiSSG1y72yHKb8BDodDafcGnYaEZGkkdrrBef+GTYvh2ufgczmYacREUkaqXvmvmsrREZDr7Mh78Kw04iIJJXULfc3fwW7t8MF94NZ2GlERJJKapb75hXw3p9gwNfg8L5hpxERSTqpWe5T/hcyW8HQu8NOIiKSlFLuDdVx45/i+hWvMr7dzbw0fiWwssafWVpUQn6ubo4tIk1Hyp25R7M682arc3m1zeVx/0x+bjbD+3drwFQiIskl5c7cv3n1pcClnB12EBGRJJZyZ+4iIlIzlbuISBqqsdzN7M9mttHMFh/kcTOz35vZh2a20MwGJj6miIjURjxn7k8Bww7x+IVA79jHCOCR+scSEZH6qLHc3f0tYOshDhkOjPPAbKCDmeUmKqCIiNReIubcuwEf7/P12tj3REQkJI26FNLMRhBM3ZCTk0MkEmnMl28w0WhUY0ky6TIO0FiSVbKPJRHlvg7osc/X3WPfO4C7PwY8BpCXl+cFBQUJePnwRSIRNJbkki7jAI0lWSX7WBJR7i8B3zKzvwODgO3uXlTTDy1fvjxqZoUJeP1k0AXYHHaIBEmXsaTLOEBjSVZhjeXIeA6qsdzN7BmgAOhiZmuBe4EsAHcfC7wCXAR8COwCboozYKG7nxznsUnNzOZqLMklXcYBGkuySvax1Fju7n5dDY87cFvCEomISL3pClURkTQUZrk/FuJrJ5rGknzSZRygsSSrpB6LBbMqIiKSTjQtIyKShkIpdzMbZmaFsc3GfhxGhvoysx5mNt3MlprZEjO7PexM9WVmGWb2vpm9HHaW+jCzDmY20cw+MLNlZnZ62JnqyszuiP35Wmxmz5hZy7Azxau6TQfNrJOZTTWzFbH/dgwzYzwOMo5fx/58LTSzF8ysQ5gZq9Po5W5mGcAfCDYcyweuM7P8xs6RABXA9909HzgNuC1Fx7Gv24FlYYdIgIeB19z9OOBEUnRMZtYN+A5wsrv3BTKAa8NNVStPceCmgz8G/uXuvYF/xb5Odk9x4DimAn3dvR+wHLizsUPVJIwz91OBD919lbvvAf5OsPlYSnH3InefH/t8B0GBpOyeOmbWHbgYeDzsLPVhZu2Bs4AnANx9j7sXh5uqXjKBVmaWCbQG1oecJ24H2XRwOPB07POngfjvlxmS6sbh7lPcvSL25WyCK/OTShjlnnYbjZlZT2AAMCfcJPXyO+CHQFXYQeqpF7CweFQQAAAB6ElEQVQJeDI2xfS4mbUJO1RduPs64EFgDVBEcPX3lHBT1dth+1zB/glwWJhhEuRm4NWwQ+xPb6jWk5m1BSYB33X3krDz1IWZXQJsdPd5YWdJgExgIPCIuw8AdpIa//Q/QGw+ejjB/2F1BdqY2VfDTZU4sQsgU3q5npn9hGCKdnzYWfYXRrnHvdFYsjOzLIJiH+/uz4edpx7OAC4zs9UE02RDzeyv4Uaqs7XAWnff+6+oiQRln4rOBf7j7pvcvRx4HvivkDPV14a993uI/XdjyHnqzMxuBC4BvuJJuKY8jHJ/D+htZr3MrDnBG0QvhZCjXszMCOZ1l7n7b8POUx/ufqe7d3f3ngS/j2nunpJniO7+CfCxmeXFvnUOsDTESPWxBjjNzFrH/rydQ4q+ObyPl4AbYp/fAEwOMUudmdkwgmnMy9x9V9h5qtPo5R57E+JbwOsEf1AnuPuSxs6RAGcAXyM4y10Q+7go7FACwLeB8Wa2EOgP3B9ynjqJ/etjIjAfWETw9zWpr4rcV2zTwVlAnpmtNbNbgF8C55nZCoJ/mfwyzIzxOMg4xgDtgKmxv/tjQw1ZDV2hKiKShvSGqohIGlK5i4ikIZW7iEgaUrmLiKQhlbuISBpSuYuIpCGVu4hIGlK5i4ikof8Hpg2Tf44vRKsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "  # Define model function to be used to fit to the data above:\n",
    "  def p2(x, *p):\n",
    "    a0,a1,a2 = p\n",
    "    a3=0\n",
    "    return a0+x*(a1+x*(a2+a3*x))\n",
    "\n",
    "  def p3(x, *p):\n",
    "    a0,a1,a2,a3 = p\n",
    "    return a0+x*(a1+x*(a2+a3*x))\n",
    "\n",
    "  # p0 is the initial guess for the fitting coefficients\n",
    "  p0 = [1., 0.5, 0.,0.]\n",
    "\n",
    "  coeff, var_matrix = curve_fit(p3, bin_centres, chi2, p0=p0, sigma=err)\n",
    "  print coeff\n",
    "  # Get the fitted curve\n",
    "  hist_fit = p3(bin_centres, *coeff)\n",
    "  print hist_fit\n",
    "  plt.step(bin_centres, chi2, label='', where='mid')\n",
    "  plt.plot(bin_centres, hist_fit, label='')\n",
    "  plt.grid()\n",
    "  plt.show()\n",
    "\n"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}