keras_g3_cae_true-table.ana.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import uproot4\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as patches\n",
    "from matplotlib.collections import PatchCollection\n",
    "from matplotlib import cm\n",
    "from matplotlib.colors import LogNorm, Normalize\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import itertools\n",
    "\n",
    "from event_display import  table_display, table_display_compare, get_bin_centers\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "pics_prefix = os.path.join(\"plots\", \"v1_\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time open 2.5428295135498047\n",
      "Time split 3.19761061668396\n",
      "Inputs shape new = (50000, 12, 12, 1)\n"
     ]
    }
   ],
   "source": [
    "from geant3_parser import build_true_answers_train_set\n",
    "from geant3_parser import Geant3DataFile\n",
    "\n",
    "def norm_func(e):\n",
    "    return np.float64(np.log(e) / 11)\n",
    "\n",
    "file_name = os.path.join('data', 'shower_geant3_new.dat')\n",
    "data_file = Geant3DataFile(file_name, skip_lines=3)\n",
    "input_data, answers, values = build_true_answers_train_set(data_file, 50000, norm_func=norm_func, rnd_shift=((-2,2), (-2,2)))\n",
    "\n",
    "# Remove x,y\n",
    "\n",
    "# Reshape to events x 11 x 11 x 1 form\n",
    "input_data = np.reshape(input_data, (len(input_data), 11, 11, 1))  # -1 => autodetermine\n",
    "answers = np.reshape(answers, (len(answers), 11, 11, 1))  # -1 => autodetermine\n",
    "\n",
    "# Pad with 1 row and column of zeroes, so it divides by 2, events x 12 x 12 x 1 now\n",
    "input_data = np.pad(input_data, ((0,0), (0,1), (0,1), (0,0)), mode='constant', constant_values=0)\n",
    "answers = np.pad(answers, ((0,0), (0,1), (0,1), (0,0)), mode='constant', constant_values=0)\n",
    "\n",
    "# printout shape\n",
    "print(f\"Inputs shape new = {np.shape(input_data)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Load NN model\n",
    "from keras.models import load_model\n",
    "model = load_model(os.path.join('trained_models', 'g3_true.hd5'))\n",
    "# model = load_model(os.path.join('trained_models', 'g3_true_table_cae.hd5'))\n",
    "\n",
    "# Run ML on data\n",
    "ml_events = model.predict(input_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[1.0004983]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.0001383]]\n",
      "\n",
      " [[1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.0000004]]\n",
      "\n",
      " [[1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]]\n",
      "\n",
      " [[1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]]\n",
      "\n",
      " [[1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]]\n",
      "\n",
      " [[1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]]\n",
      "\n",
      " [[1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]]\n",
      "\n",
      " [[1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]]\n",
      "\n",
      " [[1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]]\n",
      "\n",
      " [[1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]]\n",
      "\n",
      " [[1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.0000015]]\n",
      "\n",
      " [[1.0005186]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.       ]\n",
      "  [1.0010935]]]\n",
      "true event\n",
      "[[[  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]]\n",
      "\n",
      " [[  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]]\n",
      "\n",
      " [[  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]]\n",
      "\n",
      " [[  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]]\n",
      "\n",
      " [[  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [410.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]]\n",
      "\n",
      " [[  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]]\n",
      "\n",
      " [[  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]]\n",
      "\n",
      " [[  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]]\n",
      "\n",
      " [[  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]]\n",
      "\n",
      " [[  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]]\n",
      "\n",
      " [[  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]]\n",
      "\n",
      " [[  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]\n",
      "  [  1.]]]\n",
      "max true 410.00000000000017 max eval 0.0\n",
      "sum true 410.00000000000017 sum eval 0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsEAAAG5CAYAAABvKOlFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhV0lEQVR4nO3dfZBsaV0f8O/v3mF9g4BmpUhYLhdETBDdgBeiQRGEWLDykkqlSqSwRFNupCK1IBF5KWNZiVUGjC5GKmYL1koCQhQBlaAC5Uu0LBZ3V96WRQsWEDaQhSSK0cR1dp/80T1356W7p+f26ek583w+VbfuzDnP+Z1fP33mme/0nJmp1loAAKAnZzbdAAAAHDchGACA7gjBAAB0RwgGAKA7QjAAAN3Z2sRJL7/88nb+/PlNnBpgJTfddNPnWmtfvuk+jpM1Gxizeev2RkLw+fPnc+ONN27i1AArqapPbLqH42bNBsZs3rrtdggAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6M0gIrqr7VdWbqurDVXVrVX3DEHUBWA/rNtC7of5s8quS/Hpr7Z9U1WVJvnigugCsh3Ub6NrKIbiq7pvk8UmemySttTuT3LlqXQDWw7oNMMwrwQ9J8tkkP1dVVya5Kck1rbW/2D2oqq5OcnWSnDt37pJOVFWrdQqwS2tt0y1syqHrtjUbOImGXLdr1WJVdSHJu5M8rrV2Q1W9KsnnW2s/PO+YCxcutBtvvPFSzpXh7uDYb1tttdXurPalrH9VdVNr7cIaGjo2R123rdlqq33ctdddf7y1h1y3h/jBuE8l+VRr7Ybp+29K8ugB6gKwHtZtoHsrh+DW2meSfLKqvmq66UlJPrRqXQDWw7oNMNzr1c9P8vrpTxjfluS7B6oLwHpYt4GuDRKCW2vvTTLqe+QAemLdBnrnL8YBANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0J2tTTdwdNtqq6222ozGWK8NtdXedO111x9r7eGMMASvq+VttdVWu7varN9Yrw211d5k7XXXH3Pt4bgdAgCA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRnsBBcVWer6g+r6m1D1QRgPazZQO+GfCX4miS3DlgPgPWxZgNdGyQEV9UVSb4tyWuGqAfA+lizAZKtgepcm+TFSe4zb0BVXZ3k6iQ5d+7cCqfaXuFYtdVWW21izVZb7RHUXnf9sdYezsohuKqeluSO1tpNVfWEeeNaa9cluS5JLly40C79jEPl9v221VZb7e5q98earbbaY6i97vpjrj2cIW6HeFySZ1TVx5O8Mcm3VNXrBqgLwPCs2QAZIAS31l7aWruitXY+ybOS/GZr7TkrdwbA4KzZABN+TzAAAN0Z9KaN1tpvJ/ntIWsCsB7WbKBnXgkGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtbm27g6LbVVltttRmNsV4baqu96drrrj/W2sMZYQheV8vbaqutdne1Wb+xXhtqq73J2uuuP+baw3E7BAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOjOyiG4qh5UVb9VVR+qqluq6pohGgNgPazbAMnWADW2k7yotXZzVd0nyU1V9c7W2ocGqA3A8KzbQPdWfiW4tfbp1trN07f/PMmtSR64al0A1sO6DTDMK8EXVdX5JI9KcsOMfVcnuTpJzp07t8JZtlc4Vm211Vab3eat29ZstdXedO111x9r7eEMFoKr6t5JfinJC1prn9+/v7V2XZLrkuTChQvt0s80aG7fZVtttdXurnbfFq3b1my11d5k7XXXH3Pt4Qzy2yGq6l6ZLKSvb629eYiaAKyPdRvo3RC/HaKSvDbJra21n1y9JQDWyboNMMwrwY9L8p1JvqWq3jv9d9UAdQFYD+s20L2Vb9porf1ekhqgFwCOgXUbwF+MAwCgQ0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN3Z2nQDR7etttpqq81ojPXaUFvtTdded/2x1h7OCEPwulreVltttburzfqN9dpQW+1N1l53/THXHo7bIQAA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDoziAhuKqeUlV/VFUfqaqXDFETgPWxbgO9WzkEV9XZJK9O8tQkj0jyHVX1iFXrArAe1m2AYV4JfmySj7TWbmut3ZnkjUmeOUBdANbDug10b4gQ/MAkn9z1/qem2wA4mazbQPe2jutEVXV1kquT5Ny5cytU2h6mIbXVVltt5rJmq632pmuvu/5Yaw9niBB8e5IH7Xr/ium2PVpr1yW5LkkuXLjQLv1068rt22qrrXZ3tbt16LptzVZb7U3WXnf9MdcezhC3Q/xBkq+sqodU1WVJnpXkVwaoC8B6WLeB7q0c1Vtr21X1/Ul+I8nZJNe31m5ZuTMA1sK6DTDQ69WttbcnefsQtQBYP+s20Dt/MQ4AgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0Z2vTDRzdttpqq602ozHWa0NttTdde931x1p7OCMMwetqeVtttdXurjbrN9ZrQ221N1l73fXHXHs4bocAAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDsrheCqemVVfbiq3l9Vb6mq+w3UFwBrYN0GmFj1leB3Jnlka+1rk/xxkpeu3hIAa2TdBsiKIbi19o7W2vb03XcnuWL1lgBYF+s2wMSQ9wR/T5JfG7AeAOtl3Qa6tXXYgKp6V5IHzNj18tbaL0/HvDzJdpLXL6hzdZKrk+TcuXOX1OzE9uFD1FZbbbU7NsS6bc1WW+1N1153/bHWHs6hIbi19uRF+6vquUmeluRJrbW2oM51Sa5LkgsXLswdd7hDW75E22qrrXZ3tU+nIdZta7baam+y9rrrj7n2cFbqsqqekuTFSb65tfaXw7QEwLpYtwEmVr0n+GeS3CfJO6vqvVX1swP0BMD6WLcBsuIrwa21hw3VCADrZ90GmPAX4wAA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3tjbdwNFtq6222mozGmO9NtRWe9O1111/rLWHM8IQvK6Wt9VWW+3uarN+Y7021FZ7k7XXXX/MtYfjdggAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQnUFCcFW9qKpaVV0+RD0A1su6DfRu5RBcVQ9K8q1J/mT1dgBYN+s2wDCvBP9UkhcnaQPUAmD9rNtA97ZWObiqnpnk9tba+6rqsLFXJ7k6Sc6dO7fCWbdXOFZttdVWu2/LrtvWbLXV3nTtddcfa+3hHBqCq+pdSR4wY9fLk7wsk2+pHaq1dl2S65LkwoULK7z6sFJuX2BbbbXV7q726TTEum3NVlvtTdZed/0x1x7OoV221p48a3tVfU2ShyTZeTXhiiQ3V9VjW2ufGbRLAJZm3QY43CVH9dbaB5Lcf+f9qvp4kguttc8N0BcAA7NuA9zD7wkGAKA7g9200Vo7P1QtANbPug30zCvBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRna9MNHN222mqrrTajMdZrQ221N1173fXHWns4IwzB62p5W2211e6uNus31mtDbbU3WXvd9cdcezhuhwAAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdWTkEV9Xzq+rDVXVLVb1iiKYAWB/rNkCytcrBVfXEJM9McmVr7a+q6v7DtAXAOli3ASZWfSX4eUl+vLX2V0nSWrtj9ZYAWCPrNkBWfCU4ycOTfFNV/ViS/5fkX7TW/mDWwKq6OsnVSXLu3LkVTrm9wrFqq6222t1bat22Zqut9qZrr7v+WGsP59AQXFXvSvKAGbtePj3+y5J8fZLHJPmFqnpoa63tH9xauy7JdUly4cKFA/uXMaMsAPsMsW5bs4HT7tAQ3Fp78rx9VfW8JG+eLp7vqaq7k1ye5LPDtQjAUVi3AQ636j3Bb03yxCSpqocnuSzJ51asCcD6vDXWbYCV7wm+Psn1VfXBJHcm+a5Zt0IAcGJYtwGyYghurd2Z5DkD9QLAmlm3ASb8xTgAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd2oTfyioqj6b5BOXcOjlGeef99T38dL38eqt7we31r586GZOMmv2aOj7eI2172S8vQ+6bm8kBF+qqrqxtXZh030clb6Pl76Pl76ZZ6xzrO/jpe/jN9beh+7b7RAAAHRHCAYAoDtjC8HXbbqBS6Tv46Xv46Vv5hnrHOv7eOn7+I2190H7HtU9wQAAMISxvRIMAAArE4IBAOjOKENwVT2/qj5cVbdU1Ss23c9RVNWLqqpV1eWb7mUZVfXK6Vy/v6reUlX323RPi1TVU6rqj6rqI1X1kk33s4yqelBV/VZVfWh6TV+z6Z6OoqrOVtUfVtXbNt3LsqrqflX1pum1fWtVfcOmezrNrNnHx5q9ftbs47euNXt0IbiqnpjkmUmubK19dZKf2HBLS6uqByX51iR/sulejuCdSR7ZWvvaJH+c5KUb7meuqjqb5NVJnprkEUm+o6oesdmulrKd5EWttUck+fok/3wkfe+4Jsmtm27iiF6V5Ndba38nyZUZX/+jYc0+dtbs9bNmH7+1rNmjC8FJnpfkx1trf5UkrbU7NtzPUfxUkhcnGc1PI7bW3tFa256+++4kV2yyn0M8NslHWmu3tdbuTPLGTD75nmittU+31m6evv3nmXxwP3CzXS2nqq5I8m1JXrPpXpZVVfdN8vgkr02S1tqdrbU/3WhTp5s1+xhZs9fPmn281rlmjzEEPzzJN1XVDVX1O1X1mE03tIyqemaS21tr79t0Lyv4niS/tukmFnhgkk/uev9TGcnCtKOqzid5VJIbNtzKsq7NJCTcveE+juIhST6b5Oem3xJ8TVV9yaabOsWs2ZtjzV4za/axWNuavTVEkaFV1buSPGDGrpdn0vOXZfItiMck+YWqemg7Ab/r7ZC+X5bJt9VOnEV9t9Z+eTrm5Zl8C+j1x9lbT6rq3kl+KckLWmuf33Q/h6mqpyW5o7V2U1U9YcPtHMVWkkcneX5r7YaqelWSlyT54c22NV7W7ONlzT4ZrNnHZm1r9okMwa21J8/bV1XPS/Lm6QL6nqq6O8nlmXyVsFHz+q6qr8nkK5n3VVUy+fbUzVX12NbaZ46xxZkWzXeSVNVzkzwtyZNOwieuBW5P8qBd718x3XbiVdW9MllMX99ae/Om+1nS45I8o6quSvKFSf5GVb2utfacDfd1mE8l+VRrbeeVmzdlsqByiazZx8uavXnW7GO1tjV7jLdDvDXJE5Okqh6e5LIkn9tkQ4dprX2gtXb/1tr51tr5TJ7QR5+ExfQwVfWUTL518ozW2l9uup9D/EGSr6yqh1TVZUmeleRXNtzToWryWfa1SW5trf3kpvtZVmvtpa21K6bX9LOS/OYIFtNMP+4+WVVfNd30pCQf2mBLp91bY80+Ntbs9bNmH691rtkn8pXgQ1yf5Pqq+mCSO5N81wn/SnfsfibJFyR55/QVkXe31r5vsy3N1lrbrqrvT/IbSc4mub61dsuG21rG45J8Z5IPVNV7p9te1lp7++ZaOvWen+T100+8tyX57g33c5pZs4+XNXv9rNnHby1rtj+bDABAd8Z4OwQAAKxECAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0B0hGACA7gjBAAB0RwgGAKA7QjAAAN0RggEA6I4QDABAd4RgAAC6IwQDANAdIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuiMEAwDQHSEYAIDuCMEAAHRHCAYAoDtCMAAA3RGCAQDojhAMAEB3hGAAALojBAMA0J2tTTdwXKrqbJLHJLnXpnsBAOjMba212zfdxG7VWtvEeY/tpHfddVe2ti5LcneS2rWnZoyuOf/f8/bFLTVr38EX1qsW1Jqxr2ace/+2mnnuRT2cWTB+wXmW7KvarH37j5t1zhnj27553tPP7i37n4/5+w6tVfO62VVjxuWyTP2qg+P3n3fWvtl9zTj3EmP2vD1j/D372oFtMx/jMn3t+hDfP3ezzl11cEmY/dzurb+or72nPTj+YF8He64ZS9Xia2HW+LZvzO6dsx7H/PGz56ntGbjncex7Y09/C8cf3tfe52X++JnnuTiw7Xt/yW0zrtWZx13sYdk52V9j8XxlXz8zzz3rcc96Hmf0Ou88e8YfON+cfvaPqxmPY+Yczu/h4PG7x9Xu/+bUmrFz5uI864Ou5o5ps447sG3RvoPb2szxZ+Ye12bVysHxbWGtBduyc54z88fv+dx/5uC2A33ds6/2b9vzHJyZbpo1vva+n1zMLDfc8NFceeW5vOyl1+fBD35wjtmsj4bT+0rwwfB7Nrvn4GIQm/XEzrpY9u2rWfvq7K76e8ftvSD2jp+5L2cPbDuz4DxnZh13cd+MbTl47ov1p/vOLNi3e9uZdvAxXty3e1vb2Xbwg2Rn2/4xk3F1YNs9Nab76uC+M/uO3z1uz5cF+7bVnlo7x82qtbMtu7btr5W548/s3rfi+N37lhq/q+d7ahwMCmcWjZ/xCfPMvn2zzrn7k/c983swRF18/DPG72zb81G4r8aZGYHh6ONnPcad43b1VYv27dQ42MPFbQvPPb/WrG0Lzz3jPAtrXeK+PePO7Jz77rnn3nPcmVm93j2z5sxaZ2Ycd4nj9wTNi9t2PY4z82vlzN6e9+xfcd+e/WeWOPesfbu3Lay1sy8zjsuefXvenn5Q7P1UWbP/3z1wz8J1Zt//i/bdc6J2cdvZJcefPbAt+2rsCZsXx5/de749287O37Zr3/5te4/bmjF+77ZWWzPGb+0dmyS1d9/ut+vivnu+UX5x2/T/2nXczr7ade4zF7edPbhvp9fW8rZfvTnf/qwn58orz+WHXvqMPPTB12STTt09wXfddVeqzmZrayuTr2rPZn8ABgDg+FRVnv6Mr8vv/f6P5KlXXZlnP+tn8s++70m57ROv2lhPpyYEC78AACfbSQrDow/Bwi8AwLichDA82hAs/AIAjNsmw/DoQrDwCwBwumwiDI8mBAu/AACn23GG4RMfgoVfAIC+HEcYPrEhWPgFAOjbOsPwiQvBwi8AALutIwyfmBAs/AIAsMiQYXjjIVj4BQDgKIYIwxsLwcIvAACrWCUMH3sIrqqzwi8AAEO5lDC8iVeCH5ncPT31mQi/AAAMYScMv/mtL8zNN388P/SDPz937LGH4Nba+z72sY9l8irwXZkE4nbcbQAAcMrcccef5cU/+PN5+lWvzAtf8K/yxjf8/tyxG7kn+Pz582nt7gjDAACsaif8Pu2pr8yjHnU+73nPR/LsZz87Z8+enXvMRn87hDAMAMClmhV+n/PsVy8Mvzs2/ivSEmEYAIDlrRJ+d5yIELxDGAYAYJ4hwu+OExWCdwjDAADsGDL87jiRIXiHMAwA0K91hN8dJzoE7xCGAQD6sc7wu2NrsErHYCcMV21lEoJ3/uhGxR/dAAAYtzvu+LP8xCv/a377Nz+UF77oqrzix39l0OC726hC8I7WtpNEGAYAOAWOM/zuGGUI3iEMAwCM1ybC745Rh+AdwjAAwHhsMvzuOBUheIcwDABwcp2E8LvjVIXgHcIwAMDJcZLC745qbSO/ZuxYT3pPGJ512kWheNa+5cbPemtXQ/uOWu48y9TcX3le3UW1FvWzqP5h5642f+7293PYlyqze5x35l376hKPO6Sfe+pfWo2Z++oSj1ux5pHrz7ocj1yzHdw36FzOqH+EmnvG1fzla3H9Gcct9Xzcc9xyc76gv5knmlF/4XkWLN81681lHvf8MUcZd/DdNnPMnvELns+9Aw+OO/h8HN7f7nELn49DF8G955r9sXO0fgZdbBZd20ddnBYurAf3tZnjj1p///jDP5fv3bZkXwvGz9x3oJ+aMX56vpn7do+pfYcd3HfPPCzOKbVgXF08T+UFP/DUfMe3//Qmwu+cTw2bCcEnVlU9vbX2q5vuY6zM36Uzd0dnzpZnrg5njuYzN7OZl73GNh+j+GMZx+zpm25g5MzfpTN3R2fOlmeuDmeO5jM3s5mXvUY1H0LwQaP5CuaEMn+XztwdnTlbnrk6nDmaz9zMZl72GtV8uB0CAIDueCUYAIDuCMEAAHRHCAYAoDtCMAAA3TmVfzHuJKiqf5Tk25L8jSSvba29Y7MdjYv5W435W8z8rMb8LWZ+FjM/s5mXvY5lPlpro/yX5AuTvCfJ+5LckuRH9+3/qiTv3fXv80lesML5rk9yR5IP7tv+lCR/lOQjSV4y47gvnT55G5+zo8zfrnEvnO7/YJI3JPnCIefvsDk8ifN3hLm7Zjpvt6xy7Z3G629fj2eT/GGStw09P2O8vpadl0WPuZfra9EcJXlQkt9K8qHpx+A1PV1fh8zNUmvYabx+lllvhliTxjIvS6wzH0/ygUxy1I2nbT6O7cJbwxNXSe49ffteSW5I8vULnuTPJHnwjH33T3KffdseNmPc45M8eveTN6370SQPTXLZdEF5xL7j/m2SR296vi5l/pI8MMnHknzR9P1fSPLcoeZvmTk8ifO35Nw9MpMA/MWZfMflXfvnZdm5O43X374efyDJz89ahHu8vpaZl0WPuafr65Br52/t9JbkPkn+eEb/p/b6OmRullnDTuX1c9jH1RJzd6rmZYl15uNJLl9w/KjnY7T3BLeJ/zN9917Tf/N+6fGTkny0tfaJGfu+Oclbq+oLkqSqvjfJv5txvv+W5H/t2/zYJB9prd3WWrszyRuTPHNap6rq3yT5tdbazUd7dOt3hPnbSvJFVbWVSaD77/v2rzJ/yZw5PMnzt+Tc/d0kN7TW/rK1tp3kd5L8431jlpq76TlP1fW3o6quyOTbXa+ZM6S76ytZal4WPeYdp/r6OmyOWmuf3umttfbnSW7N5Av73U7l9bXE3Cyzhp2662eZj6uh1qTk5M/LMvOxhFHPx6jvCa6qs0luSvKwJK9urd0wZ+izMvlW/gGttV+sqock+S9V9YtJvifJP1yyhQcm+eSu9z+V5O9P335+kicnuW9VPay19rNL1jw2h81fa+32qvqJJH+S5P8meUfbd0/OivOXzJ/DEz1/S1x7H0zyY1X1NzOZu6uS3Lh7wBrnLjnh87fLtUlenMkrdQf0en3lkHlZRgfX17VZco6q6nySR2XyiudFp/j6ujaHzM0S6/9pvH6uzeHXzMIxp2xers3h89GSvKOqWpL/0Fq7bs/Okc/HqENwa+2uJH+vqu6X5C1V9cjW2gd3j6mqy5I8I8lLF9R5RVW9Mcm/T/IVu75CXqW3n07y06vWWafD5q+qvjSTr8gekuRPk/xiVT2ntfa6fXW6m7/D5q61duv0K9h3JPmLTO6numtGncHnblr3RM9fklTV05Lc0Vq7qaqeMG9cb9fXsvOyjNN6fR1ljqrq3kl+KZP78j+/f/9pu76O8HF16OfP03T9LDMvm1yTpnWPbV6O8DH0jdMXxO6f5J1V9eHpK7oXjXk+Rns7xG6ttT/N5AcgnjJj91OT3Nxa+x/zjq+qb8rkHs63JPmRI5z69kx++GLHFdNto7Jg/p6c5GOttc+21v46yZuT/IP9x68wf8nI53DRtddae21r7etaa49P8r8zuSdxj57nLsnjkjyjqj6eybfAvqWqXrd/UIdztNS8LOMUz92y1869MgnAr2+tvXlWoVM4R0e6fhatYadsbpaZl57WpKUea2vt9un/d2TyeB+7f8yo56Md4w3pQ/5L8uVJ7jd9+4uS/G6Sp80Y98Yk372gzqMyuVfsKzL5ouANSf71nLHns/eG7q0kt2XySunODd1fvem5GWr+MvmWxC2Z3AtcSf5jkucPNX9jncMjXHv3n/5/LsmHd465lLk7bdffjMf2hMz+IZTurq9l5mXRY+7t+lpw7VSS/5Tk2gXHnurra8HcLLP+n9rr57CPq0Pm7tTNy4LH+iWZ/tDb9O3fT/KU0zQfG78YV3jSvjaTX+vx/kzuv/yXu/a9Pcnfnj5p/zPJfRfUeVySr9n1/r2SfO+McW9I8ukkf53JPSv/dLr9qkxe4ftokpdvel6GnL/p2z+aSYD7YJL/nOQLhpy/Mc7hEebudzP59UzvS/KkS732TuP1N+PxzVuEu7u+Fs3Lvutr7mPu6fqaN0dJvjGT+xnfn3t+VeZVPV1fC+Zm7hrWw/Wz6ONq3pjTPC8LrpOHZvL5a+dX6R3odezzUdMGAACgG6finmAAADgKIRgAgO4IwQAAdEcIBgCgO0IwAADdEYIBAOiOEAwAQHeEYAAAuvP/AQ9JYT4Hs8xGAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 864x648 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "true_event = np.exp(answers[5]*11)\n",
    "ml_event = np.exp(ml_events[5]*11)\n",
    "print(ml_event)\n",
    "print(\"true event\")\n",
    "print(true_event)\n",
    "\n",
    "ml_event[ml_event<2] = 0\n",
    "true_event[true_event<1.5] = 0\n",
    "\n",
    "table_display_compare(true_event, ml_event, cal_size=11)\n",
    "print(f\"max true {np.max(true_event)} max eval {np.max(ml_event)}\")\n",
    "print(f\"sum true {np.sum(true_event)} sum eval {np.sum(ml_event)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def data_to_energies(data):\n",
    "    # We took log(hit) to make hits linear. /11 as a scale. Now reverting it    \n",
    "    denorm_data = np.exp(data*11)\n",
    "\n",
    "    # Get rid of array form of cells\n",
    "    # For ML you have arrays like \n",
    "    # [num_pictures][x_pixels][y_pixels][colors]  where colors=3 for RGB,  colors=1 for BW\n",
    "    # We use colors=1 (we just have an energy deposit in a cell) and we have num_pictures=num_events, [xy]_pixels=[xy]_cells\n",
    "    # No we want to remove this useless [1]\n",
    "    # (len(data), 12, 12, 1) => (len(data), 12, 12)\n",
    "    reshaped = np.reshape(denorm_data, (len(data), 12, 12))\n",
    "\n",
    "    # Sum along (x_pixels, y_pixels) axis => basicly summ everything in one event\n",
    "    sums = np.sum(reshaped, axis=(1,2))\n",
    "    return sums\n",
    "\n",
    "sr_e = data_to_energies(input_data)             # SR = simple reco = energy from original data, summed together\n",
    "ml_e = data_to_energies(ml_events)              # Energy from ML (sum everything)\n",
    "true_e = 1000*np.reshape(values[:,0], len(input_data))    # True energy, 1000 to convert from GeV to\n",
    "\n",
    "# create pandas table... for fun\n",
    "pd = pd.DataFrame({\"sr_e\": sr_e, \"ml_e\": ml_e, \"true_e\":true_e})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3MAAAF0CAYAAACJ02DrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAujElEQVR4nO3df7TddX3n++eLhED9CUhkaBJMrPH2xtuK3iPg0julIBDERWz9cWE6bXToZHUKHTt2poZpF05R10LbVdS10DZXWEVbjRF1yJW0yCBOZ24rJAiiCSLHGIekKim/2i4rEHjfP/YnuDmewznJOWfv/U2ej7X2Ot/v5/v5fvZ7n73P93zf+/P5fr6pKiRJkiRJ3XLEsAOQJEmSJB04kzlJkiRJ6iCTOUmSJEnqIJM5SZIkSeogkzlJkiRJ6iCTOUmSJEnqIJM5SZKkIUiyOsk9ScaTrJ9k+4uS3JzkriRfTrJ0GHFKGl3xPnOSJEmDlWQB8C3gLGA3sBW4sKp29NX5DPCFqro2yRnA26vqV4cSsKSRtHDYATyT448/vpYvXz7sMCTNsdtvv/3vq2rxsOOYDY9P0qFnwMemU4DxqtoJkGQjsAbY0VdnFfDOtnwL8F+na3RRjqqjefbcRippqP6Rh6Y8No10Mrd8+XK2bds27DAkzbEk3x12DLPl8Uk69Az42LQEuK9vfTdw6oQ6XwN+GfgQ8EvAc5O8oKoemKrRo3k2p+bMuY5V0hD9t7puymOT18xJkiSNpv8I/EKSO4BfAPYAT0yslGRdkm1Jtj3Oo4OOUdIQjXTPnCRJ0iFqD7Csb31pK3tKVf0dvZ45kjwHeFNVPTyxoaraAGwAeF6OczIE6TBiz5wkSdLgbQVWJlmRZBFwAbC5v0KS45PsP1e7FLhmwDFKGnEmc5IkSQNWVfuAS4AbgbuBTVW1PcnlSc5v1U4H7knyLeAE4H1DCVbSyHKYpSRJ0hBU1RZgy4Syy/qWrwOuG3RckrrDnjlJkiRJ6iCTOUmSJEnqIJM5SZIkSeogkzlJkiRJ6iCTOUmSJEnqIJM5SZIkSeogkzlJkiRJ6iCTOUmSJEnqIG8arkPe8vU3zEu7u644b17alXR42H9s8lgiSTpYJnOHKBOY+Tdfv2Pw9yxJkqTpOcxSkiRJkjrIZE6SJEmSOshhlhoJ8zlkUZIkSToU2TMnSZIkSR1kz9wQ2RslSZIk6WCZzEkjyNlIJUmSNB2HWUqSJElSB5nMSZIkSVIHmcxJkiRJUgeZzEmSJElSBzkBig6IM3BKkiRJo8FkbgZMYKTRlGQX8I/AE8C+qhpLchzwaWA5sAt4a1U9lCTAh4DXAz8E3lZVX23trAV+vzX73qq6dpCvQ5Ik6WA4zFJS1/1iVZ1cVWNtfT1wc1WtBG5u6wDnAivbYx3wUYCW/L0bOBU4BXh3kmMHGL8kSdJBMZmTdKhZA+zvWbsWeGNf+cer5yvAMUlOBM4BbqqqB6vqIeAmYPWAY5YkSTpgJnOSuqyALya5Pcm6VnZCVX2vLX8fOKEtLwHu69t3dyubqvwnJFmXZFuSbXv37p2r1yBJknRQvGZOUpe9tqr2JHkhcFOSb/ZvrKpKUnP1ZFW1AdgAMDY2NmftSpIkHQx75iR1VlXtaT/vBz5P75q3H7Thk7Sf97fqe4BlfbsvbWVTlUuSJI00kzlJnZTk2Umeu38ZOBv4BrAZWNuqrQWub8ubgV9Lz2nAI2045o3A2UmObROfnN3KJEmSRprDLCV11QnA53t3HGAh8Mmq+qskW4FNSS4Cvgu8tdXfQu+2BOP0bk3wdoCqejDJe4Ctrd7lVfXg4F6GpMNVktX0bpmyAPhYVV0xYftJ9CZyOqbVWV9VWwYdp6TRZTInqZOqaifw8knKHwDOnKS8gIunaOsa4Jq5jlGSppJkAXAVcBa9iZe2JtlcVTv6qv0+sKmqPppkFb0vpZYPPFhJI8thlpIkSYN3CjBeVTur6jFgI71bqPQr4Hlt+fnA3w0wPkkdYM+cJEnS4E12W5RTJ9T5L/Ruv/JbwLOB1w0mNEldYc+cJEnSaLoQ+LOqWkrvmt9PJPmJc7f+e2A+zqMDD1LS8JjMSZIkDd5MbotyEbAJoKr+FjgaOH5iQ1W1oarGqmrsSI6ap3AljSKTOUmSpMHbCqxMsiLJIuACerdQ6fe/aBM6Jfnf6SVzewcapaSRZjInSZI0YFW1D7iE3n0t76Y3a+X2JJcnOb9V+x3g3yb5GvAp4G1tZl5JApwARZIkaSjaPeO2TCi7rG95B/CaQcclqTvsmZMkSZKkDppRMpdkV5KvJ7kzybZWdlySm5Lc234e28qT5MNJxpPcleSVfe2sbfXvTbJ2fl6SJEmSJB36DqRn7her6uSqGmvr64Gbq2olcHNbBzgXWNke64CPQi/5A95N7x4qpwDv3p8ASpIkSZIOzGyGWa4Brm3L1wJv7Cv/ePV8BTgmyYnAOcBNVfVgVT0E3ASsnsXzS5IkSdJha6bJXAFfTHJ7knWt7ISq+l5b/j5wQlteAtzXt+/uVjZVuSRJkiTpAM10NsvXVtWeJC8Ebkryzf6NVVVJ5mSq3JYsrgM46aST5qJJSZIkSTrkzKhnrqr2tJ/3A5+nd83bD9rwSdrP+1v1PcCyvt2XtrKpyic+14aqGquqscWLFx/Yq5EkSZKkw8S0yVySZyd57v5l4GzgG8BmYP+MlGuB69vyZuDX2qyWpwGPtOGYNwJnJzm2TXxydiuTJEmSJB2gmQyzPAH4fJL99T9ZVX+VZCuwKclFwHeBt7b6W4DXA+PAD4G3A1TVg0neA2xt9S6vqgfn7JVIkiRJ0mFk2mSuqnYCL5+k/AHgzEnKC7h4irauAa458DAlSZIkSf1mc2sCSZIkSdKQmMxJkiRJUgeZzEmSJElSB5nMSZIkSVIHmcxJkiRJUgeZzEmSJElSB5nMSZIkSVIHmcxJkiRJUgeZzEmSJElSB5nMSZIkSVIHmcxJkiRJUgeZzEmSJElSB5nMSZIkSVIHmcxJkiRJUgeZzEmSJElSB5nMSZIkSVIHmcxJkiRJUgeZzEmSJElSB5nMSZIkDUGS1UnuSTKeZP0k269Mcmd7fCvJw0MIU9IIWzjsACRJkg43SRYAVwFnAbuBrUk2V9WO/XWq6j/01f8t4BUDD1TSSLNnTpIkafBOAcaramdVPQZsBNY8Q/0LgU8NJDJJnWEyJ0mSNHhLgPv61ne3sp+Q5EXACuBLA4hLUoc4zFKSJGm0XQBcV1VPTLYxyTpgHcDRPGuQcUkaMnvmJEmSBm8PsKxvfWkrm8wFPMMQy6raUFVjVTV2JEfNYYiSRp3JnCRJ0uBtBVYmWZFkEb2EbfPESkl+FjgW+NsBxyepA0zmJEmSBqyq9gGXADcCdwObqmp7ksuTnN9X9QJgY1XVMOKUNNq8Zk6SJGkIqmoLsGVC2WUT1v/LIGOS1C32zEmSJElSB5nMSZIkSVIHmcxJkiRJUgeZzEmSJElSB5nMSZIkSVIHmcxJkiRJUgeZzEnqtCQLktyR5AttfUWSW5OMJ/l0uxkvSY5q6+Nt+/K+Ni5t5fckOWdIL0WSJOmAmMxJ6rp30Lvh7n7vB66sqpcADwEXtfKLgIda+ZWtHklW0bsp78uA1cBHkiwYUOySJEkHzWROUmclWQqcB3ysrQc4A7iuVbkWeGNbXtPWadvPbPXXABur6tGq+g4wDpwykBcgSZI0CyZzkrrsg8DvAk+29RcAD1fVvra+G1jSlpcA9wG07Y+0+k+VT7LP0yRZl2Rbkm179+6dw5chSZJ04EzmJHVSkjcA91fV7YN6zqraUFVjVTW2ePHiQT2tJEnSpBYOOwBJOkivAc5P8nrgaOB5wIeAY5IsbL1vS4E9rf4eYBmwO8lC4PnAA33l+/XvI0mSNLLsmZPUSVV1aVUtrarl9CYw+VJV/QpwC/DmVm0tcH1b3tzWadu/VFXVyi9os12uAFYCtw3oZUiSJB00e+YkHWreBWxM8l7gDuDqVn418Ikk48CD9BJAqmp7kk3ADmAfcHFVPTH4sCVJkg6MyZykzquqLwNfbss7mWQ2yqr6EfCWKfZ/H/C++YtQkiRp7jnMUpIkSZI6yGROkiRJkjrIZE6SJEmSOshkTpIkSZI6aMbJXJIFSe5I8oW2viLJrUnGk3w6yaJWflRbH2/bl/e1cWkrvyfJOXP+aiRJkiTpMHEgPXPvAO7uW38/cGVVvQR4CLiolV8EPNTKr2z1SLKK3lTgLwNWAx9JsmB24UuSJEnS4WlGyVySpcB5wMfaeoAzgOtalWuBN7blNW2dtv3MVn8NsLGqHq2q7wDjTDJ9uCRJkiRpejPtmfsg8LvAk239BcDDVbWvre8GlrTlJcB9AG37I63+U+WT7CNJkiRJOgDTJnNJ3gDcX1W3DyAekqxLsi3Jtr179w7iKSVJkiSpc2bSM/ca4Pwku4CN9IZXfgg4JsnCVmcpsKct7wGWAbTtzwce6C+fZJ+nVNWGqhqrqrHFixcf8AuSJEmSpMPBtMlcVV1aVUurajm9CUy+VFW/AtwCvLlVWwtc35Y3t3Xa9i9VVbXyC9pslyuAlcBtc/ZKJEmSJOkwsnD6KlN6F7AxyXuBO4CrW/nVwCeSjAMP0ksAqartSTYBO4B9wMVV9cQsnl+SJEmSDlsHlMxV1ZeBL7flnUwyG2VV/Qh4yxT7vw9434EGKUmSJEl6ugO5z5wkSZLmSJLVSe5JMp5k/RR13ppkR5LtST456BgljbbZDLOUJEnSQUiyALgKOIve7Zq2JtlcVTv66qwELgVeU1UPJXnhcKKVNKrsmZMkSRq8U4DxqtpZVY/RmzF8zYQ6/xa4qqoeAqiq+wcco6QRZzInSZI0eEuA+/rWd7eyfi8FXprk/0vylSSrBxadpE5wmKUkSdJoWkjvVk6n07s/718n+bmqeri/UpJ1wDqAo3nWgEOUNEz2zEmSJA3eHmBZ3/rSVtZvN7C5qh6vqu8A36KX3D1NVW2oqrGqGjuSo+YtYEmjx2ROkiRp8LYCK5OsSLKI3n15N0+o81/p9cqR5Hh6wy53DjBGSSPOZE6SJGnAqmofcAlwI3A3sKmqtie5PMn5rdqNwANJdgC3AP+pqh4YTsSSRpHXzEmSJA1BVW0Btkwou6xvuYB3tock/QR75iRJkiSpg0zmJEmSJKmDTOYkSZIkqYNM5iRJkiSpg0zmJEmSJKmDTOYkSZIkqYNM5iRJkiSpg0zmJEmSJKmDTOYkSZIkqYNM5iRJkiSpg0zmJEmSJKmDTOYkSZIkqYNM5iRJkiSpg0zmJEmSJKmDTOYkSZIkqYNM5iRJkiSpg0zmJEmSJKmDTOYkSZIkqYNM5iRJkiSpg0zmJEmSJKmDTOYkSZIkqYNM5iRJkiSpg0zmJEmSJKmDTOYkSZIkqYNM5iRJkiSpg0zmJEmSJKmDTOYkSZIkqYNM5iR1UpKjk9yW5GtJtif5g1a+IsmtScaTfDrJolZ+VFsfb9uX97V1aSu/J8k5Q3pJkiRJB8RkTlJXPQqcUVUvB04GVic5DXg/cGVVvQR4CLio1b8IeKiVX9nqkWQVcAHwMmA18JEkCwb5QiQdnpKsbl8ijSdZP8n2tyXZm+TO9vj1YcQpaXSZzEnqpOr5p7Z6ZHsUcAZwXSu/FnhjW17T1mnbz0ySVr6xqh6tqu8A48Ap8/8KJB3O2pdGVwHnAquAC9uXSxN9uqpObo+PDTRISSPPZE5SZyVZkORO4H7gJuDbwMNVta9V2Q0sactLgPsA2vZHgBf0l0+yjyTNl1OA8araWVWPARvpfbkkSTNmMieps6rqiao6GVhK78ToZ+fz+ZKsS7Ityba9e/fO51NJOvTN9IukNyW5K8l1SZYNJjRJXWEyJ6nzquph4Bbg1cAxSRa2TUuBPW15D7AMoG1/PvBAf/kk+0x8ng1VNVZVY4sXL57rlyFJE/2/wPKq+nl6ow+unaxS/xdNj/PoQAOUNFwLp68i6VCxfP0N89LurivOm5d2n0mSxcDjVfVwkp8CzqI3qcktwJvpDVlaC1zfdtnc1v+2bf9SVVWSzcAnk/wx8NPASuC2gb4YSYejab9IqqoH+lY/BnxgsoaqagOwAeB5Oa7mNkxJo8xkTlJXnQhc2yYROALYVFVfSLID2JjkvcAdwNWt/tXAJ5KMAw/Sm8GSqtqeZBOwA9gHXFxVTwz4tUg6/GwFViZZQS+JuwD4V/0VkpxYVd9rq+cDdw82REmjzmROUidV1V3AKyYp38kks1FW1Y+At0zR1vuA9811jJI0laral+QS4EZgAXBN+3LpcmBbVW0G/n2S8+l90fQg8LahBSxpJJnMSZIkDUFVbQG2TCi7rG/5UuDSQcclqTucAEWSJEmSOmjaZC7J0UluS/K1JNuT/EErX5Hk1iTjST6dZFErP6qtj7fty/vaurSV35PknHl7VZIkSZJ0iJtJz9yjwBlV9XLgZGB1ktPozRp3ZVW9BHgIuKjVvwh4qJVf2eqRZBW9i3tfBqwGPtImLpAkSZIkHaBpk7nq+ae2emR7FHAGcF0rvxZ4Y1tew4/vg3IdcGaStPKNVfVoVX0HGGeSSQokSZIkSdOb0TVzSRYkuRO4n95NK78NPFxV+1qV3cCStrwEuA96MzUBjwAv6C+fZJ/+53rqxpd79+494BckSZIkSYeDGSVzVfVEVZ1M74aWpwA/O18BVdWGqhqrqrHFixfP19NIkiRJUqcd0GyWVfUwcAvwauCYJPtvbbCU3g0vaT+XAbTtzwce6C+fZB9JkiRJ0gGYyWyWi5Mc05Z/CjgLuJteUvfmVm0tcH1b3tzWadu/VFXVyi9os12uAFYCt83R65AkSZKkw8pMbhp+InBtm3nyCGBTVX0hyQ5gY5L3AncAV7f6VwOfSDIOPEhvBkuqanuSTcAOYB9wcVU9MbcvR5IkSZIOD9Mmc1V1F/CKScp3MslslFX1I+AtU7T1PuB9Bx6mJEmSJKnfAV0zJ0mSJEkaDSZzkiRJktRBJnOSJEmS1EEmc5IkSZLUQSZzkiRJktRBJnOSJEmS1EEzuc9cZyxff8OwQ5AkSZKkgbBnTpIkSZI6yGROkiRJkjrIZE6SJEmSOshkTpIkSZI66JCaAEWSJB1e+ic/23XFeUOMRJIGz545SZIkSeogkzlJkiRJ6iCTOUmSJEnqIJM5SZIkSeogkzlJkiRJ6iCTOUmSJEnqIJM5SZIkSeogkzlJkqQhSLI6yT1JxpOsf4Z6b0pSScYGGZ+k0WcyJ0mSNGBJFgBXAecCq4ALk6yapN5zgXcAtw42QkldYDInSZI0eKcA41W1s6oeAzYCayap9x7g/cCPBhmcpG4wmZMkSRq8JcB9feu7W9lTkrwSWFZVNzxTQ0nWJdmWZNvjPDr3kUoaWQuHHYAkSZKeLskRwB8Db5uublVtADYAPC/H1fxGJmmU2DMnSZI0eHuAZX3rS1vZfs8F/g/gy0l2AacBm50ERVI/kzlJkqTB2wqsTLIiySLgAmDz/o1V9UhVHV9Vy6tqOfAV4Pyq2jaccCWNIpM5SZKkAauqfcAlwI3A3cCmqtqe5PIk5w83Okld4TVzkiRJQ1BVW4AtE8oum6Lu6YOISVK32DMnSZIkSR1kMidJkiRJHWQyJ0mSJEkdZDInSZIkSR1kMidJkiRJHWQyJ0mSJEkdZDInSZI6Yfn6G1i+/oZhhyFJI8NkTpIkSZI6yGROkiRJkjrIZE6SJEmSOshkTpIkSZI6yGROUiclWZbkliQ7kmxP8o5WflySm5Lc234e28qT5MNJxpPcleSVfW2tbfXvTbJ2WK9JkiTpQJjMSeqqfcDvVNUq4DTg4iSrgPXAzVW1Eri5rQOcC6xsj3XAR6GX/AHvBk4FTgHevT8BlKROyhGze0jqDP9iJXVSVX2vqr7alv8RuBtYAqwBrm3VrgXe2JbXAB+vnq8AxyQ5ETgHuKmqHqyqh4CbgNWDeyWSJEkHx2ROUuclWQ68ArgVOKGqvtc2fR84oS0vAe7r2213K5uqXJIkaaQtHHYAkjQbSZ4DfBb47ar6hyRPbauqSlJz+Fzr6A3R5KSTTpqrZiVp7iRkwYJZNVH7npyjYCTNN3vmJHVWkiPpJXJ/UVWfa8U/aMMnaT/vb+V7gGV9uy9tZVOV/4Sq2lBVY1U1tnjx4rl7IZIkSQfBZE5SJ6XXBXc1cHdV/XHfps3A/hkp1wLX95X/WpvV8jTgkTYc80bg7CTHtolPzm5lkiRJI23aZM7pvyWNqNcAvwqckeTO9ng9cAVwVpJ7gde1dYAtwE5gHPh/gN8EqKoHgfcAW9vj8lYmSZI00mZyzdz+6b+/muS5wO1JbgLeRm/67yuSrKc3/fe7ePr036fSm/771L7pv8eAau1sbrPHSdIBqar/CWSKzWdOUr+Ai6do6xrgmrmLTpKGIwlZdOSs2qh9j89RNJLm27Q9c07/LUmSJEmj54CumXP6b0mSJEkaDTO+NcGgpv926m9JkqSDlyOc3046XMzor32Q03879bckSVq+/gaWr79h2GFI0kibyWyWTv8tSZLmjYmbJB2cmQyz3D/999eT3NnK/jO96b43JbkI+C7w1rZtC/B6etN//xB4O/Sm/06yf/pvcPpvSZIkSTpo0yZzTv8tSdJg7O+d2nXFeUOORJLUBV4hK0mSJEkdNOPZLCVJkmaj/7o4ex8hyWrgQ8AC4GNVdcWE7b9Bb7TTE8A/Aeuqasc0jcKCBfMTsKSRY8+cJEnSgCVZAFwFnAusAi5MsmpCtU9W1c9V1cnAB4A/RpL6mMxJkiQN3inAeFXtrKrHgI3Amv4KVfUPfavPBubknr6SDh0Os5QkSXPOIZXTWgLc17e+Gzh1YqUkFwPvBBYBZ0zWUJJ1wDqAo/PsOQ9U0uiyZ06SJGlEVdVVVfUzwLuA35+izoaqGquqsUU5Gp58cnYPSZ1hMidJkjR4e4BlfetLW9lUNgJvnM+AJHWPyZwkSdLgbQVWJlmRZBFwAbC5v0KSlX2r5wH3DjA+SR3gNXOSJEkDVlX7klwC3Ejv1gTXVNX2JJcD26pqM3BJktcBjwMPAWunbTiBI/yuXjpcmMxJkiQNQVVtAbZMKLusb/kdAw9KUqf41Y0kSZIkdZA9c5IkSYeKepJ67LFZNXHEUUfPav8nH/3RrPaXNHP2zEmSJElSB5nMSZIkSVIHOcxSkiTpEFFPFk/+8z/Pqo0FxxwzuyAcZikNjD1zkiRJktRBJnOSJEmS1EEOs5QkaYQtX3/DU8u7rjhviJFIkkaNyZwkSZKekiM9PZS6wmGWkiRJktRBJnOSJEmS1EH2o0uSJOnHnqxhRyBphuyZkyRJkqQOMpmTJEmSpA5ymKUkSZJ+7IgMOwJJM2TPnCRJkiR1kMmcJEmSJHWQyZwkSZIkdZDXzEmSJOnHjvC7fqkrTOYkSeqY5etveGp51xXnDTGSnlGLR5IOF371IkmSJEkdZM+cJEmSfsxhllJn+NcqSZIkSR1kMidJkiRJHWQyJ0mSJEkd5DVzkiQdApxRUpIOPyZzkiQdouYjwTNplKTRYTInSZKkH3M2S6kz/GuVJEmTWr7+hqf1xGluJVmd5J4k40nWT7L9nUl2JLkryc1JXjSMOCWNLpM5SZIOMyZpw5dkAXAVcC6wCrgwyaoJ1e4Axqrq54HrgA8MNkpJo85kTpIkmeAN3inAeFXtrKrHgI3Amv4KVXVLVf2wrX4FWDrgGCWNOK+ZkyRJGrwlwH1967uBU5+h/kXAX062Ick6YB3A0Txr9pEtXDD7NiQNhMmcJEnSCEvyr4Ex4Bcm215VG4ANAM/LcTXA0CQNmcMsJXVWkmuS3J/kG31lxyW5Kcm97eexrTxJPtwmGrgrySv79lnb6t+bZO0wXoukw84eYFnf+tJW9jRJXgf8HnB+VT06oNgkdYTJnKQu+zNg9YSy9cDNVbUSuLmtQ2+SgZXtsQ74KPSSP+Dd9IY3nQK8e38CKM0Vr0fTJLYCK5OsSLIIuADY3F8hySuAP6WXyN0/sMgWLZrdQ9LAmMxJ6qyq+mvgwQnFa4Br2/K1wBv7yj9ePV8BjklyInAOcFNVPVhVDwE38ZMJoiTNqaraB1wC3AjcDWyqqu1JLk9yfqv2h8BzgM8kuTPJ5imak3SY8po5SYeaE6rqe235+8AJbXmyyQaWPEO5JM2rqtoCbJlQdlnf8usGHpSkTpm2Z85rUiR1VVUVMGeTASRZl2Rbkm179+6dq2YlabRUze4haWBmMszyz/CaFEnd8YM2fJL2c/91JlNNNjCjSQigN2NcVY1V1djixYvnPHBJkqQDMW0y5zUpkjpmM7C/938tcH1f+a+1EQSnAY+04Zg3AmcnObZ9yXR2K5MkSRppB3vN3Lxdk9J/48uTTjrpIMOTdDhI8ingdOD4JLvpjQC4AtiU5CLgu8BbW/UtwOuBceCHwNsBqurBJO+hN7McwOVVNfELLEk6bNRRTqkgdcWs/1qrqpLM2QDp/htfjo2NOfBa0pSq6sIpNp05Sd0CLp6inWuAa+YwNEmSpHl3sLcmmLdrUiRJkiRJ0zvYZM5rUiRJkiRpiKYdZuk1KZIkSYeRPfdPX+cZLFw6u1t17tvt4C1ppqZN5rwmRZIkSZJGz8EOs5QkSZIkDZFzz0qSJOkpTzz88Kz2P+JFPz27ABxmKc2YPXOSJEmS1EEmc5IkSZLUQQ6zlCRJ0pypoxYMOwTpsGHPnCRJkiR1kMmcJEmSJHWQwywlSZojy9ff8NTyrivOG2Ik0vA8uWh2wywzR3FIhwN75iRJkiSpg0zmJEmSJKmDTOYkSZIkqYO8Zk6SJElz5omjZ3d66cmpNHP2zEmSJElSB5nMSZIkSVIH2ZMtSZKkOXPEY0/Mav+Fi4+fdQz79v79rNuQusCeOUmSJEnqIJM5SZIkSeogh1lKkiRpzhzx13fMav+/X/vqWcdwzLUOs9ThwZ45SZIkSeogkzlJkg7C8vU3PPWQDkaS1UnuSTKeZP0k2/9lkq8m2ZfkzcOIUdJoc5ilJEnSgCVZAFwFnAXsBrYm2VxVO/qq/S/gbcB/HHyEw7PvWRl2CFJnmMxJkiQN3inAeFXtBEiyEVgDPJXMVdWutu3JYQQoafQ5zFKSJGnwlgD39a3vbmUHLMm6JNuSbHucR+ckOEndYDInSZLUYVW1oarGqmrsSI4adjiSBshhlpIkSYO3B1jWt760lR32Fv1jzbqNBS976az2f2L7t2YdgzQIJnOSJOmQ0j/D6K4rzhtiJM9oK7AyyQp6SdwFwL8abkiSusZhlpIkSQNWVfuAS4AbgbuBTVW1PcnlSc4HSPKqJLuBtwB/mmT78CKWNIrsmZMkSRqCqtoCbJlQdlnf8lZ6wy8PK8/787+ddRvf+uBps9r/Z3571iFIA2HPnCRJkiR1kMmcJEmSJHWQwywlSZpGRybUkNQ8+7uz66/44ZtmN0wT4Fmf/cqs25CmY8+cJEmSJHWQyZwkSZIkdZDDLCVJ6uOQSqn7/sWVfzOr/Z/3PxfPOoYf3rliVvvv+/Z3Zh2DDn32zEmSJElSB5nMSZIkSVIHmcxJkiRJUgd5zZwk6bC1//o4r42T1O8fXrt31m2su/fWWe3/J7/+5lntf8R//+qs9lc32DMnSTosLF9/w9MmN5EkqetM5iRJkiSpgxxmKUnqrKluI+DwSUnDtmHli2e1/6pt35jV/g8+fsys9ge4/9UPz7oNzS+TOUlSJ3j/N0mSns5hlpIkSZLUQfbMSZIGZrLhj5P1uNkLJ+lw982xx2a1/w9/+RWzjuF//N2GWbcxW+e96txZ7b9vz9/NUSSjyZ45SdKsOVOkJEmDZ8+cJGnGDqTHzElIJEmaXwPvmUuyOsk9ScaTrB/080vSZDw2SZKkrhloz1ySBcBVwFnAbmBrks1VtWOQcUhSP49NXrcmSYeaZ33u1lm3cc7nXj4HkczOH+363Kz2f+mRs093dj/x6Kz2/6UP/afZBfBH1025adDDLE8BxqtqJ0CSjcAa4LA5YZI0kg6rY5NJmiRJh4ZBD7NcAtzXt767lUnSMB2yxyYnJpEk6dCVqhrckyVvBlZX1a+39V8FTq2qS/rqrAPWtdX/DbhnQjPHA38/gHBnowsxQjfiNMa5M0pxvqiqFg87iP1mcmxq5dMdn2ZqVN6LUYhjFGIA4xi1GGA4cYzUselgJNkLfHfYcUxiVD5XB6qLcXcxZuhm3IOKecpj06CHWe4BlvWtL21lT6mqDcCUN7VIsq2qxuYnvLnRhRihG3Ea49zpSpxDMu2xCaY/Ps3UqLwXoxDHKMRgHKMXwyjF0TWjmox29f3sYtxdjBm6GfcoxDzoYZZbgZVJViRZBFwAbB5wDJI0kccmSZLUOQPtmauqfUkuAW4EFgDXVNX2QcYgSRN5bJIkSV008JuGV9UWYMssmpj1EKcB6EKM0I04jXHudCXOoZiDY9OBGJX3YhTiGIUYwDj6jUIMMDpxaG509f3sYtxdjBm6GffQYx7oBCiSJEmSpLkx6GvmJEmSJElzoDPJXJLVSe5JMp5k/RCe/5ok9yf5Rl/ZcUluSnJv+3lsK0+SD7dY70ryyr591rb69yZZO8cxLktyS5IdSbYneceoxZnk6CS3Jflai/EPWvmKJLe2WD7dJqEgyVFtfbxtX97X1qWt/J4k58xVjH3tL0hyR5IvjHCMu5J8PcmdSba1spF5vw9V0x2PDuYzMdkxppVP+n4OIY4/TPLN9tn5fJJjhhFH3/bfSVJJjh9GDEl+q/0+tif5wDB+F0lOTvKV/X//SU6Zrzgyxf+Xtm2qY84gY5jy86n5M9ln8xk+D8kI/A+a6nPUgbg7c/40Sewjfz41Id5unVtV1cg/6E1I8G3gxcAi4GvAqgHH8C+BVwLf6Cv7ALC+La8H3t+WXw/8JRDgNODWVn4csLP9PLYtHzuHMZ4IvLItPxf4FrBqlOJsz/WctnwkcGt77k3ABa38T4B/15Z/E/iTtnwB8Om2vKp9Do4CVrTPx4I5fs/fCXwS+EJbH8UYdwHHTygbmff7UHzM5Hh0MJ8JJjnGTPN+DjqOs4GFbfn9w4qjbVtGb7Ka79K7x8+gfxe/CPw34Ki2/sIhvSdfBM7t+/v+8nzFwRT/X6b6jA4hhkk/nz7m/Xg48udGk8Q88udKU8TdmfOnSWIf+fOpCfHuokPnVvP2xs3xL/XVwI1965cClw4hjuUTDlj3ACe25ROBe9rynwIXTqwHXAj8aV/50+rNQ7zXA2eNapzAs4CvAqfSu+Hi/n/ET73f9E7YXt2WF7Z6mfgZ6K83R7EtBW4GzgC+0J5zpGJsbU52wBnJ9/tQeczkeHSwn4mJx5hp3s+BxjGh3V8C/mJYcQDXAS/f//kfwnuyCXjdCHw2bgT+77Z8IfDJ+Y6jr/x64KypPqODjmGqz6eP+X9M/GxO9nloyyP5P4gRP1eaIuaRPX+aJNZOnE9NiHkXHTq36sowyyXAfX3ru1vZsJ1QVd9ry98HTmjLU8U7sNfRuqZfQe+bm5GKs3W33wncD9xE7xuWh6tq3yTP91QsbfsjwAvmO0bgg8DvAk+29ReMYIwABXwxye1J1rWykXq/D0Ez+X3N5WfiQN/P+Yqj37+h903kwONIsgbYU1Vfm6z9QcQAvBT4v9owoP+e5FVDiuO3gT9Mch/wR/ROeOY9jgn/X2Dyz+igY+jX//nU4HXmf9AonytNEW8Xzp8m+iDdOJ/q16lzq4HfmuBQVVWVpIYdB0CS5wCfBX67qv4hyVPbRiHOqnoCOLld0/B54GeHGc9ESd4A3F9Vtyc5fcjhTOe1VbUnyQuBm5J8s3/jKLzfmjuj8H4m+T1gH/AXQ3juZwH/md6QumFaSG/ozGnAq4BNSV48hDj+HfAfquqzSd4KXA28bj6fcOL/l4nbB/EZfaYYhvn51E8ahWPWVEb9XGkyo37+NFHHzqf6dercqis9c3voXSOx39JWNmw/SHIiQPt5fyufKt55fx1JjqR3cPqLqvrcqMYJUFUPA7fQ62I/Jsn+Lxf6n++pWNr25wMPzHOMrwHOT7IL2EhvaMCHRixGAKpqT/t5P70D+ymM6Pt9CJnJ72suPxMH+n7OVxwkeRvwBuBXqo0bGXAcP0Pveomvtb/PpfSGGv3zAGOA3jesn6ue2+h943z8DNuayzjWAvuP85+h9/f/tOeYyzim+P8Ck39GBx3DVJ9PDd7I/w/q0rnSZEb4/GmizpxP9evcudV8jTedywe9b0F30vsnvv9C6pcNIY7lPH1c+B/y9IshP9CWz+PpF0Pe1sqPA75D70LIY9vycXMYX4CPAx+cUD4ycQKLgWPa8k8B/4PeP9/P8PSLYX+zLV/M0y+G3dSWX8bTL4bdyTxcDAuczo8v2B2pGIFnA8/tW/4bYPUovd+H4mMmx6OD/Uww+XVRU72fg45jNbADWDzM38eEdnfRS6IG/bv4DeDytvxSekNpMoQ47gZOb8tnArfP13vCFP9fpvqMDiGGST+fPub/MfGzOdnnoS2PxP+gqT5HHYi7U+dPk8R/OiN6PjUhzs6dW83rGzfHv9zX05tx6NvA7w3h+T8FfA94nN63shfRG8d7M3AvvZnNjmt1A1zVYv06MNbXzr8Bxtvj7XMc42vpjfO9C7izPV4/SnECPw/c0WL8BnBZK38xcFt7vs/w41nijm7r4237i/va+r0W+z20Gd3m4X3vP/iMVIwtnq+1x/b9fxej9H4fqo/JjkfA5cD5B/uZYJJjzDO9n0OIY5xe0nJne/zJMOKY8D7sol2kPuDfxSLgz+kdw74KnDGk9+S1wO30jgG3Av/nfMXBFP9fpjnmDDKGKT+fPub1WDjy50aTxDzy50pTxN2p86dJ4j+dET2fmhBn586t0p5MkiRJktQhXblmTpIkSZLUx2ROkiRJkjrIZE6SJEmSOshkTpIkSZI6yGROkqQOSvKWJNuTPJlkbJq6C5LckeQLk2z7cJJ/6ls/Ksmnk4wnuTXJ8r5tl7bye5Kc01e+upWNJ1nfV76itTHe2lw0qOd4ht/FryS5K8nXk/xNkpc/U31JGmUmc5Ikjbgkpyf5swnF3wB+GfjrGTTxDnr3pZvY7hi9eyD1uwh4qKpeAlwJvL/VXUXv3k8vo3ffpY+0JHEBvam5zwVWARe2urR9r2xtPdTaHtRzTOU7wC9U1c8B7wE2TFNfkkaWyZwkSR1UVXdX1T3T1UuylN6NbT82oXwBvRvh/u6EXdYA17bl64Azk6SVb6yqR6vqO/TunXRKe4xX1c6qegzYCKxp+5zR2qC1+cZBPUeSxUk+m2Rre7ym/d7+pqoeavW/Aiyd7ncoSaPKZE6SpEPbB+klbE9OKL8E2FxV35tQvoTeDbipqn3AI/RumPtUebO7lU1V/gLg4dZGf/mgnuND9HrsXgW8iQnJbHMR8JeTlEtSJywcdgCSJGlySW4FjgKeAxyX5M626V1VdeMM9n8DcH9V3Z7k9L7ynwbeApw++Z6HhNcBq3qddwA8L8lzquqfAJL8Ir1k7rVDik+SZs1kTpKkEVVVp0LvmjngbVX1tgNs4jXA+UleDxxNL6H5c+BTwEuA8ZbsPCvJeLvubA+wDNidZCHwfOCBvvL9lrYypih/ADgmycLWc9ZffxDPcQRwWlX9aOIvJcnP0+upO7eqHpjqlydJo85hlpIkHaKq6tKqWlpVy+lNLPKlqvrXVXVDVf2Lqlretv2wJXIAm4G1bfnNbZ9q5Re0mShXACuB24CtwMo2q+Si9jyb2z63tDZobV4/wOf4IvBb+38XSU5uP08CPgf8alV962B+r5I0KkzmJEnqoCS/lGQ38GrghiQ3tvKfTrJlFk1fDbwgyTjwTmA9QFVtBzYBO4C/Ai6uqidaj9glwI30Zszc1OoCvAt4Z2vrBa3tQT3HvwfG2m0IdgC/0cova/U+kuTOJNtm8buSpKFK70stSZIkSVKX2DMnSZIkSR1kMidJkiRJHWQyJ0mSJEkdZDInSZIkSR1kMidJkiRJHWQyJ0mSJEkdZDInSZIkSR1kMidJkiRJHfT/AxhmJp5c1u/zAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1080x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# using tuple unpacking for multiple Axes\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15,6))\n",
    "\n",
    "ax1.hist(sr_e, bins=10)\n",
    "ax2.hist(ml_e, bins=100)\n",
    "ax3.hist2d(sr_e, ml_e/sr_e, bins=(15, 100))\n",
    "plt.savefig(pics_prefix + \"e_2dhist.png\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(11, 19, 19)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Build numpy 2d histogram\n",
    "err_bins = np.arange(0.4, 1.4, 0.05)\n",
    "e_bins = np.arange(0, 6000, 500)\n",
    "err_bin_centers = get_bin_centers(err_bins)\n",
    "e_bin_centers = get_bin_centers(e_bins)\n",
    "de = ml_e/sr_e\n",
    "dee_hist, de_bin_edges, e_bin_edges = np.histogram2d(sr_e, de, bins=[e_bins, err_bins])\n",
    "e_bins_num, dee_bins_num = dee_hist.shape\n",
    "e_bins_num, dee_bins_num, len(err_bin_centers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total subplots: 11 plot rows=3 columns=4\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-34-a694d72342d5>:16: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  return H + A * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))\n",
      "<ipython-input-34-a694d72342d5>:16: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return H + A * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))\n",
      "/home/name/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py:807: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "<ipython-input-34-a694d72342d5>:19: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  mean = sum(x * y) / sum(y)\n",
      "/home/name/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py:807: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "<ipython-input-34-a694d72342d5>:19: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  mean = sum(x * y) / sum(y)\n",
      "/home/name/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py:807: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "<ipython-input-34-a694d72342d5>:19: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  mean = sum(x * y) / sum(y)\n",
      "/home/name/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py:807: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "<ipython-input-34-a694d72342d5>:19: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  mean = sum(x * y) / sum(y)\n",
      "/home/name/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py:807: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "<ipython-input-34-a694d72342d5>:19: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  mean = sum(x * y) / sum(y)\n",
      "/home/name/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py:807: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "<ipython-input-34-a694d72342d5>:19: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  mean = sum(x * y) / sum(y)\n",
      "/home/name/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py:807: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "<ipython-input-34-a694d72342d5>:19: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  mean = sum(x * y) / sum(y)\n",
      "/home/name/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py:807: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "<ipython-input-34-a694d72342d5>:19: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  mean = sum(x * y) / sum(y)\n",
      "/home/name/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py:807: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "<ipython-input-34-a694d72342d5>:19: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  mean = sum(x * y) / sum(y)\n",
      "/home/name/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py:807: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "<ipython-input-34-a694d72342d5>:19: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  mean = sum(x * y) / sum(y)\n",
      "/home/name/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py:807: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHAAAANWCAYAAACWJHt/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAABr+ElEQVR4nOz9fbxldV03/r/eMqAlICZjKXejCSka3o1oWUlmCZRil10GeYNmknphN1pXVIaEZVlf08skFa8MtRSpvBkNw36KlbcxXCo6KDShyIDIiIIgKmKf3x9rHdwczpmzz8ze56zZ5/l8PNaDvdft+7P3OS/2vM9aa1drLQAAAAAM1x1WuwAAAAAAdkwDBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DZw2qqrOq6o8mtK8fqqpPVNUNVfVrk9gnMCwyA5gEWQJMi3xhrdDA2U1U1eer6htVdePI9KoJ7v93q+olVXVUVf33vOPcWFU/ssim/zvJ+a21fVprr5xAHadV1bfnHfu6ndjPhqo6v6puqqrPVtVjlli/VdXXR475f0eWVVW9tKqu7aeXVlWNLH9QVV3YH+vCqnrQcuuFSZMZy97P/NfrvfOW/2ZVXV1VX6uq11fVHUeWLTdvzqyqS/rX7ekLLN/pY+1oW9gZsmTZ+3lxVX2qqm6pqtMWWP5LVXV5/5njHVX1fSPLvq+q3t4vu7yqfmncbRc4zgOq6ryq+nJVtQWWf6Cqvjky1ksmVSeMS74sez+L5ssiYzxxZPkk8+XE6v7N87Wq2lZVf1ZV6yZxLPmyNA2c3cvjWmt7j0wnT3DfP5vk3P7xVfOOs3dr7SOLbHdIki07c8DRX/R53jrv2PvtxO7fkuTjSe6W5PeT/ENVrV9imweOHPNXRuaflOQJSR6Y5Igkj0vyq/0Y9kryziR/m+SuSd6Q5J39fFhtMmN5Rl+vnxk57mOTnJLkp9LVf+8kfziy3XLz5pNJnpvk/81fsCvHGmNb2FmyZHxb0/3j758WOO79k7w2yVOTfH+Sm5L81cgqZyS5uV/25CSv7rcZZ9v5vp3knCTP3ME6J4+M9YcmUSfsBPkyvkXzpTd/jG8YWTbJfPneJL+RZP8kD0/3ueO3JnQs+bIEDZwZUFVPr6oPVdXLq+q6qrqsqn60n39FVV0z2oFdYPu7JjksyWIhtth270/yk0le1Xd5D6uqu1TVG6tqe981fWFV3WGBOq9NctrOj3qHdR2W5CFJXtRa+0Zr7R+TfCrJE3dylycmeVlrbVtr7cokL0vy9H7ZUUnWJXlFa+1bfZe+kjx6F4YAUyUzlu3EJH/dWtvSWvtqkhenz4CdyZvW2hmttfcl+eaEj7XotjANsuT2WmtvaK29J8kNCyx+cpJ3tdb+rbV2Y5I/SPI/qmqfqrpzut/lP2it3dha+2CSTen+kbPDbRep45LW2l9n5/4Ruit1wkTIl9tbIl8WNYV8eXVr7d9bazf3/zb6uySP3NVjyZfxaODMjocnuSjdX2XfnOTsJA9Lcp8kT0kXQnsvsu1jk7yvtfad5RywtfboJP+e7/4F59Ikf5nkLun+8vuoJE9L8ox5dV6Wrqv6x8s5XpJU1UV9iC80zXVv75/kstbaaLh9sp+/I/9W3aUHb6uqDSPz799vv9C+7p/kotba6OnJF41xLFhtMuO7mTHn7/oPZ++tqgeOzF8oA76/qu6Wnc+bxezKsXa0LUyLLLl9lizmNr+jrbX/SveX5sP66ZZ+LHMW/f2et+3O+pPqLrH6UFUdNaE6YZLky/j5kiR3r6ovVdXn+obSnfv5086Xn8h3m8W7ciz5MgYNnN3LO+b9Qj9rZNnnWmt/04fUW5MclOT0/qyQ96b7xbjPIvsdPb0wSe65QHjceZFtb1VVeyQ5PsnvttZuaK19Pt3ZKqNd06taa3/ZWrultfaNRXb1pHnHPn9uQWvtiNbafotMz+1X2zvJ9fP2eX2SBbvIvUcl2ZDkvkmuSvLu+u4pkPP3d32SvauqdvJYsFJkxniZkXR/EdqQ7rTp85OcV1X79csWyoCk+z2fdAbsyrF2tC3sClkyfpbsyI5+h/dO8rVFli217c74nXT/GD0gyZlJ3lVVPziBOmG55Mtk8uWzSR6U5B7prgR4aJK/6JdNLV+q6peTbEzy/03gWPJlDBo4u5cnzPuFft3Isi+NPP5GkrTW5s+7XYe6P/3vp5P888jsqxYIj6+PUd/+SfZMcvnIvMvTfTiYc8UY+zln3rF/coxtRt2YZN958/ZNf7ph3fbmXgcnSX8a382tteuS/HqSeyW53yL72zfJjf1ZNzs8FqwymTGm1tqHWndZ0k2ttT9Jcl2SH+8XL5QBSfd7vuy8WcJOH2uJbWFXyJLJ2NHv8HJ/v29dXlVPHsmZ94xTSGvtY/0/Rr/VuvtkfCjJsROoE5ZLvkxAa+3q1trFrbX/bq19Lt29cuYusZ5KvlTVE5L8SZJjWmtf3tVjjbEt0cChOw3x8tba9gns68vpbpp3yMi8g5NcOfL8dt+EsBxVtaVufwf5uek1/Wpbkty7bnvd5gP7+Wm3vbnXFxY5VEt3L5u5/Y1eTnHrvvr/HtGfjTPniOzkjc9gNzCLmbGQpTLgS621azOZvBm108daYlsYmrWSJaNu8ztaVfdOcsckl/bTuqo6dGT9RX+/R7dtrf3dSM4cs5NDXDTzllknDMFazJf5Wr77b/2J50tVHZ3kdeluQv2pkf3u9LHG2JZo4ND9tWWxO5kvS3964zlJ/ri6G1EdkuT56b6haSJaa/dvt7+D/Nz07H6dS5N8IsmLqupOVfXz6Zoq/7jQPqvq/tV9Ffge1V1H+7J0ofyZfpU3Jnl+VR1QVfdM8oIkZ/XLPpDkO0l+raruWFVzd85//6TGDAMzc5lRVQdX1SOraq8+M3473V/cPtTv5o1JnllVh/eXVb0wfQYsN2/64+1VVXdK94+lPfvt5v5/vCvHWnRbGKCZy5Ikqao9+9/vO6T7h8id+kswku5Gn4+rqh/vL904Pcnb+jNhvp7kbUlOr6o7V9UjkxyX5E1LbbtQvdW5U5K9+ud3qqo79o/3q6rH9vPWVdWT093DYu5shV2pE4ZgzeVLVf1kVR3S/+4flORP031TbqaQL4/ut3lia+0/5o1lp48lX8ajgbN7ede8juzbJ7DP+deHJt01ovO7v+N+g9Pzknw93Y28PpjuhmOvX2ZNv7jA8e++zH0cn+56zK+mC7Bf2EEX/vvTXVf7tb7uDUl+rrX27X75a5O8K923vXw63f8QXpskrbWb033F+NPSXXLxy+lOBb15mfXCNMiM8eyT5NXp8uLKJEenOx342iRprf1zkj9Ld2+cL6Q7dfpFI9svJ2+S5L3pTvv+0XT3nvhGun887dKxxtgWdpYsGd/r0v1On5Dk9/vHT02S1tqWJM9O9w+Ya9Jlz+j9LZ6b5Hv6ZW9J8px+m3G2ne+Q/thzf7n+RpJL+sd7JvmjJNvTnV3wvHSfXS7d1TphJ8iX8S2aL0kenOTDfZ0fTvfvll8b2XaS+fIH6W7qfO7IWEYvr9qVY8mXJVRru3TGF7uxqvr+JB9PckDzgwAsQWYAkyBLgGmRL8w6Z+CsbXdJ8gLhBoxJZgCTIEuAaZEvzDRn4AAAAAAMnDNwAAAAAAZOA4eZU1UbqqpV1brVrgWYbfIGmAbZAkyLfNm9aeCsEVX1U1X12aq6qarO778Cb7F1N/Tr3NRv85h5y3+zqq6uqq9V1evnvpZynG13J9V9tfiF/VgurKoH7WDd+1XV+6vq+qraWt3X+84te0RV/UtVfaWqtlfV31fVPUaWn1ZV3553V/p770wdAAAAzCYNnDWgqvZP8rZ0X/n2fUk2p/va7MW8Jd3d2++W7ivq/qGq1vf7emySU5L8VLqvqLx3kj8cZ9vdSVXtleSdSf42yV2TvCHJO/v589dd16/77nSv70lJ/raqDutXuWu6rwrekO41uyHJ38zbzVtba3uPTJcttw4AAABmlwbOFFTn5VV1TX+Wyqeq6gH9srtV1bv6+RdU1R9V1QdHtm1V9dyq+s+quqGqXlxVP1hVH+63OWcn/vH+P5Jsaa39fWvtm0lOS/LAqrrvArUfluQhSV7UWvtGa+0fk3wqyRP7VU5M8tettS2tta8meXGSp4+57fxjnVVVr+nPTrmhqv519MygqrrvyJkrl1TVk0aW/WxVfbx/Ta6oqtMWG3xVPbGqPj/3HozpqCTrkryitfat1tork1SSRy+w7n2T3DPJy1tr32mtvT/Jh5I8NUlaa+/pX/uvtdZuSvKqJI+cQh2wQ0PLpqo6qqq2VdUL+pq+WFXPGFm+6O95fff03xOr6gtV9eWq+v0dHGuwebPUWKrqyKr6SFVd179Grxp9rfttn92/N9dV1RlVVeMeH3aVbJEtMC3yRb4wT2vNNOEpyWOTXJhkv3T/2L5fknv0y87up+9NcniSK5J8cGTblu6Mi32T3D/Jt5K8L92ZLndJcnGSE/t1D05y3Q6mX+rX+z9JXj2vxk8neeICtf98ks/Mm/eqJH/ZP/5kkl8cWbZ/X/Pdltp2gWOdle5slJ9Icse+zg/2y+7cvzbPSNfAeHCSLyc5vF9+VJIfTteEPCLJl5I8oV+2oa9pXb/91iT3GTnujl6zU/p1fjPJe+bV++50X0s4fxwPSHJj+m916+f9S5K3LzLu30jy0ZHnpyW5PslXkmxJ8pyRZWPXYTItNWV42XRUkluSnJ5kzyTHJrkpyV1Hli/1e/66JN+T5IF9TfdbZOxnZbh5s8OxJHlokkf0x9iQ5DNJfmPee/PudO/rwUm2Jzl6tX/eTGtnimyRLSbTlKbIF/liuu3PxWoXMItTurMjLu1/aO8wMn+PJN9O8kMj8/4otw+aR448vzDJ74w8f1m6szGWU89fJ/nTefM+lOTpC6z71Iw0F/p5f5zkrP7xf43+cqULrtb/Yu5w2wWOdVaSs0ee753kO0kOSvKLSf593vqvTXd2z0L7ekW6M2BGA+W30gXzgTvxHv7BaG39vL9LctoC6+6Z5LIk/7t//DNJbk5y3gLrHpGuUfPjI/MOT3cGzx5JfjTJF5OcsNw6TKalpgFm01FJvpFk3ci8a5I8YpH1F/o9P3Bk+X8kOX6RbYecN8sdy29kpEHcb/tjI8/PSf8By2RaiUm2yBaTaVqTfJEvpttOLqGagtZdQvOqJGckuaaqzqyqfZOsT9eFvGJk9SsW2MWXRh5/Y4Hney+zpBvTdZ5H7Zuuo7vcdecvn3t8wzKPM+fW8bfWbkzX3LhnunvFPLw/pe66qrouyZOT/ECSVNXDq7tZ8vaquj7Js9OdDTTqt5Oc0VrbtoPjL2bssbTWvp3kCUl+NsnVSV6QLoRuc9yquk+S9yT59dbav49sf3Fr7arWXX714XTd9V9Ybh2wlAFmU5Jc21q7ZeT5TXP7GfP3/OqFtl3EUPNmh2OpqsOq6t3V3zw+yUsWOP5yXgeYKNkiW2Ba5It84bY0cKaktfbK1tpD051dcVi6X5Dt6U65O3Bk1YN29hhVdXDd9puL5k9P7lfdku60trnt7pzkB/v5821Jcu+q2mdk3gNH1r3NvvrHX2qtXTvGtgu5dfxVtXe6mwBflS6s/rW1tt/ItHdr7Tn96m9OsinJQa21uyR5TbrTKkf9TJIXVtVt7sGzxGv2eyPjPGLetZhHLDaW1tpFrbVHtdbu1lp7bLpTM/9j5JiHJPn/JXlxa+1NO3g9kq4jPXfcZdUBSxlYNi1lnN/z5Rhq3izl1Uk+m+TQ1tq+SX5vgePDqpItt9YoW2DC5MutNcoXNHCmoaoe1nc990zy9STfTPLfrbXvpPs2qNOq6nuru4nw03b2OK21L7TbfnPR/Onv+lXfnuQB1d2g6k5JTk1yUWvtswvs89Ikn0jyoqq6U3Vfh31Ekn/sV3ljkmdW1eFVtV+SF6Y7vW+cbRdybFX9WHU3tXpxukuwrkh3TeRhVfXUqtqznx5WVffrt9snyVdaa9+sqiOT/NIC+96S5OgkZ1TV40fGuKPX7CX9ah9Id4rir1XVHavq5H7++xcaRFUd0Y/5e6vqt5LcY+51qaoD+u1e1Vp7zQLbHldVd63OkUl+Ld31usuuA3ZkgNm0lHF+z5djqHmzlH2SfC3Jjf1785wl1ocVJVtkC0yLfJEv3JYGznTsm+6GTl9NcnmSa5P8eb/s5HQ3zbo6yZvSfe32t6ZZTGtte7pvgvrjvqaHJzl+bnl1dzcfbSwcn2Rjv+6fJvmFfh9prf1zkj9Lcn6SL6Qb34vG2XYRb+63/0q6m109pT/ODem6wsen6zJfneSl6W7glSTPTXJ6Vd2QriF1ziJj/2SSn0vyuqo6Zgd1zN/u5nSXRT0t3Q27fjndjb9uTpKq+r2qes/IJk9Nd++aa9J9xfpPt9bm3tdfSXdGzmmjHeyRbY9Pd/OwG9I1yF7aWnvDOHXAMg0qm8Yw1u/5Mgwyb8bwW+k+eN2Q7v176wT3DZMgW2QLTIt8kS+MqNbdNIhVUlUvTfIDrbUTV7uWlVZVZyXZ1lp74WrXAtzWrGWTvIFhkC3AtMgX1gJn4Kywqrpvf7nN3OUyz0x3iRPAqpFNwDTIFmBa5Atr0brVLmAN2ifd6X33THcX9Jflu/c7AVgtsgmYBtkCTIt8Yc1ZtUuo9t9//7Zhw4ZVOTYwPBdeeOGXW2vrJ7Ev+QKMmlS+yBZglM8uwLQsli+rdgbOhg0bsnnz5tU6PDAwVXX5pPYlX4BRk8oX2QKM8tkFmJbF8sU9cAAAAAAGTgMHAAAAYOA0cAAAAAAGTgMHAAAAYOA0cAAAAAAGbskGTlW9vqquqapPL7K8quqVVbW1qi6qqodMvkwAAACAtWucM3DOSnL0DpYfk+TQfjopyat3vSwAAAAA5izZwGmt/VuSr+xgleOSvLF1Pppkv6q6x6QKBAAAAFjr1k1gHwckuWLk+bZ+3hcnsG/WqDd/7At55yeuXNY2xz3ogPzSww+eUkUAAACwelb0JsZVdVJVba6qzdu3b1/JQ7ObeecnrszFX/za2Otf/MWvLbvhw2yRL8A0yBZgWuQLsFyTOAPnyiQHjTw/sJ93O621M5OcmSQbN25sEzg2M+zwe+ybt/7qj4y17i++9iNTroahky/ANMgWYFrkC7BckzgDZ1OSp/XfRvWIJNe31lw+BQAAADAhS56BU1VvSXJUkv2raluSFyXZM0laa69Jcm6SY5NsTXJTkmdMq1gAAACAtWjJBk5r7YQllrck/2tiFQEAAABwGyt6E2MAAAAAlk8DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgNHAAAAAABk4DBwAAAGDgxmrgVNXRVXVJVW2tqlMWWH5wVZ1fVR+vqouq6tjJlwoAAACwNi3ZwKmqPZKckeSYJIcnOaGqDp+32guTnNNae3CS45P81aQLBQAAAFirxjkD58gkW1trl7XWbk5ydpLj5q3TkuzbP75LkqsmVyIAAADA2rZujHUOSHLFyPNtSR4+b53Tkry3qp6X5M5JHjOR6gAAAACY2E2MT0hyVmvtwCTHJnlTVd1u31V1UlVtrqrN27dvn9ChAeQLMB2yBZgW+QIs1zgNnCuTHDTy/MB+3qhnJjknSVprH0lypyT7z99Ra+3M1trG1trG9evX71zFAAuQL8A0yBZgWuQLsFzjNHAuSHJoVd2rqvZKd5PiTfPW+UKSn0qSqrpfugaONjIAAADABCzZwGmt3ZLk5CTnJflMum+b2lJVp1fV4/vVXpDkWVX1ySRvSfL01lqbVtEAAAAAa8k4NzFOa+3cJOfOm3fqyOOLkzxysqUBAAAAkEzuJsYAAAAATIkGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAjdXAqaqjq+qSqtpaVacsss6TquriqtpSVW+ebJkAAAAAa9e6pVaoqj2SnJHkp5NsS3JBVW1qrV08ss6hSX43ySNba1+tqrtPq2AAAACAtWacM3COTLK1tXZZa+3mJGcnOW7eOs9KckZr7atJ0lq7ZrJlAgAAAKxd4zRwDkhyxcjzbf28UYclOayqPlRVH62qoydVIAAAAMBaN6mbGK9LcmiSo5KckOR1VbXf/JWq6qSq2lxVm7dv3z6hQwPIF2A6ZAswLfIFWK5xGjhXJjlo5PmB/bxR25Jsaq19u7X2uSSXpmvo3EZr7czW2sbW2sb169fvbM0AtyNfgGmQLcC0yBdgucZp4FyQ5NCquldV7ZXk+CSb5q3zjnRn36Sq9k93SdVlkysTAAAAYO1asoHTWrslyclJzkvymSTntNa2VNXpVfX4frXzklxbVRcnOT/Jb7fWrp1W0QAAAABryZJfI54krbVzk5w7b96pI49bkuf3EwAAAAATNKmbGAMAAAAwJRo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcBo4AAAAAAOngQMAAAAwcGM1cKrq6Kq6pKq2VtUpO1jviVXVqmrj5EoEAAAAWNuWbOBU1R5JzkhyTJLDk5xQVYcvsN4+SX49yccmXSQAAADAWjbOGThHJtnaWrustXZzkrOTHLfAei9O8tIk35xgfQAAAABr3jgNnAOSXDHyfFs/71ZV9ZAkB7XW/mlHO6qqk6pqc1Vt3r59+7KLBViMfAGmQbYA0yJfgOXa5ZsYV9UdkvxFkhcstW5r7czW2sbW2sb169fv6qEBbiVfgGmQLcC0yBdgucZp4FyZ5KCR5wf28+bsk+QBST5QVZ9P8ogkm9zIGAAAAGAyxmngXJDk0Kq6V1XtleT4JJvmFrbWrm+t7d9a29Ba25Dko0ke31rbPJWKAQAAANaYJRs4rbVbkpyc5Lwkn0lyTmttS1WdXlWPn3aBAAAAAGvdunFWaq2dm+TcefNOXWTdo3a9LAAAAADm7PJNjAEAAACYLg0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIEbq4FTVUdX1SVVtbWqTllg+fOr6uKquqiq3ldVh0y+VAAAAIC1ackGTlXtkeSMJMckOTzJCVV1+LzVPp5kY2vtiCT/kOTPJl0oAAAAwFo1zhk4RybZ2lq7rLV2c5Kzkxw3ukJr7fzW2k39048mOXCyZQIAAACsXeM0cA5IcsXI8239vMU8M8l7FlpQVSdV1eaq2rx9+/bxqwRYgnwBpkG2ANMiX4DlmuhNjKvqKUk2JvnzhZa31s5srW1srW1cv379JA8NrHHyBZgG2QJMi3wBlmvdGOtcmeSgkecH9vNuo6oek+T3kzyqtfatyZQHAAAAwDhn4FyQ5NCquldV7ZXk+CSbRleoqgcneW2Sx7fWrpl8mQAAAABr15INnNbaLUlOTnJeks8kOae1tqWqTq+qx/er/XmSvZP8fVV9oqo2LbI7AAAAAJZpnEuo0lo7N8m58+adOvL4MROuCwAAAIDeRG9iDAAAAMDkaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADNxYDZyqOrqqLqmqrVV1ygLL71hVb+2Xf6yqNky8UgAAAIA1askGTlXtkeSMJMckOTzJCVV1+LzVnpnkq621+yR5eZKXTrpQAAAAgLVq3RjrHJlka2vtsiSpqrOTHJfk4pF1jktyWv/4H5K8qqqqtdYmVulRR91+3pOelDz3uclNNyXHHnvr7M9f+/V8/Vvfyb/+yLH51x/92exz43X5zdf+/u02/5dH/Xw+svExudtXvpT/9Ten3275u3/6hPy/I34s97j68jzr7/7sdsvffuzT86n7PSyHXHFpTjzn/9xu+dlPeHYu/cEfzmH/9akc/47X3G75G57067n8oMPyw5+5ID9/7lm3W/66J//vfPEHDslDLvpgfu5f3nK75Wc849Rc+33fnx/Z/P/LT//r22+3/OW/+se5Ye/98qgP/1Me9ZFzb7f8T5/3sty8153yMx/4xzziwvffbvnpLzgjSfJz731zHvKpD91m2c173TF/+ry/SJL8j3/6mzzgs5tvs/zGve+Sv/jVlyRJTnj7q3PoZZ++zfKv3HV9XvXLpyVJTjznFTnkiv+8zfL/uc8P5O9P+oPuyUknJZdeetviHvSg5BWv6B4/5Sk59T8uzk0335ItL+t+pP/z3g/IW37+OUmS57/297L3jdffZvNP33dj3vazz0iSnPKXz89eN3/rNsv/3w8/Mu/+mV9Kkpz6sv81/6XJRx/66Lz3qCdmr5u/mVP+8gW3W+5n749z0GEH50WPu//tlg3JH75rSy6+6murXQawTIffc1/5Akzc7pAtiXyB3dGk8mWcS6gOSHLFyPNt/bwF12mt3ZLk+iR3m7+jqjqpqjZX1ebt27fvXMWsCft975457kHzf8wWt//ed8z37jVOP5JZJV+AaZAtwLTIF2C5aqmTZKrqF5Ic3Vr7lf75U5M8vLV28sg6n+7X2dY//69+nS8vtt+NGze2zZs3L7YYWGOq6sLW2sZJ7Eu+AKMmlS+yBRjlswswLYvlyzhn4FyZ5KCR5wf28xZcp6rWJblLkmt3rlQAAAAARo3TwLkgyaFVda+q2ivJ8Uk2zVtnU5IT+8e/kOT9E73/DQAAAMAatuRNQ1prt1TVyUnOS7JHkte31rZU1elJNrfWNiX56yRvqqqtSb6SrskDAAAAwASMddfX1tq5Sc6dN+/UkcffTPI/J1saAAAAAMl4l1ABAAAAsIo0cAAAAAAGTgMHAAAAYOA0cAAAAAAGrlbr276ranuSy5exyf5JvjylclaScQyLcQzHIa219ZPY0U7ky5DNwns736yNadbGk8zemCaSLzOWLcnsvc+zNp7EmIbOZ5eFzdJ7PMeYhm/WxrNgvqxaA2e5qmpza23jatexq4xjWIyDoZvF93bWxjRr40lmc0zc3qy9z7M2nsSY2D3N4ntsTMM3a+NZjEuoAAAAAAZOAwcAAABg4HanBs6Zq13AhBjHsBgHQzeL7+2sjWnWxpPM5pi4vVl7n2dtPIkxsXuaxffYmIZv1sazoN3mHjgAAAAAa9XudAYOAAAAwJqkgQMAAAAwcINr4FTV0VV1SVVtrapTdrDeE6uqVdUgvypsnHFU1ZOq6uKq2lJVb17pGsex1Diq6uCqOr+qPl5VF1XVsatR545U1eur6pqq+vQiy6uqXtmP8aKqeshK1ziOMcbx5L7+T1XVh6vqgStdIztvVjJjzixkx3yzkiVzZMraMGvZksxevsxatiTyZa2QL/JlpcmWJK21wUxJ9kjyX0nunWSvJJ9McvgC6+2T5N+SfDTJxtWue2fGkeTQJB9Pctf++d1Xu+6dHMeZSZ7TPz48yedXu+4FxvETSR6S5NOLLD82yXuSVJJHJPnYate8k+P40ZGfp2OGOg7Tgu/dTGTGMscz+OxYYFwzkSXLGI9M2c2nWcuWZYxpt8qXWcuWMcckX3bzSb7Il4GOZ+azZWhn4ByZZGtr7bLW2s1Jzk5y3ALrvTjJS5N8cyWLW4ZxxvGsJGe01r6aJK21a1a4xnGMM46WZN/+8V2SXLWC9Y2ltfZvSb6yg1WOS/LG1vlokv2q6h4rU934lhpHa+3Dcz9P6ZqbB65IYUzCrGTGnJnIjvlmJUvmyJQ1YdayJZnBfJm1bEnkyxohX+TLipMtw7uE6oAkV4w839bPu1V/WtdBrbV/WsnClmnJcSQ5LMlhVfWhqvpoVR29YtWNb5xxnJbkKVW1Lcm5SZ63MqVN1Djj3N08M103nd3DrGTGnLWSHfPNYpbMkSm7p1nLlmRt5sssZ0siX3ZX8kW+DN1MZsvQGjg7VFV3SPIXSV6w2rVMwLp0pxUeleSEJK+rqv1Ws6CddEKSs1prB6Y7Be9N/fvEKqmqn0wXWL+z2rUwUbOSGXNkx25Cpsy8WcuWRL7sNuTLzJMvrIpZzpah/bBdmeSgkecH9vPm7JPkAUk+UFWfT3ed3qYa3o2MlxpH0nU3N7XWvt1a+1ySS9MF3JCMM45nJjknSVprH0lypyT7r0h1kzPOOHcLVXVEkv+b5LjW2rWrXQ9jm5XMmLNWsmO+mcmSOTJltzdr2ZKszXyZuWxJ5MsMkC/yZZBmPVuG1sC5IMmhVXWvqtoryfFJNs0tbK1d31rbv7W2obW2Id11bY9vrW1enXIXtcNx9N6Rrhudqto/3SmGl61gjeMYZxxfSPJTSVJV90sXYttXtMpdtynJ0/q7sD8iyfWttS+udlHLVVUHJ3lbkqe21i5d7XpYllnJjDlrJTvmm4ksmSNTZsKsZUuyNvNlprIlkS8zQr7Il8FZC9mybrULGNVau6WqTk5yXrq7gL++tbalqk5Psrm1Nv8XaJDGHMd5SX6mqi5O8p0kvz20DuGY43hButMhfzPdTb2e3lp32++hqKq3pPufx/799aovSrJnkrTWXpPu+tVjk2xNclOSZ6xOpTs2xjhOTXK3JH9VVUlyS2ttaGensYBZyYw5s5Id881KlsyRKbNv1rIlmc18mbVsSeTLWiBf5MtqkC1JDfjnDQAAAIAM7xIqAAAAAObRwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA2c3VRVnVVVfzShff1QVX2iqm6oql+bxD6B3Zd8AaZBtgDTIl9YKzRwVlBVfb6qvlFVN45Mr5rg/n+3ql5SVUdV1X/PO86NVfUji2z6v5Oc31rbp7X2ygnUcVpVfXvesa9b5j7uXlVvqaqrqur6qvpQVT183jq/VFWXV9XXq+odVfV9I8u+r6re3i+7vKp+adxtF6jlxKq6sKq+VlXbqurPqmrdyPIPVNU3R8Z6yaTqhHHJl2XtY4f5ssgYTxxZPsl8Ob6qLunruKaq3lBV+07iWPKFSZAty97P+VW1vf/M8MmqOm7e8pX67LJUtvjswqqTL8vez6L5ssgYp/LZZd5276uqVrf9t9GGvtabquqzVfWYedv8ZlVd3Y/j9VV1x3G3RQNnNTyutbb3yHTyBPf9s0nO7R9fNe84e7fWPrLIdock2bIzBxz9ZZ3nrfOOvd8yd713kguSPDTJ9yV5Q5J/qqq9++PeP8lrkzw1yfcnuSnJX41sf0aSm/tlT07y6n6bcbad73uT/EaS/ZM8PMlPJfmteeucPDLWH5qbuSt1wk6QL+PZYb705o/xDSPLJpkvH0ryyNbaXZLcO8m6JKN/QdyVY8kXJkW2jO/Xk9yjtbZvkpOS/G1V3aM/7kp+dlkqWxKfXRgG+TK+RfOlt1KfXdJv9+Qkey6w6C1JPp7kbkl+P8k/VNX6fpvHJjkl3b+nDkmXT384zrb0WmumFZqSfD7JYxZZ9vR0/7N9eZLrklyW5Ef7+VckuSbJiSPrn5Xkj0ae37VfZ48kRyXZNmZN70/ynSTfTHJjksOS3CXJG5NsT3J5khcmucMCdV47WsPIPk9L8rdTeP2+luSh/eOXJHnzyLIfTBdK+yS5c//4sJHlb0ryp0ttO2Ydz0/yrpHnH0jyK4usu9N1mkzLmeTLLr9+o/my6BinmS/pGktvTHLurh5LvpgmNcmWXXrtjuxrPLJ/viqfXeZnSz/vA/HZxbTKk3zZpddufr4sOsZp5Ev/mlya5BFJWpJ1/fzDknxrdNsk/57k2f3jNyd5yciyn0py9TjbmrrJGTjD8vAkF6XrOL45ydlJHpbkPkmekuRV8/5CPOqxSd7XWvvOcg7YWnt0ul+Mub/CXJrkL9P9Ut47yaOSPC3JM+bVeVm6Du0fL+d4SVJVF1XVdYtMC3Z7q+pBSfZKsrWfdf8knxwZx3+lD6Z+uqUfy5xP9tsste04fiK378r/SVV9ubpLMY4amb8rdcIkyZfx8yVJ7l5VX6qqz1XVy6vqzv38iedLVf1YVV2f5IYkT0zyigkcS76wUmTLvGypqndX1TeTfCxdo2Rzv2hFP7vsIFvm+OzC0MmX8fMlWcHPLumaPq9OcvW8+fdPcllr7YZxjtU//v6qutsY2xKXUK2Gd8z7pXzWyLLPtdb+pg+atyY5KMnprbVvtdbem+4X6T6L7Hf0FMEkuecCAXDnRba9VVXtkeT4JL/bWruhtfb5JC9Ld0rdnKtaa3/ZWrultfaNRXb1pHnHPn9uQWvtiNbafotMz12gpn3TdYn/sLV2fT977yTXz1v1+nR/Hdo73V/TF1q21LY7VFW/nGRjkv9vZPbvpAv0A5KcmeRdVfWDE6gTlku+TCZfPpvkQUnukeTR6S61+ot+2cTzpbX2wdZd5nBgkj9P9xfJXT2WfGGSZMsysqW19nPpfteOTfLe1tp/94tW9LPLDrIl8dmF4ZAvk8mXFfvsUlUbkzwyXWNrvqX2NX/53ON9Fli2wzrWKg2clfeEeb+UrxtZ9qWRx99Iktba/Hm36zJX1R2S/HSSfx6ZfdUCAfD1MerbP921jJePzLs83f/g51wxxn7OmXfsnxxjm9upqu9J8q4kH22t/cnIohuT7Dtv9X3T/ZVpR8t2uG1VPbm+e+Ov98yr5QlJ/iTJMa21L8/Nb619rA/0b7XuWtMPpQvVXa0Tlku+LMNi+dJau7q1dnFr7b9ba59LdzPDJ/aLp5Iv/XGvTPc6n72rxxpjW1gO2bJMrbVvt9bek+Rnqurx/ewV/+zS1zI/W3x2YUjkyzItlC8r9dmlf23/Ksmvt9ZuWaC85R5r7rF8GZMGzmx4WJLLW2vbJ7CvLyf5drqbSs05OMmVI8/brhygqrbU7e8CPze9ZmS9OyZ5R5JtSX513m62JHngyLr3TnLHdNdiXppkXVUdOrL+A/Pdy54W3ba19nftuzf+OmZknaOTvC7djdY+tcQQW5KaQJ0wBGsxX+Zr+e7/LyeeL/OsS3ft+S4da4xtYbXNZLYsYPR3ekU/u+ygjoX47MIsWYv5Mt+0Prvsm+5qhLdW1dXpvhQiSbZV1Y/3+7p3VY2eNbPosfrHX2qtXTvGtkQDZ1Ycm+SfJrGj/hTFc5L8cVXtU1WHpLtp799OYv/9Me7fbn8X+Lnp2UlSVXsm+Yd0nfUTR04PnPN3SR5XVT/en/54epK39X9N+nqStyU5varuXFWPTHJcusskdrjtQvVW1aP7bZ7YWvuPecv2q6rHVtWdqmpddXdj/4l8t+O/K3XCEKy5fKmqn6yqQ6pzUJI/TfLOfv+TzpcnV9XB/eND0l07/75dPZZ8YTcwi9ly36o6pqq+p6r2rKqnpPtM8K/9blbys8ui2eKzC2vAmsuXFfzscn2Se6a7XOtB+e6Zew9N8rHW3WfnE0le1GfMzyc5Isk/9uu9Mckzq+rwqtov3Q2hz+rrXGpbooGzGt41r6v69gnsc/41nkl3nef8Du4TF9p4Ac9L8vV0N+P6YLqbhr1+mTX94gLHv/sytv/RJD+X5GeSXDeyjx9PktbaliTPThc416S7NnL0GtHnJvmeftlbkjyn32acbef7g3Q3Ljt3pI65U5T3TPe1nNvTdeifl+5U0Et3tU7YCfJlPDvMlyQPTvLhvs4PJ/lUkl8b2X6S+XJ4kg9X1dfTXcJwSZLR6/935VjyhUmRLeOpdN82c026zwW/nuQXW2v/L1nxzy47yhafXRgS+TKeHeZLVuizS+tcPTf1tSTdWTQ394+PT3eWzlfTNZJ+Ye5sqNbaPyf5syTnJ/lCusvRXjRyiEW3pVOt7dIZX6yyqvr+JB9PckDzZgITJF+AaZAtwLTIF2adM3B2f3dJ8gIBBUyBfAGmQbYA0yJfmGnOwAEAAAAYOGfgAAAAAAycBg4AAADAwGngMPOqakNVtapat9q1ALNFvgDTIFuAaZEvuzcNHAAAAICB08ABAAAAGDgNnIGpzsur6pqq+lpVfaqqHtAvu1tVvauff0FV/VFVfXBk21ZVz62q/6yqG6rqxVX1g1X14X6bc6pqr2XWc1RVbauqF/Q1fbGqnjGy/Ger6uP9/q+oqtNGls2dnndiVX2hqr5cVb+/g2OdVVWvqap/6ev/16o6ZGT5fftlX6mqS6rqSePUscBxnlhVn597Xcd8HXY4lqo6sqo+UlXX9a/Rq0Zf637bZ/fvzXVVdUZV1bjHh0mQL/IFpkG2yBaYFvkiX5intWYa0JTksUkuTLJfkkpyvyT36Jed3U/fm+TwJFck+eDIti3JO5Psm+T+Sb6V5H1J7p3kLkkuTnJiv+7BSa7bwfRL/XpHJbklyelJ9kxybJKbktx1ZPkPp2sGHpHkS0me0C/b0Nf0uiTfk+SBfU33W2TsZyW5IclPJLljkv8zN74kd+7H+4wk65I8OMmXkxy+jDrW9dtvTXKfkePu6HU4ZZyxJHlokkf0x9iQ5DNJfmPee/PudO/rwUm2Jzl6tX/eTGtrinyRLybTFKbIFtliMk1pinyRL6bb/lysdgGmeW9I8ugkl/Y/8HcYmb9Hkm8n+aGReX+U24fUI0eeX5jkd0aevyzJK5ZZz1FJvpFk3ci8a5I8YpH1X5Hk5f3juV/sA0eW/0eS4xfZ9qwkZ4883zvJd5IclOQXk/z7vPVfm+RFy6jjt9IF9YELbbPE67DcsfxGkrfPe29+bOT5OXMBaDKt1CRf5IvJNI1JtsgWk2lak3yRL6bbTi6hGpjW2vuTvCrJGUmuqaozq2rfJOvTdTCvGFn9igV28aWRx99Y4PneO1HWta21W0ae3zS3n6p6eFWdX1Xbq+r6JM9Osv+87a9eaNtF3Dqm1tqNSb6S5J5JDkny8P4Uu+uq6rokT07yA8uo47eTnNFa2zbWqBe24Fiq6rCqendVXV1VX0vykgWOv5zXASZOvsgXmAbZIltgWuSLfOG2NHAGqLX2ytbaQ9OdCnhYul+u7elO1ztwZNWDdvYYVXVwVd24g+nJY+7qzUk2JTmotXaXJK9Jd3rjzrp1TFW1d5LvS3JVuvD619bafiPT3q215yyjjp9J8sKqeuLozCVeh98bs+5XJ/lskkNba/sm+b0Fjg+rTr7cWqN8gQmSLbfWKFtgwuTLrTXKFzRwhqaqHtZ3TPdM8vUk30zy36217yR5W5LTqup7q+q+SZ62s8dprX2h/yVfbPq7MXe1T5KvtNa+WVVHJvmlna2pd2xV/Vh1N7l6cZKPttauSHeN5GFV9dSq2rOfHlZV91tGHVuSHJ3kjKp6/NzMJV6Hl4xZ9z5Jvpbkxv69ec4S68OKky/yBaZBtsgWmBb5Il+4LQ2c4dk33c2gvprk8iTXJvnzftnJ6W64dXWSNyV5S7qbRa2m5yY5vapuSHJquusXd8Wbk7wo3emBD03ylCRprd2Qrkt8fLqu89VJXpruhl5j19Fa+2SSn0vyuqo6ZhdrHfVb6YLxhnTv31snuG+YFPkiX2AaZItsgWmRL/KFEdW6mwaxG6qqlyb5gdbaiatdyyRU1VlJtrXWXrjatcBaJ1+AaZAtwLTIF9YCZ+DsRqrqvlV1RHWOTPLMJG9f7bqA3Z98AaZBtgDTIl9Yi9atdgEsyz7pTg28Z7o7qL8syTtXtSJgVsgXYBpkCzAt8oU1Z9Uuodp///3bhg0bVuXYwPBceOGFX26trZ/EvuQLMGpS+SJbgFE+uwDTsli+rNoZOBs2bMjmzZtX6/DAwFTV5ZPal3wBRk0qX2QLMMpnF2BaFssX98ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGLglGzhV9fqquqaqPr3I8qqqV1bV1qq6qKoeMvkyAQAAANaucc7AOSvJ0TtYfkySQ/vppCSv3vWyAAAAAJizZAOntfZvSb6yg1WOS/LG1vlokv2q6h6TKhAAAABgrZvEPXAOSHLFyPNt/bzbqaqTqmpzVW3evn37BA4N0JEvwDTIFmBa5AuwXCt6E+PW2pmttY2ttY3r169fyUMDM06+ANMgW4BpkS/Ack2igXNlkoNGnh/YzwMAAABgAibRwNmU5Gn9t1E9Isn1rbUvTmC/AAAAACRZt9QKVfWWJEcl2b+qtiV5UZI9k6S19pok5yY5NsnWJDcleca0igUAAABYi5Zs4LTWTlhieUvyvyZWEQAAAAC3saI3MQYAAABg+TRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZurAZOVR1dVZdU1daqOmWB5QdX1flV9fGquqiqjp18qQAAAABr05INnKraI8kZSY5JcniSE6rq8HmrvTDJOa21Byc5PslfTbpQAAAAgLVqnDNwjkyytbV2WWvt5iRnJzlu3jotyb7947skuWpyJQIAAACsbeM0cA5IcsXI8239vFGnJXlKVW1Lcm6S5y20o6o6qao2V9Xm7du370S5AAuTL8A0yBZgWuQLsFyTuonxCUnOaq0dmOTYJG+qqtvtu7V2ZmttY2tt4/r16yd0aAD5AkyHbAGmRb4AyzVOA+fKJAeNPD+wnzfqmUnOSZLW2keS3CnJ/pMoEAAAAGCtG6eBc0GSQ6vqXlW1V7qbFG+at84XkvxUklTV/dI1cJwHCAAAADABSzZwWmu3JDk5yXlJPpPu26a2VNXpVfX4frUXJHlWVX0yyVuSPL211qZVNAAAAMBasm6clVpr56a7OfHovFNHHl+c5JGTLQ0AAACAZHI3MQYAAABgSjRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZurAZOVR1dVZdU1daqOmWRdZ5UVRdX1ZaqevNkywQAAABYu9YttUJV7ZHkjCQ/nWRbkguqalNr7eKRdQ5N8rtJHtla+2pV3X1aBQMAAACsNeOcgXNkkq2ttctaazcnOTvJcfPWeVaSM1prX02S1to1ky0TAAAAYO0ap4FzQJIrRp5v6+eNOizJYVX1oar6aFUdvdCOquqkqtpcVZu3b9++cxUDLEC+ANMgW4BpkS/Ack3qJsbrkhya5KgkJyR5XVXtN3+l1tqZrbWNrbWN69evn9ChAeQLMB2yBZgW+QIs1zgNnCuTHDTy/MB+3qhtSTa11r7dWvtckkvTNXQAAAAA2EXjNHAuSHJoVd2rqvZKcnySTfPWeUe6s29SVfunu6TqssmVCQAAALB2LdnAaa3dkuTkJOcl+UySc1prW6rq9Kp6fL/aeUmuraqLk5yf5Ldba9dOq2gAAACAtWTJrxFPktbauUnOnTfv1JHHLcnz+wkAAACACZrUTYwBAAAAmBINHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICBG6uBU1VHV9UlVbW1qk7ZwXpPrKpWVRsnVyIAAADA2rZkA6eq9khyRpJjkhye5ISqOnyB9fZJ8utJPjbpIgEAAADWsnHOwDkyydbW2mWttZuTnJ3kuAXWe3GSlyb55gTrAwAAAFjzxmngHJDkipHn2/p5t6qqhyQ5qLX2TzvaUVWdVFWbq2rz9u3bl10swGLkCzANsgWYFvkCLNcu38S4qu6Q5C+SvGCpdVtrZ7bWNrbWNq5fv35XDw1wK/kCTINsAaZFvgDLNU4D58okB408P7CfN2efJA9I8oGq+nySRyTZ5EbGAAAAAJMxTgPngiSHVtW9qmqvJMcn2TS3sLV2fWtt/9bahtbahiQfTfL41trmqVQMAAAAsMYs2cBprd2S5OQk5yX5TJJzWmtbqur0qnr8tAsEAAAAWOvWjbNSa+3cJOfOm3fqIusetetlAQAAADBnl29iDAAAAMB0aeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADJwGDgAAAMDAaeAAAAAADNxYDZyqOrqqLqmqrVV1ygLLn19VF1fVRVX1vqo6ZPKlAgAAAKxNSzZwqmqPJGckOSbJ4UlOqKrD56328SQbW2tHJPmHJH826UIBAAAA1qpxzsA5MsnW1tplrbWbk5yd5LjRFVpr57fWbuqffjTJgZMtEwAAAGDtGqeBc0CSK0aeb+vnLeaZSd6z0IKqOqmqNlfV5u3bt49fJcAS5AswDbIFmBb5AizXRG9iXFVPSbIxyZ8vtLy1dmZrbWNrbeP69esneWhgjZMvwDTIFmBa5AuwXOvGWOfKJAeNPD+wn3cbVfWYJL+f5FGttW9NpjwAAAAAxjkD54Ikh1bVvapqryTHJ9k0ukJVPTjJa5M8vrV2zeTLBAAAAFi7lmzgtNZuSXJykvOSfCbJOa21LVV1elU9vl/tz5PsneTvq+oTVbVpkd0BAAAAsEzjXEKV1tq5Sc6dN+/UkcePmXBdAAAAAPQmehNjAAAAACZPAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4DRwAAAAAAZOAwcAAABg4MZq4FTV0VV1SVVtrapTFlh+x6p6a7/8Y1W1YeKVAgAAAKxRSzZwqmqPJGckOSbJ4UlOqKrD5632zCRfba3dJ8nLk7x00oUCAAAArFXrxljnyCRbW2uXJUlVnZ3kuCQXj6xzXJLT+sf/kORVVVWttTaJIv/wXVty8VVfm8SugBVy+D33zYsed//VLmOHZAvsnuQLMA27Q7Yk8gV2R5PKl3EuoTogyRUjz7f18xZcp7V2S5Lrk9xt/o6q6qSq2lxVm7dv375zFQMsQL4A0yBbgGmRL8ByjXMGzsS01s5McmaSbNy4ceyzc3aHTjiwunYmX2QLsBSfXYBpkS/Aco1zBs6VSQ4aeX5gP2/BdapqXZK7JLl2EgUCAAAArHXjNHAuSHJoVd2rqvZKcnySTfPW2ZTkxP7xLyR5/6TufwMAAACw1i15CVVr7ZaqOjnJeUn2SPL61tqWqjo9yebW2qYkf53kTVW1NclX0jV5AAAAAJiAse6B01o7N8m58+adOvL4m0n+52RLAwAAACAZ7xIqAAAAAFaRBg4AAADAwGngAAAAAAycBg4AAADAwNVqfdt3VW1PcvmqHHzy9k/y5dUuYsKMafcwS2M6pLW2fhI7ki+DN2tjmrXxJLM3ponky4xlSzJ77/OsjScxpqHz2WVhs/QezzGm4Zu18SyYL6vWwJklVbW5tbZxteuYJGPaPczimLitWXyPZ21MszaeZDbHxO3N2vs8a+NJjInd0yy+x8Y0fLM2nsW4hAoAAABg4DRwAAAAAAZOA2cyzlztAqbAmHYPszgmbmsW3+NZG9OsjSeZzTFxe7P2Ps/aeBJjYvc0i++xMQ3frI1nQe6BAwAAADBwzsABAAAAGDgNHAAAAICB08BZhqo6uqouqaqtVXXKIus8qaourqotVfXmla5xuZYaU1UdXFXnV9XHq+qiqjp2NeocV1W9vqquqapPL7K8quqV/XgvqqqHrHSNyzXGmJ7cj+VTVfXhqnrgStfIrpu1fJm1bElmL19ky9owa9mSzF6+zFq2JPJlrZAv8mWlyZYkrTXTGFOSPZL8V5J7J9krySeTHD5vnUOTfDzJXfvnd1/tuicwpjOTPKd/fHiSz6923UuM6SeSPCTJpxdZfmyS9ySpJI9I8rHVrnkCY/rRkZ+5Y3aHMZlu9x7OVL7MYrb0dc5UvsiW2Z9mLVuWMabdKl9mLVvGHJN82c0n+SJfBjqemc8WZ+CM78gkW1trl7XWbk5ydpLj5q3zrCRntNa+miSttWtWuMblGmdMLcm+/eO7JLlqBetbttbavyX5yg5WOS7JG1vno0n2q6p7rEx1O2epMbXWPjz3M5fko0kOXJHCmKRZy5eZy5Zk9vJFtqwJs5YtyQzmy6xlSyJf1gj5Il9WnGxxCdVyHJDkipHn2/p5ow5LclhVfaiqPlpVR69YdTtnnDGdluQpVbUtyblJnrcypU3NOGPenT0zXRed3cus5ctazJZktvNFtuyeZi1bkrWZL7OcLYl82V3JF/kydDOZLRo4k7Uu3amCRyU5Icnrqmq/1SxoAk5IclZr7cB0p9i9qar83AxQVf1kuqD6ndWuhamYtXyRLbsJ2TLzZi1bEvmy25AvM0++sCpmOVv8sI3vyiQHjTw/sJ83aluSTa21b7fWPpfk0nShNVTjjOmZSc5JktbaR5LcKcn+K1LddIwz5t1OVR2R5P8mOa61du1q18OyzVq+rMVsSWYwX2TLbm/WsiVZm/kyc9mSyJcZIF/kyyDNerZo4IzvgiSHVtW9qmqvJMcn2TRvnXek6zCnqvZPd9rgZStY43KNM6YvJPmpJKmq+6ULqe0rWuVkbUrytP6O649Icn1r7YurXdSuqKqDk7wtyVNba5eudj3slFnLl7WYLcmM5YtsmQmzli3J2syXmcqWRL7MCPkiXwZnLWTLutUuYHfRWrulqk5Ocl66O5S/vrW2papOT7K5tbapX/YzVXVxku8k+e0hd/3GHNML0p3u+Jvpbtr19Na623oPUVW9Jd3/KPbvr019UZI9k6S19pp016oem2RrkpuSPGN1Kh3fGGM6NcndkvxVVSXJLa21jatTLTtj1vJlFrMlmb18kS2zb9ayJZnNfJm1bEnky1ogX+TLapAtSQ345w0AAACAuIQKAAAAYPA0cAAAAAAGTgMHAAAAYOA0cAAAAAAGTgMHAAAAYOA0cAAAAAAGTgMHAAAAYOA0cAAAAAAGTgMHAAAAYOA0cAAAAAAGTgMHAAAAYOA0cAAAAAAGTgMHAAAAYOA0cAAAAAAGTgMHAAAAYOA0cHZTVXVWVf3RhPb1Q1X1iaq6oap+bRL7BHZf8gWYFvkCADtPA2cFVdXnq+obVXXjyPSqCe7/d6vqJVV1VFX997zj3FhVP7LIpv87yfmttX1aa6+cQB2nVdW35x37ul3Y36Oqqs3/wFdVv1lVV1fV16rq9VV1x5FlG6rq/Kq6qao+W1WPGXfbHdSxV1V9pqq2zZvfqurrI2P9vyPLqqpeWlXX9tNLq6pGlj+oqi7s67ywqh60Ey8RyJed39/t8qWqnl5V35l3jKNGlk8sXxYZz71Hli+aEfKFlSJflr2f+a/Xe+ctXzQjVjhffH4B2M1o4Ky8x7XW9h6ZTp7gvn82ybn946vmHWfv1tpHFtnukCRbduaAVbVukUVvnXfs/XZy/3sm+T9JPjZv/mOTnJLkp9LVf+8kfziyyluSfDzJ3ZL8fpJ/qKr1Y267mN9Osn2RZQ8cGeuvjMw/KckTkjwwyRFJHpfkV/s69kryziR/m+SuSd6Q5J39fNgZ8mV5+18wX3ofmXeMD4wsm3S+zB/PZf2+lsoI+cJKki/LM/p6/czIcVf688uC+TLC5xeA3YgGzkBU9xffD1XVy6vquqq6rKp+tJ9/RVVdU1Un7mD7uyY5LMliH3IW2+79SX4yyav6v74cVlV3qao3VtX2qrq8ql5YVXdYoM5rk5y286MeywuSvDfJZ+fNPzHJX7fWtrTWvprkxUme3td4WJKHJHlRa+0brbV/TPKpJE9catvFVNW9kjwlyZ8ss/4Tk7ystbattXZlkpeNHOuoJOuSvKK19q3+r4eV5NHLPAbskHxZ1GL5sqhp5MsOHJUdZ4R8YdXJl2Vb0c8vu1infAEYGA2cYXl4kovS/dXlzUnOTvKwJPdJ1zx4VVXtvci2j03yvtbad5ZzwNbao5P8e5KT+7++XJrkL5PcJd1fdh6V5GlJnjGvzsuSfH+SP17O8ZKkqi7qP+QtNP3VyHqHJPnlJKcvsJv7J/nkyPNPJvn+qrpbv+yy1toN85bff4xtF/OXSX4vyTcWWf5v1Z3S/Laq2rBEnaN1XNRaayPLLxpZDpMkX8bPlyR5cFV9uaourao/GPlr/TTy5XFV9ZWq2lJVzxmZv1RGyBeGQr6M5Evv7/pG0nur6oEj81f688ti+TLH5xeA3YgGzsp7x7z/4T9rZNnnWmt/03+IeWuSg5Kc3v91471Jbk73YWgho6cfJ8k9F/hwceeliquqPZIcn+R3W2s3tNY+n+6vLk8dWe2q1tpfttZuaa0t1tB40rxjnz+3oLV2RGttv0Wm547s45VJ/qC1duMC+987yfUjz+ce77PAsrnl+4yx7e1U1c8n2aO19vZFxvqoJBuS3DfJVUnePfKPvYWOtXdV1Rh1wnLJl8nky78leUCSu6f7y/cJ6S6hTCacL0nOSXK/JOuTPCvJqVV1wi4cS74wLfJl/Hx5crrPBYckOT/JeVW1X79sxT6/ZMf5kvj8ArDb0cBZeU+Y9z/8140s+9LI428kSWtt/rzb/QWrPz34p5P888jsqxb4cPH1MerbP8meSS4fmXd5kgNGnl8xxn7OmXfsnxxjm1tV1eOS7NNae+siq9yYZN+R53OPb1hg2dzyub9oLbptVf1effdmfq/pPzT+WZJFv92itfZvrbWbW2vXJfn1JPdK94FpsWPd2P/Vaqk6YbnkyxiWypfW2mWttc+11v67tfapdGfp/EK/eGL50h/r4tbaVa2177TWPpzunjy7ciz5wrTIlzG11j7Uukugbmqt/UmS65L8eL94RT6/9HXsKF98fgHYDWngzIaHJbm8tbbYDXaX48tJvp3ur0ZzDk5y5cjzll3Qn8Y7/xsmbvOhI90N+jb2p/VeneQXk/xGVb2zX74l3Y315jwwyZdaa9f2y+5dVfvMW75lqW1bay9p372Z37OTHJrur1P/3tfxtiT36OvasMgQW7prwRc71mgdR/R/zZpzRHbyhowwJWsxX+ab/zs9qXwZ51g7ygj5wu5uFvNlIUt9LpjG55el6tiZOuULwCrTwJkNxyb5p0nsqD/9+Zwkf1xV+1R3n4jnp/umgYlord2/3f4bJuZ/6PiDdDc1fFA/bUryunz3WvY3JnlmVR3en5b8wiRn9fu/NMknkryoqu7UXwJ1RJJ/XGrbBXw63angc3X8Srq/ND4oyRVVdf/qvkpzj+qu739Zug+Lnxk51vOr6oCqume6m6bOHesDSb6T5Neq6o5VNfeNHu9f+lWEFbPm8qWqjqmq7+8f37df/539/ieZL6mq46rqrtU5Mt3ZfnONpA9kxxkhX9jdzVy+VNXBVfXIqtqrz4jfTnd20If63azU55cd5ovPLwC7Jw2clfeueX+xWey+Kssx//rxpLuGfP5fh5640MYLeF6Sr6e70d8H092Q8PXLrOkXFzj+3cfduHXXr189N6U7/frrrbWv9Mv/Od2lTecn+UK606RfNLKL45NsTPLVJH+a5Bfm/sI3xrajddwyr46vJPnv/vl30t0I8a1Jvpbu9dqQ5Odaa9/ud/HaJO9K9y0Sn073QfW1/b5vTvcVnU9Ld3r1L6c7Rf3mcV8nmEe+jGGpfEl3hs5FVfX1dGN/W5KXjOxiIvkysq+t6S49eGOSl7bW3tDva6mMkC+sJPkynn2SvDpdPlyZ5Ogkx/Rn2KzY55eRfS2YL/H5BWC3VK3t0tmkrLL+r8QfT3JA82YCEyRfgGmRLwCwfM7A2f3dJckLfPgBpkC+ANMiXwBgmZyBAwAAADBwzsABAAAAGDgNHGZeVW2oqlZV61a7FmC2yBdgGmQLAAvRwAEAAAAYOA0cAAAAgIHTwBmY6ry8qq6pqq9V1aeq6gH9srtV1bv6+RdU1R9V1QdHtm1V9dyq+s+quqGqXlxVP1hVH+63Oaeq9lpmPUdV1baqekFf0xer6hkjy3+2qj7e7/+KqjptZNnc6b8nVtUXqurLVfX7OzjWWVX1mqr6l77+f62qQ0aW37df9pWquqSqnjROHQsc54lV9fm513XM12GHY6mqI6vqI1V1Xf8avWr0te63fXb/3lxXVWdUVY17fJgE+SJfYBpki2wBYIW01kwDmpI8NsmFSfZLUknul+Qe/bKz++l7kxye5IokHxzZtiV5Z5J9k9w/ybeSvC/JvdN9XefFSU7s1z04yXU7mH6pX++oJLckOT3JnkmOTXJTkruOLP/hdM3AI5J8KckT+mUb+ppel+R7kjywr+l+i4z9rCQ3JPmJJHdM8n/mxpfkzv14n5FkXZIHJ/lyksOXUce6fvutSe4zctwdvQ6njDOWJA9N8oj+GBuSfCbJb8x7b96d7n09OMn2JEev9s+baW1NkS/yxWSawhTZIltMJpPJtCLTqhdgmveGJI9Ocmn/P9Q7jMzfI8m3k/zQyLw/yu0/BD1y5PmFSX5n5PnLkrximfUcleQbSdaNzLsmySMWWf8VSV7eP5774HDgyPL/SHL8ItueleTsked7J/lOkoOS/GKSf5+3/muTvGgZdfxWug+CBy60zRKvw3LH8htJ3j7vvfmxkefnzH3AMplWapIv8sVkmsYkW2SLyWQymVZmcgnVwLTW3p/kVUnOSHJNVZ1ZVfsmWZ/uLyRXjKx+xQK7+NLI428s8HzvnSjr2tbaLSPPb5rbT1U9vKrOr6rtVXV9kmcn2X/e9lcvtO0ibh1Ta+3GJF9Jcs8khyR5eH8K73VVdV2SJyf5gWXU8dtJzmitbRtr1AtbcCxVdVhVvbuqrq6qryV5yQLHX87rABMnX+QLTINskS0ArAwNnAFqrb2ytfbQdKcaH5buf97b050OfODIqgft7DGq6uCqunEH05PH3NWbk2xKclBr7S5JXpPu9OmddeuYqmrvJN+X5Kp0H47+tbW238i0d2vtOcuo42eSvLCqnjg6c4nX4ffGrPvVST6b5NDW2r5Jfm+B48Oqky+31ihfYIJky601yhYApkYDZ2Cq6mH9X2T2TPL1JN9M8t+tte8keVuS06rqe6vqvkmetrPHaa19of8Qsdj0d2Puap8kX2mtfbOqjkzySztbU+/Yqvqx/iZ6L07y0dbaFemuwT6sqp5aVXv208Oq6n7LqGNLkqOTnFFVj5+bucTr8JIx694nydeS3Ni/N89ZYn1YcfJFvsA0yBbZAsDK0MAZnn3T3Wzuq0kuT3Jtkj/vl52c7oZ+Vyd5U5K3pLsZ3Wp6bpLTq+qGJKemuz56V7w5yYvSnX780CRPSZLW2g3p/gp1fLq/al2d5KXpbhg4dh2ttU8m+bkkr6uqY3ax1lG/le6D1w3p3r+3TnDfMCnyRb7ANMgW2QLACqjW2mrXwE6qqpcm+YHW2omrXcskVNVZSba11l642rXAWidfgGmQLQCw85yBsxupqvtW1RHVOTLJM5O8fbXrAnZ/8gWYBtkCAJOzbrULYFn2SXfq8T3TfUPDy5K8c1UrAmaFfAGmQbYAwISs2iVU+++/f9uwYcOqHBsYngsvvPDLrbX1k9iXfAFGTSpfZAswymcXYFoWy5dVOwNnw4YN2bx582odHhiYqrp8UvuSL8CoSeWLbAFG+ewCTMti+eIeOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADt2QDp6peX1XXVNWnF1leVfXKqtpaVRdV1UMmXyYAAADA2jXOGThnJTl6B8uPSXJoP52U5NW7XhYAAAAAc5Zs4LTW/i3JV3awynFJ3tg6H02yX1XdY1IFAgAAAKx1k7gHzgFJrhh5vq2fdztVdVJVba6qzdu3b5/AoQE68gWYBtkCTIt8AZZrRW9i3Fo7s7W2sbW2cf369St5aGDGyRdgGmQLMC3yBViuSTRwrkxy0MjzA/t5AAAAAEzAJBo4m5I8rf82qkckub619sUJ7BcAAACAJOuWWqGq3pLkqCT7V9W2JC9KsmeStNZek+TcJMcm2ZrkpiTPmFaxAAAAAGvRkg2c1toJSyxvSf7XxCoCAAAA4DZW9CbGAAAAACyfBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwI3VwKmqo6vqkqraWlWnLLD84Ko6v6o+XlUXVdWxky8VAAAAYG1asoFTVXskOSPJMUkOT3JCVR0+b7UXJjmntfbgJMcn+atJFwoAAACwVo1zBs6RSba21i5rrd2c5Owkx81bpyXZt398lyRXTa5EAAAAgLVtnAbOAUmuGHm+rZ836rQkT6mqbUnOTfK8hXZUVSdV1eaq2rx9+/adKBdgYfIFmAbZAkyLfAGWa1I3MT4hyVmttQOTHJvkTVV1u3231s5srW1srW1cv379hA4NIF+A6ZAtwLTIF2C5xmngXJnkoJHnB/bzRj0zyTlJ0lr7SJI7Jdl/EgUCAAAArHXjNHAuSHJoVd2rqvZKd5PiTfPW+UKSn0qSqrpfugaO8wABAAAAJmDJBk5r7ZYkJyc5L8ln0n3b1JaqOr2qHt+v9oIkz6qqTyZ5S5Knt9batIoGAAAAWEvWjbNSa+3cdDcnHp136sjji5M8crKlAQAAAJBM7ibGAAAAAEyJBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwGngAAAAAAycBg4AAADAwI3VwKmqo6vqkqraWlWnLLLOk6rq4qraUlVvnmyZAAAAAGvXuqVWqKo9kpyR5KeTbEtyQVVtaq1dPLLOoUl+N8kjW2tfraq7T6tgAAAAgLVmnDNwjkyytbV2WWvt5iRnJzlu3jrPSnJGa+2rSdJau2ayZQIAAACsXeM0cA5IcsXI8239vFGHJTmsqj5UVR+tqqMX2lFVnVRVm6tq8/bt23euYoAFyBdgGmQLMC3yBViuSd3EeF2SQ5McleSEJK+rqv3mr9RaO7O1trG1tnH9+vUTOjSAfAGmQ7YA0yJfgOUap4FzZZKDRp4f2M8btS3Jptbat1trn0tyabqGDgAAAAC7aJwGzgVJDq2qe1XVXkmOT7Jp3jrvSHf2Tapq/3SXVF02uTIBAAAA1q4lGzittVuSnJzkvCSfSXJOa21LVZ1eVY/vVzsvybVVdXGS85P8dmvt2mkVDQAAALCWLPk14knSWjs3ybnz5p068rgleX4/AQAAADBBk7qJMQAAAABTooEDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADp4EDAAAAMHBjNXCq6uiquqSqtlbVKTtY74lV1apq4+RKBAAAAFjblmzgVNUeSc5IckySw5OcUFWHL7DePkl+PcnHJl0kAAAAwFo2zhk4RybZ2lq7rLV2c5Kzkxy3wHovTvLSJN+cYH0AAAAAa944DZwDklwx8nxbP+9WVfWQJAe11v5pgrUBAAAAkAncxLiq7pDkL5K8YIx1T6qqzVW1efv27bt6aIBbyRdgGmQLMC3yBViucRo4VyY5aOT5gf28OfskeUCSD1TV55M8IsmmhW5k3Fo7s7W2sbW2cf369TtfNcA88gWYBtkCTIt8AZZrnAbOBUkOrap7VdVeSY5PsmluYWvt+tba/q21Da21DUk+muTxrbXNU6kYAAAAYI1ZsoHTWrslyclJzkvymSTntNa2VNXpVfX4aRcIAAAAsNatG2el1tq5Sc6dN+/URdY9atfLAgAAAGDOLt/EGAAAAIDp0sABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGDgNHAAAAICB08ABAAAAGLixGjhVdXRVXVJVW6vqlAWWP7+qLq6qi6rqfVV1yORLBQAAAFiblmzgVNUeSc5IckySw5OcUFWHz1vt40k2ttaOSPIPSf5s0oUCAAAArFXjnIFzZJKtrbXLWms3Jzk7yXGjK7TWzm+t3dQ//WiSAydbJgAAAMDaNU4D54AkV4w839bPW8wzk7xnV4oCAAAA4LsmehPjqnpKko1J/nyR5SdV1eaq2rx9+/ZJHhpY4+QLMA2yBZgW+QIs1zgNnCuTHDTy/MB+3m1U1WOS/H6Sx7fWvrXQjlprZ7bWNrbWNq5fv35n6gVYkHwBpkG2ANMiX4DlGqeBc0GSQ6vqXlW1V5Ljk2waXaGqHpzktemaN9dMvkwAAACAtWvJBk5r7ZYkJyc5L8lnkpzTWttSVadX1eP71f48yd5J/r6qPlFVmxbZHQAAAADLtG6clVpr5yY5d968U0ceP2bCdQEAAADQm+hNjAEAAACYPA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIHTwAEAAAAYOA0cAAAAgIEbq4FTVUdX1SVVtbWqTllg+R2r6q398o9V1YaJVwoAAACwRi3ZwKmqPZKckeSYJIcnOaGqDp+32jOTfLW1dp8kL0/y0kkXCgAAALBWjXMGzpFJtrbWLmut3Zzk7CTHzVvnuCRv6B//Q5KfqqqaXJkAAAAAa9e6MdY5IMkVI8+3JXn4Yuu01m6pquuT3C3JlydR5B++a0suvuprk9gVsEIOv+e+edHj7r/aZeyQbIHdk3wBpmF3yBZgbVvRmxhX1UlVtbmqNm/fvn0lDw3MOPkCTINsAaZFvgDLNc4ZOFcmOWjk+YH9vIXW2VZV65LcJcm183fUWjszyZlJsnHjxjZukTrhwFJ2Jl9kC7AUn12AadnZfAHWrnHOwLkgyaFVda+q2ivJ8Uk2zVtnU5IT+8e/kOT9rTUhBAAAADABS56B09/T5uQk5yXZI8nrW2tbqur0JJtba5uS/HWSN1XV1iRfSdfkAQAAAGACxrmEKq21c5OcO2/eqSOPv5nkf062NAAAAACSFb6JMQAAAADLp4EDAAAAMHAaOAAAAAADp4EDAAAAMHAaOAAAAAADV6211Tlw1fYkl6/KwSdv/yRfXu0iJsyYdg+zNKZDWmvrJ7Ej+TJ4szamWRtPMntjmki+zFi2JLP3Ps/aeBJjGrpZ/eyyu79Hu3P9u3Ptye5d/9BqXzBfVq2BM0uqanNrbeNq1zFJxrR7mMUxcVuz+B7P2phmbTzJbI6J25u193nWxpMYE6tjd3+Pduf6d+fak927/t2ldpdQAQAAAAycBg4AAADAwGngTMaZq13AFBjT7mEWx8RtzeJ7PGtjmrXxJLM5Jm5v1t7nWRtPYkysjt39Pdqd69+da0927/p3i9rdAwcAAABg4JyBAwAAADBwGjgAAAAAA6eBswxVdXRVXVJVW6vqlEXWeVJVXVxVW6rqzStd43ItNaaqOriqzq+qj1fVRVV17GrUOa6qen1VXVNVn15keVXVK/vxXlRVD1npGpdrjDE9uR/Lp6rqw1X1wJWukV03a/kya9mSzF6+yJa1YdayJZm9fJm1bEnky+6gqr6vqv6lqv6z/+9dF1nvxH6d/6yqExdYvmmx93madqX+qvreqvqnqvpsn3t/ukI1L5Vdd6yqt/bLP1ZVG0aW/W4//5KqeuxK1Duvtp2qvap+uqou7H/XL6yqR6907X0dO/3a98sPrqobq+q3VqzoxbTWTGNMSfZI8l9J7p1krySfTHL4vHUOTfLxJHftn999teuewJjOTPKc/vHhST6/2nUvMaafSPKQJJ9eZPmxSd6TpJI8IsnHVrvmCYzpR0d+5o7ZHcZkut17OFP5MovZ0tc5U/kiW2Z/mrVsWcaYdqt8mbVsGXNM8mX136M/S3JK//iUJC9dYJ3vS3JZ/9+79o/vOrL8fyR582Lv81Dr//+3d3chUpVxHMe/fxK1iLSUTNtCi4TMokAqoUAxNrsohbwwKKMskCCoCLqwCysJCiyJpAgvsi40kl6EAtFKDEsyTJCK1k2j1uwFNaEkM/x38TxLp23WPfNyzpl55veBgTMzZ4bfM3P2x+6z5wU4C5gb1xkNfAzcUnDePN31APByXF4MvBGXZ8T1xwDT4vucUeJn3Uz2a4ApcXkmcLCCbaXh/JnnNwJvAo+WnX/oTXvg5Hct0O/u+939L2ADsGDIOvcDa9z9KIC7/1JyxnrlGZMD58TlccCPJearm7tvB46cZpUFwGse7ATGm9nkctI1ZqQxufsng9scsBPoKSWYtFJq/ZJct0B6/aJu6QqpdQsk2C+pdQuoXzrEAmBdXF4HLKyxzs3AFnc/Er+vLcB8ADM7G3gEWFl81Joazu/ux939I4DYI7spfhvM013ZMW0E5pmZxcc3uPsJdz8A9Mf3K0vD2d39C3cf7OAvgTPNbEwpqf/VzGePmS0EDhDyV04TOPldCPyQuT8QH8uaDkw3sx1mttPM5peWrjF5xrQCuNPMBoD3gQfLiVaYPGPuZEsJ/6WTzpJav3Rjt0Da/aJu6UypdQt0Z7+k3C2gfqnKJHc/FJd/AibVWOd0295TwCrgeGEJT6/Z/ACY2XjgVuCDAjLWlSW7jrv/DRwDJuR8bZGayZ51O7Db3U8UlHM4DeePE5WPAU+UkDOXUVUHSMwowq7IcwizuNvN7Ep3/63KUE26A3jV3VeZ2WzgdTOb6e6nqg4m/2Vmcwm/BN1QdRYpRGr9om7pEOqW5KXWLaB+6Rjql2KZ2VbgghpPLc/ecXc3M6/jfa8GLnX3h4eeK6SVisqfef9RwHrgBXff31hKycPMrgCeAXqrzlKnFcDz7v573CGncprAye8gcFHmfk98LGuAcAzvSeCAmfURfinaVU7EuuUZ01LirpLu/qmZjQUmAu2+i/Vw8oy545jZVcBawvG7h6vOI3VLrV+6sVsgwX5Rt3S81LoFurNfkusWUL+Uwd1vGu45M/vZzCa7+6F4SF6tn4+DhMndQT3ANmA2MMvMviP8PXm+mW1z9zm0UIH5B70C7HP31c2nHVGen+PBdQbi5NI44HDO1xapmeyYWQ/wNrDE3b8tPu7/NJP/OmCRmT0LjAdOmdmf7v5i4amHoUOo8tsFXGZm08xsNOHkRpuGrPMOsSTMbCJht+R2ns3NM6bvgXkAZnY5MBb4tdSUrbUJWGLB9cCxzO6XHcnMLgbeAu5y976q80hDUuuXbuwWSKxf1C1JSK1boDv7JaluAfVLm9gEDF5V6m7g3RrrbAZ6zexcC1d56gU2u/tL7j7F3acS9p7qa/XkTQ4N5wcws5WEP9IfKj4qkK+7smNaBHzo7h4fX2zhSknTCJPsn5WUG5rIHg9Re49wwukdZQUeouH87n6ju0+N2/pq4OkqJ28AXYWqnhvhKgB9hLNYL4+PPQncFpcNeA74CtgLLK46cwvGNAPYQThb9x6gt+rMI4xnPXAIOEn4r+JSYBmwLPMdrYnj3QvMqjpzC8a0Fjgav589wOdVZ9atoe85qX5JrVti5qT6Rd3SHbfUuiXnmDqqX1LrlpxjUr9U/x1NIJz3ZR+wFTgvPj4LWJtZ717CSXP7gXtqvM9UqrkKVcP5CXtgOPB1Zhu8r4TMI3XXWMKVjvoJEzSXZF67PL7uGwq+YlYrswOPA39kPuc9VHC1w2Y++8x7rKANrkJlMYyIiIiIiIiIiLQpHUIlIiIiIiIiItLmNIEjIiIiIiIiItLmNIEjIiIiIiIiItLmNIEjIiIiIiIiItLmNIEjIiIiIiIiItLmNIEjIiIiIiIiItLmNIEjIiIiIiIiItLm/gE/Q1FscHzFNQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1152x864 with 12 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_columns = 4\n",
    "plot_rows = int(e_bins_num / plot_columns) + 1\n",
    "(int(e_bins_num/plot_columns) + 1)*plot_columns\n",
    "print(f\"Total subplots: {e_bins_num} plot rows={plot_rows} columns={plot_columns}\")\n",
    "fig, axes = plt.subplots(plot_rows, plot_columns, figsize=(plot_columns*4, plot_rows*4), sharey=True)\n",
    "\n",
    "from statistics import NormalDist\n",
    "from __future__ import print_function\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "# Fit fragment is taken from here\n",
    "# https://gist.github.com/cpascual/a03d0d49ddd2c87d7e84b9f4ad2df466\n",
    "def gauss(x, H, A, x0, sigma):\n",
    "    return H + A * np.exp(-(x - x0) ** 2 / (2 * sigma ** 2))\n",
    "\n",
    "def gauss_fit(x, y):\n",
    "    mean = sum(x * y) / sum(y)\n",
    "    sigma = np.sqrt(sum(y * (x - mean) ** 2) / sum(y))\n",
    "    popt, pcov = curve_fit(gauss, x, y, p0=[min(y), max(y), mean, sigma])\n",
    "    return popt\n",
    "\n",
    "# axes is a 2d array consisting of rows and colums [[...], [...], ...]\n",
    "# flatten it as we need to feel it one by one\n",
    "axes = list(itertools.chain.from_iterable(axes))\n",
    "sigmas = []\n",
    "dee_peaks = []\n",
    "for i in range(e_bins_num):\n",
    "    hist1d = dee_hist[i]\n",
    "    ax = axes[i]\n",
    "    # (!) using ax.step() for histograms\n",
    "    # But commented ax.plot() to check ax.step() centers correctness\n",
    "    # ax.plot(err_bin_centers, hist1d, \"-o\")    \n",
    "    ax.step(err_bins[1:], hist1d)\n",
    "\n",
    "    # Fit using Gaussian function\n",
    "    H, A, x0, sigma = gauss_fit(err_bin_centers, hist1d)\n",
    "    std = np.std(hist1d)    \n",
    "\n",
    "    # Plot gauss\n",
    "    smooth_gauss_x = np.arange(err_bins[0], err_bins[-1], (err_bins[-1]-err_bins[0])/1000)\n",
    "    smooth_gauss_y = gauss(smooth_gauss_x, H, A, x0, sigma)\n",
    "    ax.plot(smooth_gauss_x, smooth_gauss_y, '--r', label='fit')\n",
    "\n",
    "    # titles\n",
    "    ax.set_title(f\"Eml/E for E={e_bins[i]}-{e_bins[i+1]}\\n sgm={sigma:.3f} peak={x0:.4f}\", wrap=True)  \n",
    "\n",
    "    # Add our error\n",
    "    sigmas.append(sigma)\n",
    "    dee_peaks.append(x0)\n",
    "    \n",
    "\n",
    "plt.tight_layout()\n",
    "plt.savefig(pics_prefix + \"ml_div_true_fit\")\n",
    "plt.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEICAYAAABS0fM3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUEklEQVR4nO3de7TdZX3n8ffHRFJaWsIlXExIkwriRJkB1xG0lBk6KoHewlBmFtTRtEuHlpFZ09GqYNcIpV2raC84TrUOrbTRYpGhSLNqXRREu6aOIicgchOJoCUxQCBcRCIY/M4f+5fO5nSHnJOzc3Z2nvdrrb3O7/c8z2/v73P2yf7s3/Pb5yRVhSSpXS8adQGSpNEyCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQaI+X5I1J/m4PqOPkJBtGXYc0bAaB9ghJfirJ/03yRJItSb6Q5NUAVXVFVZ0y6hqlvdX8URcgJfkx4G+Ac4GrgH2Ak4BnRlnX3iDJ/KraNuo6tGfzjEB7gpcBVNVfVtVzVbW1qv6uqr4KkOSXk/zD9sFJTklyT3f28OEkf5/krX1jv5Dk0iSPJ7kvyU927Q8keTjJ6r77+tkktyZ5suu/aGfFJnlPkkeSfDPJG7u2Vyd5KMm8vnFnJLltB/exIMnvJ/nH7riPJNm36zs5yYYk7+jq3ZTkV2Z47LuTPAj8WZJ9k6xJ8liSu5O8a/sSV5J3JvmrKbV9MMn/2Nn3QXsPg0B7gq8Dz3UvVqclOWBHA5McDFwNXAAcBNwD/OSUYScAX+36PwFcCbwaOBL4j8AfJdmvG/td4M3AQuBngXOTnP4CtR4GHAwsBlYDlyU5uqpuBh4F+pew3gR8bAf3cwm9ADy2q2sx8N4pj7N/1/4W4EN935fpHHsg8OPAOcCFwDLgJ4A3dN+D7f4CODXJQuidQQBnvUDd2htVlTdvI78B/wL4c2ADsA1YCxza9f0y8A/d9puBL/YdF+AB4K19Y+/t6z8GqO331bU9Chy7gzo+AFy6g76Tu9p+pK/tKuC/d9vvBq7otg8EngYOH3A/oRdAL+1rey1wf9/jbAXm9/U/DLxmmsc+C/xQX/99wMq+/bcCG/r2PwP8p27754C7Rv3z4G1ub14j0B6hqu6m9yJOkpfTe6f6AeDsKUNfQu+Ff/txNeCTPA/1bW/txk1t2697rBPovcN+Jb1rEwuA//0CpT5WVd/t2/9WVxNdzXcn+RHgPwD/p6o2DbiPRcAPA+uSbG8LMK9vzKP1/LX9p7uap3Ps5qr6Xt/+875nU7YB1tC7PvMn9M4WPj6gZu3FXBrSHqeqvkbv7OCVA7o3AUu276T3arhkwLjp+gS9s48jqmp/4CP0Xlh35IDuhX67pcC3u7o3Al8EzqC3LLSjF9RH6IXRK6pqYXfbv6r228H4mR479U8KP+97Bhwxpf9a4F8meSW9M4IrplGH9iIGgUYuycu7C6NLuv0j6J0JfGnA8E8DxyQ5vVvPfhu9NfFd9aPAlqr6XpLjgV+axjG/lWSfJCfRe+HsP4P4GPAuektS1ww6uKp+QO/d96VJDgFIsjjJyp098C4eexVwQZIDkiwGzptyn9+jd93lE8CXq+ofd1aH9i4GgfYE36F3gfemJN+lFwB3AO+YOrCqHgH+PfB+emv9K4BJdv2jpv8ZuDjJd+hdcL1qJ+MfBB6jdxZwBfBr3RnMdp+id5H2U1X19Avcz7uB9cCXkjwJ3AAcPc2aZ3rsxfSuvdzfjb2af/79WkMvvFwWalCq/I9pNL6SvIjei9wbq+pzo64HIMk3gF+tqhtGXcsgSc4Fzqqqf9PXthT4GnBYVT05suI0Ep4RaOwkWZlkYZIFwHvorekPWkaac0l+kd4a/Y2jrmW7JIcnOTHJi5IcTe9M61N9/S8C3g5caQi0yU8NaRy9lt569j7AXcDpVbV1tCVBks/TW6p6U7eWv6fYB/hfwHLgcXq/V/FhgO7C90P0Pv106ojq04i5NCRJjXNpSJIaN5ZLQwcffHAtW7Zs1GVI0lhZt27dI1W1aGr7WAbBsmXLmJycHHUZkjRWknxrULtLQ5LUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUuKEEQZJTk9yTZH2S8wf0L0jyya7/piTLpvQvTfJUkt8YRj2SpOmbdRAkmQd8CDgNWAGcnWTFlGFvAR6rqiOBS4H3Ten/Q+Azs61FkjRzwzgjOB5YX1X3VdWzwJXAqiljVgFruu2rgdclCUCS04H7gTuHUIskaYaGEQSLgQf69jd0bQPHVNU24AngoCT7Ae8GfmtnD5LknCSTSSY3b948hLIlSTD6i8UXAZdW1VM7G1hVl1XVRFVNLFq0aPdXJkmNmD+E+9gIHNG3v6RrGzRmQ5L5wP7Ao8AJwJlJ3g8sBH6Q5HtV9UdDqEuSNA3DCIKbgaOSLKf3gn8W8EtTxqwFVgNfBM4EbqyqAk7aPiDJRcBThoAkza1ZB0FVbUtyHnAdMA+4vKruTHIxMFlVa4GPAh9Psh7YQi8sJEl7gPTemI+XiYmJmpycHHUZkjRWkqyrqomp7aO+WCxJGjGDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcUMJgiSnJrknyfok5w/oX5Dkk13/TUmWde1vSLIuye3d1387jHokSdM36yBIMg/4EHAasAI4O8mKKcPeAjxWVUcClwLv69ofAX6+qo4BVgMfn209kqSZGcYZwfHA+qq6r6qeBa4EVk0ZswpY021fDbwuSarq1qr6dtd+J7BvkgVDqEmSNE3DCILFwAN9+xu6toFjqmob8ARw0JQxvwjcUlXPDKEmSdI0zR91AQBJXkFvueiUFxhzDnAOwNKlS+eoMkna+w3jjGAjcETf/pKubeCYJPOB/YFHu/0lwKeAN1fVN3b0IFV1WVVNVNXEokWLhlC2JAmGEwQ3A0clWZ5kH+AsYO2UMWvpXQwGOBO4saoqyULg08D5VfWFIdQiSZqhWQdBt+Z/HnAdcDdwVVXdmeTiJL/QDfsocFCS9cDbge0fMT0POBJ4b5KvdLdDZluTJGn6UlWjrmHGJiYmanJyctRlSNJYSbKuqiamtvubxZLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNW4oQZDk1CT3JFmf5PwB/QuSfLLrvynJsr6+C7r2e5KsHEY90ly79taNnHjJjSw//9OceMmNXHvrxlGXJE3brIMgyTzgQ8BpwArg7CQrpgx7C/BYVR0JXAq8rzt2BXAW8ArgVODD3f1JY+PaWzdywTW3s/HxrRSw8fGtXHDN7YaBxsYwzgiOB9ZX1X1V9SxwJbBqyphVwJpu+2rgdUnStV9ZVc9U1f3A+u7+pLHxe9fdw9bvP/e8tq3ff47fu+6eEVUkzcwwgmAx8EDf/oaubeCYqtoGPAEcNM1jAUhyTpLJJJObN28eQtnScHz78a0zapf2NGNzsbiqLquqiaqaWLRo0ajLkf7JSxbuO6N2aU8zjCDYCBzRt7+kaxs4Jsl8YH/g0WkeK+3R3rnyaPZ98fMvbe374nm8c+XRI6pImplhBMHNwFFJlifZh97F37VTxqwFVnfbZwI3VlV17Wd1nypaDhwFfHkINUlz5vTjFvO7ZxzD4oX7EmDxwn353TOO4fTjBq5ySnuc+bO9g6raluQ84DpgHnB5Vd2Z5GJgsqrWAh8FPp5kPbCFXljQjbsKuAvYBrytqp4b+EDSHuz04xb7wq+xld4b8/EyMTFRk5OToy5DksZKknVVNTG1fWwuFkuSdg+DQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcbMKgiQHJrk+yb3d1wN2MG51N+beJKu7th9O8ukkX0tyZ5JLZlOLJGnXzPaM4Hzgs1V1FPDZbv95khwIXAicABwPXNgXGL9fVS8HjgNOTHLaLOuRJM3QbINgFbCm214DnD5gzErg+qraUlWPAdcDp1bV01X1OYCqeha4BVgyy3okSTM02yA4tKo2ddsPAocOGLMYeKBvf0PX9k+SLAR+nt5ZhSRpDs3f2YAkNwCHDej6zf6dqqokNdMCkswH/hL4YFXd9wLjzgHOAVi6dOlMH0aStAM7DYKqev2O+pI8lOTwqtqU5HDg4QHDNgIn9+0vAT7ft38ZcG9VfWAndVzWjWViYmLGgSNJGmy2S0NrgdXd9mrgrweMuQ44JckB3UXiU7o2kvwOsD/w67OsQ5K0i2YbBJcAb0hyL/D6bp8kE0n+FKCqtgC/Ddzc3S6uqi1JltBbXloB3JLkK0neOst6JEkzlKrxW2WZmJioycnJUZchSWMlybqqmpja7m8WS1LjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUuFkFQZIDk1yf5N7u6wE7GLe6G3NvktUD+tcmuWM2tUiSds1szwjOBz5bVUcBn+32nyfJgcCFwAnA8cCF/YGR5AzgqVnWIUnaRbMNglXAmm57DXD6gDErgeuraktVPQZcD5wKkGQ/4O3A78yyDknSLpptEBxaVZu67QeBQweMWQw80Le/oWsD+G3gD4Cnd/ZASc5JMplkcvPmzbMoWZLUb/7OBiS5AThsQNdv9u9UVSWp6T5wkmOBl1bVf0uybGfjq+oy4DKAiYmJaT+OJOmF7TQIqur1O+pL8lCSw6tqU5LDgYcHDNsInNy3vwT4PPBaYCLJN7s6Dkny+ao6GUnSnJnt0tBaYPungFYDfz1gzHXAKUkO6C4SnwJcV1V/XFUvqaplwE8BXzcEJGnuzTYILgHekORe4PXdPkkmkvwpQFVtoXct4ObudnHXJknaA6Rq/JbbJyYmanJyctRlSNJYSbKuqiamtvubxZLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMalqkZdw4wl2Qx8a9R1zNDBwCOjLmKOOec2OOfx8eNVtWhq41gGwThKMllVE6OuYy455zY45/Hn0pAkNc4gkKTGGQRz57JRFzACzrkNznnMeY1AkhrnGYEkNc4gkKTGGQRDkuSIJJ9LcleSO5P81yn970hSSQ7u9pPkg0nWJ/lqkleNpvJdswvzPTnJE0m+0t3eO5rKd92O5pzkoiQb++b2M33HXNA9x/ckWTm66nfNTOecZFmSrX3tHxntDGbuhX62k/yXJF/r2t/f1z7WzzNV5W0IN+Bw4FXd9o8CXwdWdPtHANfR+yW4g7u2nwE+AwR4DXDTqOewm+d7MvA3o657d8wZuAj4jQHjVwC3AQuA5cA3gHmjnsdunvMy4I5R172b5vzTwA3Agq7vkL3lefaMYEiqalNV3dJtfwe4G1jcdV8KvAvovzK/CvhY9XwJWJjk8LmseTZ2Yb5jbydzHmQVcGVVPVNV9wPrgeN3f6XDswtzHnsvMOdzgUuq6pmu7+HukLF/ng2C3SDJMuA44KYkq4CNVXXblGGLgQf69jcwpv/ApjlfgNcmuS3JZ5K8Yk6LHLL+OXdN53VLfJcnOaBr22ueY5j2nAGWJ7k1yd8nOWnOCx2iKXN+GXBSkpu6ub26Gzb2z7NBMGRJ9gP+Cvh1YBvwHmDs1sOnawbzvYXe3zn5V8D/BK6doxKHrn/OVfUk8MfAS4FjgU3AH4yuut1jBnPeBCytquOAtwOfSPJjc1/x7A2Y83zgQHpLue8ErkqSEZY4NAbBECV5Mb0fnCuq6hp6/1CWA7cl+SawBLglyWHARnpr6dst6drGxkzmW1VPVtVTAFX1t8CLt19IHicD5kxVPVRVz1XVD4A/4f8vC4z9cwwzm3O3PPJot72O3nr5y0ZT+a4bNGd67/Sv6ZZzvwz8gN4fnxv759kgGJLuncFHgbur6g8Bqur2qjqkqpZV1TJ6P0ivqqoHgbXAm7tPD70GeKKqNo2q/pma6XyTHLb93VOS4+n97D06ovJ3yaA5d+3913b+HXBHt70WOCvJgiTLgaOAL89VvcMw0zknWZRkXrf9E/TmfN/cVTx7O5ozvbPYn+7GvAzYh95fIB3753n+qAvYi5wIvAm4PclXurb3dO9+B/lbep8cWg88DfzKbq9wuGY63zOBc5NsA7YCZ1X3kYsxMnDOwNlJjqV3cfybwK8CVNWdSa4C7qK3bPa2qnpujmuerRnNGfjXwMVJvk/vHfOvVdWWuSx4CHY058uBy5PcATwLrO5+hsf+efZPTEhS41wakqTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcf8PqlukNBXOYOsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(e_bin_centers, sigmas, \"--o\")\n",
    "plt.title(\"Sigma by energy\")\n",
    "plt.savefig(pics_prefix + \"errors\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAASPElEQVR4nO3df7TkdV3H8edLFmnNEpBFcZftboLWmifAadXjjzCVH5UtGdXaqba06BeezDRRT/5AzwGzwiy1SOmQp0SOIe7JbEPQfp0C7vJDWBXZkIIVdQXENglcfPfHfDeHe2bZOzuzd5j9PB/n3HO/38/3MzPvz87d+5rv5zvzuakqJEntesS0C5AkTZdBIEmNMwgkqXEGgSQ1ziCQpMYtm3YB++KII46oubm5aZchSTNly5YtX6mqFQvbZzII5ubmmJ+fn3YZkjRTkvznsHanhiSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcRMJgiSnJLkpybYkZw05fkiSD3bHr0wyt+D46iQ7k7xqEvVIkhZv7CBIchDwLuBUYC3wkiRrF3R7GXB3VR0DnAe8bcHxPwQ+Nm4tkqTRTeKMYB2wrapuqar7gYuA9Qv6rAcu7LY/BDw/SQCSnAZ8Htg6gVokSSOaRBCsBG4b2L+9axvap6p2AfcAj03yaOA1wJv39iBJzkgyn2R+x44dEyhbkgTTv1j8JuC8qtq5t45VdX5V9aqqt2LFiv1fmSQ1YtkE7mM7cPTA/qqubVif25MsAx4D3Ak8HTg9ye8BhwLfTPK/VfUnE6hLkrQIkwiCq4Fjk6yh/wt/A/AzC/psAjYC/wacDlxRVQU8Z3eHJG8CdhoCkrS0xg6CqtqV5ExgM3AQcEFVbU1yNjBfVZuA9wHvT7INuIt+WEiSHgbSf2E+W3q9Xs3Pz0+7DEmaKUm2VFVvYfu0LxZLkqbMIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJatxEgiDJKUluSrItyVlDjh+S5IPd8SuTzHXtL0yyJckN3fcfmkQ9kqTFGzsIkhwEvAs4FVgLvCTJ2gXdXgbcXVXHAOcBb+vavwK8qKqeCmwE3j9uPZKk0UzijGAdsK2qbqmq+4GLgPUL+qwHLuy2PwQ8P0mq6tqq+kLXvhVYnuSQCdQkSVqkSQTBSuC2gf3bu7ahfapqF3AP8NgFfX4CuKaq7ptATZKkRVo27QIAkjyF/nTRSQ/R5wzgDIDVq1cvUWWSdOCbxBnBduDogf1VXdvQPkmWAY8B7uz2VwEfBn6+qv5jTw9SVedXVa+qeitWrJhA2ZIkmEwQXA0cm2RNkkcCG4BNC/pson8xGOB04IqqqiSHAh8Fzqqqf51ALZKkEY0dBN2c/5nAZuAzwMVVtTXJ2Ul+rOv2PuCxSbYBrwR2v8X0TOAY4A1Jruu+jhy3JknS4qWqpl3DyHq9Xs3Pz0+7DEmaKUm2VFVvYbufLJakxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXETCYIkpyS5Kcm2JGcNOX5Ikg92x69MMjdw7LVd+01JTp5EPdJSu/Ta7Tzr3CtYc9ZHeda5V3DptdunXZK0aGMHQZKDgHcBpwJrgZckWbug28uAu6vqGOA84G3dbdcCG4CnAKcA7+7uT5oZl167nddecgPbv3ovBWz/6r289pIbDAPNjEmcEawDtlXVLVV1P3ARsH5Bn/XAhd32h4DnJ0nXflFV3VdVnwe2dfcnzYy3b76Je7/xwIPa7v3GA7x9801TqkgazSSCYCVw28D+7V3b0D5VtQu4B3jsIm8LQJIzkswnmd+xY8cEypYm4wtfvXekdunhZmYuFlfV+VXVq6reihUrpl2O9P+ecOjykdqlh5tJBMF24OiB/VVd29A+SZYBjwHuXORtpYe1V5/8ZJYf/OBLW8sPPohXn/zkKVUkjWYSQXA1cGySNUkeSf/i76YFfTYBG7vt04Erqqq69g3du4rWAMcCV02gJmnJnHb8Ss558VNZeehyAqw8dDnnvPipnHb80FlO6WFn2bh3UFW7kpwJbAYOAi6oqq1Jzgbmq2oT8D7g/Um2AXfRDwu6fhcDnwZ2Ab9RVQ8MfSDpYey041f6i18zK/0X5rOl1+vV/Pz8tMuQpJmSZEtV9Ra2z8zFYknS/mEQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1bqwgSHJ4ksuS3Nx9P2wP/TZ2fW5OsrFre1SSjyb5bJKtSc4dpxZJ0r4Z94zgLODyqjoWuLzbf5AkhwNvBJ4OrAPeOBAYv19V3wMcDzwryalj1iNJGtG4QbAeuLDbvhA4bUifk4HLququqrobuAw4paq+XlWfAKiq+4FrgFVj1iNJGtG4QfC4qrqj2/4i8LghfVYCtw3s3961/b8khwIvon9WIUlaQsv21iHJx4HHDzn0+sGdqqokNWoBSZYBHwDeWVW3PES/M4AzAFavXj3qw0iS9mCvQVBVL9jTsSRfSnJUVd2R5Cjgy0O6bQdOHNhfBXxyYP984Oaqesde6ji/60uv1xs5cCRJw407NbQJ2NhtbwQ+MqTPZuCkJId1F4lP6tpI8lbgMcArxqxDkrSPxg2Cc4EXJrkZeEG3T5JekvcCVNVdwFuAq7uvs6vqriSr6E8vrQWuSXJdkl8asx5J0ohSNXuzLL1er+bn56ddhiTNlCRbqqq3sN1PFktS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1LixgiDJ4UkuS3Jz9/2wPfTb2PW5OcnGIcc3JblxnFokSftm3DOCs4DLq+pY4PJu/0GSHA68EXg6sA5442BgJHkxsHPMOiRJ+2jcIFgPXNhtXwicNqTPycBlVXVXVd0NXAacApDk0cArgbeOWYckaR+NGwSPq6o7uu0vAo8b0mclcNvA/u1dG8BbgD8Avr63B0pyRpL5JPM7duwYo2RJ0qBle+uQ5OPA44ccev3gTlVVklrsAyc5DnhiVf1Wkrm99a+q84HzAXq93qIfR5L00PYaBFX1gj0dS/KlJEdV1R1JjgK+PKTbduDEgf1VwCeBZwK9JLd2dRyZ5JNVdSKSpCUz7tTQJmD3u4A2Ah8Z0mczcFKSw7qLxCcBm6vqPVX1hKqaA54NfM4QkKSlN24QnAu8MMnNwAu6fZL0krwXoKruon8t4Oru6+yuTZL0MJCq2Ztu7/V6NT8/P+0yJGmmJNlSVb2F7X6yWJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1LhU1bRrGFmSHcB/TruOER0BfGXaRSwxx9wGxzw7vquqVixsnMkgmEVJ5quqN+06lpJjboNjnn1ODUlS4wwCSWqcQbB0zp92AVPgmNvgmGec1wgkqXGeEUhS4wwCSWqcQTAhSY5O8okkn06yNclvLjj+20kqyRHdfpK8M8m2JJ9KcsJ0Kt83+zDeE5Pck+S67usN06l83+1pzEnelGT7wNh+eOA2r+2e45uSnDy96vfNqGNOMpfk3oH2P53uCEb3UD/bSV6e5LNd++8NtM/080xV+TWBL+Ao4IRu+zuAzwFru/2jgc30PwR3RNf2w8DHgADPAK6c9hj283hPBP522nXvjzEDbwJeNaT/WuB64BBgDfAfwEHTHsd+HvMccOO0695PY34e8HHgkO7YkQfK8+wZwYRU1R1VdU23/d/AZ4CV3eHzgN8BBq/Mrwf+svr+HTg0yVFLWfM49mG8M28vYx5mPXBRVd1XVZ8HtgHr9n+lk7MPY555DzHmXwPOrar7umNf7m4y88+zQbAfJJkDjgeuTLIe2F5V1y/othK4bWD/dmb0P9gixwvwzCTXJ/lYkqcsaZETNjjmrunMborvgiSHdW0HzHMMix4zwJok1yb5xyTPWfJCJ2jBmJ8EPCfJld3YfqDrNvPPs0EwYUkeDfwN8ApgF/A6YObmwxdrhPFeQ3+dk+8H/hi4dIlKnLjBMVfV14D3AE8EjgPuAP5getXtHyOM+Q5gdVUdD7wS+Osk37n0FY9vyJiXAYfTn8p9NXBxkkyxxIkxCCYoycH0f3D+qqouof8fZQ1wfZJbgVXANUkeD2ynP5e+26qubWaMMt6q+lpV7QSoqr8DDt59IXmWDBkzVfWlqnqgqr4J/DnfmhaY+ecYRhtzNz1yZ7e9hf58+ZOmU/m+GzZm+q/0L+mmc68Cvkl/8bmZf54NggnpXhm8D/hMVf0hQFXdUFVHVtVcVc3R/0E6oaq+CGwCfr5799AzgHuq6o5p1T+qUceb5PG7Xz0lWUf/Z+/OKZW/T4aNuWsfvLbz48CN3fYmYEOSQ5KsAY4Frlqqeidh1DEnWZHkoG77u+mP+Zalq3h8exoz/bPY53V9ngQ8kv4KpDP/PC+bdgEHkGcBPwfckOS6ru113avfYf6O/juHtgFfB35xv1c4WaOO93Tg15LsAu4FNlT3losZMnTMwEuSHEf/4vitwK8AVNXWJBcDn6Y/bfYbVfXAEtc8rpHGDDwXODvJN+i/Yv7VqrprKQuegD2N+QLggiQ3AvcDG7uf4Zl/nl1iQpIa59SQJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAIdkJI8MLAC5nVJzpp2TQsl+WSS/xr8dGqSS5Ps3MvtPrFwhcskr0jynv1Vqw5sfo5AB6p7q+q4Sd5hkmVVtWuS9wl8lf771v8lyaH0V77cmw8AG+iv8LrbBvoL/Ukj84xATUlya5I3J7kmyQ1Jvqdr//Zu8bSrugXT1nftv5BkU5IrgMuTPCrJxemvVf/hbgGyXpKXJnnHwOP8cpLzFlHSRfR/iQO8GLhk8GCSVye5ulvc7c1d84eAH0nyyK7PHPAE4J/3+R9GTTMIdKBavmBq6KcHjn2lqk6gv3Daq7q21wNXVNU6+ssIvD3Jt3fHTgBOr6ofBH4duLuq1gK/Czyt63Mx8KJujRrof1L8gkXUeTnw3G5Zhg3AB3cfSHIS/eUK1tFf3O1pSZ7bfVL3KuDUrusG4OIZ/KS2HiacGtKB6qGmhna/6t5C/1U4wEnAjyXZHQzfBqzuti8bWCbh2cAfAVTVjUk+1W3v7M4afjTJZ4CDq+qGRdT5APAv9H+ZL6+qWwcuGZzUfV3b7T+afjD8E9+aHvpI9/1li3gsaSiDQC26r/v+AN/6PxDgJ6rqpsGOSZ4O/M8i7/e99Nek+SzwFyPUcxHwYfp/9etBDw+cU1V/NuQ2HwHOS/9PnD6qW+lT2idODUl9m4GXD6yQevwe+v0r8FNdn7XAU3cfqKor6S9H/DP0X7HT9bs8yUP9oZJ/Bs4ZvM1ATS/t1sUnycokR3aPtRP4BP3pp4W3k0biGYEOVMsHVo4E+Puqeqi3kL4FeAfwqSSPAD4P/OiQfu8GLkzyafqv/LcC9wwcvxg4rqruBuju6xhgjytwdnP7vz+k/R+SfC/wb10+7QR+Ftj9JxI/QP9MYsPC20qjcPVRaQTdRd2Dq+p/kzyR/h8zf3JV3d8d/1vgvKq6vNv/PuClVfXKqRUt7YVBII0gyXfQn5I5mP4c/muq6mPdZwCuAq6vqp+cYonSyAwCSWqcF4slqXEGgSQ1ziCQpMYZBJLUOINAkhr3f5NGQynmuprCAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sigma_div_e = sigmas/e_bin_centers\n",
    "\n",
    "plt.plot(e_bin_centers, sigma_div_e, \"--o\")\n",
    "plt.xlabel(\"Energy, MeV\")\n",
    "plt.savefig(pics_prefix + \"errors_div_energy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "array must not contain infs or NaNs",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-37-077e30400007>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mx\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mpopt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpcov\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcurve_fit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcalo_error\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0me_bin_centers\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msigma_div_e\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mp0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      5\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpopt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"a={a} b={b} c={c}\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/scipy/optimize/minpack.py\u001b[0m in \u001b[0;36mcurve_fit\u001b[0;34m(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)\u001b[0m\n\u001b[1;32m    711\u001b[0m     \u001b[0;31m# NaNs can not be handled\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    712\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mcheck_finite\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 713\u001b[0;31m         \u001b[0mydata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray_chkfinite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mydata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfloat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    714\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    715\u001b[0m         \u001b[0mydata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mydata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfloat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.8/site-packages/numpy/lib/function_base.py\u001b[0m in \u001b[0;36masarray_chkfinite\u001b[0;34m(a, dtype, order)\u001b[0m\n\u001b[1;32m    483\u001b[0m     \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0morder\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0morder\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    484\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mchar\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mtypecodes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'AllFloat'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misfinite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 485\u001b[0;31m         raise ValueError(\n\u001b[0m\u001b[1;32m    486\u001b[0m             \"array must not contain infs or NaNs\")\n\u001b[1;32m    487\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: array must not contain infs or NaNs"
     ]
    }
   ],
   "source": [
    "def calo_error(x, a, b, c):\n",
    "    return a/np.sqrt(x) + b/x + c\n",
    "\n",
    "popt, pcov = curve_fit(calo_error, e_bin_centers, sigma_div_e, p0=[1, 1, 0])\n",
    "a, b, c = popt\n",
    "print(f\"a={a} b={b} c={c}\")\n",
    "\n",
    "# Plot gauss\n",
    "smooth_cal_x = np.arange(e_bin_centers[0], e_bin_centers[-1], (e_bin_centers[-1]-e_bin_centers[0])/1000)\n",
    "smooth_cal_y = calo_error(smooth_cal_x, a, b, c)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(e_bin_centers, sigma_div_e, \"o\")\n",
    "ax.plot(smooth_cal_x, smooth_cal_y, '--r', label='fit')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAbf0lEQVR4nO3de5QU5b3u8e/DcBcEgYHAgIJh1OANcLxtNRKN1ySinmRvycXLzjqsZWJiTnZMdLlWTMzZuZnsZJsYjeYYReMtXpAYdvAeohF1CIoQgg7eYEBAARFBB4bf+aNrTNsyMw3TM9VT/XzW6kXV+1ZV/962fbrmrZoeRQRmZpZdPdIuwMzMOpeD3sws4xz0ZmYZ56A3M8s4B72ZWcY56M3MMs5Bb1ZCkr4j6ea06zDL56C3Nkl6WVKTpGEF7QskhaSxXVzPlOR57yloPzhpfzSvLSSN38ExzpXULGlTwWNU3jajJK1Ill+WtKVg21924jDNSspBb8V4CZjWsiLpQKB/euWwFjhS0tC8tnOA53fiGE9ExICCx8q8/lOBP+Wtf6pg2ws6UH9ZkdQz7RqscznorRg3AWfnrZ8DzMjfQFIfST+R9Kqk1ZKukdQv6dtD0n2S1kpanyyPztv3UUnfk/S4pLck3V/4E0SBJmAmcFayfxXwb8DvSjLanFOB2bu4b19Jtydj+ZukgwEkXSTprvwNJV0p6b93dJDkp4q7ktftJUlfzev7jqQ7JM1InmexpLqd2PdOSTdL2gicK2mcpLnJsR6UdFXLFJSkP0r6SkFtCyWdsYuvj3UxB70VYx6wu6SPJKF6FlA4D/1DYB9gIjAeqAG+nfT1AH4L7AXsCWwBCqc+PgucBwwHegPfaKemGfzzw+ckYBGwsvXNiyepF/BR4IFdPMRU4PfAEOAWYGZyzJuBkyUNTp6nJ7nXckbhAST1AP4APEvutTwe+Jqkk/I2Ow24DRgMzCJ5TYvcdypwZ7Lv75I6nwKGAt8BvpC37Y3A5/NqOzg57h+Lf0ksTQ56K1bLWf0JwBKgsaVDkoDpwP+JiHUR8RbwfZIz7oh4IyLuiojNSd9/AscWHP+3EfF8RGwB7iD3gdGqiPgrMETSvkldHwjLdhwhaUPeY1le30eBZ5NaW8ws2P5/t3Hs+RFxZ0RsBf4L6AscERGrgLnAZ5LtTgZej4j5OzjGoUB1RFweEU0R8SJwHclrmngsImZHRDO5/z4H78S+T0TEzIjYDlQn+3w72f4xch8cLWYB+0iqTda/ANweEU1tvAZWRjw3Z8W6iVxIjeODoVpNbs5+fi7zARBQBSCpP/AzcsG2R9I/UFJVElIAr+UdbzMwoMiaLgA+Bvw7uZ8KijUvIo5upW9H0zanR8SDRR57ectCRGxPLuq2XOi9ETifXPB+ntwYdmQvYJSkDXltVcBf8tYLX7O+yU8Jxey7PG95FLAuIjYX9I9JxvCOpNuBz0v6LrnrNZ9upW4rQw56K0pEvCLpJXIh+MWC7tfJTcfsHxGNH9gZ/gPYFzg8Il6TNBFYQO7DoCNuAhqAGRGxOe9DpqNOBc7swP5jWhaSaZTR/HNaaSZwtaQDgE8C32zlGMuBlyKitpX+thSzb/7X1q4i99NR/7ywH1Ow/Y3kXu/HgM0R8cQu1GUp8dSN7YwvAsdFxNv5jcmP/9cBP5M0HEBSTd6c8EByHwQbJA0BLitFMRHxErkpoEvb2Ky3pL55j6q2jilpHNAnIpZ0oLRDJJ2ZnF1/DXiX3HUOIuIdcnPjtwBPRcSrrRzjKeAtSd+S1E9SlaQDJB1axPPv1L4R8QpQD3xHUm9JRwKfKtjmCWA78FNa/ynEypSD3ooWEcsior6V7m+RO7uel9zJ8SC5s3iAnwP9yJ35z+P9ty12tKbHCm6LLLSY3IdMy+O8pP1IffA++kOBT7Dju23+ULDtPTvYpsW95O4CWk9uPvvMZL6+xY3AgbQRmMmU1ifJXat4idxr9xtgUBvP25F9PwccCbwB/F/gdnIfUPlmJHX7F8K6GfkPj5j9k6TZwC8jYldvrSzmOfYE/gF8KCI2dtbzdEQyJ/+PiLgsr+1sYHob1zasTPmM3uz9HgUe6ayDJ3P2XwduK6eQl3SopA9L6iHpZHK3X87M6+8PfAm4NqUSrQN8MdYsT0T8uLOOLWk3YDXwCrk7kMrJh4C7yd1HvwI4PyIWACTXWu4mNx13S2oV2i7z1I2ZWcZ56sbMLOPKbupm2LBhMXbs2LTLMDPrVubPn/96RFTvqK/sgn7s2LHU17d2B5+Zme2IpFda6/PUjZlZxjnozcwyzkFvZpZxDnozs4xz0JuZZZyD3sws4xz0ZmYZ56A3M8s4B72ZWcY56M3MMs5Bb2aWcQ56M7OMc9CbmWWcg97MLOMc9GZmGeegNzPLOAe9mVnGtRv0kq6XtEbSolb6JelKSQ2SFkqaXNC/u6QVkn5ZqqLNzKx4xZzR3wCc3Eb/KUBt8pgOXF3Q/z1g7q4UZ2ZmHddu0EfEXGBdG5tMBWZEzjxgsKSRAJIOAUYA95eiWDMz23mlmKOvAZbnra8AaiT1AH4KfKO9A0iaLqleUv3atWtLUJKZmbXozIuxXwJmR8SK9jaMiGsjoi4i6qqrqzuxJDOzytOzBMdoBMbkrY9O2o4EjpH0JWAA0FvSpoi4uATPaWZmRSpF0M8CLpB0G3A48GZErAI+17KBpHOBOoe8mVnXazfoJd0KTAGGSVoBXAb0AoiIa4DZwKlAA7AZOK+zijUzs53XbtBHxLR2+gP4cjvb3EDuNk0zM+ti/s1YM7OMc9CbmWWcg97MLOMc9GZmGeegNzPLOAe9mVnGOejNzDLOQW9mlnEOejOzjHPQm5llnIPezCzjHPRmZhnnoDczyzgHvZlZxjnozcwyzkFvZpZxDnozs4xz0JuZZZyD3sws4xz0ZmYZ56A3M8s4B72ZWcY56M3MMs5Bb2aWcQ56M7OMc9CbmWWcg97MLOMc9GZmGeegNzPLOAe9mVnGOejNzDLOQW9mlnEOejOzjHPQm5llXLtBL+l6SWskLWqlX5KulNQgaaGkyUn7RElPSFqctP9bqYs3M7P2FXNGfwNwchv9pwC1yWM6cHXSvhk4OyL2T/b/uaTBu1ypmZntkp7tbRARcyWNbWOTqcCMiAhgnqTBkkZGxPN5x1gpaQ1QDWzoYM1mZrYTSjFHXwMsz1tfkbS9R9JhQG9gWQmez8zMdkKnX4yVNBK4CTgvIra3ss10SfWS6teuXdvZJZmZVZRSBH0jMCZvfXTShqTdgT8Cl0bEvNYOEBHXRkRdRNRVV1eXoCQzM2tRiqCfBZyd3H1zBPBmRKyS1Bu4h9z8/Z0leB4zM9sF7V6MlXQrMAUYJmkFcBnQCyAirgFmA6cCDeTutDkv2fVfgY8CQyWdm7SdGxHPlK58MzNrTzF33Uxrpz+AL++g/Wbg5l0vzczMSsG/GWtmlnEOejOzjHPQm5llnIPezCzjHPRmZhnnoDczyzgHvZlZxjnozcwyzkFvZpZxDnozs4xz0JuZZZyD3sws4xz0ZmYZ56A3M8s4B72ZWcY56M3MMs5Bb2aWcQ56M7OMc9CbmWWcg97MLOMc9GZmGeegNzPLOAe9mVnGOejNzDLOQW9mlnEOejOzjOuZdgGlMnNBI1fMWcrKDVsYNbgfF520L6dPqkm7LDOz1GUi6GcuaOSSu59jy9ZmABo3bOGSu58DcNibWcXLxNTNFXOWvhfyLbZsbeaKOUtTqsjMrHxkIuhXbtiyU+1mZpUkE0E/anC/nWo3M6skmQj6i07al369qt7X1q9XFRedtG9KFZmZlY9MXIxtueB6xZylNCbTNV8/YR9fiDUzIyNBD7mwP31SDWveeodjfvQIL76+Ke2SzMzKQmaCvsXwgX255vOHMGnPwWmXYmZWFtqdo5d0vaQ1kha10i9JV0pqkLRQ0uS8vnMkvZA8zill4W352H7DGdy/d1c9nZlZWSvmYuwNwMlt9J8C1CaP6cDVAJKGAJcBhwOHAZdJ2qMjxe6MxSvfZOpVj7N83eauekozs7LUbtBHxFxgXRubTAVmRM48YLCkkcBJwAMRsS4i1gMP0PYHRkkN3a0PS1Zu5FePNnTVU5qZlaVS3F5ZAyzPW1+RtLXW/gGSpkuql1S/du3aEpQEHxrUl7MOG8Pv61ewYr3P6s2scpXFffQRcW1E1EVEXXV1dcmOe/6UD9ND4lePLivZMc3MuptSBH0jMCZvfXTS1lp7lxk5qB//euhofl+//L37683MKk0pgn4WcHZy980RwJsRsQqYA5woaY/kIuyJSVuXOn/KeL572gFUD+jT1U9tZlYW2r2PXtKtwBRgmKQV5O6k6QUQEdcAs4FTgQZgM3Be0rdO0veAp5NDXR4RbV3U7RQ1g/vx2cP37OqnNTMrG+0GfURMa6c/gC+30nc9cP2ulVZav3vyFdZtauIrx9emXYqZWZcqi4uxXWHh8jf5xSMNrN74TtqlmJl1qYoJ+i9/bDzN24Nr/uw7cMysslRM0O85tD9nTqrhlidfZY3P6s2sglRM0ANccNx4tm0Pfj33xbRLMTPrMpn79sq27DV0N/7jxH3Yf9SgtEsxM+syFRX0AF+aMj7tEszMulRFTd20WP92Ez/+0z94fdO7aZdiZtbpKjLo33i7iav/vIzr/uK5ejPLvooM+vHDB/Cpg0Zx0xOvsO7tprTLMTPrVBUZ9ABfPX48W7Y28xuf1ZtZxlVs0I8fPpBPHDiSG//6Mut9Vm9mGVZxd93k++rxtWxrDjZvbabL/sahmVkXq+ig32fEQK75wiFpl2Fm1qkqduomX8OaTfxp0aq0yzAz6xQVfUbf4qf3L+WxF17nyL2HMah/r7TLMTMrKZ/Rk/sOnLfe3cZv//pS2qWYmZWcgx7Yf9QgTpgwgusfe4mN72xNuxwzs5Jy0CcuPL6Wje9s44bHX067FDOzknLQJw6oGcQnDhrJtubtaZdiZlZSvhib55fTJiEp7TLMzErKZ/R5JBER/HXZ62x6d1va5ZiZlYSDvsALazbx2eueZMYTL6ddiplZSTjoC+wzYiDH7lPNdXNf5G2f1ZtZBjjod+DCj9eyfvNWbp73StqlmJl1mIN+BybvuQfH1A7j2rkvsrnJZ/Vm1r056FvxtY/XAvDC6k0pV2Jm1jG+vbIVh+w1hMcvPo6+varSLsXMrEN8Rt+Gvr2q2Na8nYY1Pqs3s+7LQd+Ob931HNOum8c7W5vTLsXMbJc46NvxmbrRrH3rXW596tW0SzEz2yUO+nYcsfdQDh83hGv+vMxn9WbWLTnoi3Dh8bWs3vgud9QvT7sUM7Od5qAvwpEfHsqhY/fggb+vTrsUM7Od5tsriyCJX33uEIbs1jvtUszMdlpRZ/SSTpa0VFKDpIt30L+XpIckLZT0qKTReX0/lrRY0hJJV6qbfg9w9cA+VPUQm97dRtM2f2e9mXUf7Qa9pCrgKuAUYAIwTdKEgs1+AsyIiIOAy4EfJPv+C3AUcBBwAHAocGzJqu9iK9Zv5ugfPcyd81ekXYqZWdGKOaM/DGiIiBcjogm4DZhasM0E4OFk+ZG8/gD6Ar2BPkAvoNtOdNcM7sfYobtx1SMNPqs3s26jmKCvAfJvN1mRtOV7FjgzWT4DGChpaEQ8QS74VyWPORGxpGMlp0cSF368lsYNW7j7bz6rN7PuoVR33XwDOFbSAnJTM41As6TxwEeA0eQ+HI6TdEzhzpKmS6qXVL927doSldQ5puxTzUGjB3HVow1s9d+XNbNuoJigbwTG5K2PTtreExErI+LMiJgEXJq0bSB3dj8vIjZFxCbgf4AjC58gIq6NiLqIqKuurt61kXQRSVx4fC3L123hLy+U94eSmRkUF/RPA7WSxknqDZwFzMrfQNIwSS3HugS4Pll+ldyZfk9Jvcid7XfbqZsWx+03nFkXHMVx+41IuxQzs3a1G/QRsQ24AJhDLqTviIjFki6XdFqy2RRgqaTngRHAfybtdwLLgOfIzeM/GxF/KO0Qup4kDho9GMDTN2ZW9hQRadfwPnV1dVFfX592GUW5du4y7py/gtlfPYaeVf4lYzNLj6T5EVG3oz6nUwfsOWQ3nl+9iVnPrky7FDOzVjnoO+DECSMYuXsfvnnnQsZd/EeO+uHDzFzQ2P6OZmZdyN910wGznl3JG283sW17bvqrccMWLrn7OQBOn1T4qwZmZunwGX0HXDFnKU3N77/GsWVrM1fMWZpSRWZmH+Sg74CVG7bsVLuZWRoc9B0wanC/nWo3M0uDg74DLjppX/r1qnpfW79ePbjopH1TqsjM7IN8MbYDWi64XjFnKY3JdM0Xjx7nC7FmVlYc9B10+qQaTp9Uw+ambRzyvQdZt3lr2iWZmb2Pp25KpH/vnpy4/whmP7fK31VvZmXFQV9CUyeOYsPmrf5WSzMrKw76EjqmtpqDRw9iy9bmtEsxM3uP5+hLqFdVD+694Oi0yzAzex+f0XeCrc3bWbPxnbTLMDMDHPSd4oxfPc637lqYdhlmZoCDvlMcPb6auS+8zhub3k27FDMzB31nmDpxFM3bg9mLXku7FDMzB31n2O9DA9lnxABmPePvpjez9DnoO4Ekpk6s4emX17Ni/ea0yzGzCufbKzvJZw4ZzaFjhzBqkL/J0szS5aDvJMN378vw3fumXYaZmaduOtOaje/wnVmLeX71W2mXYmYVzEHfiap6iJvmvcI9/oPhZpYiB30nGjqgD8fUDmPWMyvZvj3a38HMrBM46DvZ6RNraNywhb+9uj7tUsysQjnoO9kJE0bQt1cPZvqeejNLiYO+k+3WpydnTh5N/96+wcnM0uH06QLfP+PAtEswswrmM/ouEhH+LVkzS4WDvov89P7nOfFnc9nctC3tUsyswjjou8jRtcPY3NTMA39fnXYpZlZhHPRd5LCxQxg5qC+znlmZdilmVmEc9F2kRw9x2sGj+PPza1n/dlPa5ZhZBXHQd6HTJo5i2/Zg9qJVaZdiZhXEt1d2oQkjd+c3Z9dx1PhhaZdiZhWkqDN6SSdLWiqpQdLFO+jfS9JDkhZKelTS6Ly+PSXdL2mJpL9LGlvC+rsVSXx8wgj69a5KuxQzqyDtBr2kKuAq4BRgAjBN0oSCzX4CzIiIg4DLgR/k9c0AroiIjwCHAWtKUXh3tX17cNUjDdw1f0XapZhZhSjmjP4woCEiXoyIJuA2YGrBNhOAh5PlR1r6kw+EnhHxAEBEbIqIiv6toR49xINLVvObx15KuxQzqxDFBH0NsDxvfUXSlu9Z4Mxk+QxgoKShwD7ABkl3S1og6YrkJ4T3kTRdUr2k+rVr1+78KLqZqQePYsmqjf6DJGbWJUp11803gGMlLQCOBRqBZnIXe49J+g8F9gbOLdw5Iq6NiLqIqKuuri5RSeXrEweNoqqHuNffaGlmXaCYoG8ExuStj07a3hMRKyPizIiYBFyatG0gd/b/TDLtsw2YCUwuQd3dWvXAPhw1fhj3PrOSCP9BEjPrXMUE/dNAraRxknoDZwGz8jeQNExSy7EuAa7P23ewpJbT9OOAv3e87O7vf02uoXb4AN7csjXtUsws49q9jz4itkm6AJgDVAHXR8RiSZcD9RExC5gC/EBSAHOBLyf7Nkv6BvCQJAHzges6Zyjdy9SJNUydWHipw8ys9FRuUwd1dXVRX1+fdhldZtWbWxg2oA+9qvxLyma26yTNj4i6HfU5XVL05Itv8C8/fJjHG15PuxQzyzAHfYom7jmYgX16+hstzaxTOehT1KdnFaceOJI5i19jS1Nz2uWYWUY56FN22sRRvN3UzEP/8B8kMbPO4aBP2eHjhjJi9z7c6+kbM+sk/prilFX1EL+YNpmxw/qnXYqZZZSDvgwcNm5I2iWYWYZ56qZM3L/4Nb4/e0naZZhZBjnoy8SSVW9x3V9eZOWGLWmXYmYZ46AvE1MnjiIC7lvoi7JmVloO+jIxdthuHDxmsO++MbOSc9CXkakHj2Lxyo00rPEfJDGz0nHQl5FPHjySA2sGsX6zv7rYzErHt1eWkeED+/KHrxyddhlmljE+oy9Db7+7jTc2vZt2GWaWEQ76MtO0bTtH/ehhfvFwQ9qlmFlGlN0fHpH0FrA07Tq62DCg0r6U3mOuDB5z19krIqp31FGOc/RLW/srKVklqd5jzj6PuTKU45g9dWNmlnEOejOzjCvHoL827QJS4DFXBo+5MpTdmMvuYqyZmZVWOZ7Rm5lZCTnozcwyrqyCXtLJkpZKapB0cdr1dISk6yWtkbQor22IpAckvZD8u0fSLklXJuNeKGly3j7nJNu/IOmcNMZSDEljJD0i6e+SFku6MGnP8pj7SnpK0rPJmL+btI+T9GQyttsl9U7a+yTrDUn/2LxjXZK0L5V0UkpDKpqkKkkLJN2XrGd6zJJelvScpGck1Sdt3ee9HRFl8QCqgGXA3kBv4FlgQtp1dWA8HwUmA4vy2n4MXJwsXwz8KFk+FfgfQMARwJNJ+xDgxeTfPZLlPdIeWyvjHQlMTpYHAs8DEzI+ZgEDkuVewJPJWO4AzkrarwHOT5a/BFyTLJ8F3J4sT0je732Accn/B1Vpj6+dsX8duAW4L1nP9JiBl4FhBW3d5r2d+guY96IdCczJW78EuCTtujo4prEFQb8UGJksjyT3y2EAvwamFW4HTAN+ndf+vu3K+QHcC5xQKWMG+gN/Aw4n91uRPZP2997XwBzgyGS5Z7KdCt/r+duV4wMYDTwEHAfcl4wh62PeUdB3m/d2OU3d1ADL89ZXJG1ZMiIiViXLrwEjkuXWxt4tX5Pkx/NJ5M5wMz3mZArjGWAN8AC5M9MNEbEt2SS//vfGlvS/CQylm40Z+DnwTWB7sj6U7I85gPslzZc0PWnrNu/tcvwKhIoQESEpc/e2ShoA3AV8LSI2SnqvL4tjjohmYKKkwcA9wH7pVtS5JH0SWBMR8yVNSbmcrnR0RDRKGg48IOkf+Z3l/t4upzP6RmBM3vropC1LVksaCZD8uyZpb23s3eo1kdSLXMj/LiLuTpozPeYWEbEBeITctMVgSS0nUfn1vze2pH8Q8Abda8xHAadJehm4jdz0zX+T7TETEY3Jv2vIfaAfRjd6b5dT0D8N1CZX73uTu3AzK+WaSm0W0HKl/Rxy89gt7WcnV+uPAN5MfiScA5woaY/kiv6JSVvZUe7U/f8BSyLiv/K6sjzm6uRMHkn9yF2TWEIu8D+dbFY45pbX4tPAw5GbrJ0FnJXcoTIOqAWe6pJB7KSIuCQiRkfEWHL/jz4cEZ8jw2OWtJukgS3L5N6Ti+hO7+20L3IUXNw4ldzdGsuAS9Oup4NjuRVYBWwlNxf3RXJzkw8BLwAPAkOSbQVclYz7OaAu7zj/DjQkj/PSHlcb4z2a3DzmQuCZ5HFqxsd8ELAgGfMi4NtJ+97kQqsB+D3QJ2nvm6w3JP175x3r0uS1WAqckvbYihz/FP55101mx5yM7dnksbglm7rTe9tfgWBmlnHlNHVjZmadwEFvZpZxDnozs4xz0JuZZZyD3sws4xz0ZmYZ56A3M8u4/w9qxHAEsAZwfAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(e_bin_centers, dee_peaks , \"--o\")\n",
    "plt.title(\"Mean MLE/E by energy\")\n",
    "plt.ylim(0.4,1.4)\n",
    "plt.savefig(pics_prefix + \"ml_div_true_peaks\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
  },
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "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.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}