{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d563033c",
   "metadata": {},
   "source": [
    "# Import tools/models be used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "fc9e6e26",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from math import pi\n",
    "from sklearn.gaussian_process import GaussianProcessRegressor   # import model to be used\n",
    "from sklearn.gaussian_process.kernels import RBF                # import kernel to be used\n",
    "import joblib                                                   # to store model\n",
    "import os\n",
    "import errno\n",
    "from sklearn.model_selection import train_test_split  # this is very important"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "313e157e",
   "metadata": {},
   "source": [
    "# Read data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "bdc08ed0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>$M_{inv}$</th>\n",
       "      <th>$rapidity$</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>121.079631</td>\n",
       "      <td>-0.304619</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>196.219429</td>\n",
       "      <td>-0.401002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>478.422633</td>\n",
       "      <td>-1.010334</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>99.018601</td>\n",
       "      <td>0.097534</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>423.950224</td>\n",
       "      <td>2.223391</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9996</th>\n",
       "      <td>128.145991</td>\n",
       "      <td>0.246759</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9997</th>\n",
       "      <td>576.782956</td>\n",
       "      <td>-0.554320</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9998</th>\n",
       "      <td>35.437849</td>\n",
       "      <td>0.628855</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9999</th>\n",
       "      <td>189.229229</td>\n",
       "      <td>0.462469</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10000</th>\n",
       "      <td>443.127754</td>\n",
       "      <td>-1.858528</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>10001 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        $M_{inv}$  $rapidity$\n",
       "0      121.079631   -0.304619\n",
       "1      196.219429   -0.401002\n",
       "2      478.422633   -1.010334\n",
       "3       99.018601    0.097534\n",
       "4      423.950224    2.223391\n",
       "...           ...         ...\n",
       "9996   128.145991    0.246759\n",
       "9997   576.782956   -0.554320\n",
       "9998    35.437849    0.628855\n",
       "9999   189.229229    0.462469\n",
       "10000  443.127754   -1.858528\n",
       "\n",
       "[10001 rows x 2 columns]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "directory = '../data/'\n",
    "filein    = 'checkLO.txt'\n",
    "file      = directory+filein\n",
    "data = np.genfromtxt(file)\n",
    "\n",
    "m        = data[:,0:1]\n",
    "rap      = data[:,1:2]\n",
    "x2       = data[:,2:3]\n",
    "\n",
    "m = np.concatenate((m),axis=0) \n",
    "rap = np.concatenate((rap),axis=0) \n",
    "x2 = np.concatenate((x2),axis=0) \n",
    "\n",
    "Data_ave = np.stack((m,rap),axis=-1)\n",
    "df = pd.DataFrame(Data_ave)\n",
    "df.columns = ['$M_{inv}$','$rapidity$']\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0248b96b",
   "metadata": {},
   "source": [
    "# Split data into training and testing sub-sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "de4cf784",
   "metadata": {},
   "outputs": [],
   "source": [
    "# for reproducibility\n",
    "#X_train, X_test, y_train, y_test = train_test_split(Data_ave, x2, test_size=0.4, train_size=0.6 , random_state=0)\n",
    "X_train, X_test, y_train, y_test = train_test_split(Data_ave, x2, test_size=0.4, train_size=0.6)\n",
    "\n",
    "\n",
    "# rename data for convenience\n",
    "x2_train   = y_train\n",
    "m_train    = X_train[:,0:1]\n",
    "rap_train  = X_train[:,1:2]\n",
    "\n",
    "x2_test   = y_test\n",
    "m_test    = X_test[:,0:1]\n",
    "rap_test  = X_test[:,1:2]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2624e86c",
   "metadata": {},
   "source": [
    "# Now let's do some training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "b8ef444a",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''Training and prediction'''\n",
    "#-----------------------------------------------------------------------------#\n",
    "# Input:                                                                      #\n",
    "# X_train  := training data set                                               #\n",
    "# y_train  := targets for training                                            #\n",
    "# X_trest  := testing data set                                                #\n",
    "# y_test   := targets for testing                                             #\n",
    "# filename := name for training storage                                       #\n",
    "# model    := name of the mode we want to use                                 #\n",
    "# ----------------------------------------------------------------------------#\n",
    "# Output:                                                                     #  \n",
    "# loaded_model := trained model                                               #\n",
    "# R2 := quality of the training on the testing set                            #\n",
    "# y_fit := predicted targets                                                  #\n",
    "# y_est := estimated targets                                                  #\n",
    "# ----------------------------------------------------------------------------#\n",
    "# Output labelling :                                                          #\n",
    "# Training(X,y,filename)[0] -> R2                                             #\n",
    "# Training(X,y,filename)[1] -> y_new                                          #\n",
    "# ----------------------------------------------------------------------------#\n",
    "# Note: not all models have the same characteristics\n",
    "\n",
    "def training (X_train, y_train, X_test, y_test,filename,model):\n",
    "\n",
    "    # train \n",
    "    model.fit(X_train, y_train)\n",
    "\n",
    "    # save trained model\n",
    "    joblib.dump(model, filename)\n",
    "\n",
    "    # load model \n",
    "    loaded_model = joblib.load(filename)\n",
    "\n",
    "    # R^2\n",
    "    R2   = loaded_model.score(X_test, y_test)\n",
    "    \n",
    "    # predicted targets \n",
    "    y_fit = loaded_model.predict(X_test)\n",
    "    # estimated targets\n",
    "    y_est = loaded_model.predict(X_train)\n",
    "\n",
    "    \n",
    "    return R2, y_fit,y_est, model.get_params(deep=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "id": "06b6bed9",
   "metadata": {},
   "outputs": [],
   "source": [
    "x2_train,x2_test;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "ca19e5ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "m_train,m_test;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "be16d388",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/piazurita/opt/anaconda3/lib/python3.9/site-packages/sklearn/gaussian_process/_gpr.py:506: ConvergenceWarning: lbfgs failed to converge (status=2):\n",
      "ABNORMAL_TERMINATION_IN_LNSRCH.\n",
      "\n",
      "Increase the number of iterations (max_iter) or scale the data as shown in:\n",
      "    https://scikit-learn.org/stable/modules/preprocessing.html\n",
      "  _check_optimize_result(\"lbfgs\", opt_res)\n"
     ]
    }
   ],
   "source": [
    "# name of the model for storage\n",
    "filename = 'storage/x2_LO_Gaussian_RBF_model.sav'\n",
    "\n",
    "# train \n",
    "Training_machine_learning = training(X_train,y_train,X_test,y_test,filename,GaussianProcessRegressor(kernel=RBF(1.0),random_state=0))\n",
    "\n",
    "# predict\n",
    "R2    = Training_machine_learning[0] \n",
    "y_fit = Training_machine_learning[1]\n",
    "y_est = Training_machine_learning[2]\n",
    "coefs = Training_machine_learning[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "3b225d13",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_est,y_train;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "1cf2c418",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_fit,y_test;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "49492be3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999992731542482"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93efe92e",
   "metadata": {},
   "outputs": [],
   "source": [
    "coefs;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56253913",
   "metadata": {},
   "source": [
    "# Let's plot and see what we have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "c97b86ca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAELCAYAAAAybErdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABJKElEQVR4nO2deXyU1bn4v887E8gCIQkEQYlCq5VSAojgcl2LLS54cUGr9mrFttZeauX216tiKwUUKxaV4tYiWqi0VgVcUrGKmhbqdhtAxBWLkLBLIAlLEiAz8/z+mJkwSWYms887k/P9fAIz75z3vM9558x53nPOs4iqYjAYDAZDIFa6BTAYDAaD/TDKwWAwGAwdMMrBYDAYDB0wysFgMBgMHTDKwWAwGAwdcKZbgETRp08fHThwYLrFMGQpq1ev3q2qpem4tunbhmQSqm9njXIYOHAgq1atSrcYhixFRGrSdW3Ttw3JJFTfNstKBoPBYOiAUQ4Gg8Fg6IBRDgaDwWDoQNbsOQSjpaWFrVu3cvDgwXSLkjXk5uYyYMAAcnJy0i2KwWBIIlmtHLZu3UrPnj0ZOHAgIpJucTIeVWXPnj1s3bqVQYMGpVscg8GQRLJaORw8eNAohgQiIvTu3Zva2tp0i5IQ3n1oIqP3vIQDD24sqnpfwum3LEy3WBGT6fIb7E3W7zkYxZBYsuV+vvvQRE7b8wJO8SACTvFw2p4XePehiekWLSIyXX6D/cl65WAwBGP0npdor+dEvMczgUyX32B/jHJIIdOnT+f+++8P+fmLL77IJ598kkKJspeqinnUTx+ATuuFTutF/fRjqKqY1/q5A0/Q80IdtxuZLr/B/hjlYCOMckgMVRXzGL76FxSzn1d65HN+2dGcM7CYX+z6LVMXTALAHaLrhzpuJ959aCLLCvIZO+Bohg0sY+yAo1lWkA9khvyGzCCrN6Sj5cX3tzH7tfVsb2jm6KI8bj3/RC496Zi46rznnnt46qmnKCsro7S0lJNPPpn58+fz+OOPc/jwYY4//ngWLVrE2rVrqaioYMWKFcycOZOlS5dSWVnZoVx+fn6CWpudzKhcxMpdD1E7qD+Fbg/7HRYe3/rLzhwnf/WsxFm5iLG9L+G0PS+0WZpRxbupmybZI6X24HKm9u3d2q4dOU5+UdobVeide77t5TdkBuYxw8eL72/jjuc/ZFtDMwpsa2jmjuc/5MX3t8Vc5+rVq3nmmWd4//33ef7556mqqgLg8ssvp6qqig8++ICvf/3rPPnkk/zHf/wH48ePZ/bs2axdu5avfvWrQcsZQjOjchGLa+awK8dCRdjrdLQOoH7clrCk+hFOv2Uh7/W+DJdaqIJLLd7rfVlGWPvMLC3u0C6PCHeVlmSE/IbMwCgHH7NfW09zi7vNseYWN7NfWx9znf/85z+57LLLyM/Pp7CwkPHjxwPw0UcfcdZZZ1FeXs6f//xnPv7446DnR1rO4GXppvmI1dJpObWamFG5iJsKNjJi0ACGHzuce855KGMG1uYQFmPNIsyoXJRiaQzZSlqUg4hcICLrRWSDiEwJ8nkvEfmriHwgIh+LyA3Jlml7Q3NUxyMlmOnnxIkTeeSRR/jwww+ZNm1aSA/uSMsZvHgc9RGXXVwzB3XWIwLqrGdxzZy4B9bU9esQ5sQiCWmHwQBpUA4i4gAeBS4EhgDXiMiQdsV+AnyiqsOBc4EHRKRbMuU6uigvquORcPbZZ/PCCy/Q3NzM/v37+etf/wrA/v376d+/Py0tLfz5z39uLd+zZ0/279/f+j5UOUNwLHdxhCWlwwxDrBaWbpof87VT2a81lHIg/nYYDH7SMXM4BdigqhtV9TDwDHBJuzIK9BTvY3cPoA5wJVOoW88/kbwcR5tjeTkObj3/xJjrHDlyJFdddRUjRoxgwoQJnHXWWQDcfffdnHrqqXz7299m8ODBreWvvvpqZs+ezUknncQXX3wRspwhOBMG3Yh62sZ8Um1bRj0OvN2rI9HMPIKQsn49sNt5HdoVSJztMBgAEA3Xy5JxQZErgAtU9Ye+99cBp6rqzQFlegIVwGCgJ3CVqi4LUtePgB8BHHvssSfX1LTNWfHpp5/y9a9/PWLZkmGtlI1Ee19TycVP/4zqw28iKH4lYAEeBHEXM2HQjSzdNB91dhxA81ry+NcP/xW0XhFZraqjQl03kf3aVzZs3z5twZU0ymdBV5jCtcNgaE+ovp2OmUOwOXF7DXU+sBY4GhgBPCIihR1OUn1cVUep6qjS0vgzOF560jG8PWUMm2aN4+0pY4xiyDBmVC6i+tAKRNTby0RABI8IuerhEmsoFx1o4ra6neR62jmLqZIjB1r9IGIgYf3aK07ovj2jchGNrA9xReXCpt1tHP4MhlhIh5/DVqAs4P0AYHu7MjcAs9Q7rdkgIpvwPm2ZxyFDKzMqF7F003w8jnosdzEqhxBHcGulg5bFq54VVNStxFPaDQFyPB5afAoEEfY5HfzN5wcxbcx10YqTsn69pPoRcISY8YvwSs88Tlj/AKO5KZpqDYY2pGPmUAWcICKDfJtxV+OdageyGTgPQESOAk4ENqZUSoOtmbpgEn+tua+NxZFaTWHPOeibRSCCihxRDAEcsiTWDd2U9OsZlYugs3ZaFk8VO8KWMRg6I+UzB1V1icjNwGuAA/iDqn4sIj/2ff574G5goYh8iHfyfLuq7k61rAZ7UlUxj/fcf+dQTtvu22nA2GCR6oIQy4Zuqvr1so0PQwR5lnY6jXIwxEdawmeo6ivAK+2O/T7g9XZgbKrlMmQGZWtm8+XA3KCfqUagJDohcpPY9tdOfr8+6GwipJ9DALkuE2bFEB/GQ9qQcfTVWvq53EE/E08+4ioOZa3akXbWeuqxmDDoxjglTB6h2t0Ghb75ZyRfGENWY5RDBvGPf/yDiy++GICKigpmzZoVsmxDQwOPPfZY6/vt27dzxRVXJF3GVLBLSplc39DB4qi7R7li4M2s+8FK8lyROi+2fwq3dzKj79W7O1patUeg+tAK4yltiAujHGyA2x3B02A7xo8fz5QpHSI0tNJeORx99NEsWbIkJvnsxpaRtzLmgIvpu+vo3+JCVOnX4uJCzmbamOuoqpjHhU27g3jAtXurdNAFYrlt7WF84ok/547affRvcXkbEMJPyXhKG+LFKIdA1j0Hc4bC9CLv/+uei7vK6upqBg8ezPXXX8+wYcO44ooraGpqYuDAgdx1112ceeaZLF68mOXLl3P66aczcuRIrrzySg4cOADAq6++yuDBgznzzDN5/vnnW+tduHAhN9/s9a/68ssvueyyyxg+fDjDhw/nnXfeYcqUKXzxxReMGDGCW2+9lerqaoYOHQp4c2vfcMMNlJeXc9JJJ/H3v/+9tc7LL7+cCy64gBNOOIHbbrst7vYnkhmVixj25NncUPcI55QNop6evLplB8urD/Lrvv/D3Td4lWHZmtm8m989yAY0qIp3X8IVel/Bzh7Go8ffRNngO/levS9zQ5gNFju3w2B/TD4HP+ueg7/eAi2+QHt7t3jfAwz7TlxVr1+/nieffJIzzjiD73//+61P9Lm5ubz11lvs3r2byy+/nDfeeIOCggLuu+8+HnzwQW677TZuvPFGKisrOf7447nqqquC1n/LLbdwzjnn8MILL+B2uzlw4ACzZs3io48+Yu3atYBXSfl59NFHAfjwww/57LPPGDt2LJ9//jkAa9eu5f3336d79+6ceOKJ/PSnP6WsrKz9JVOOPxy3OFsQoDmnmVml3XnmmBt4+btz6BdQtq/WstMZSmblo4kfAnDKk6fQ7OwYWDHWDelU8UqPfBb3yUMk/PKS3dthsDdm5uDnzbuOKAY/Lc3e43FSVlbGGWd4NwivvfZa3nrrLYDWwf69997jk08+4YwzzmDEiBH88Y9/pKamhs8++4xBgwZxwgknICJce+21QeuvrKzkv//7vwFwOBz06tUrrDxvvfUW113ndfIaPHgwxx13XKtyOO+88+jVqxe5ubkMGTKE9mEb0kWwcNwiUH34jQ5r67ukNOTG7dEuFzunH88XC27il/UdPaUtj8PWG9LgdYITK7xiyPV4+GFhyGgfBkOnmJmDn71bozseBe3DdvvfFxQUAKCqfPvb3+Yvf/lLm3Jr164NGvI7XsLF0+revXvra4fDgcuV1HiHEeNx1AfdKhbxDpbTOOLR/OhxY2lgRUe7VlXObmqiHw14qp/hqwJObWFucRE7nQ76udz8cF8L37khau/olKJWU/htc1Uu2X+AH+56BXgwRVIZsg0zc/DTa0B0x6Ng8+bNvPvuuwD85S9/4cwzz2zz+Wmnncbbb7/Nhg0bAGhqauLzzz9n8ODBbNq0iS+++KL13GCcd955/O53vwO8m9v79u3rEP47kLPPPrs1BPjnn3/O5s2bOfHE2KPPpoJwSyT+5D3gXX56gfdotqygTm8rfWlWLd9H4xqbWL51O+uqt7B863au2LczKfKnFF87c5uzoC2GtGGUg5/zfgU57cwfc/K8x+Pk61//On/84x8ZNmwYdXV1rUtAfkpLS1m4cCHXXHMNw4YN47TTTuOzzz4jNzeXxx9/nHHjxnHmmWdy3HHHBa1/7ty5/P3vf6e8vJyTTz6Zjz/+mN69e3PGGWcwdOhQbr311jblJ02ahNvtpry8nKuuuoqFCxe2mTHYkQmDbgwZplqEVsuczrLBdeY5vN3TO2YZU4V4Ondw2+l00KAFKZDGkK2kPGR3shg1apSuWrWqzbGoQ0uve867x7B3q3fGcN6v4t6Mrq6u5uKLL+ajjz6Kqx47ka6Q3Rc//TNqDr8R1BVBFT6a+CFDF5aH9ZDu3+Ji+dbtree0W3XiWTmfq6d3tFLrLGR3Mmnft2dULmLx5t902s6/bvmSdSffy+jxJgCfITR2CtltX4Z9B372EUxv8P4fp2IwJJZpPQZTGMYnZNiTZ4d9qs71eJhc3wDAXSVFjBhURvnAMoYPLGNmSREicDarEy12wpk25jpy3GECLKmy0+lgdu+elK2ZnTrBDFmFUQ5JZuDAgVk1a0gnZWtm84u6+qAewq2RWeWgL9vbEVTB6cph8u4WLjzQzB0l/Vlc2LM1QqtHhGcLezKzpIh+uidVzYmLn9cdQEJ5Svuizj5b2JPHS+xhUGDIPIxyMNgWv9Pb0IXlDHvybKryGxnX2MT03XX0crmDegeL5UG0O+Iq9n7szkc8+bQ4WvhNST/uPudh/lrYLehm9eLCnuySPqlpXJz00v2dRxgUYWlhj9QIZMg6jHIw2JKLn/4Zizf/pk2+hhmlvVlWkM+4xibyw4RfVauJdT9YyZXH3oZKCziaWutYXDPHl0K0Ix68oTnszozKRfyqtDcagZmzx+axogz2xSgHg+2YUbmI6sNvdBj7D1nCb4uLgPBWR36z16COc1YLoUO2Cq/0sH+o6yXVj+CyzKBvSC5GORhsx9JN80OumOx0OtlJaA9oVVo9nMPGFgpmpSewuGaO7aOZdpbxLhBRj+3bY7AnRjnYhIULF7J9+5GUwz/84Q/55JNP4q63urqap59+OurzJk6cmLYorhpmUO/vcrGp5ExOdXyz48a0KgMajmV48RggtONcf5ebq/btD7FnkV3RTNWysqo9htRhlINNaK8cnnjiCYYMGRJ3vbEqh3Qwo3IRpzxxCiGXfVSZXN/A6D0v8ZOa5Vyy/wCWP2y1KvmqbC+qYebqy7j+ud95Hec8bbu402fOemddQ8jVeLtHMw1rxhoEu7fHYE+Mcghg2cZljF0ylmF/HMbYJWNZtnFZ3HX+6U9/4pRTTmHEiBHcdNNNuN1uJk6cyNChQykvL2fOnDksWbKEVatW8V//9V+MGDGC5uZmzj33XPyOTz169OD222/n5JNP5lvf+hb/+te/OPfcc/nKV75CRYU3h311dTVnnXUWI0eOZOTIkbzzzjsATJkyhX/+85+MGDGCOXPm4Ha7ufXWWxk9ejTDhg1j3rx5gDfe0s0338yQIUMYN24cu3btirvt0eCPutqc0xx8o1mV05qbGdfYhAMPVfmNvNSzR6s5KiI0WRYqwqGcZtY2zqNq5xrae8wpVuuEITeEA6jdo5le0NgUMo9DMOzeHoM9MYH3fCzbuIzp70znoPsgADsadzD9nekAjPvKuJjq/PTTT3n22Wd5++23ycnJYdKkScycOZNt27a1+j40NDRQVFTEI488wv3338+oUR2dcBsbGzn33HO57777uOyyy7jzzjt5/fXX+eSTT7j++usZP348ffv25fXXXyc3N5d///vfXHPNNaxatYpZs2Zx//338/LLLwPw+OOP06tXL6qqqjh06BBnnHEGY8eO5f3332f9+vV8+OGHfPnllwwZMoTvf//7MbU7FpZumo84Q4e9QISanG4AuLH4bUkJB63QzzYey0314TcRq+0g6rbgF337cEeI2Yl6hCtsHpV1db5EnChbFQqso5IskSEbMcrBx9w1c1sVg5+D7oPMXTM3ZuXw5ptvsnr1akaPHg1Ac3MzF1xwARs3buSnP/0p48aNY+zYzvPNd+vWjQsuuACA8vJyunfvTk5ODuXl5a15GlpaWrj55ptZu3YtDoejNQR3e5YvX866deta9xP27t3Lv//9b1auXMk111yDw+Hg6KOPZsyYMTG1ORZmVC4KGXU1kJ1OB6pQ1fsSdjkj8WQOoQAEQqYDFa8Hsp3pLD5UICJwgPVJlMaQrRjl4GNnY/AIlqGOR4Kqcv3113Pvvfe2OX7PPffw2muv8eijj/Lcc8/xhz/8IWw9OTk5raG7LctqDZJnWVZrSO05c+Zw1FFH8cEHH+DxeMjNzQ0p08MPP8z555/f5vgrr7ySlPDgndGaxCeCBc5+LjcfdhvBl2fdg2f11Vg5DWHLx9Yc+8cac3q60eIIM8tqRyi/DoMhHGbPwUe/gn5RHY+E8847jyVLlrSu39fV1VFTU4PH42HChAncfffdrFmzBiBsiO1I2Lt3L/3798eyLBYtWtSal7p9veeffz6/+93vaGnxDi6ff/45jY2NnH322TzzzDO43W527NjRmjo02XQWRdWPPy5S+eG1fOvFk/hV3dagYTTaE31cSfv7D7QQ3YY0KO8+NDEZohiyGDNz8DF55OQ2ew4AuY5cJo+cHHOdQ4YMYebMmYwdOxaPx0NOTg4PPvggl112GR7fwOafVUycOJEf//jH5OXlteZ+iIZJkyYxYcIEFi9ezDe/+c3WRELDhg3D6XQyfPhwJk6cyOTJk6murmbkyJGoKqWlpbz44otcdtllVFZWUl5ezte+9jXOOeecmNsdDZ0uJ6nS3+Vmcn0D4xqbQGBFgYMni3M5KIKlSquKCDJVENQ7GYhkGqHKwG7fir4RKabTZD9BGL3npaTIYsheTMjuAJZtXMbcNXPZ2biTfgX9mDxycsz7DdlMIkJ2z6hcxNJN873KIcxIJ6qsq97S+n5ZQT7T+7TbjPb34ViXxVSxgCv27WfqLcHTotopZPewJ89GnVGYp6qybtMWZMbeJEhnyHRC9W0zcwhg3FfGGWWQAlr3GZwtnT4Bt/eEnltc1NFKKc69kv4uN8u3bselmbHKemzeyUHDi4RC8Fp4mR+7IRoy49dgyCpC7TO0n8R29yi31DW0ORaNpU5Q2l3Ev5fht4LKBDY3r45KHyrwSbdhSZPHkJ1kvXLIlmUzu5CI+xnKY1eAQrfHu8yjyiGBh0qKWFaQ3zqmh4qp1E7IkDvRRR4P/VtciCr9W1xMq63j/AMHea/3ZZx+y8LYGpRiovV47u9y07dlW5KkMWQrWT3TzM3NZc+ePfTu3TstZprZhqqyZ8+ekGaykWK5i4OvmQs0CeQALb7va0eOkymlvbm3pJjb9jTw07q9/Kq0JHxU0lDftSoNlsU+C76zbz+/2LMXa0YDAKfH1aLUEvL+BUOVAwJ9tTa5QhmyjqxWDgMGDGDr1q3U1pofRqLIzc1lwIABcdUxYdCNLK65H7E6mqK6gnk9i7DX6eDu0mKO6j6WlkOVCJ2bsQarB7x5G54t7Ekjedwb/gxbMmHQjZ3mkG5FhP0OBzN7F/OrpEtmyCay2lrJYF+G/mEkEoUjlx9VQSQxfVYU1k38MLKyNrJWAhi6sDyqfQdLlQ8mmnS1ho4YayVD2qmqmEfZmtn01VpkUBmxOZwl7mFGM9hzONo7F8M8y9DFyfoNaYM9qKqYx9DVd9KPWsImMVMlJ8xsNpFbRxbAuucSV2GKqKqYR0xKMgPbakgfRjkYUkLZmtnkyeHW90UhQl/kq3J37R6K3O5YYl9EjipX7ttP098ybyW+bM3skPcvHJnYVkP6MMrBkDRmVC5i2JNnU76wnO8NyGFZwZH8zFP21HeYIeSo8svaei480MzizYe5VM8OrR8CkvxEjK+8pcpV+/ZzZ10DeU078Ezrxc7px/ueyO1PX63l/AON0eV0AHKbdyRPKEPWYZSDISn4vaDVWQ/iNUmd3qekVUGMa2zi7to9bXwOZuzaQ//B07BmNNBv+gbG9x7OUa7gT8j9XW76u9wRrzN19yizavfwYfUWPqjewp0+5zoRsAT6UcvQ1XdmhIJYUlDCSz17RLXG5gFQMqJ9BntglIMhoVRVzGPn9ON5e+OvO3hBH7Qs7i05kpVsXGMTy7duZ131FpZv3c64xiZGj7+ptZ6hq+/k/9XXdYi+2t3jTRcasbe0CnfU7uOiA03MLCli+MAyygP+xg44mmUF+eTJYcrWzI7vBqSAx0vywiY6CsVN/frw9dVTkyCRIRsxysGQMAI3nUMN3HsdFssK8llWkM/YAUczLGBw3iWlreX8exTjGpuYvruudYZxVIubqbUNQGQWO+IRZu6qY0LTXu7pXcSzhT3bpBZFpNXRrnxgGd8bkMOMykUJuBvJoapiHl86Q/xswy0zifBeXh7XHlOUFLkM2UdalIOIXCAi60Vkg4hMCVHmXBFZKyIfi8iKVMtoiJ7ATedCd4gNUxFm9S5mep8SduQ4Ud/gPL1PCQuKv8aL72/jjFmVrR69ywrymVtcxE6ng34uN/9TV49T3EzvU+Id5Duh0OPikqb9LCvI59nCnqGXYgIUxeKaOTEpiFT067I1s0MrRZFOFcQX3brZWvkZ7EPK/RxExAE8Cnwb2ApUiUiFqn4SUKYIeAy4QFU3i0jfVMtpiJ6+Wtv6OB9u3G6wrA4FDloWfy/YxI9eOIWzxVtN+/DcO3KczCgtIVc14mWVvQ6Ls449hgP+mUIEiNXC0k3zmUbk6UJT1a/7ai1KWThBOhM06rYZuibpmDmcAmxQ1Y2qehh4BmgfDvO7wPOquhlAVXelWEZDDPiXhZYV5HsVQJTsdDrobR2gRA4gEjw890HLiq5uERocjuBhOcIQbXA7UtSvA5feYiWGthm6IOlQDscAWwLeb/UdC+RrQLGI/ENEVovI91ImnSFmtoy8lefzC5nepyRs8LteIZacBNrsQUQbnjuRbhHiye+8UFtS0q+3jLyVvAQ01CwtGTojHcoh2KjRvrc7gZOBccD5wFQR+VqHikR+JCKrRGSVCa6XfkaPv4lZJf1DL/mocuW+A9yyp7lj/mdVPCJt9iBC7VvkeTSoBdM3Dn2VvJa8yLVEtH4S4UlYv4bQfXv0+Jto8XSPT1Dx5tQwGMKRDuWwFdosmg4Atgcp86qqNqrqbmAlMLx9Rar6uKqOUtVRpaXxT7cN8dPkbA56XBWuOPZ2fnVLNd+5bRMXcA49XZ4jA3SQPYgWS3AG8QR2WcIl+w+09ZHYvYeF29/mta3bvP4PnRFmBgPePM1RkrB+7RUvdN9uiSFgYXvM0pKhM9KhHKqAE0RkkIh0A64GKtqVeQk4S0ScIpIPnAp8mmI5DVESbqnCchczbcx1rX4Qp9S+wGErwKQ0CE0idAtyvEWElfn5HXwk8uQwRXqAyfUNHWcmQdjrdIS8tuUuDno8DCnr1zHI1hFPTvx1GLKalCsHVXUBNwOv4f1hPKeqH4vIj0Xkx74ynwKvAuuAfwFPqKqJN2xzlm6aH3SsVfXmIAj0g3i4pIhDYSPwASI0hRi8w+1HBPpGBF06CjJTCSZvNKSyX08YdGP8q2GWK84KDNmOyedgiIuqinmsX/8ATxU72BHiSVzVt8FrNWIRED46jhCr/VtcLN/aftWm47g/s6Soo39DBMrho3Z5HuyWz6F8YXlsEc99BGujoWti8jkYEk5VxTwe23Y/q0q7hx1sRQBHE8SWv60DuR4Pk+sbOhz3KG3CgYd0fOtEKSVk2SaJHAnZHU/8cpM21xAeoxwMMfN89RxWFYZXDIlGVJm+u45xjUc2jFXhSyllVUEjD5UUscPpODJDiVK27h7lQmtoQmVONM9Xz4HC+H66lscoB0N4TGwlQ8y8UuhMnGKIcHlTgfe7d2t3TJh39jSml/ZmR45XJk+Yje72FVq+fQm/1dMvqhfbOnppQu675bZ1Gw3pxygHQ8wkNPVkpIOdCIsLe7Y5tEv6sHTT/M43uIOiHOVyM6t2TxurJztHZ03EfS90e2zdRkP6McrBEDV+c9SEdZ4ojSICB8dD6mDLyFtjt9sPiMo6s6So9XBf3R1bfSkh/ju/12HxeImxWDKExigHQ1RUVcxj2Oo76EctV+7bH7+HcQznB3ZaQajY80Hn26uq9HK5Q6cfFeHZwp6tyYh2SZ+o5UoVQw59Jf77LsLiwh4mjIYhJEY5GCLCP1sYtfo2uovXA/nOugZOa26OKwxFRN7MgfhyP/t5vUc3XmVFWOObXI+HWbV7eGvLNqbsCTPD8IUTb9ZubBl5a3RypZC5Oz7iqkQoZhEW1/w2ITIZsg+jHAydEui8Ju1MRVfl5ob1cvajHgv1tHVcC2WSGo481dYUnxA8cmvbC7e1brq3pDisrA2WxUcnz2zNSGdHjtJab/5rTwJ8lHz5NwyG9hjlYOiUwCQ+gcwtLgobCrs1bJKrmCuP+1+uPO7n9G3xtMZD8g/a4WIcBZLr8TBtd11rFrnygWVex7swBEq3rCCfvY7wXV4RWysGALevVQdj2oBvizFoNYTC+DkYOiUwiU9gZrbOnltFvIphwqAbWVL9CGo1IU6LQrebyfUNjGtsYllBPk2dDNgAqGLhNWN9qWePiJP9eESYUtqbGX1KaI5ghhNDqO6UY/m25BMT2yA7IiQYEo9RDoZOebrHUUfCY0BUNvYeRz2Lax5AHO7Wp9R9Tge/7NMbjwoPF/eiJZL6fHGWwqb6DHNucwTniEeZMPDm6OpOA7uklCdK4o/MajCEwygHQ1BmVC5i6ab5eBz1SGn3ONYfBLE6bjq7LeGOvr2RaJ9ck+GNrUqeR5m6u57/vMH+6TMfPW4sL8rKhN2LGZWLmDbG/u02pBaz52DowIzKRSyumYM6673jT4xjkHoswi1beOu2waq3CEUeD2c22UCWCHiB9xJ333w5pQ2G9hjlYOjA0k3zESuOZQtVcjweLuOMmILY5agGTfITUH3C2el0UKAHMyOkhCR2SUlN4h9DEIxyMHQg3ixh+aqsqdnKT2qWe3MPeKLIBa3K3bV7mLm7jl6ujg5r6smhhw5GVRKa5bOfy003cXH8mrsTU2EGYZlNaUMQjHIwdCDekNX+BD19dTeumnfJwxXxKN7L7WFucRF3lPYmX5Wr9u1vTQd6VIubbxw+FhefIqKRuFdERKC/RZHuz4DZQ2KXvzyQAW02pBqjHAwdODbv5LgsHP2d6oWCQl5lhdceP4JR3Onx0OSw2JHjRH0xj17q2YPJ9Q2sq97CT+sa+Kz7hvAB9qKcTliqXLL/AHOLixg2sIzzy45m/foHIj4/HfTQExO6tmYBGz+blbD6DNmBUQ6GDmxuXh3Xw6l/t2BeSUFYfwRVyHHntDrF9VDtYNZ60LKYW1zEsoJ8ZpYWe0Nxd0KknTrX4+HKfft5qWePNgppbp8cW8ccavJ8mdCNfI8I95fm27rNhtRjlIMBOBI7Saf1Akdd3PXN9CXdCYcALdZhBNjhdNAQQpHsdDo6D5PRWqkvl0O4J2tfEL7pu+tYmZ/fod6DlmVrC55494SCYfc2G1KP8XMwtMZOypPDIN7N2R05QbpGJ7mXW4nUWU28/3QWPKOfy83OThRNMBlaFUSQNKH7HRZ3lPYOuXqWjAE4UeS78mjOaU54vXZusyH1mJlDF2dG5SKm7HqIUwcdxdgBR7OsIJ/J9Q0dTUlVcUSznp+gZY8cVSbXN9AvVPTWMDKFc9HwiKBhdrTtnEf6vKZ9SbHnzYTQIYbUYWYOXRi/s5vkeJ8RduQ4mVraG6dHcQV52nanwWHNpd5IqnsdVseZiyp5qiFDYwixZU1TTw5XDLoxJnmTzYzKRbzR05EU50FNhgOJIWMxyqELs3TTfMTZ1qGqRYQWh308hdUS9lohlpTCxUxSjWjzOqC4v1IGdj/HtuEklm6ajzqTM+EXR+KXqgyZi1lW6sIkZY053U+fqliR7o0E4F9hElGqD62wreVOcvcF1LbtNqQeoxy6MPmuvMRXmuZYSbEuJbWpw2qxreVOUvdCRFhS/Ujy6jdkFEY5dFGqKuYxpW4HuWFiGMVMZ7OHJM4uCiNJHBTB9e1quTNh0I1Y0YQjiRarKXl1GzIKoxy6KGVrZnN50z6m765rDU9R5HZ7l2QCiWUgDzd7SKJiyPV4Oo/0GuH17WqtNG3Mddy6u4kid8e4U4nBbEobvBjl0EXpq7UAjGtsYvnW7dxbuwfVKJZkYo16l6iASO1ksVQ5KBLSkc5XjJN1cKedXj05TLCptRLAdw98yT83b+OqffvTv8djyFqMtVIXw5/ERweVYeFVBkUeDwdEgueDFvE+TAaO5zFs+CabziyTRBWnuxurHOsBCR0dRNXW1krgzQS3uqAxtqx4EXDx0z/j5e/OSXi9hszCKIcuRKtfg7OFQM/kBkcna9i+8UfV79ScuAEplBNzVHRyco4qqorL2dJ5yCgRb2wpG/PocWN5mRXJUdAiVB9+M/H1GjIOs6zUhYg3iY93LErsMob4QmgkDVWcqsFnRSGw62a0n5c8H0XVnmgRY9JqwCiHLoVdM34J2mHpPGFL6eEc5cJIZGeSrrwE25ryGlKHUQ5dhKqKefRzueKvqDNLoJg2qjVYbLw0Yu9N3lRYUtl99mRIPlEpBxEZJCJvishGEXlQRHIDPvtX4sUzxMuMykUMe/JsflD3MPv88YmSRE4sGwiJ2tzuJABfe38Op8cT5l7Ye+YwYdCN4Em2ArP3PTAkn2hnDo8BzwNXAiXAmyLS0/dZTiIFM8SPfwNanfWoCE2WFXwgTtA444rSTLWDT0UchOvIClyy/4D3ej6z1wn7D3Ryhs1J+tidAffAkFSiVQ5HqeqjqrpaVScCL+NVEL0wvcl2RLoBnagvLvrFJG/O6EQQzpRVgGcLe3rL+JIBvdSzR8hr29UBzs/STfOTvu6W4zbPel2daE1Zuwe+UdV7RaQFeBPoGfwUQ7rwOOojfMAUkqrbQywd9fJ4aHQkcNvLn+CnXVhvDXLtg5ZFrrrJ9XjaZIKzc7huP5F/r7GTYx1O8hUMdifaX+bnIvLtwAOqej/wNPDVhEllSAiRPAHnejz8577DdG+/hh1syaezZaBQ6TmDDM5e3wM65IxOBK2hJTrZz9hrWUze3UJeS563qKuYK4/7ma0d4CA1M5tms+XQ5el05iAiI1V1je/t1cHKqOqDIvJsQiUzxExVxTyOX3MXv873MLW0pI1NvHqEHI8Tt+Mw/VxuJtc3MK6xiTMO5nNXaYnX7DPY3kEU+wMS4mk9kBbwJvDpBEvV66wXhRI5GOHeR6lLufbWz7k24prtwYRBN/J8zQN4rIDseAnxJjQYjhDJzOHvIvJNAFU9pKqHghVS1W2RXlRELhCR9SKyQUSmhCk3WkTcInJFpHV3daoq5jFs9R0Uc8A3xrcdLCxgyp561lVvYfnW7Yxr9EbhfKln/hHFEIxIN5tFIlugimLzuijKyLEHI3EQU6WP+wSqKuaxc/rxeKb1Yuf046mqmBfVtQJJVb++6EATl+9vaKuwkxCzauqCSQmtz5BZRKIcngZeEZEJ7T8QkTNF5K1oLigiDuBR4EJgCHCNiAwJUe4+4LVo6u/qlK2ZTXfxPlHOLS7qsGyjljKzbwHlA8va/L2Xl2fLp06PP5heonNXi/BZt2q2fDaTftRiCfSjlqGr74xJQaSyX5etmc3b+Un+vkR4lRVGQXRhOlUOqvrfwL3AMyLyYwARKReRvwIrgWgXQE8BNqjqRlU9DDwDXBKk3E+BpcCuKOvv0lTlNzJ2wNEMG1jGDmfo9JpB/xKFP1hfIutLwkDosdz8vqRHm2N5cpiyNbNjqS5l/bqv1ob+bhPIQcviHfeKpF/HYE8islZS1btEZBvwOxG5BjgD2AJ8H3gqymse4zvXz1bg1MACInIMcBkwBhgdqiIR+RHwI4Bjjz02SjGyB3+kVY+jHunb2x4zgIC9aYeCR5IYyTUOR7qdTgfLCvKZW1zETqeDfi43P63by39GX1XC+rWvbMi+vbigJHrpYmSX0wZ9yZAWIrJWEpES4GuAGzgLeA84QVUXqmq0hurBelv758zfArerqjtI2SMnqT6uqqNUdVRpaWmUYmQHgY5unSa6STSdLPX4H/iTqhj8F4oRpyrT+5SwI8eJirAjx8mM0pJYAs8lrF9D+L79UO/8lH7PJghf16RT5SAi04CNwE+AB/DOFkYBD8Z4za1AWcD7AcD2dmVG4V3GqgauAB4TkUtjvF5WsyTOSKspwQ4zmRC0iHTYwD5kSSyB51LWr/c5Ung/TV7pLksky0q/BJ4AZqjqlwAishl4QUSOAq5V1WhGpyrgBBEZBGzDax773cACqjrI/1pEFgIvq+qLUVyj65DOAGk2HvQjJkQbYgg8l7X9Wk1e6S5JJMtKX1fVSX7FAKCqlcA3gXOAV6O5oKq6gJvxWmt8Cjynqh+LyI/9G96GyOnrijP8RMyRVOMkkddURT2J3aCN1tEslf26V9KD7rVFMEtLXZFOZw6q+kWI42tE5ExiMMlT1VeAV9od+32IshOjrb+rUFUxj5/V1XFHPJvQIvRvcaXE+gU4snnslzchUVkFxJ2wCK+xhtBIVb++fU89vygtSd3MTWBxzf1Mw96e44bEEldgG1XdAPxHgmQxREFVxTyGrr6Ti5uaInISs8I8qe90OpKT2CNYBp9kJG7w65pEKAbF9iE0vtmcIkUegEhiAiQaMoe4x4TA5SZD6ihbM5vKHk7GDji6cycxVa7ct5/+ruBGMgLeEBVJWFoqcrsRVfq3JCDRUAoQsLViAHB5lLwULy2BctW8S1N8TUM6MZngMpSq/MZWE8zWZZpQg7sIK/PzmVzf0CHpDaqtoaxb64hWSYS57l7L4t7aPUyub4iuzkQRQ1Y6u6+vF+p+xh84kNp9IhE+6b7B9vfGkDiMcshQfltS0jGGUJhlFf+ewvTddfRvcR3RAcGWeaJdnglTXkWYWtqbO/qkcI08kBjaYvf8ybuklJX5qfV1ADLi3hgSh1EOGciMykXsckb51YkwvY/Xs/aWugZytWPe5rgI8xTbIoJGEgzPJtg9f/KWkbeyM1UGBO2w+70xJI7M+cUaAK9ieL7mgZjSRB60LOYWF/FQSRGHLJv6KASa1qbavNaHePLTct1IGT3+Jpyebmm7vlla6hoY5ZBhLN/4YNs4/lGy0+lIvNmqKn1drsQM5iL0d7n5sHqLSXEfhjTpTUQwS0tdBKMcMoiqinnsd8YXKqNXlLkRIkKEXU5nwtbAdzgdDBtY1nnBJJEJHsFuR/rSeJqlpa6BUQ4ZRNma2fQLYY4aEaocSlL460SH/FbfX1jieXwOs2yVijSc8RJXP4gTuy+7GRKDUQ4ZRF+tDW6OGgXNdomHpHok13MM58btDR1CSarCsXknx15virilrgFHmtaW1Goy+w5dAKMcMohdUsq4xqZWc9QOm7edbOKGcoJLF//cHHFm2Y5EqxgiHEhFoPrQCtsPfmc0WdxTu4d8jyflm/cisLhmju3vkSE+jHKwOVMXTOK8J77BsIVDuXZAN2aUlDC3uOjIpnK7TG6nNTfTy9XxiTzX42FyfUP8nrXRDkShlm6AmSVFqdl0jnLgFKvF9puuG0b+Co9Ct3RZdGXAPTLEh1EONmbqgkm8ygp25VioCF/mOFhSWNDWKzoQEf6Vl8fKzdv49a499G9xtYaumL67jnGNTXSPJ3+nfyCK9Kldle/s20/3IMtgHhGeLezZ+b5CMKI8p7/LTS93dEtxdt90rdjzAdNKe9PgcCQtjWpnqM3vkSE+IkoTakgP77r/wcGcdmannQwCHuC0shM53zGKp2qWc5RuRwSWFeQficPUCZaqN9ZS4LViWOO38Lpj5KlyKJqgewmKruqva4fTEfUMxe6b0u+6/0FL+76RYvLSZU9rSAlm5mBDZlQu4pQnTuHLaL2gAURozmnmBfkn3z62N88VlLCsIL9tHKZwhJodxDBY+2cHrU+3kRJN2XDLXAHhwaOZoVgeBxNiCNmdKqoq5sXWNxJMsxiHuGwm/T3M0AZ/TujmnOa4np5FAEcT95fmc29Jccc4TGFO9CTYLDWZCGE22mO4tqXKtNo9XHTAvr4O69c/kG4RvIjwUvWcdEthSBJGOdiMpQnOCX3QstjrsOnXnIBlCcUXVDBBSxwKXN60j7I1sxNSXzJ4qjjKmVgSabHS54xnSC5mz8FmdKlNvoQk+olvX6Q9fg/yvro7rnqSScqy9hm6NDZ9pOyaVFXMw4rHmigUIZy9LHcnSYIyjQQomwMiLCvIZ5f0SYBAyUFsFnXK7DtkJ0Y52IAZlYsY9uTZ/KDuYSIyuEyAn4EA79VsZ1btnrg8rrMNl2Xx65Jitoy8Nd2ihEST8QARKyIsqX4k3VIYkoBRDmnGH4JbnfVei5pInn6j8DMIbW6oVPZwdvS47gKoAu78kM3d53DwSg/7xg+KJGd4KjHhNLIToxzSzLKND8cVgjssInQPZeopwtziIgDGNTYxub7B6yiWCAVhlzpCICh9PQdCL87YPCy1R+21rGTCeGcnZkM6jVRVzOOgs4mYMvdESDintx1OB6ccO4Bmf+KfRFnAxFtPCmYwu3LCPxfZ2UN6v8NeygG6mCFFF8HMHNJEVcU8tnw2M/lbiyJoqKuI0OywUh5+odPwTMmUJ0KLJjt7SKczXHdo1CwtZRlGOaSJ9esf4N7SwsQ6nIXAJibxraRFHlUk0hmJYmsP6TObmu23PyTCa5tmUVUxL92SGBKEUQ5p4qliR+Rey1lIOhRExBv+qK09pN/Kz7Ofxse73DV09Z1GQWQJXXd0SjM7E+HIlOI4/h2unUlEMZj2d7lt7SGdkL6TJB7onW/re2eIHKMcUozfpyHkPkAUCHDVvv3pGaijDY6XbFSRRJh4qnJcy2Fbe0jnumxqZusLtPh4iSvdkhgSgFEOKcQfVE+d9QkxUCp0ezjpUAbEtknFEogIvXy5K+JSRiK8l5fH0z36Jk62BNM3/wzs5AfXBhEWF/Ywm9NZgFEOKWLqgkks3XxfQoPq7XVY/KK0ty3Xn9NBg2WxfGsCvL5FeMjnA2JHNjevTqb1c/yIsHzjg+mWwhAnxs8hBUxdMIllrIwt61k4RCILt9FF8N/dcY3ezeSpfUpoidEstsnZnEDJEksm+BTscybuIciQHoxySAHvuFfQ0onTlSF+FDj1uAE0BSqDGBWyeGy6rg94W2rnqYOX0xZcyXs3LE63GIYYMSNWCtjltP8POWYiWd9XJc+ToNAc4RChyQpw6jPLbelDhAPymdl7yGCMckgBGTlERTqQRxgosLtq+iyrYkAt+/o5ZMqPVgQW19yfbjEMMZIp/SyjSYTZaspJ8FN3g2XxfM+eGfM0b+fwGYMPHZ8xSlbE7IplKkY5pICk/YzbDxDRDhgJyAsRMSK0WJmhGNTm4TOevelFMmc+qkxdMCndQhhiwCiHJFFVMY+d04/HM60X+Zq8p6deLjeiihVLikyRxC4fZQOqDOz2LaaNuS7dkoQlx5OTbhEiQ4RKXZFuKQwxYJRDEqiqmMfw1b+gH7VYgjenQjIQIV+VVZt24In1SbKrDPpR8PJ356RbhLBUVcxj+u6dGbO0tM8SszGdgaRFOYjIBSKyXkQ2iMiUIJ//l4is8/29IyLD0yFnrBy/5m66yZEQAnuTGGBvh9PhSxtpBvlE0NcV+4Cbqn59/Jq7Gd/UGLOcKUeEVzb9Jt1SGKIk5cpBRBzAo8CFwBDgGhEZ0q7YJuAcVR0G3A08nlop46NI97e+XlaQn/Rh+40e3cG4w8VNrsfDfzjOiencVPZrf//qb8u8DsFpcij3zr0o3WIYoiAdTnCnABtUdSOAiDwDXAJ84i+gqu8ElH8PGJBSCWPk4qd/RvXhN5FBZVjAKc3NrM3NTW7OBhHuKO2dvPoziXgmUKqM33+Aqbc8FuvVU9avlxXkc1+fYm+Wv1j2mtKBCM8UbuZw5SLb7+cYvKRjWekYYEvA+62+Y6H4AfC3YB+IyI9EZJWIrKqtrU2giNFz8dM/o+bwG4h4f6weXwC3VORsiDxPQZYTzy0QoaJnXAHjEtavveIE79szKhcxtbQ3DQ5Hxjn6eSxhWbXxe8gU0qEcgvXmoAu9IvJNvD+i24N9rqqPq+ooVR1VWlqaQBGjp+bwGx1blkE/XAMctCyWbXw41tMT1q8hdN9eumk+rgwxCQ5Gs+Xm3YcmplsMQwSkQzlsBcoC3g8AtrcvJCLDgCeAS1R1T4pkiwnv02YSLEdUIRE5CgwR0+xsinX2kJJ+7cmAoHud8ZOeVcZ6KQNIh3KoAk4QkUEi0g24GqgILCAixwLPA9ep6udpkDEqlm6an/hZgs9MsZdHyfHHJbKD6WIqHefSgQjP1zwQy+CVkn5tZ8/tiBChxbJYUvMboyBsTsqVg6q6gJuB14BPgedU9WMR+bGI/NhX7FdAb+AxEVkrIqtSLWekzKhclJynOd968l6n40jYaTssU0Ujgx3kjQGP5Y56eSlV/XrCoBszU+m2x4KlNca81c6IZkNHA0aNGqWrVqVWh0xdMIm/sZJDGbwGbAiOqLJu4kdH3ousVtVR6ZClfd8e9ofhqCMLlhtVuVTP5u4bYrYQMySAUH3beEjHSFXFPN5z/71rKYYseZCIhH429iHIzZavQYRKXZluKQwhMMohRsrWzGan05FuMQxJINfj4Xv19lUOBx32lS1a9lmYwHw2xSiHGHm8xKRBzEZElTtq93HiiT9PtyghsfOsJmpEqGAFVRXz0i2JoR1GOcTAjMpFLC60cW6CJAb6y2pUmbZrH2WD72T0+JvSLU1Ibq7bm1VLfB4RvvxsulEQNsMohxhIiulqosiUcAo25fclBVTs+SDdYoSlNHdsukVILCLc0bc3f6u+1ygIG2GUQ5RUVcwDR126xQiNUQyxI8LOHCd/Y6WtbfBPv2UhWReFV4TFhT14dNuD6ZbE4MMohyiYumASv9j12+RldjPYgkOWeGeHhtQiwuq8HC5++mfplsSAUQ4RM3XBJF5lBTtznObpvAtg5zAVMyoXIVm059AGEWoOv85V8y5NtyRdHqMcIsDv05CKCKsGe2DnMBVLN81Hs9m/RoRPum8wJq5pxox2EVC++g7j09CFUPWFqbApauNZTcIQ4TVM7ul0YpRDJ6y75xyuPqZPusUwpBg7J6Q5Kpv8HMLQLMKwBcNtbRyQzRjlEIbn/vAA5YfX8kW3bmafoaugirgL0y1FWP6nrp7crhDKXQS1PCzZbCK4pgOjHELw4vvbOK/6QaMTuhoiqGOfrS1mxjU2MX13XVY5woVFYMnm+8weRIoxyiEEa5c9TokcYFlBfrpFMaQYEag+/Ga6xQjLuMamdIuQWkR4kRVGQaQQZ7oFsCMzKhfxdr/FDHf6EnuZ6UMXxL5P5S8X5PNwSVG6xUg9lsVLnhXcnW45ughm5tCOGZWLeLFmNjtzHKgIahRDF8We3/uMykXcVVrCji7qb6Mi3PXQcSYPdQowyqEdS6ofwWXZ96nRkAJUGdjtvHRLEZSlm+Z3bX8bERYX9uRHhav434e+kW5pspou3Ms6MqNyEWp1sbVcQxCE0f1GpluIoNjZcztl+FLmvlYojFww1FgyJQmjHAJYuml+V5ypG9oj3hmkHbGz53bKEaHFEpZs/o2trcsyFaMcAjBPZQY/dp1BThh0Y5exYI0YgZrDr5sZRIIxysHHVfMuRWxsoWJILXadQE4bcx3iMebVHRBhyeb7GPmHcqMkEoRRDngVwyfdN3RJ6w9DcAo99n1QuGLgzagnJ91i2A8RWhxeh7nTHz893dJkPF1eOcyoXGQUg6EtqkzZbd+ETtPGXMfA7ueY5aVQiHCg237KF5abvYg46NLKYUblIhbX3G8Ug6EDe6VHukUIyYzKRVQfWmG6bThEWvcifvDIWemWJiPp0sph+cYHEasLBDAzRIcIDxUXpVuKkCzdNB+xWtItRmYgwr961FO+sJzyBSPNfkQUdFnlcO/ci9jnND8wQ3CanM3pFiEkxqouSnyzCKwWlmy+z8wkIqRLKod7517EXwprzHKSISR2tggyvg5x4JtJnLLgGzx0/ynplsbWdDnlsO6ec3ihZw3alUMQGDKaCYNuxNkV8jkkCxGaLYv5fZooXziUEX8YYZabgtClRsiqinm81LOa5mzOv2tICHZ1ggO46EATPYypUvz4wnC4HW6WbL7PWDe1o0sph42fzeK5wp5mOckQAfbtI2VrZrPXzHwTS4B1U/mCb3Dv3IvSLVHa6TL5HKoq5vGb0nyjGAwRYt8n875aSz/X0d6w3YbE4ptNPN1rM08vHOo9psqlnMPdNzyWXtlSTJd5/Hi5+jccNorBECG9XfbtK3ulB5PrG7pOmtB04FMSiIBl8aKspHzhUMoXDmXYgq4RoqNLKIepCybxfGGemTUYIkOVW+tq0y1FGIRxjU1ctW+/URCpIkBZqIVvj8KrLL6RpcqiS8xLX2UFSJfQg4YuQC/dDwJ31jVQ0bMHzeahJ/UE3HNLvMpiycL7/B8i7n6s+8Hr6ZEtQWS9cpi6YBIHzY/HEA0i3NenmIvTLUcIdkkp/fDObEzftgntvgd17KDcv2fhI8fTjUsG/oxpY65LpWQxk9XKwURbNcRKg+VItwgh2TLyVnqtvpM8OUyuqpk52JEg30mLo6XdDCMAG256Z61ymLpgklEMhjiw71r+6PE3UQVs+myWUQyZRqjvS4QXdSUvtpttBOJQZeLBr/A///3XJAnXlqxVDu+6/wE59n36MxjiYfT4m/hB7Z9BTJylrKETRe8W4cm8TTwZRoGE4/hDPXjhpvciLp+1u7RfOrO2aQYDYALwdUkCTWyj/NvQ/QDlC74R8aXSMoKKyAUisl5ENojIlCCfi4g85Pt8nYiMTIecBkM0pLpfhwzAZ8xbDcHwKYkRESqIlCsHEXEAjwIXAkOAa0RkSLtiFwIn+P5+BPwupUIaDFGSjn49YdCNWJ52S6eqfPXwYaMgDMERwR3hPlU6Zg6nABtUdaOqHgaeAS5pV+YS4Cn18h5QJCL9Uy2ooesSwzZvyvv1tDHXcflxPyevJQ9RpX+Lizt3NfKnbXu5d9ceerncXiURiaIwysTQjnRsSB8DbAl4vxU4NYIyxwA7AguJyI/wPoFx7LHHJlxQQxdFlUv2RZ3sJ2H9GiLv29PGXBfUbv5i3x9404ou3TTfu0fhyQGrpe3ep2+28UW3bm02RR2qdPOoN4qxsYrqcqRj5hCsl7V/bImkDKr6uKqOUtVRpaWlCRHOkMH4n5KDPQVH+vSsyqX7mrn7li+ivXrC+rVXlMT17WljrmPdD1by0cQP+ej7a7jy2NsQVzGqIK5irjj2dn7Z7+fcXnuY/i2u1lnIXbvq+L+ardxeexjxtBM98F57PBS53W2PBXxmZiWZSTpmDluBsoD3A4DtMZQJj79D2v2Jx/xwEoCQ5/Hwq911jGtsYllBPvf1KabBF9Y6123xrcYWVucLO5wOwEJR1JODWC6847MwsNu3ePm7c2IVIjX9OgFMG3Md0+g42xjNTVwbpPy1wBcBsw/LXcyEQTcybcx1VFXMo2zNbPrqbnZJH7aMvJXR429qPbeqYh4Pbfsta/P8eyNCW33Y/n2chPq9q9p/LEgFqvQ43DOioqIpHpxExAl8DpwHbAOqgO+q6scBZcYBNwMX4Z2aP6SqYXP6jRo1SletWtXmWPmCb9i4QwgK9NDBvHfD4nQLY+gEEVmtqqPCfJ6Ufg3B+3ZXozMlBO2Wz3wIUOhRzj+wn9d7FLQ+MHRVehzuybs/erfNsVB9O+XKwSfMRcBvAQfwB1W9R0R+DKCqvxcRAR4BLgCagBtUNeyvw/yADMmkM+XgK5Pwfg2mbxuSS6i+nRYPaVV9BXil3bHfB7xW4CeplstgiAfTrw3ZRNeeYxkMBoMhKEY5GAwGg6EDRjkYDAaDoQNGORgMBoOhA2mxVkoGIlIL1AQc6gPsTpM48WJkTw/hZD9OVdPiaRmkb8dLJn9HkZDN7UtG24L27axRDu0RkVWdmR7aFSN7eshk2aMh29uZze1LZdvMspLBYDAYOmCUg8FgMBg6kM3K4fF0CxAHRvb0kMmyR0O2tzOb25eytmXtnoPBYDAYYiebZw4Gg8FgiBGjHAwGg8HQgaxUDp0lek83IlImIn8XkU9F5GMRmew7XiIir4vIv33/Fwecc4evPetF5Pz0Se/Nlywi74vIy773mSJ3kYgsEZHPfPf+9EyRPdGIyGzffVgnIi+ISFG6ZYoXu//u4yHUmJFUVDWr/vCGS/4C+ArQDfgAGJJuudrJ2B8Y6XvdE28egCHAb4ApvuNTgPt8r4f42tEdGORrnyON8v8/4GngZd/7TJH7j8APfa+7AUWZInsS7sVYwOl7fZ+/3Zn6lwm/+zjbF3TMSOY1s3HmEEmi97SiqjtUdY3v9X7gU7y5hC/BO4Dh+/9S3+tLgGdU9ZCqbgI24G1nyhGRAcA44ImAw5kgdyFwNvAkgKoeVtUGMkD2ZKCqy1XV5Xv7Ht6sdJmM7X/38RBmzEga2agcQiVxtyUiMhA4Cfg/4ChV3QHezgD09RWzU5t+C9wGeAKOZYLcXwFqgQW+JbEnRKSAzJA92Xwf+Fu6hYiTLvN9tRszkkY2KoeIk7inGxHpASwF/kdV94UrGuRYytskIhcDu1R1daSnBDmWru/CCYwEfqeqJwGNeJeRQmEn2WNCRN4QkY+C/F0SUOaXgAv4c/okTQgZ/31FQhRjRtykJRNckrFFEvfOEJEcvF/yn1X1ed/hL0Wkv6ruEJH+wC7fcbu06QxgvC8dZi5QKCJ/wv5y+2XZqqr+p60leJVDJsgeE6r6rXCfi8j1wMXAeepbzM5gMv776owQY0bSyMaZQxVwgogMEpFuwNVARZplaoMvl/CTwKeq+mDARxXA9b7X1wMvBRy/WkS6i8gg4ATgX6mS14+q3qGqA1R1IN77Wqmq12JzuQFUdSewRURO9B06D/iEDJA9GYjIBcDtwHhVbUq3PAnA9r/7eAgzZiSPdO/CJ2ln/yK8u/lfAL9MtzxB5DsT75R3HbDW93cR0Bt4E/i37/+SgHN+6WvPeuBCG7ThXI5YK2WE3MAIYJXvvr8IFGeK7Em4FxvwrtH7+9/v0y1TAtpk6999nG0LOmYk85omfIbBYDAYOpCNy0oGg8FgiBOjHAwGg8HQAaMcDAaDwdABoxwMBoPB0AGjHAwGg8HQAaMcDAZDRiMiC0WkOoJyE0VEfeEn/MeqRWRhwPtzRWS6iHT5sbHL3wCDwZDx3A1cFuO5l/nO93MuMA0zNmZl+AyDwWBTRKS7qh5KZJ2q+kUc576fSFmyiS6vHbMd3xRZRWSwiLwmIo0isllEbvB9fp0v6csBXzKRr6ZbZkN2END3hvr63gHgOREZKyKviMgOEWnyBQP8uYg42p1fLSJ/EpEbfQl8DorIGhH5ZrtyHZaVROQrIrLMV3+tiMzFm5ujvYyty0oiMh3vrAGgxSe7+kKo1IrInCDn+5eqBsdxq2yJmTl0HRYD84H7gUnAH0TkBLzT6ClADjAXbxKfU9MkoyE7eQlvXKD78IZ6H4w3VMnDwEFgFDAdKKVjpNxzgJPxhjI5hDce1N9EZLiqrg92MV9spdeBPOAneIMp3gRc3omcT+AN2PcDvOEq3ACqekhEFgA/FJE7VPVgwDk3AStU9bNO6s44jHLoOsxW1acARGQV8J94O/Yg9YX+9UUlnSsix6lqTfpENWQZD6nq3ID3//C/8AWU+yfe7G3/KyK/UNU2uUKAM1R1s6/8m0ANcCdwXYjrXY83f8fpqvqe77y/AR+GE1JVt4rIVt/b/9MjyZAAfgf8HLgSWOSrcxhwGnBNuHozFbOs1HVoTeaiqvV4n6be07Yx4f1PP4Ghjw2GeHkh8I2I9BeReSJSAxwGWoCZeNO29m137nt+xQCtWdCWAaeHud7pwBa/YvCd5wGei7UB6s0G+BreByo/N+FNIJX08NnpwCiHrkN9u/eHQxwDb64GgyFR7PC/8JmIVuDNIzETGAOMBu7xFWnf974MUt+XhM/y1j/MefHwGHCGbw+lALgWWKDetKRZh1lWMhgMySYw9PNX8e4xXKeqf/IfFJH/DHHuUSGObQtzvR3ANyKsKxpeAarxzhg+AHoCj8dZp20xMweDwZBK8n3/t/gP+DKc/VeI8qeJSFlA2Z7AOODdMNd4FygTkdMCzrOA70Qgn9/MNq/9B76lqXl49zpuBt6Ix4zW7hjlYDAYUsmneDeU7xGRK3z5rF8PU/5LYLmIXCUilwLLgQLaOq6154/ARuB5n6npRXiTOxVGIN8nvv9/LiKnisiodp8/iXfpazjw+wjqy1iMcjAYDCnDtz5/KbATeAp4FFgJzApxygrgAeDXwLN4B+YLVfXzTq7xbbzZ0h7Dqyw24d3j6IyXfedMwjsDqWpXd61Pph1kURrSYJhMcAaDwZb4HNveUm+eclsgIsXAZuC3qjo13fIkE7MhbTAYDJ0gIqXAicBkvCsuj6VXouRjlpUMBoOhc8bhddY7BbheVXd0Uj7jMctKBoPBYOiAmTkYDAaDoQNGORgMBoOhA0Y5GAwGg6EDRjkYDAaDoQNGORgMBoOhA/8fIxwR+4EWR2QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "\n",
    "ax1.scatter(m, x2, label='data')\n",
    "ax1.scatter(m_test, y_fit,label='prediction')\n",
    "ax1.scatter(m_train, y_est,label='estimated')\n",
    "ax1.legend()\n",
    "\n",
    "ax1.set_xlabel(xlabel='m',fontsize=16)\n",
    "ax1.set_ylabel(ylabel=r'$x_2$',fontsize=16)\n",
    "\n",
    "ax2.scatter(rap, x2, label='data')\n",
    "ax2.scatter(rap_test, y_fit,label='prediction')\n",
    "ax2.scatter(rap_train, y_est,label='estimated')\n",
    "\n",
    "ax2.set_xlabel(xlabel='rapidity',fontsize=16)\n",
    "\n",
    "plt.savefig('Plots/x2_LO_Gaussian.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "1292c0f8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAELCAYAAAD3HtBMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtY0lEQVR4nO3de5wcVZn/8c8zQycZEDOIgSWTQKJiQiQhhAHhlyg3IZEgJIAigr8gq4gsqLtuNKz+BBGWICyILhezyAYE5SZmI6ABCa5yW5kwgZCQQIAAM2EhAkEgQzKZeX5/dPfQ0+lrdXV3dff3/XrlxXR1dfXp5px6us459Rxzd0RERIrVVO0CiIhIbVIAERGRQBRAREQkEAUQEREJRAFEREQC2a7aBaiUD37wgz5mzJhqF0Pq1LJly/7q7iOq8d6q21JOuep2wwSQMWPG0NHRUe1iSJ0ysxeq9d6q21JOueq2urBERCQQBRAREQlEAURERAJpmDGQRtTb20tXVxfvvvtutYtSN4YNG8aoUaOIxWLVLopI1SmA1LGuri523HFHxowZg5lVuzg1z9157bXX6OrqYuzYsdUujkjVNXQAWdTZzSVL1rB+Yw8jW1uYO30cs/Ztq3axQvPuu+8qeITIzNh5553ZsGFDtYsiEliY572GDSCLOrs5544V9PT2AdC9sYdz7lgBUFdBRMEjXPo+pZaFfd5r2EH0S5asGfgSk3p6+7hkyZoqlUhEpLzCPu81bABZv7GnqO1SfX/84x85+uijAVi8eDHz58/Puu/GjRu56qqrBh6vX7+eE044oexlFImysM97DRtARra2FLVdyqevry//TmmOOeYY5s2bl/X59AAycuRIbr/99kDlE6kXYZ/3GjaAzJ0+jpZY86BtLbFm5k4fV6USVd+izm6mzl/K2Hl3MXX+UhZ1dpd8zHXr1jF+/HjmzJnDpEmTOOGEE9i0aRNjxozh/PPPZ9q0adx2223cc889HHTQQUyZMoXPfvazvP322wD8/ve/Z/z48UybNo077rhj4LgLFy7krLPOAuCVV15h9uzZ7LPPPuyzzz489NBDzJs3j2effZbJkyczd+5c1q1bx9577w3EJxd86UtfYuLEiey7777cf//9A8c87rjjmDFjBnvuuSff/va3S/78IlES9nmvYQfRkwNG9TwLqxjlnFSwZs0afv7znzN16lROO+20gSuDYcOG8cADD/DXv/6V4447jj/84Q/ssMMOXHzxxVx22WV8+9vf5itf+QpLly7lIx/5CCeeeGLG43/961/n4IMP5je/+Q19fX28/fbbzJ8/nyeffJLly5cD8UCWdOWVVwKwYsUKVq9ezZFHHsnTTz8NwPLly+ns7GTo0KGMGzeOs88+m9GjR5f0+UWiIuzzXsMGEIh/mY0aMNLlGlwr9TsaPXo0U6dOBeCUU07hJz/5CcBAQHjkkUdYtWrVwD5btmzhoIMOYvXq1YwdO5Y999xz4LULFizY5vhLly7lhhtuAKC5uZnhw4fzxhtvZC3PAw88wNlnnw3A+PHj2WOPPQYCyOGHH87w4cMBmDBhAi+88IICiNSVMM97DR1A5D3lnFSQPvU1+XiHHXYA4jfoHXHEEfzqV78atN/y5cvLMm3W3bM+N3To0IG/m5ub2bp1a+jvL1IvGnYMRAYr56SCF198kYcffhiAX/3qV0ybNm3Q8wceeCAPPvgga9euBWDTpk08/fTTjB8/nueff55nn3124LWZHH744Vx99dVAfED+b3/7GzvuuCNvvfVWxv0/+clPctNNNwHw9NNP8+KLLzJuXOOOfYkEpQAiQHknFey1115cf/31TJo0iddff52vfe1rg54fMWIECxcu5KSTTmLSpEkceOCBrF69mmHDhrFgwQJmzpzJtGnT2GOPPTIe/4orruD+++9n4sSJ7LfffqxcuZKdd96ZqVOnsvfeezN37txB+5955pn09fUxceJETjzxRBYuXDjoykNECmO5LufrSXt7uzfaojtPPfUUe+21V8H7lyO1y7p16zj66KN58sknSzpOlGT6Xs1smbu3V6M8jVi3pXJy1W2NgcgATSoQkWKoC0vKasyYMXV19SEi71EAERGRQBRAREQkkEgGEDObYWZrzGytmWVNeGRm+5tZn5kpS57UBNVtqSeRCyBm1gxcCXwamACcZGYTsux3MbCksiUUCUZ1W+pN5AIIcACw1t2fc/ctwM3AsRn2Oxv4NfBqJQsn5bNw4ULWr18/8PjLX/4yq1atKvm469at45e//GXRrzv11FPDzuCrui11JYoBpA14KeVxV2LbADNrA2YD1+Q6kJmdbmYdZtahZUijLz2AXHvttUyYsM0P9KIFDSBloLotdSWKASRT8qP0ux1/DHzH3XMuJOHuC9y93d3bR4wYEVb56tcTt8Lle8N5rfH/PnFrKIe98cYbOeCAA5g8eTJf/epX6evr49RTT2Xvvfdm4sSJXH755dx+++10dHRw8sknM3nyZHp6ejjkkENI3iD3vve9j+985zvst99+fOpTn+Ivf/kLhxxyCB/60IdYvHgxEA8Un/jEJ5gyZQpTpkzhoYceAmDevHn8+c9/ZvLkyVx++eX09fUxd+5c9t9/fyZNmsTPfvYzIJ4j66yzzmLChAnMnDmTV18N/QJAdVvqShRvJOwCUtOfjgLWp+3TDtycSLT3QeAoM9vq7osqUsJ69MSt8NuvQ28ieeKbL8UfA0z6XODDPvXUU9xyyy08+OCDxGIxzjzzTC644AK6u7sH7g/ZuHEjra2t/Pu//zuXXnop7e3b3vT6zjvvcMghh3DxxRcze/Zsvve973HvvfeyatUq5syZwzHHHMMuu+zCvffey7Bhw3jmmWc46aST6OjoYP78+Vx66aXceeedACxYsIDhw4fz6KOPsnnzZqZOncqRRx5JZ2cna9asYcWKFbzyyitMmDCB0047LfBnz0B1W+pKFAPIo8CeZjYW6AY+D3whdQd3H5v828wWAneqgZXovvPfCx5JvT3x7SUEkPvuu49ly5ax//77A9DT08OMGTN47rnnOPvss5k5cyZHHnlk3uMMGTKEGTNmADBx4kSGDh1KLBZj4sSJA2t99Pb2ctZZZ7F8+XKam5sHUrSnu+eee3jiiScGxjfefPNNnnnmGf70pz9x0kkn0dzczMiRIznssMMCf+4sVLelrkQugLj7VjM7i/gMlGbgOndfaWZnJJ7P2TcsAb3ZVdz2Ark7c+bM4aKLLhq0/cILL2TJkiVceeWV3HrrrVx33XU5jxOLxQZSuzc1NQ0kP2xqahpIuX755Zez66678vjjj9Pf38+wYcOylumnP/0p06dPH7T97rvvLkv6+JT3Vd2WuhLFMRDc/W53/6i7f9jdL0xsuyZTA3P3U91di12Xavio4rYX6PDDD+f2228fGE94/fXXeeGFF+jv7+f444/nhz/8IY899hhAzhTshXjzzTfZbbfdaGpq4he/+MXAWuvpx50+fTpXX301vb29QDyl+zvvvMMnP/lJbr75Zvr6+nj55ZcHlroNk+q21JPIXYFIlRz+/cFjIACxlvj2EkyYMIELLriAI488kv7+fmKxGJdddhmzZ8+mv78fYODq5NRTT+WMM86gpaVlYP2QYpx55pkcf/zx3HbbbRx66KEDC1ZNmjSJ7bbbjn322YdTTz2Vb3zjG6xbt44pU6bg7owYMYJFixYxe/Zsli5dysSJE/noRz/KwQcfXNJnF6l3Sudex4pN584Tt8bHPN7sil95HP79ksY/6pXSuUsjUTp3KcykzylgiEjBIjkGIiIi0acAUucapYuyUvR9irxHAaSODRs2jNdee00nvZC4O6+99lrW6cEijUZjIHVs1KhRdHV1oVxJ4Rk2bBijRpU2tVmkXiiA1LFYLMbYsWPz7ygiEoC6sEREJBBdgUhgizq7uWTJGtZv7GFkawtzp49j1r5t+V8oIlUXRvtVAJFAFnV2c84dK+jpjacL6d7Ywzl3rABQEBGJuLDar7qwJJBLlqwZqHxJPb19XLJkTZVKJCKFCqv9KoBIIOs39hS1XUSiI6z2qwAigYxsbSlqu4hER1jtVwFEApk7fRwtseZB21pizcydPq5KJRKRQoXVfjWILoEkB9o0C0uk9oTVfhVAJLBZ+7YpYIjUqDDar7qwREQkEAUQEREJRAFEREQCUQAREZFAFEBERCSQSAYQM5thZmvMbK2Zzcvw/LFm9oSZLTezDjObVo1yihRLdVvqSeSm8ZpZM3AlcATQBTxqZovdfVXKbvcBi93dzWwScCswvvKlFSmc6rbUm8gFEOAAYK27PwdgZjcDxwIDjczd307ZfwdAa7ZWgdK5F011W6qmHO01igGkDXgp5XEX8PH0ncxsNnARsAswszJFkySlcw9EdVuqolztNYpjIJZh2za/wtz9N+4+HpgF/DDjgcxOT/Qjd2hd8HApnXsgqttSFeVqr1EMIF3A6JTHo4D12XZ29z8BHzazD2Z4boG7t7t7+4gRI8IvaQNTOvdAVLelKsrVXqMYQB4F9jSzsWY2BPg8sDh1BzP7iJlZ4u8pwBDgtYqXtIEpnXsgqttSFeVqr5ELIO6+FTgLWAI8Bdzq7ivN7AwzOyOx2/HAk2a2nPislhPdXYONFaR07sVT3ZZqKVd7tUapm+3t7d7R0VHtYtQVzcJ6j5ktc/f2ary36rYUImh7zVW3ozgLS2qE0rmL1I5ytNfIdWGJiEhtUAAREZFAFEBERCQQBRAREQlEAURERAJRABERkUAUQEREJBAFEBERCUQBREREAlEAERGRQBRAREQkEAUQEREJRAFEREQCUQAREZFAFEBERCQQBRAREQlEAURERAJRABERkUAUQEREJBAFEBERCUQBREREAlEAERGRQCIZQMxshpmtMbO1ZjYvw/Mnm9kTiX8Pmdk+1SinSLFUt6WeRC6AmFkzcCXwaWACcJKZTUjb7XngYHefBPwQWFDZUooUT3Vb6k3kAghwALDW3Z9z9y3AzcCxqTu4+0Pu/kbi4SPAqAqXUSQI1W2pK1EMIG3ASymPuxLbsvl74HeZnjCz082sw8w6NmzYEGIRRQJR3Za6sl21C5CBZdjmGXc0O5R4I5uW6Xl3X0CiC6C9vT3jMSQ6FnV2c8mSNazf2MPI1hbmTh/HrH1znV9rjuq2VF2Y7SyKAaQLGJ3yeBSwPn0nM5sEXAt82t1fq1DZJE1YlXFRZzfn3LGCnt4+ALo39nDOHSsA6imIqG5LxaW20eEtMd7ZspXevvhvjlLbWRS7sB4F9jSzsWY2BPg8sDh1BzPbHbgD+KK7P12FMgrvnfS7N/bgvFcZF3V2F32sS5asGQgeST29fVyyZE1IpY0E1W2pqPQ2urGndyB4JJXSziIXQNx9K3AWsAR4CrjV3Vea2RlmdkZit+8DOwNXmdlyM+uoUnEbWpgn/fUbe4raXotUt6XSMrXRTIK2syh2YeHudwN3p227JuXvLwNfrnS5ZLAwT/ojW1vozvC6ka0tRR8rylS3pZIKbYtB21nkrkCkdmSrdEEq49zp42iJNQ/a1hJrZu70cYHKJiKFtcVS2pkCiAQW5kl/1r5tXHTcRNpaWzCgrbWFi46bWE8D6CIVl6mNxpqMnbaPhdLOItmFJbUhWenCmhI4a982BQyREIXdRtMpgEhJdNIXibZytlF1YYmISCAKICIiEogCiIiIBKIAIiIigSiAiIhIIEUFkEQOn/vM7Dkzu8zMhqU895fwiyciIlFV7BXIVcQTvX0W+ABwn5ntmHguFmbBREQk2ooNILu6+5XuvszdTwXuJB5EhpNlXQMREalPxd5IODT1gbtfZGa9wH3AjplfIvWqARaAEqkZ1WiPxQaQp83sCHe/N7nB3S81s37g0nCLJlHWIAtAidSEarXHvF1YZjYl5eHngT+l7+PulzF4pTWpcw2yAJRITSi2PS7q7Gbq/KWMnXcXU+cvDbQIHBQ2BnJ/Yn1m3H2zu2/OtJO7ByuB1KRMa3fk2i4i5VPM2jxhriRaSAD5JXC3mR2f/oSZTTOzB4p+V6l5zWZFbReRcKVeRTRlaXeZ1gMJs/cgbwBx968BFwE3J5fdNLOJZvZb4t1ZOxX9rlLz+jzzpLts20UkPOlXEZnaXba1ecJcSbSgabzufj5wBvATM/tvoBPYGzgNmFj0u0rN0xWISPVkW+u82SzvQlFhriRa0CwsM/sA8FGgD/gE8BBwiLtvLfodpS7oCkSkerJdLfS78/z8mTlfO3f6uEEztiD4SqKFzMI6F3gO+Afg34hfdbQDlxX9blI32rL8Wsm2XUTCU8pVRJjLRxdyBfJd4FrgB+7+CoCZvQj8xsx2BU5x996i31lqWpi/YkSkOKW2v7BWKSxkDGQvdz8zGTwA3H0pcChwMPD7kkuRxsxmmNkaM1trZvMyPD/ezB42s81m9s9hv7/kF+avmEaiui1hiEr7y3sF4u7PZtn+mJlNA5aEWSAzawauBI4AuoBHzWyxu69K2e114OvArDDfW4qj9dCLo7otYYpC+ytpPRB3Xwv8n5DKknQAsNbdn3P3LcDNwLFp7/uquz8KqOtMaonqttSVkheUSu3aCkkb8FLK467EtqKZ2elm1mFmHRs2bAilcCIlUN2WuhLFFQkz3UgQaG6ouy9w93Z3bx8xYkSJxRIpmeq21JVis/FWQheDEzOOAtZXqSySh1K6F0V1W8qq0u0xigHkUWBPMxsLdBPPAPyF6hZJMlFK96KpbkvZVKM9Ri6AuPtWMzuL+OyuZuA6d1+ZzMPl7teY2d8BHcD7gX4z+yYwwd3/Vq1yN6JcSdkUQLalui3lVGx7DONqJXIBBMDd7wbuTtt2Tcrf/0v88l+qSCndi6e6LUHlO+EHSele6tVKFAfRpQbkWjtACRVFwlXIGh7FpDcJK6V7JK9AKkUDwMHlqmhKqCgSrmwn/G/d+jj/eMtyRra2cOj4Efx6WXdB6U3CSunesFcgYa7K1YhyVTQlVBQJV7b21uc+cP769bJujt+vraD0JmGldG/YKxANAJdmZGtLxrEOAyVUFAlZtvaWqqe3j/tXb+DBeYflPV5YyVAb9gokzFW5GtHc6eNoiTUP2mbAyQfurgAsErJM7S2TQs9fYSVjbNgrkGwRPciqXI0oWdE0hiRSfuntrcks41hjMeevMJIxNmwA0XoWpYtCNlCRRpHa3tKn4UJ1zl8NG0D0C1pEalVUzl8NG0BAv6DrnaZpSz0r5vxVrrbQ0AFE6pfydInElbMtNOwsLKlvYd1pK1LrytkWFECqZFFnN1PnL2XsvLuYOn+pbmAMmaZpi8Rlq/PdG3tKPveoC6uMsvU71kP3yqLObs5bvJKNPfGVV3faPsa5n/lYQeWvxNhE0GnaGjeRWpCvniafz3fzYannHl2BlEmuVCm13r2yqLObubc9PhA8AN7Y1Mvc2x/P+2umUilkDh2feZW+bNsrWTaRUuSrp6nPF6KUc48CSJnkChK13r1y3uKV9PZvexNTb5/nrYiVCp73r868Tni27ZUsm0gp8tXTTM/nE/Tcoy6sMskVJGrxLvhCL4nzVcRKBc8g71PrgV0aQ756GqS+Dm+JBSqLrkCKUMzAd7b/IcNbYhnz2kT5LvhiLonzBcGwsoAGLUeu96lU2USCWtTZTVOW9XaazBg7766sz+fy1uatgbpqFUAKVGz/+JatmS8hzcJLZFZM2UuZ8VXoJXGT5c/EO3f6OGLNgyt4rNlCD55BgnStBXZpLMlzULb1dpKp3YOsx9PX7/zgtyuLfl1Dd2EVM+OmmPTvizq72dTbn/E4GzfFB54rdRd8GDO+Cr0kLrjapu+Y44WFzjZJfz5Iqof012w/pJl3tvTxzVuW881blrN9rIl/PW6SZmUVQbPaipftOwsytlGMNzb1sqizu6j/Pw0bQBZ1djP39sfp7Yufvbo39jD39seBzCfWXHOp0+UadE3tDqlE4yok8OUrR+v2Md7Y1Es+7mwTnFKPPbwlxps9vdvEi95+57zFK7cpA5A1+AH84LcrB5Wre2MP37xlOR0vvM4FsyaWFKQdeGfL4O9tU28//3Tr8kGfT7Krh+nqmZSz3eb6zioxFlfsekjmDbL8aHt7u3d0dAw83vf8e7KeFNvS7tlIP1Fl0toSw4yC9jvvmI8BZMymGaQrK1sZW1tig6bapjLg+fkzs2b1TJbj5P94mAeffb2o8iTtuuMQXn1rS+FXJilaYs0M3a4pY/lbW2Js3tqf89eYEQ8CzWac9PHRXDBrYt73zPRdZNLW2rLNoj1mtszd2/O+SRmk1+1sKn01MHX+0ow/sDJ9f+WQ6fNCaQkI87WXUmX7zlpbYuwwdLuCp+aWYt38mYMe56rbDRtAxsy7K+f+yRNQpbW2xFh+7pEAfG/RCm565MWBcuwwpJkLZ8dPhKm/6t/avJW+DNNqc0k24sk/uCdrkGkEqTdA5vpRka6YRlZu2QLI9xat4Kb/eZFsTTzWZFzy2X3KFkTGzrsrYxtK/njJJKwgl+lEH2s2+vqd1KYSazYuOaHw7yBfUEy/wRbin/f/fPgDrFz/1jZtLf0HTr7zUrk1m/HsRUcN2lZzAcTMZgBXAM3Ate4+P+15Szx/FLAJONXdH8t1zEwB5OkhXyBW/IQFEXqG7sL2//LMwONCA0gl6jZA//eHE2AyjjQ4d2g6/81B23LV7ciNgZhZM3AlcATQBTxqZovdfVXKbp8G9kz8+zhwdeK/BUsGDzUyCaJl86tw6Xj459UFv6ZSdTsZPFS3JYj+84bTdN6befeDaE7jPQBY6+7PufsW4Gbg2LR9jgVu8LhHgFYz262YN1HwkFKYAW+/XOzLKlK3FTwkKDOwIjqlohhA2oCXUh53JbYVuw9mdrqZdZhZx4YN2VNYiFSI6rbUlSgGkEy/ndJjYiH74O4L3L3d3dtHjMieRE+kQlS3pa5EMYB0AaNTHo8C1gfYJ6deJ+vsFJF83IH3FdWzBBWq2666LRUSxQDyKLCnmY01syHA54HFafssBv6vxR0IvOnuRXVInz/lgYEgUq7GVkuN2Mv8XRRblrCPl/6v1GOt99aiBtATKlK3F89elfEz6190/0WFO9zQ96mC94/cLCx332pmZwFLiE91vM7dV5rZGYnnrwHuJj7NcS3xqY5fKvZ9Lpg1kTGP/BKAH2x3HSc3L6WZ/oG+gnxjkG/7UMyMHXgXgM1sxxC2DnrdG97CTgyeM+4FHLsUm4kxlN5t3qeXJmLE06u8w1CafSvDLD5Hvh/jxr7DOXfrafxuyFzGs22urF43Yimja4V8jj6aBr7T1H1zvbaPJh7q34sDm1YTI/MNff0YW7xpoPyZ3rOPJproZ71/kB9t/RyL+6cN7HNM0wOcu90NfMDezvMJBnPgF32f4tytpwGwrqhXV65uz9q3jUWs2mbBrwm77cgjz70RKFfS9rEmenr7B/WltcSaGBZrZuOm3m0WTEtmbm42o8990M256Uq94S/T++20fQx3eLOnN+frC71RuFAtsSbeTfue8glaH8vhhkT9nlPg/pG8D6QcMs2Vz3ZTULGaDPLdx1fIPtWQvGGyLUuK+UxOOXD3gRufir0RsdgbNHfaPsbMSbsNuqEyKlJvJozijYRSe9JvHq6GIc3G0xe+dzNhrrodxS6sismUfTWIXIGhrbWFdfNnctnnJhNrquzcylMO3J1182fy4xMn05qSXn6HIc3x1CvEc3P9+MTJPDjvMNoKSFve2hIblBbkvGM+VtDnaok18+MTJ3P5iZMHshA355lrasDMSbtx/+oNgRqUEf9FWI6vfaftg62fIJLLBbMmDmojba0tnHLg7hWtbz86YZ+C941cF1YlpWZfLUeOmdRU4Mn3Sk9zkEsp+W922j7GTY+8yP2rNzB3+riB9CjJFA/JMqQma5s7fVzOXFAtseaBPF5J2bLeZtqW3Dc10eLc2x7PuLohxK9Ufr2sO1AG0vR8S8WkKckn1myc+5mP5d9RJIBMSUAvmDUxY5qUXNpaW3hn89aypipq6C6sVNm6s3YY0symLX1F/wLO1edbSNdZaoK2UvPjpB4r13u3tbYwZueWgX5yI36y3JLIWJyaNyoshTSKZL92oTIlt8vW1ZavSy3WbJy4/2juX70hZ3+8urCkktKzXJuxzVgUxLvEbnzkxaKOnf7jq6ZSmVRLpl/fLbFmYs1NeJbB3GzyZRvNlZY52a2UWgmyjU/stH2Md3tzZ6WFwesl5wpc3Rt7Bj3vMBA8IJ5pOOx03Km/trIl3ysmeDSbZcyM+maWAOXEkwpmugoqR8AUCUMhSxUs6uzm18uKX2WwmLTxCiAJ6V0xrYlZHMVe/hWygl22NdGzBZ5swS3ZjZLsgsv1azp9LY2gsi2iFYZc38sb72zOukhXUq602tmODYDFuwvzzdgRqSU/+O3KQN2/2w8pfFy4oQfR083at40H5x3G5SdO5t3e/qKDR7Zfv+mKXTo11xK4yTK3tbbk7WYLazWzci1sk+t76ckRPApZFjjXhInePmeHodvx/PyZPDjvMAUPqXmLOrsDj/lt2lL4eUJXIBkEXTqy3z3vXPVkP/rx+7Xl7VdPle+StRKrlSWNLGC2VhC5lqHNNtGh0MWJksf+5i3LMz5fye9PpNxyrYqaTzHjvQogGeQ6mey0fYyNm7ZdlhWyn1gzLVP562Xdoa1itqizm6YiB5qDKqSLLlmmIAsDZQuU2brxCilL6rGzBaJyBUWRaijlB1Exs97VhZVBtpNJctA602k618ks17rkpUoGpyDBo5CKcsqBu2fsOiukTN0be3DeG39Z1Fn8gF5Srm68YhTbfQjxzzN1/lLGzruLqfOXlvQ5RMopWVdL+Sk5ZLvCw4KuQDLI9mvXPfM4Qr6xj2y/BsLoNsnW3dZsxvtbtsu57nsh95cUspZ4IWUKY/C9kJknhRwjWcZC02SkXz2GPRNNJAyZlvENYvPW3JNVUimAZJDtJPOPWfrP+9wHriaKmQEURrdJtiDU7865n/lYxkBYyD0hQEF3pmeS7ZjrN/Zs07V16PgRRY0FhaGYQFSuYCgStlxjt8XeS1UodWFlkZzdlDozJ9cJP1c3TZBuk0JlK9PI1pa83T65ZiYFLd+izu6sXWPDW2LbdG3d+MiLoXZ1ha2cV48iYcpWJw34t8/tU3DaplgRUUEBpAj5cmdlG9cIq/++0DKlp1BJD4SZygXv5aYqpXyXLFmTsf/ViC+XWcxNj+VQ7HhGrgAtEiX5fkwev19b3vxzAH1FXKioCyuLXLOIcuXOyvYrIIz++2zHTZYpU1nzzYYKu1zZPr8TT7VQyjFKFWQ8I4zZXyKVkKuuFpPlt5is4QogGeQ70eQaP6jGL9NsQaAaA8C57iZPlqGQY5RDkPGMYgfdRaolV2LTYlLEF3KVkqQAkkEhJ5pa+GVajQHgfN9Lvlki5fwOg45nlOvqUSRsmepqsdN6T/r46Pw7JSiAZFDIiaYWfplWYwC4kO+lkrOwUrvwst1sqfEMqWeFtvcmgy98fPeipu4rgGRQ6LTbqP8yLef04VxyfS+V/M7Su/AyBY+oXTWKhC1nIlHiE1wuP3FyoHapWVgZlHPabSXVy+cIKtdNlmHPhhOJqnyzR53gubN0BZJBlLqnguaUgmh9jmrIdZPl8ynrmYvUs0Jmj3YnbvIt9tygAJJFFLqnwphFFYXPUS3V6sITiZp8s0eBQDM01YUVYeVMwtgIGr0LTyRdru6sIOeWSAUQM/uAmd1rZs8k/rtTlv2uM7NXzezJSpexkpRGozTlzABQLNVtiYJkm8im2HNLpAIIMA+4z933BO5LPM5kITCjUoWqFqXRKF2uVC4VprotkTBr37asiVKLPbdELYAcC1yf+Pt6YFamndz9T8DrFSpT1agLpq6obktkhHVuidog+q7u/jKAu79sZrtUu0DV1OizqOqM6rZERljnlooHEDP7A/B3GZ76bhne63TgdIDdd9897MNXRJBZVKVM/ZXgVLelVoR1jqh4AHH3T2V7zsxeMbPdEr/QdgNeLfG9FgALANrb28u/YHgEaAW96lHdlloQ5jkiamMgi4E5ib/nAP9VxbLUJE39DV9Ia6Krbksk/OC3K0M7R0QtgMwHjjCzZ4AjEo8xs5FmdndyJzP7FfAwMM7Muszs76tS2gjS1N9wJX+thbBqouq2VN2izm7eyLIuT5BzRKQG0d39NeDwDNvXA0elPD6pkuWqJbr7OlxhpcRX3ZYoyHWVEeQcEbUrEClRpul5RvyXcwndLw1LV3RST3LV2yC3ByiA1Jn0dc4NBhaTKaH7pWHpZk6pJ2HXWwWQOpS8+7qttWWblcg0oF4c3cwp9WTu9HFkW7C2HgbRJUTqfildlPJpiZRq1r5tWZe3rflBdAmXBtTD0cgp8aX+tIV4XtAVSB1T94uIpAvzvKArkDqmXFoiki7M84ICSJ1T94uIpAvrvKAuLBERCUQBREREAlEAERGRQBRAREQkEAUQEREJRAFEREQCUQAREZFAFEBERCQQBRAREQlEd6JLXVvU2a1ULiIpwmwTCiBSt5LrmSeXpE0uqAUoiEhDCrtNqAtL6lau9cxFGlHYbUIBROqWFtQSGSzsNqEAInVL65mLDBZ2m1AAkbqlBbVEBgu7TUQqgJjZB8zsXjN7JvHfnTLsM9rM7jezp8xspZl9oxplleiL0nrmqtsSBWG3CXPPtsR65ZnZj4DX3X2+mc0DdnL376Ttsxuwm7s/ZmY7AsuAWe6+Ktex29vbvaOjo2xll8ZmZsvcvT3H86rbUpNy1e1IXYEAxwLXJ/6+HpiVvoO7v+zujyX+fgt4CtCcTIk61W2pO1ELILu6+8sQb0zALrl2NrMxwL7A/2R5/nQz6zCzjg0bNoRdVpFiqG5L3an4jYRm9gfg7zI89d0ij/M+4NfAN939b5n2cfcFwAKIX+YXWVSRoqhuS6OpeABx909le87MXjGz3dz95UR/8KtZ9osRb2A3ufsdZSqqSFFUt6XRRK0LazEwJ/H3HOC/0ncwMwN+Djzl7pdVsGwipVDdlroTtVlYOwO3ArsDLwKfdffXzWwkcK27H2Vm04A/AyuA/sRL/8Xd785z7A3AC2mbPwj8NczPUEEqe3VkK/se7j4i24uqULdLVcv/j/LRZytO1rodqQBSaWbWkWvqZZSp7NVRy2UvRj1/Tn228EStC0tERGqEAoiIiATS6AFkQbULUAKVvTpquezFqOfPqc8WkoYeAxERkeAa/QpEREQCUgAREZFAGjKAmNkMM1tjZmsTmVEjJVta71wpwc3snMTnWWNm06tX+oHyNJtZp5ndmXhcE2U3s1Yzu93MVie+/4NqpexhMrNLEt/BE2b2GzNrrXaZShX1dl+Kqi0F4O4N9Q9oBp4FPgQMAR4HJlS7XGll3A2Ykvh7R+BpYALwI2BeYvs84OLE3xMSn2MoMDbx+Zqr/Bn+CfglcGficU2UnXim3C8n/h4CtNZK2UP+Ho4Etkv8fXHyM9fqv1po9yV+voznjHK/byNegRwArHX359x9C3Az8VTbkeHZ03pnSwl+LHCzu2929+eBtcQ/Z1WY2ShgJnBtyubIl93M3g98kng6Edx9i7tvpAbKHjZ3v8fdtyYePgKMqmZ5QhD5dl+KHOeMsmrEANIGvJTyuIsIr7mQltY7W0rwqH2mHwPf5r10HFAbZf8QsAH4z0T327VmtgO1UfZyOg34XbULUaJG+X+VdymAMDViALEM2yI5l7mQtN7JXTNsq8pnMrOjgVfdfVmhL8mwrVr/P7YDpgBXu/u+wDvEu6yyiVLZi2ZmfzCzJzP8OzZln+8CW4GbqlfSUNT0/6tCFXHOCEXF07lHQBcwOuXxKGB9lcqSVZa03tlSgkfpM00FjjGzo4BhwPvN7EZqo+xdQJe7J3+53U48gNRC2YvmOdLPA5jZHOBo4HBPdK7XsJr+f1WIaiwF0IhXII8Ce5rZWDMbAnyeeKrtyMiR1jtbSvDFwOfNbKiZjQX2BP5SqfKmcvdz3H2Uu48h/t0udfdTqI2y/y/wkpmNS2w6HFhFDZQ9bGY2A/gOcIy7b6p2eUIQ+XZfiqotBVDt2QPV+AccRXyWwrPAd6tdngzlm0b88voJYHni31HAzsB9wDOJ/34g5TXfTXyeNcCnq/0ZEmU6hPdmYdVE2YHJQEfiu18E7FQrZQ/5e1hLfMwgWf+uqXaZQvhMkW73JX62jOeMcr+vUpmIiEggjdiFJSIiIVAAERGRQBRAREQkEAUQEREJRAFEREQCUQARkbpnZgvNbF0B+51qZp5IB5Lcts7MFqY8PsTMzjOzhj9/NvwXICIN4YfA7ICvnZ14fdIhwLno/NmQqUxEJMLMbKi7bw7zmO7+bAmv7QyzLPWk4SNoo0tciruZjTezJWb2jpm9aGZfSjz/xcTCQm8nFqz5cLXLLPUjpf7tnah/bwO3mtmRZna3mb1sZpsSSR6/ZWbNaa9fZ2Y3mtlXEgtFvWtmj5nZoWn7bdOFZWYfMrO7EsffYGZXEF/bJb2MA11YZnYe8asPgN5E2T2RzmaDmV2e4fXJbrHxJXxVkaQrEEm6DfgP4FLgTOA6M9uT+OX6PCAGXEF8kaiPV6mMUr/+i3gup4uJLwMwnnjamJ8C7wLtwHnACLbNkHwwsB/xtDKbiefw+p2Z7ePuazK9WSIf1r1AC/APxBNkfhU4Lk85ryWeiPHviacP6QNw981m9p/Al83sHHd/N+U1XwX+291X5zl2zVEAkaRL3P0GADPrAD5DvOKP9URa6EQm2ivMbA93f6F6RZU69BN3vyLl8R+TfyQSBf6Z+EqC/2xm/+Lug9aaAaa6+4uJ/e8DXgC+B3wxy/vNIb7+y0Hu/kjidb8DVuQqpLt3mVlX4uH/+HuLbgFcDXwL+Czwi8QxJwEHAiflOm6tUheWJA0sGOTubxD/RfaID15TIPkLKjUttkgYfpP6wMx2M7OfmdkLwBagF7iA+BLDu6S99pFk8ICBFfnuAg7K8X4HAS8lg0fidf3ArUE/gMdXpVxC/IdX0leJL1JWkfTqlaYAIklvpD3ekmUbxNf5EAnTy8k/EtNjFxNfi+QC4DBgf+DCxC7p9e+VDMd7hdwrDu6W43WluAqYmhjT2QE4BfhPjy+jW3fUhSUiUZCaFvzDxMc8vujuNyY3mtlnsrx21yzbunO838vAxwo8VjHuBtYRv/J4HNgRWFDiMSNLVyAiEjXbJ/7bm9yQWG3v5Cz7H2hmo1P23RGYCTyc4z0eBkab2YEpr2sCPldA+ZJTjFvSn0h0g/2M+NjLWcAfSplCHHUKICISNU8RHwS/0MxOSKzRfm+O/V8B7jGzE81sFnAPsAODb/5Ldz3wHHBHYprtUcQXEHt/AeVblfjvt8zs42bWnvb8z4l3s+0DXFPA8WqWAoiIREpivGAW8L/ADcCVwJ+A+Vle8t/AvwH/CtxC/OT9aXd/Os97HEF85b6riAeU54mPueRzZ+I1ZxK/knk07dgbEmV6mTpaNjcTrUgoIjUrcXPgA+5+SrXLkmRmOwEvAj929/9X7fKUkwbRRURCYGYjgHHAN4j37lxV3RKVn7qwRETCMZP4DY8HAHPc/eU8+9c8dWGJiEggugIREZFAFEBERCQQBRAREQlEAURERAJRABERkUD+P71D/DYoV9CPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "\n",
    "#ax1.scatter(m, x2, label='data')\n",
    "ax1.scatter(m_test,(y_test-y_fit)/y_test,label='prediction')\n",
    "ax1.scatter(m_train, (y_train-y_est)/y_train,label='estimated')\n",
    "ax1.legend()\n",
    "\n",
    "ax1.set_xlabel(xlabel='m',fontsize=16)\n",
    "ax1.set_ylabel(ylabel=r'$x_2$',fontsize=16)\n",
    "\n",
    "#ax2.scatter(rap, x2, label='data')\n",
    "ax2.scatter(rap_test, (y_test-y_fit)/y_test,label='prediction')\n",
    "ax2.scatter(rap_train, (y_train-y_est)/y_train,label='estimated')\n",
    "\n",
    "ax2.set_xlabel(xlabel='rapidity',fontsize=16)\n",
    "\n",
    "plt.savefig('Plots/x2_LO_Gaussian_relative.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "223a95b9",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}