fig7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PhysRevA.96.022128 - Figure 7\n",
    "\n",
    "Positron mean energy\n",
    "\n",
    "Dots - data from paper (extracted with WebPlotDigitizer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import API\n",
    "import pwpp\n",
    "from pwpp import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "xdim = 10\n",
    "# tau\n",
    "E0 = 1; #[GeV]\n",
    "I21 = 5; #[]\n",
    "lbd = 1; #[\\mu m]\n",
    "taulst = np.logspace(log10(5),log10(60),xdim);\n",
    "def Ep_tau(tau):\n",
    "    return gp(E0/m, 0.855*sqrt(I21*1e21/1e18)*lbd, 2*pi*c/(lbd*1e-6)*hbar/e*1e-9, tau*c/(lbd*1e-6)*1e-15)\n",
    "Ep_tau = np.vectorize(Ep_tau)\n",
    "Eplst_tau = Ep_tau(taulst);\n",
    "\n",
    "# lbd\n",
    "E0 = 1; #[GeV]\n",
    "I21 = 5; #[]\n",
    "tau = 30; #[fs]\n",
    "lbdlst = np.logspace(log10(0.1),log10(10),xdim);\n",
    "def Ep_lbd(lbd):\n",
    "    return gp(E0/m, 0.855*sqrt(I21*1e21/1e18)*lbd, 2*pi*c/(lbd*1e-6)*hbar/e*1e-9, tau*c/(lbd*1e-6)*1e-15)\n",
    "Ep_lbd = np.vectorize(Ep_lbd)\n",
    "Eplst_lbd = Ep_lbd(lbdlst);\n",
    "\n",
    "# I21\n",
    "E0 = 1; #[GeV]\n",
    "lbd = 1.0; #[\\mu m]\n",
    "tau = 30; #[fs]\n",
    "I21lst = np.logspace(log10(1),log10(30),xdim);\n",
    "def Ep_I21(I21):\n",
    "    return gp(E0/m, 0.855*sqrt(I21*1e21/1e18)*lbd, 2*pi*c/(lbd*1e-6)*hbar/e*1e-9, tau*c/(lbd*1e-6)*1e-15)\n",
    "Ep_I21 = np.vectorize(Ep_I21)\n",
    "Eplst_I21 = Ep_I21(I21lst);\n",
    "\n",
    "# g0\n",
    "lbd = 1.0; #[\\mu m]\n",
    "I21 = 5; #[]\n",
    "tau = 30; #[fs]\n",
    "E0lst = np.logspace(log10(0.4),log10(4),xdim);\n",
    "def Ep_E0(E0):\n",
    "    return gp(E0/m, 0.855*sqrt(I21*1e21/1e18)*lbd, 2*pi*c/(lbd*1e-6)*hbar/e*1e-9, tau*c/(lbd*1e-6)*1e-15)\n",
    "Ep_E0 = np.vectorize(Ep_E0)\n",
    "Eplst_E0 = Ep_E0(E0lst);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAGoCAYAAAATsnHAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABNHUlEQVR4nO3deXxU5dn/8c+VQAJFfuKCVQFRK264AAYqKMqmYBVwq4qitVYQLa1W2z7i0iruWtvap7jgBigiuDYoj6jIIhCUiFYFqyKigAuouFQlLLl+f5xJHIaEzCQzc85kvu/Xa17JnLnnnCszc+Wac8597tvcHRERkagpCDsAERGRmqhAiYhIJKlAiYhIJKlAiYhIJKlAiYhIJKlAiYhIJIVaoMxsgJm9bWZLzezSWtqcYmZLzGyxmT2U7RhFwqQckXxmYV0HZWaFwDvAUcBKYCEwxN2XxLXpAEwB+rj7WjPbyd1XhxKwSJYpRyTfhbkH1Q1Y6u7L3H098DAwOKHNMGCMu68FUOJJnlGOSF5rEuK22wAr4u6vBH6a0GZvADObBxQCV7n7M4krMrPhwHCAFi1aHLLvvvtmJGDJH6+88spn7t465DDSliOxNsoTSZts5EiYBSoZTYAOQC+gLTDHzA509y/jG7n7WGAsQElJiZeXl2c5TGlszOyDsGNIUlI5AsoTSa9s5EiYh/hWAe3i7reNLYu3Eih19w3u/j7B8fgOWYpPJGzKEclrYRaohUAHM9vDzIqA04DShDZPEnwzxMx2JDicsSyLMYqESTkieS20AuXuG4GRwHTgLWCKuy82s9FmNijWbDrwuZktAWYCf3D3z8OJWCS7lCOS70LrZp4pOrYu6WBmr7h7SdhxZIryRBoqGzmikSRERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSVKBERCSSQi1QZjbAzN42s6VmdmkNj59tZmvM7LXY7dww4hQJi3JE8lmTsDZsZoXAGOAoYCWw0MxK3X1JQtPJ7j4y6wGKhEw5IvkuzD2obsBSd1/m7uuBh4HBIcYjEjXKEclrYRaoNsCKuPsrY8sSnWRmr5vZo2bWrqYVmdlwMys3s/I1a9ZkIlaRMKQtR0B5Irkn6p0kpgK7u/tBwHPA+JoauftYdy9x95LWrVtnNUCRkCWVI6A8kdwTZoFaBcR/22sbW1bN3T9394rY3XuAQ7IUm0gUKEckr4VZoBYCHcxsDzMrAk4DSuMbmNkucXcHAW9lMT6RsClHJK+F1ovP3Tea2UhgOlAI3Ofui81sNFDu7qXAb81sELAR+AI4O6x4RbJNOSL5ztw97BjSqqSkxMvLy8MOQ3Kcmb3i7iVhx5EpyhNpqGzkSNQ7SYiISJ5SgRIRkUhSgRIRkUhSgRIRkUhSgRIRkUhSgRIRkUhSgRIRkUhSgRIRkUhSgRIRkUhSgRIRkUhSgRIRkUhSgRIRkUhSgRIRkUiqc7oNM/uwHut14Fh3f7MezxXJKcoRkcxIZj6otsA0YE2S6ywAhgJF9Q1KJMcoR0QyINkJC0e7+8vJNDSzJsCZ9Q9JJCcpR0TSLJlzUFcCK5JdobtvjD1nVV1tzWyAmb1tZkvN7NKttDvJzNzMGu0EcpLTMpYjoDyR/FXnHpS7X5fqSpN5jpkVAmOAo4CVwEIzK3X3JQntWgIXAi+lGodINmQqR0B5Ivmtzj0oM7vPzI7IwLa7AUvdfZm7rwceBgbX0O4a4CZgXQZiEGmwDOYIKE8kjyVziO9UYKaZvW9mV5vZXmnadhs2PyyyMrasmpl1Adq5+9NbW5GZDTezcjMrX7Mm2fPUImmTqRwB5YnksWQK1I+Bc4HlwBXA22Y2z8yGmdm2mQrMzAqAvwKX1NXW3ce6e4m7l7Ru3TpTIYnUJpQcAeWJNG51Fih3/6+73+/uvYHdCU7ubgfcBXxsZg+b2TGxREnFKqBd3P22bH7SuCVwADDLzJYDhwKlOgEsUZPBHAHlieSxlBLG3Ve4+/Xuvj9BItwH9AGeAlaZ2V9SWN1CoIOZ7WFmRcBpQGnctr5y9x3dfXd33x1YAAxy9/JUYhbJpjTnCChPJI/Ve6gjd3/Z3UcSHA//G7AT8LsUnr8RGAlMB94Cprj7YjMbbWaD6huXSFQ0NEdi61CeSN5K9kLdLcROBJ9FcEX87sDXwJRU1uHu0wiuwI9f9qda2vaqT5wiYUlHjoDyRPJXSgXKzLYjOMRwFkH3VweeAy4DnnR3dXGVvKYcEUmfZAaLbQocR5BwxxCMH7YEuBR40N0/zmiEIhGnHBHJjGT2oD4FtgW+AMYC4939lYxGJZJblCMiGZBMgZoNjAeedvcNGY5HJBcpR0QyIJmx+E7IRiAiuUo5IpIZKXczN7POZva4mX1mZhtjw6xgZteb2YD0hyiSW5QjIumRUoEys8OBMmBf4KGE51cCI9IXmkjuUY6IpE+qe1A3Elww2BG4OOGxRUCXdAQlksOUIyJpkuqFul2AE93dzcwTHvsM0AiUku+UIyJpkuoe1DrgR7U8tgvwVcPCEcl5yhGRNEm1QM0FLorN8lml6lvir4AX0hKVSO5SjoikSaqH+K4E5gH/Bh4lSLxfmNlfgUOArukNTyTnKEdE0iTV6Tb+DRxBcOX85YARjLQMcKS7v53e8ERyi3JEJH1SHs3c3RcBfc2sGbA98KW7f5f2yERylHJEJD3qPd0GsB74BKqnnXZ3T+y1JJLPlCMiDVDnIT4z29nMnjazs+KWFQIbEm5fmtmPMxapSEQpR0QyI5lzUBcQXNvxSMJyA+4BRgPXAB+R4lXyZjbAzN42s6VmdmkNj48wszfM7DUzm2tm+6eyfpEsyViOgPJE8lcyh/gGAHe7+/cJyx24K3a8HTNbQzAfztXJbDj2DXMMcBSwElhoZqXuviSu2UPufmes/SDgr7F4RKIkIzkSe47yRPJWMntQ+wDza1huCfffibVNVjdgqbsvc/f1wMPA4PgG7v513N0W/HA9iUiUZCpHQHkieSyZPahmwH/jF7j7JjPbhWDolirrYm2T1QZYEXd/JfDTxEZm9muCMc2KgD41rcjMhgPDAVq1akVZWRndu3dPIRSRBslUjkCG8mS33XZLMQyR7EtmD2o1sGfiQnf/1N03xS3aA1iTrsDitjPG3X8C/A9wRS1txrp7ibuXfPnllxx++OGcfPLJ3H777cyePZvPPvuspqeJpEuoORLbVkp50rq1hgSU6EtmD2oucCYwoY52ZxFcQZ+sVUC7uPttY8tq8zBwRzIrrqys5KmnnuKxxx6rXta6dWs6duzI/vvvX/3zu+++49///je9evXSHpc0RKZyBDKYJyJRl0yB+gcw18z+Alzq7hvjHzSzJsDNQC+gZwrbXgh0MLM9CBLuNOD0hHV3cPd3Y3ePBd4lCc2bN+f5559nt912Y/HixSxZsoTFixezePFiHnzwQb7++uvN2hcUFHD88cfTp08fOnbsSKdOnWjVqlUKf4rkuUzlCGQwT0SizpK5btDMLiFIsDXAc8CHsYd2I+hdtCMwyt1vSWnjZj8D/g4UAve5+3VmNhood/dSM7sN6EdwDclaYKS7L97aOtu2beuPPPJIrXtE7s6qVav405/+xLhx46j6+4uLi6moqKhu16FDB0pKSjjkkEMoKSlhw4YNLFy4UHtbecLMXnH3khTaZyRHYutOe56UlJR4eXl5qqGIVEs1R+q1jWQvbDez3gTHt4/ghxO964A5wM3uHolRmpNNvLKyMvr27cv69espKiqq3uN68803WbRoEeXl5ZSXl7NixYrNnldYWMjIkSM5+eST6dy5My1atMjUnyIhqk/y5UqOgAqUNFykClT1E4LrMnaI3f084SRw6FJJvLKyMmbNmrXVvaLVq1fzxz/+kQkTJpD4WhUUFLDffvtRUlJSvbe1bt06FixYoD2tHNeQ5It6joAKlDRcJAtU1GUi8RL3tiZPnoyZVe9llZeX8+mnn272nMLCQq699lrOPfdcdtxxx7TGI5mXjeQLkwqUNFQkCpSZvQBc4O7/SWqFwaCYzwPnxZ24zZpMJd7W9raqzmuNGjWKiRMnbrGntd9++3HEEUfQs2dPjjjiCNq1a4dEWyrJl2s5AipQ0nDZKFDJ9OLrBbRMYZ1Wj+dEXvfu3Ws9ZGdmtG3blgsuuIDHHnuM9evX07RpU/7+97/zxRdf8OKLLzJp0iTuuusuANq3b88RRxzBrrvuyvr16zn55JPp0aNHNv8cSa9eKEdE0i7Z6TaeNLOKuptVa1zHDZPUvXt3ZsyYscWe1qhRo9i0aROvv/46L774InPmzOGpp55i7dq1APztb3+jV69eHH/88fTs2ZODDz6YwsLCrW1Kokc5IpJmyRziu7+e6/6zu39Yd7P0ypVDF9dffz1XXnkllZWVmBmtWrWqLljbb789/fr1o3///uy4444sXrxYnS6yLMVDfDmVI5A7eSLRFYlDfO7+y0wGkK969+5NcXFxdceLp59+mrZt2zJnzhyee+45pk+fzpQpU6rbN2nShFtuuYURI0bQrFmqw7lJJilHRDJDvfhCVFfHiwsvvJB//vOfm3W6aN68OUceeST9+/enf//+7LvvvpglDpotDaVefCJbF4k9KMmcujpeDBkyhHvuuae608W1117Lhx9+yPTp0/nd734HQLt27aqLVd++fdluu+2A5K7xEhGJMhWoCKut0wXA8uXLefbZZ6sPBd5zzz0UFBTw05/+lAMOOIAJEyawYcMGiouLmTFjhoqUiOSctB3iSxiwMjT5eOhi48aNvPTSS0yfPp1nnnmGhQsXVj9mZpx77rnccccd6hmYgkwcvohKjkB+5omkVzYO8SUzH1SyJqdxXZKCJk2acNhhhzF69GhefvllSktLadq0KRCcy7r77rvZeeedOeeccygtLeX77xNnJpcsUY6IpCCdh/h0pj4iBg4cyOzZs5k1axZdu3Zl7dq1PPnkkzz++OPcf//9NG/enP79+3P88cdz3HHHscMOO9S9UkkH5YhIChpUoMysat5oA5qaWbvY74R1fYcEEjtg/PznP2fDhg3Mnj2bf/3rXzz55JM8+eSTFBQU0LNnT44//ngGDx7MJ598os4VaaQcEam/Bp2DMrOZBFfEG1BCMLmaAe7ufdISYYp0bD057s6iRYuqi9Ubb7wBUN1lvaioiJkzZ+ZtkUrX8fUo5ggoT6ThIjFYbNIrMnvV3TunZWUNoMSrn2XLlvGb3/yGadOmVS/baaed+PWvf82pp57KPvvsE2J02ZehThKRyBFQnkjD5VonicZ1xW+e2XPPPbniiito3rw5BQUFNG3alF133ZWrrrqKfffdl4MPPpjrrruOd9+NRCe0XKUcEUlBOgvUO6k+wcwGmNnbZrbUzC6t4fGLzWyJmb1uZjPMrH16QpWaVF13de211zJ79mxeffVVVqxYwW233UbLli254oor2HvvvenSpQs33ngj7733HmVlZdxwww2UlZWFHX4uUI6IpCC0oY5is46+AxwFrCQ4Nj/E3ZfEtekNvOTu35nZ+UAvdz91a+vVoYvMWblyJY8++iiTJ09mwYIFQOM9ZxWFoY4ylSOgPJGGy7VDfKnqBix192Xuvh54GBgc38DdZ7r7d7G7C4C2WY5R4rRt25aLLrqIsrIyli9fzjHHHIO74+5UVFRwyimncOutt/Lhh+qclibKEclrdRYoy9xIpG2AFXH3V8aW1eZXwP/V9ICZDTezcjMrX7NmTRpDlNq0b9+eK6+8svqcVZMmTdhmm234/e9/T/v27enZsyd33nknn332WdihZlwu5AgoTyT3JLMH9bWZHZjxSLbCzIYSdNG9pabH3X2su5e4e0nr1q2zG1weiz9nNWfOHN566y2WLl3Ktddey+eff87555/PLrvswnHHHcekSZP49ttvww45UyKfI6A8kdyTzIW6LYDmVXfMrIDgWPjJ7v5+3PJmQJG7f53ktlcB7eLut40t24yZ9QMuB45091RmLJUsSLwg+Cc/+QmXX345l112Ga+//joPPfQQkyZN4umnn+ZHP/oRhx9+ODvvvDPnnnsuPXv2DDHytMq5HHn77bfp1atXkmGIhKM+56AM6Axsl7D8YOCLFNazEOhgZnuYWRFwGlC62YbMOgN3AYPcfXU9YpWQmBkHH3wwN910E8uXL2f27NkcffTRPPvss0yYMIEjjjiCE088kblz51JZWRl2uOmmHBFJg3RPt5H0sXh332hmI4HpQCFwn7svNrPRQLm7lxIcrtgGeCR2mP9Ddx+U5pglwwoKCjjiiCOYN28epaWl1dPcP/XUUzzxxBO0b9+eIUOGcPrpp3PggaEeKcuGSOTIPvvsw6xZs+oTvwjwQw/eTAp1Pih3nwZMS1j2p7jf+2U9KMmYXr16bTbN/dSpU/nkk0+YOHEit9xyCzfeeCMHHHAAp59+Ovvss0/1YajG0nW9PpQjks+SLVC6Al4arLYJGM844wzWrFnDI488wkMPPcRll11W/ZymTZsydepU+vfvH1bYyVKOiKRZsuegZphZmZndCfyaIBmbZi4saay6d+/OqFGjttgrat26NRdccAFz587l97//ffXhgw0bNnDssccyePBgpkyZEuW5rJQjImmWzB7UMIITvp2A0wmOdwPMM7NlwBvA68DGTAQo+efEE09kzJgxVFRU0LRpU0488URmz55NaWkpLVu25MQTT2To0KH07t07KrMEK0dEMqDOAuXu98bfN7MOBInYiSApuwMnVDVPb3iSj2o6FLhp0yZmzZrFxIkTeeyxxxg/fjy77LILp512GmeccQZdunRhwYIFocxlpRwRyYy0jMVnZjsBXYCD3f2mBq+wATTGWOP3/fff8/TTTzNx4kSefvppNmzYwG677cZHH31EZWUlxcXFzJgxo0FFKt3jjEUpR0B5Ig2XM2Pxuftqd38mCoknjV/z5s05+eSTeeKJJ/jkk0+46667KCgoYOPGjVRWVvL9999zww038Omnn4YdajXliEjqwhwsVqTBtt9+e4YPH85DDz1EcXExZoaZMXXqVHbddVeOPvpo7r//fr788suwQxWRFKlASaPQvXt3Zs6cyXXXXce8efN48803GTVqFO+99x7nnHMOP/7xjznhhBOYMmUK3333Xd0rFJHQhTYfVKbo2LrEc3cWLlzIpEmTmDx5Mh9//DHbbLMNgwcPplOnTqxbt46+fftucb4qCvNBZZLyRBoqGzmiAiV5Y9OmTcyZM4dJkybx8MMP88033wBQWFjIbbfdxvnnn09BQXBQQQVKZOtyppOESC4oLCykd+/ejB07lj/84Q/VxWjTpk2MHDmSdu3acfHFF6N/3CLRoAIlealfv34UFxdTWFhIs2bNuPrqqykpKeGf//wnXbt2DTs8ESHkwWJFwlLbuIBr167l8ccf59xzzw05QhHROSiRGugclMjW6RyUiIjkLRUoERGJpFALlJkNMLO3zWypmV1aw+NHmNkiM9toZieHEaNImJQjks9CK1BmVgiMAY4B9geGmNn+Cc0+BM4GHspudCLhU45IvguzF183YKm7LwMws4eBwcCSqgbuvjz2WGUYAYqETDkieS3MQ3xtgBVx91fGlqXMzIabWbmZla9ZsyYtwYlEQNpyBJQnknsaRScJdx/r7iXuXtK6deuwwxGJJOWJ5JowC9QqoF3c/baxZSISUI5IXguzQC0EOpjZHmZWBJwGlIYYj0jUKEckr4VWoNx9IzASmA68BUxx98VmNtrMBgGYWVczWwn8HLjLzBaHFa9ItilHJN+FOhafu08DpiUs+1Pc7wsJDmuI5CXliOSzRtFJQkREGh8VKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiaRQC5SZDTCzt81sqZldWsPjxWY2Ofb4S2a2ewhhioRKeSL5KrQCZWaFwBjgGGB/YIiZ7Z/Q7FfAWnffC/gbcFN2oxQJl/JE8lmYe1DdgKXuvszd1wMPA4MT2gwGxsd+fxToa2aWxRhFwqY8kbwVZoFqA6yIu78ytqzGNu6+EfgK2CEr0YlEg/JE8laTsANIBzMbDgyP3a0wszeTfOq2BMnckDZbe7y2x3YEPqszuvBkOr5kXvdMrCOVv6tDiuuOvAbkSSpSfV8amoO1PVbT8sRlW7ufqRyoz+c2qq/RPnXE1HDuHsoN6A5Mj7s/ChiV0GY60D32exOCD4zVsd7yFGIY29A2W3u8tsdSiTGk9yaj8SXzumdiHen+bGTpvQg9TzL5vjQ0B7eSY1ssT1y2tftReX3y8TWKv4V5iG8h0MHM9jCzIuA0oDShTSnwi9jvJwMveOyVSZOpaWiztceTWX8+SsfrkunXNirvXRTyJBWpvm4NzcHaHqtpeeKyuu5nQn22kW+vUTUL73MMZvYz4O9AIXCfu19nZqMJKnOpmTUDHgA6A18Ap7n7sjrWWe7uJRkOvUGiHmPU46uvXP278jVPwqTXp27ZeI1CPQfl7tOAaQnL/hT3+zrg5ymudmwaQsu0qMcY9fjqKyf/rjzOkzDp9albxl+jUPegREREaqOhjkREJJJUoEREJJJUoEREJJJUoCLAzPY0s3vN7NGwY5H8ps9izcyshZmNN7O7zeyMsOOJokx8dlSgMsTM7jOz1YlX69c0MrUH46z9KpxIpbFI5TNXm3z6LKb4ep0IPOruw4BBWQ82JGH/H1OBypxxwID4BUmOTC1SX+NI8jNnZgea2VMJt52yH3KoxpF8jrblhzERN2UxxrCNI8T/Y41iLL4ocvc5NczLUz0yNYCZVY1MvSTL4UkjlMpnzt1vAI7LcoiRkmKOriQoUq+RR1/sw/4/ljcvdETUODK1me1gZncCnc1sVDihSSOVzGjo1fRZrPX1ehw4yczuIDrDYIUla//HtAcVAe7+OTAi7DhE9Fmsmbt/C/wy7DiiLBOfHe1BZdcqoF3c/baxZSKZos9cavR61S1rr5EKVHYlMzK1SDrpM5cavV51y9prpAKVIWY2CSgD9jGzlWb2Kw9mOx1JMH/PW8AUd18cZpzSeOgzlxq9XnUL+zXSYLEiIhJJ2oMSEZFIUoESEZFIUoESEZFIUoESEZFIUoESEZFIUoESEZFIUoESEZFIUoESEZFIUoFqpMxsnJl57DYr4bFzzOxdM1tvZl8mub4r4ta3MhMxi9RHbMgdj02sV+OEeVvLh2xQ/tSPClTj9gnQHbigaoGZ7QqMBeYDfYB+Sa7r/ti6pqU5RpGG+gw4ClgLjN5Kuy3yIYuUP/Wg6TYatwp3X5CwrANQCIx397nJrsjdVwGrzGxNOgMUaSh3/wZ43szuBW4ys53cfXUNTWvKh6xQ/tSP9qBykJldFne4oKbb7bU8bxwwK3Z3RqztuNhje5vZE7HDJOvM7EMze8TM9CVGcsXbsZ8HJ/sEMysws2/M7M8Jy7eL5ccvYvevit3f18ymm9m3sRz5ZezxM83sP2b2XzObaWY/Sdtflcf0zyc3PQy8EPv9FOB3QG9gXWzZcuDGGp53DfAK8A/g18AioOob3dMEh0jOJzhk0gb4GfoSIznAzAy4OHb3IOC5JJ+6N7ANQS7E6xz7+WrC8keAu4G/EBwqvM/MOgC9gEuBpsBtwEPAT5P/C6QmKlA5yN2XAcsAzOxsYLm7z4pvE+TrFs97z8zeit1dUnW4w8x2BPYCBrt7/LwuD6U9eJHMuIDgHM+3pLAHBXSJ/aypQFUASxKW3+LuEwDMrBwYCJwH7OHuX8eW7wLcZmbt3f2DlP4K2Yy+Hee+g4DXG7iOzwkK3o1mNiz2jVAkJ5jZ7gRHDK4DXiLIiWR1AVbHzhElLl8cm/so3v9V/eLua4HVwIKq4hTzn9jP+FlnpR5UoHJY7LDGAcC/G7IeDyYFOwooB24A3jGzZWZ2fsOjFMm4e4ClwPXAa8B+ZtY0yecewpZ7TxDsQSUe3oPgMHi89bUsA2iWZAxSCxWo3NYeaEnD96Bw92XufhbQmiA5XwBuN7NjGrpukUwxs+HAEcDZ7r6BoEAVAfsmuYpOJBQiM9sJ2CdxuWSfClRu2zX2c3m6VuiB1/jhhPMB6Vq3SDqZWVvgFuBad686ivBa7Gedh/liPe1aAZsSHvoNwf/G15BQqZNEbvs29vNkM2tS32s8zOwggp5HkwkOlRQCZwMb+aG3oEjU3A28R3Bor8pbBJ0bDgYm1vH8Q2I/zzWzFQTnk44GfhFbXmJmi9z9+/SFLKnQHlRuex24ExgOPNiA9XwCfEiw11QKTCLYOzvO3V9paJAi6RbrvdqX4NBedUeG2O9vklxHiS7AFwTdw68GHiA4ZP5z4GvgVBWncFlwflwam9gFuL0Iuo+7uycexkh1fUawZ3Uv0Nfd2zY0RpFsqSkfzOxZgv+BR2Vh+8qfetAeVOPWHtgAzEjDui6PreusNKxLJAyJ+dCF4ML1bFD+1IP2oBqp2LUhO8bufuPub2+leTLr24VgdAmA9e7e4J6DItmSmA8Eo64sB05x90eysH3lTz2oQImISCTpEJ+IiESSCpSIiESSCpSIiESSCpSIiESSCpSIiESSCpSIiESSCpSIiESSCpSIiESSCpSIiESSCpSIiESSCpSIiESSCpSIiERSqAXKzAaY2dtmttTMLq2lzSlmtsTMFpvZQ9mOUSRMZnafma02szdredzM7B+xHHrdzLpkO0aRTAmtQJlZITAGOAbYHxhiZvsntOkAjAIOc/eOwEXZjlMkZOOAAVt5/BigQ+w2HLgjCzGJZEWYe1DdgKXuvszd1wMPA4MT2gwDxrj7WgB3X53lGEVC5e5zCKYlr81gYIIHFgCtYnMPieS8JiFuuw2wIu7+SuCnCW32BjCzeQTTJV/l7s8krsjMhhN8e6RFixaH7LvvvhkJWPLHK6+88pm7tw47jiTUlEdtgI8TGypPJJ2ykSNhFqhkNCE4dNELaAvMMbMD3f3L+EbuPhYYC1BSUuLl5eVZDlMaGzP7IOwY0k15IumUjRwJ8xDfKqBd3P22sWXxVgKl7r7B3d8H3iEoWCISSCaPRHJSmAVqIdDBzPYwsyLgNKA0oc2TBHtPmNmOBIf8lmUxRpGoKwXOivXmOxT4yt23OLwnkotCO8Tn7hvNbCQwneD80n3uvtjMRgPl7l4ae+xoM1sCbAL+4O6fhxWzSLaZ2SSCL2k7mtlK4M9AUwB3vxOYBvwMWAp8B/wynEhF0i/Uc1DuPo0gweKX/Snudwcujt1E8o67D6njcQd+naVwRLJKI0mIiEgkqUCJiEgkqUCJiEgkqUCJiEgkqUCJiEgkqUCJiEgkqUCJiEgkqUCJiEgkqUCJiEgkqUCJiEgkqUCJiEgkqUCJiEgkqUCJiEgkqUCJiEgkqUCJiEgkqUCJiEgkhVqgzGyAmb1tZkvN7NIaHj/bzNaY2Wux27lhxCkiItkX2oy6ZlYIjAGOAlYCC82s1N2XJDSd7O4jsx6giIiEKsw9qG7AUndf5u7rgYeBwSHGIyIiERJmgWoDrIi7vzK2LNFJZva6mT1qZu2yE5qIiIQt6p0kpgK7u/tBwHPA+JoamdlwMys3s/I1a9ZkNUAREcmMMAvUKiB+j6htbFk1d//c3Stid+8BDqlpRe4+1t1L3L2kdevWGQlWJCxJdCbazcxmmtmrsaMNPwsjTpF0C7NALQQ6mNkeZlYEnAaUxjcws13i7g4C3spifCKhi+tMdAywPzDEzPZPaHYFMMXdOxPk0e3ZjVIkM0LrxefuG81sJDAdKATuc/fFZjYaKHf3UuC3ZjYI2Ah8AZwdVrwiIanuTARgZlWdieJ7uzrw/2K/bwt8lNUIRTIktAIF4O7TgGkJy/4U9/soYFS24xKJkJo6E/00oc1VwLNm9hugBdAvO6GJZFbUO0mISN2GAOPcvS3wM+ABM9sit9WZSHKNCpRItNXZmQj4FTAFwN3LgGbAjokrUmciyTUqUCLRVmdnIuBDoC+Ame1HUKC0iyQ5TwVKJMLcfSNQ1ZnoLYLeeovNbHSsAxHAJcAwM/s3MAk42909nIhF0ifUThIiUrckOhMtAQ7LdlwimaY9KBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiSQVKBERiaQ6L9Q1sw/rsV4HjnX3N+vxXJGcohwRyYxkRpJoS3AVe7JjexUAQ4Gi+gYlkmOUIyIZkOxQR6Pd/eVkGppZE+DM+ockkpOUIyJplsw5qCvZfMK0rYoNbnklW04JINJYKUdEMqDOAuXu17n7x6msNPacT+tqZ2YDzOxtM1tqZpdupd1JZuZmVpJKHCLZkMkcEclndRYoM7vPzI5I94bNrBAYAxwD7A8MMbP9a2jXErgQeCndMYikQ6ZyRCTfJXOI71Rgppm9b2ZXm9leadp2N2Cpuy9z9/XAw8DgGtpdA9wErEvTdkXSLVM5IpLXkilQPwbOBZYDVwBvm9k8MxtmZts2YNtt2Py4/crYsmpm1gVo5+5Pb21FZjbczMrNrHzNGk0kKlmXqRwRyWvJnIP6r7vf7+69gd0JTu5uB9wFfGxmD5vZMWaW1ot+Y+v7K8FsoXXFONbdS9y9pHXr1ukMQ6ROYeWISGOXUsK4+wp3v97d9wcOBe4D+gBPAavM7C8prG4V0C7ufls279XUEjgAmGVmy2PbK1VHCYmyNOeISF6r9zc6d3/Z3UcSHJb7G7AT8LsUVrEQ6GBme5hZEXAaUBq3/q/cfUd3393ddwcWAIPcvby+MYtkUxpyRCSvJXuh7hZiJ4LPIrgifnfga2BKss93941mNhKYDhQC97n7YjMbDZS7e+nW1yASbQ3NEZF8l1KBMrPtCPZ0ziLohefAc8BlwJPunlJPO3efRjBETPyyP9XStlcq6xYJQ7pzRCSfJTNYbFPgOIKEO4Zg/LAlwKXAg6leoCjS2ChHRDIjmT2oT4FtgS+AscB4d38lo1GJ5JaM5YiZDQBuIzgMfo+731hDm1OAqwj21v7t7qenY9siYUumQM0GxgNPu/uGDMcjkosykiNxo60cRXCd4EIzK3X3JXFtOgCjgMPcfa2Z7ZSu7YuErc4C5e4nZCMQkVyVwRypHm0FwMyqRltZEtdmGDDG3dfGYlmdoVhEsi7lbuZm1tnMHjezz8xsY2y0B8zs+tjhCJG8lsYcqXO0FWBvYO/YyBULtrZ+jbgiuSalAmVmhwNlwL7AQwnPrwRGpC80kdwTQo40AToAvYAhwN1m1qqmhhpxRXJNqntQNxJct9QRuDjhsUVAl3QEJZLD0pkjdY22AsFeVam7b3D394F3CAqWSM5LtUB1Ae5wdyfoMRTvM0BfyyTfpTNHtjraSsyTBHtPmNmOBIf8lqUetkjyysrKAHbO9HZSLVDrgB/V8tguwFcNC0ck56UtR2Iz71aNtvIWMKVqtBUzGxRrNh343MyWADOBP7j75/WOXqQOZWVl9OnTB7Y8H5p2qQ51NBe4yMz+Fbes6lvir4AX0hKVSO5Ka47UNdpKbE/tYrY8nCiSEdOmTWPduuwMiJJqgboSmAf8G3iUIPF+YWZ/BQ4BuqY3PJGcoxyRRmvt2rVMmVI9nGTiIey0S3W6jX8DRxBcOX85YASHIACOdPe30xueSG5Rjkhj9dVXX9G/f3+WL1/OrbfeCvBRpreZ8mjm7r4I6GtmzYDtgS/d/bu0RyaSo5Qj0th88803HHPMMbz66qs8/vjjDBw4kEsuueSTTG+33tNtAOuBT6B69luPHQ8XkYByRHLet99+y7HHHsvLL7/MlClTGDhwYNa2XechPjPb2cyeNrOz4pYVAhsSbl+a2Y8zFqlIRClHpLH67rvvGDhwIPPmzeOhhx7ixBNPzOr2kzkHdQHBtR2PJCw34B5gNHANwfFIjSQh+Ug5Io3OrFmzOPjgg5k1axYTJkzglFNOyXoMyRziGwDc7e7fJyx34K7Y8XbMbA3BfDhXJ7vxuqYSMLMRwK+BTcB/geHxIzmLRETGckQkDLNnz6Zv375UVlbStGlT9txzz1DiSGYPah9gfg3LLeH+O7G2SYmbSuAYYH9giJntn9DsIXc/0N07ATcDf012/SJZlJEcEQnD+vXrGTFiBJWVlQBUVlYya9asUGJJZg+qGcHeSzV332RmuxAM3VJlXaxtsuqcSsDdv45r34Is9LsXqYdM5YhIVm3YsIHTTjuN//znPzRt2pTKykqKioro1atXKPEkU6BWA3sSXCFfzd0/TWi3B5DKGP41TSXw08RGZvZrgqvki4A+Na3IzIYDwwF22223FEIQSYtM5YhI1mzcuJGhQ4fyxBNPcNttt9G1a1dmzZpFr1696N69eygxJVOg5gJnAhPqaHcWwRX0aeXuY4AxZnY6cAXwixrajCWYapuSkhLtZUm2hZojIg21adMmzj77bKZMmcJf/vIXfvvb3wKEVpiqJHMO6h9AHzP7i5ltUdDMrElsGJdeBB0ekpXMVALxHgaOT2H9ItmSqRwRybjKykrOPfdcJk6cyPXXX88ll1wSdkjVkpnyvczM/kjQSWGomT0HfBh7eDfgKGBHYJS7l6Ww7eqpBAgK02nA6fENzKyDu78bu3ss8C4iEZPBHBHJqMrKSs477zzGjRvH1VdfzahRo8IOaTNJjSTh7rea2SLgf4CT+OFE7zpgDnCzu6c6SvNGM6uaSqAQuK9qKgGg3N1LgZFm1o/gIse11HB4TyQKMpEjIplSVlbGzJkzKS8v54knnuDyyy/nyiuvDDusLViqI6/EuofvELv7ubtvSntUDdC2bVt/5JFHQj92KrnNzF5x95J6PjfSOQLBudry8vKww5AQVM3nVFFRgbtzxhln8MADD2CWeFXE1jUkR5KV6oSFuPsmd18du0Uu8VatWsWRRx7JxIkT0bBnEoao54jktxdeeKG6OJkZ+++/f8rFKVuSGYvvBTPbN9kVmllB7DkdGhZa/W3YsIGhQ4ey6667cuaZZzJ+/HhWrlwZVjjSyOVijkh++uabb3j++eeri1NxcTG9e/cOO6xaJXMOqhfQMoV1Wj2ek1bFxcVcfPHFvP/++0yfPp0HH3wQgH322YcDDjiA5s2bc+aZZ3L00UeHFaI0Lr3IsRyR/LN06VKOP/54/vOf/3DhhRey00470bt370ifDkl2uo0nzawihfWGdmytTZs2xJ+Dqqys5M033+T555/n0Ucf5bHHHgPgwQcfpGPHjgwePJh+/frRo0cPiouLwwpbcl/O5Ijkn2eeeYYhQ4ZQWFjI9OnT6du3b9ghJSWZAjW+nuv+rO4m6bfzzjtv9o2goKCAgw46iIMOOoiKigpeeuklKisrMTO++eYbbrrpJq6//nqaN29Oz5496dChA02aNOGUU06hR48eYfwJkntyKkckf7g7N998M6NGjeKggw7iiSeeYI899gg7rKQlcx3UL7MRSDb06tWL4uJi1q9fT1FREQ8//DAdO3Zkzpw5PP/885SWlvLss88C8I9//INTTjmFESNGcPjhh9OkSUPmdpTGrDHliDQe3377Leeccw5Tpkzh1FNP5d5776VFixZhh5Uad29Ut0MOOcS3Zv78+X799df7/Pnzt3js+uuv94KCAic4/OJNmjRxwLfbbjs/44wzfPLkyf7VV19tdR3SOBBcixf65zkIhQHA28BS4NKttDsp9tktqWuddeWJ5Kaq/02PPvqoH3TQQW5mftNNN3llZWXat5WNHMm73YLu3bvXelIwcQ9r6tSpfP3115SWlvLUU08xceJEmjRpQmVlJe5OUVERM2fOjPRJRsltcdPSHEUwoPJCMyv1hHnRzKwlcCHwUvajlCgoKyujb9++VFRUUFlZScuWLZk2bRoDBgwIO7R6S/k6qMase/fuzJgxg2uuuYYZM2bQt29fTjjhBO6//34++eQT5s6dS48ePaoLVEVFBSeddBJ//vOfeeWVV6q+xYqkU/W0NO6+nmBMysE1tLsGuIlg5ArJQzNnzmTdunXV8zgNHz48p4sTJN+Lr04J4+blrNr2sAoLCznssMO48cYbq7+lFBQUsNNOO3HttdcyevRo2rRpw6BBg9h33335+uuv6du3r/aupFo9c6TOaWnMrAvQzt2fNrM/bGX7mpamkfr+++958cUXq78kN2vWjJNOOinkqBounYf4JgNd0ri+SKray4qfJ+Wzzz5j2rRplJaWMm7cOL7/Ppj5+6qrruK2225jxIgRFBYWhhy5REDac8TMCghmmj67rrauaWkapRUrVnDCCSfwyiuvMGzYMNq3b0+fPn0ax5fjdJ3MAl7N9AmzZG5hn/y9+uqrN+toAfjOO+/sv/nNb/zFF1/0TZs2hRqfJIcMnACuT44A3YHpcfdHEYyKXnV/W4Lu6stjt3XAR9TRUSLsPJH0mD17trdu3dpbtmzppaWlWd12JnIk8dagc1Bmtlvs1h5oambtqpY1ZL257KijjqK4uJjCwkKaNWvGNddcQ48ePbj77rvp2bMn7du35+KLL2bBggW4O2VlZdxwww2UlWkWhsYoDTlSPS2NmRURTEtTWvWgu3/l7ju6++7uvjuwABjk7hoJthFzd8aMGUPfvn3Zfvvtefnllxk4cGDYYaVdQw/xjSfYSzCgfey+xZbVOD17Y1fTIUAIxsCaOnUqkydPZsyYMfztb39j55135rPPPqOyspKioiJeeOGFxrFbLvEalCOe3LQ00siVlZVV/0858MADufDCC7nvvvs47rjjePDBB9l2223DDjEz0rUrhg7xJW3t2rU+btw432effTY7FHjkkUf6q6++mpFrFiQ1ROQQX6ZuuZAnEpg/f743b97cCwsLvaioyHfeeWcH/Iorrgj1lEEmciTxls5OEjrpmqRWrVrxi1/8gr333rt6XhYzY+7cuXTu3Jn99tuPoUOHcsYZZ9C+ffuww5X0UY5IymbNmlV9bdOmTZvYtGkTL774IocffnjYoWVcOq+DeifVJ5jZADN728yWmtmlNTx+sZktMbPXzWxG7Dh+o9G9e3deeOEFrrvuOubOncvHH3/MHXfcwQ477MDll1/O7rvvTq9evbjnnnv48ssvww5XGi7lHJH85u6bXdvUpEkTHnnkkbwoTlCPGXXTtuHgCvl3iLtCHhjicVfIm1lv4CV3/87Mzgd6ufupW1tvY5kp9P3332fixIk88MADvPPOOxQXFzNw4EC6du1KRUUF/fr10/mqDMrGbKFhaix50pgtX76cESNGMH36dDp27Ejfvn057bTTIpP32ciRMAtUd+Aqd+8fuz8KwN1vqKV9Z+Cf7n7Y1tbb2BLP3SkvL+eBBx7ggQceqN6TKiws5I477uDcc8+N7GyYuUwFSsJQVlbGCy+8wOrVq7nnnnsoKCjghhtu4Pzzz4/ctZTZyJE6z0GZmXlmqlidV8gn+BXwfxmII9LMjK5du9K1a1d22mkn/vznP1cfix4+fDg333wzQ4cOZejQofzkJz8JO9y8lMEckTxSVlZG7969qagIphXr0aMHkyZNyutRP5I5B/W1mR2Y8Ui2wsyGAiXALbU8PtzMys2sfM2aNdkNLov69u272TVWl112GbvtthtXX301e+21F4cddhh33HEH06dP17VV2RV6jkhuW7duHZdddll1cTIzjj322LwuTpDcdVAtgOZVd2JDqywETnb39+OWNwOK3P3rJLe9CmgXd79tbNlmzKwfcDlwpLvXOGOp58kQLrVdY7VixQomTZrEAw88wAUXXFDdvqioiOeff56ePXuGFXK+yFSOSB6YPXs2w4YN491336WwsBB3p7i4mN69e4cdWvjq6ocOVALd4u4XxpZ1SWj3U2Bjsv3bCYrjMmAPoAj4N9AxoU1n4D2gQ7LrzefrOyorK33kyJFuZtXXVrVo0cIvuugif/XVV8MOL6eQwjUemcqRTN7yOU+iYu3atT5s2DAHfI899vDnnnsup+aaSyVH6ntL93QbSZ+td/eNQNUV8m8BUzx2hbyZDYo1uwXYBnjEzF4zM101vxVmxumnn06zZs0oKCigqKiIrl27cvvtt9O5c2cOPvhg/vrXv/Lpp58CaJilcKhHS55zdx577DH2228/7r33Xv7whz/w5ptvVvfMHTVqVGR66oWurgpGat8ON2W6otZ10zfDLWcN/vzzz33MmDHerVs3B7ywsNB79OjhRUVFXlBQ4M2bN8+Jb2zZROb2oELPEVeehGbWrFnes2dPB7xTp05eXl4edkj1lkqO1PeW7EgSjfa8TmOUOKfV9ttvzwUXXMAFF1zAW2+9xYQJE7j99ttZv349EJygnTBhAoceeqi6rNefckQ2Ez9+HsCVV17JjBkz2GWXXfjf//1fzjvvPJo2bRpukBGXbIGaYWaLCc4TLSFIRr2yOWi//fbjhhtu4Nhjj6Vv375s2LABd+fOO+9k5syZnHXWWZx55pm0a9eu7pVJPOWIVIuffh2gsrKS1q1bc+utt3L++efTvHnzOtYgkFyBGkbQWaETcDrBOSGAeWa2DHgDeB3YmIkAJTMOP/xwZs2axaxZs+jatSsffPAB48eP5/LLL+eKK66gT58+nHXWWZx44om88cYbW/QclM0oR2QzDz/8cPXEpQD9+/fn0UcfZZttttnKsyRRyiNJmFkHgkTsxA9JuXPsYXf3UC931hXyDfPee+/x4IMPMmHCBJYtW0azZs2q97LyaUqQhlwlH/UcAeVJprz99ttcddVVTJ48GXfHzCguLm6UeZMzQx2Z2U4EU1kf7O43NXiFDaDESw93Z/78+fzud79j4cKF1cuPOOII7rjjDvbff/8Qo8u8dCdflHIElCfptmzZMkaPHs0DDzxA8+bN+e1vf8uRRx7JokWLGu2Rh6wMB5bpXhjZvql3UnrNnz/fmzVr5gUFBV5QUOCFhYUO+CGHHOK33Xabf/rpp9XtcuX6jWSQhR5KYd6UJ+nx4Ycf+vDhw71JkyberFkzv+SSS6pzorHLRo6kcz4oaYSqpgSpOgf1k5/8hEmTJjFhwgQuvPBCLrnkEg499FBefvllNm7cSHFxMTNmzGiU3xhFqnrmHXDAATz77LOMHTsWgBEjRjBq1Ch23XXXkCNsZDJdAbN90zfD7HnjjTf8j3/8o/+///f/qkeuMDM/77zzQp3pMx3QHpQkmD9/vhcXF1eP1FJYWOjDhg3zDz74IOzQQpGNHEn3SBKSRw444ABuuukmnn76aYqKijAz3J277rqLPffck1GjRvH6668TfJZFclNFRQUPPfQQQ4YMoaKiAveg88NFF13E2LFj835A10zSIT5psMQu6x9//DEPPfQQt9xyCzfeeCMdO3ZkyJAhDBkyhE8//VRd1iUnfPjhh9x1113cc889rF69mrZt29K0aVM2bdpEcXExJ510UtghNnqhTViYKeqdFB1r1qzhkUceYdKkScydOxegeqSKqHe91YSF+amyspLnn3+e22+/nalTpwIwcOBALrjgAvr168dLL72kL1gxOdPNPEqUeNH0wQcfMGLECJ555pnqZXvttRejRo3ixBNPpFWrVpsNDRN28qtA5Yeqz1yXLl1YsmQJd9xxB++++y6tW7dm2LBhDB8+nPbt24cdZiSpm7lO/jYq8+fP9+bNm3tBQYE3adLE27Rp44AXFRX5kUceGanBa4lIJwlgAPA2sBS4tIbHLyYYWul1YAbQPpn1Kk+27PQA+GGHHeYTJ070devWhR1e5GUjR9RJQrKmasLFa6+9ljlz5rBixQpeeuklLrjgAhYtWsT69euprKxk3bp13H777ZsNFZOPzKwQGAMcA+wPDDGzxCukXwVK3P0g4FHg5uxGmXu+//57JkyYwKmnnrpZp4ff/OY3zJ07l9NPP53i4uKwwxSSm/JdJG3i57sxM7p168bf/vY3pk2btllPwAcffJDWrVtzyimnMHnyZL755hsg7+aw6gYsdfdl7r4eeBgYHN/A3We6+3exuwsIZqaWBBUVFfzrX//i9NNPp3Xr1vziF7+goKCApk2bUlBQQLNmzRgyZEjYYUoC9eKTSIjvCXj44Yezbt06HnvsMZ544gkeeeQRiouLKSkp4eWXX67uRZUHFwS3AVbE3V9JMKdUbX4F/F9tD5rZcGA4kBddo9evX8/zzz/P5MmTefLJJ/n666/ZYYcdOP300xkyZAi9evViwYIFkTnvKTXI9DHEbN90bL1x2bhxo8+ZM8cvuugib9WqVfW5AsAHDRrkH3300RbPScewS0TgHBRwMnBP3P0zgX/W0nYowR5UcTLrbqx5smHDBp8+fbqfc845vt122zngrVq18l/+8pf+zDPP+Pr168MOsdHIRo6EnYB1nQA+AlhEME3Bycmss7EmnrjPmzev+qR21YltM/MePXr4rbfe6u+//351R4zCwsIGdbaISIHqDkyPuz8KGFVDu37AW8BOya471/Mk/kvIxo0bfcaMGX7eeef5jjvu6IC3bNnSzzzzTJ86dapXVFSEHW6jlI0cCe0QX9wJ4KMIDl0sNLNSd18S1+xD4Gzg99mPUKKmR48ezJw5k1mzZnHkkUey7bbb8vjjj/PYY49xySWXcMkll7Drrruybt063J2KigpmzZqVy4duFgIdzGwPYBVwGsF8U9XMrDNwFzDA3VdnP8TsKysro0+fPqxfvx4zY9ttt+WLL76gRYsWDBw4kFNPPZUBAwbQrFmzsEOVBgrzHFT1CWAAM6s6AVxdoNx9eeyxyjAClOhJnM6+Y8eOXHnllbz33ns8/vjjjB8/no8++ggILrqcP38+jz32GP369WPbbbeN1LVWdXH3jWY2EpgOFAL3uftiMxtN8O21FLiFYILER2IXQX/o7oNCCzqDvvnmG2bOnMn111/PunXrqpfvuuuu3HnnnRx77LH86Ec/CjFCSbfQLtQ1s5MJvvWdG7t/JvBTdx9ZQ9txwFPu/mgt64o/+XvIBx98kLG4JfqmTp3K+PHj+eKLL1i0aBFfffUVTZo04cADD+SNN96gsrKyzskXdaFu+CorK3n99dd55plnmD59OvPmzWPDhg00b968unt41EckacyykSONohefu48FxkKQeCGHIyEbOHAgAwcOBGDDhg0sWLCAadOmMX78eDZuDGZdX7duHRdeeCGjRo2iX79+tGzZcrO9KwnHmjVreO6553jmmWd49tln+fTTTwE4+OCD+d3vfseAAQPo0aMHixYtypk9Yam/MAvUKqBd3P22sWUiadO0aVN69uxJz549GTRoEH369KGiogIzY/HixZx44ok0bdqUAw88kNdff7167wpoEXbsjVHiIdYNGzbw0ksvVe8lvfLKK7g7O+ywA0cffTT9+/fn6KOPZpdddtlsPYmHeqVxCrNA1XkCWCSdEidfPOSQQ5g/fz7/93//t8XeFdAy1GAbofjODYWFhfTo0YNXX32Vr7/+msLCQg499FBGjx5N//796dKlC4WFhWGHLCELdbBYM/sZ8Hd+OAF8XfwJYDPrCjwBbAesAz5x945bW2cuHFuX6CkrK6N3797VPcMqKyv/4+77hR1XpmQrTzZs2MCrr77KvHnzuPfee1m8eHH1Y61ateLnP/85/fv3p2/fvrRq1Srj8Uj6NPpzUO4+DZiWsOxPcb8vREO3SBZ07969ugt7r1696NGjx7dhx5SL1q5dS1lZGfPmzWPevHm8/PLL1WMq7rLLLhQWFlYfRn366afp0aNHyBFLlDWKThIi6aDzGsmpOo905JFH8uMf/7i6GM2bN696D6mwsJDOnTszfPhwDjvsMA477DB23XXXnOrmL+FTgRKRpHzxxRdMnDiRiy++uPp8XZVtt92W7t27c9ppp3HYYYfRrVs3WrTYsp+JvgRIKlSgRGQz7s7KlSt59dVXN7t9+OGHm7UzMwYNGsQ111xDx44dKSjQ5AiSXipQInnok08+oaysjG7duvHuu+9uUYw+//xzIChCe++9Nz169ODXv/41xcXFXHrppWzYsIGioiL+53/+hwMPPDDkv0YaKxUokTy0atUqDj/8cIqLi6s7MRQVFXHAAQdw/PHH07lzZzp37sxBBx3ENttss9lzu3XrpvNIkhUqUCJ5qrKykk6dOjF8+HA6d+7MfvvtV3WR8lbpPJJkiwqUSJ5q3rw5t956q4qNRJbOaorkoTZt2uTDjMSS41SgRPLQzjvvrOIkkacCJSIikaQCJSIikaQCJSIikaQCJSIikaQCJSIikaQCJSIikaQCJSIikRRqgTKzAWb2tpktNbNLa3i82Mwmxx5/ycx2DyFMkVApTyRfhVagzKwQGAMcA+wPDDGz/ROa/QpY6+57AX8DbspulCLhUp5IPgtzD6obsNTdl7n7euBhYHBCm8HA+NjvjwJ9zcyyGKNI2JQnkrfCHCy2DbAi7v5K4Ke1tXH3jWb2FbAD8Fl8IzMbDgyP3a0wszczEvHWbQt8FcI6knlOXW229nhNjyW7bEcS3qssSMf7ANAhDetIh0zlyTozW1yPeFJ5fZNt25DP59YeD+Pz1xDp+uxmazuZzxF3D+UGnAzcE3f/TOCfCW3eBNrG3X8P2LGO9ZaH9PeMDWMdyTynrjZbe7ymx1JYlvX3Ih3vQzrXk4Y4MpUn9fr7Unlesm0b8vnc2uNh/S9owHudlc9cLuVImIf4VgHt4u63jS2rsY2ZNSGo/J9nJbrUTQ1pHck8p642W3u8pseSXRaGdMURlb8nU3lS378vlecl27Yhn89UthN12fo7ciZHLFYJsy6WSO8AfQkSbCFwursvjmvza+BAdx9hZqcBJ7r7KXWst9zdSzIYuiRJ70XDZSpP8oE+f7kvtHNQHhwrHwlMBwqB+9x9sZmNJtg1LwXuBR4ws6XAF8BpSax6bMaCllTpvWigDOZJPtDnL8eFtgclIiKyNRpJQkREIkkFSkREIinM66BEJEeZ2Z7A5cC27n5yrm8nV5lZC+B2YD0wy90nRnk7qa5He1AiOcLMCs3sVTN7qgHruM/MVtd0MXtszL93zey/sTZvmVn3mtbjwcgWv2rAdmodWzCV7USBmbUzs5lmtsTMFpvZhQ1Y19Zet4Fm9r2ZVZjZp2Z2NXAi8Ki7DwMG1WM73yd+nmp5f2rcTj3e55TibfQFysxamNl4M7vbzM4IO558ZWZ7mtm9ZvZo2LHksAuBt2p6wMx2MrOWCcv2qqHpOGBADc+vGvPvdeAS4GOC3oBmZk/F3Z41sz3SsJ3NxhY0swMTtvOUme1U099aXxn8X7ARuMTd9wcOBX6dOF5imt6fvwNdgZbAJ8AJQHd+GGlkU4rbeQCoqGE7Y4ChBMNsVY392Dau2aYk461pDMm28fHWENdmcrJA1Va101GxJXmpvA+58E04ysysLXAscE8tTY4EnjSz4lj7YcD/JjZy9zkEXdETdQOWA10Iumc/DPzM3ee7+3FVN+Bu4J6q7QB71GM7W4wt6O5vxG8ndltdy99aLQr/C9z9Y3dfFPv9G4IvEW0SmqXj/Vnq7m/GXrfHCIazWs0PxaMg2e0AywiK6dqatgPsDjwSuw0G9gL+HLedZONNHENyZUK8W5WTBYoaqna6KrakZBzJvw/SMH8H/ghU1vSguz9CcK3U5NjewTnAz1NYfxuC8dnWAPcD5wFDY+cMatrO42b2PNARWJDidhLHFkz8Z17NzHYwszuBzmY2qoYm44jQ/wILpjrpDLwUvzxN788KCw7zvgZcRrAXdQtwkpndAUxNYTt/B26obTtx6zkB6A0cALxVtZ1k4427X/U+Px4fb10ryclOEu4+x7ac86a6YgOYWWLFfo3cLciRlOL7sCTL4TUaZnYcsNrdXzGzXrW1c/ebY6/3HcBP3P2/KW6qgGAP6jfADOD3wKXAlTVspwvQE2jn7mtS3E7S3P1zYMRWHo/M/wIz24Zgz+Yid/+6hlgb+v7g7puATmZ2HsH7soe7/zKV7VR9ngjGcNzatm42sxOAXkCbdLzP7v4t8Ms6G8Y0pn/YaanY0mA1vg9JfBOW2h0GDDKz5QSHSvqY2YOJjcysJ8E33Sf44XBMslYRjOG30t1fIvhHvoCgYKV7O3WNLdhQWf9fYGZNCYrTRHd/vJY26Xzdtic4lFjT+Z+6tnMYwSHOuQSvf/znKX5cx54Eh3DfSEO89XqfG1OBqpG7f+vuv3T38zPVBVPq5u6fu/sId/+Ju9d0aEFq4e6j3L2tu+9O0HHhBXcfGt/GzDoTnDsaTPANdQczuzaFzSwE2gOrzeyA2HY2krDnm6btdDCzPcysKLad0hSeX2+Z+l9gZkYw3NRb7v7XWtqk43Xbx8wOir1uQwiK1H9S3U7V5wk4nKB4x3+eqt6fY2PrWQv8op7xNvh9bkwFKhvfzKRueh/C8SPgFHd/z90rgbOADxIbmdkkoIzgn91KM/sVBGP+ASOBHwPlsZ9tgOsztJ3pBHsAU+IHvk2TbH8GDyOYBqWPmb0Wu/0soU06XrebgZeBbwjen3+5e+IlB6lupx3BBJeJ788dwDbAg+7+Zj3jbfD7nLNj8cWOOz/l7gfE7tc56rOkn94HCZs+g41XTu5B1VS1s/TNTOLofZCw6TPYuOXsHpSIiDRuObkHJSIijZ8KlIiIRJIKlIiIRJIKlIiIRJIKlIiIRJIKlIiIRJIKlIiIRJIKlIhEmpmda2Zey+3b2PQaW3v+uLj2s2p4vLuZPRy70He9mX1tZgvN7Boz26Ue8T5pZmvthzmzEh9vGYt7XOz+FXHxrUx1e42ZClQjFZfUh4Qdi0gDdQa+JZg9NvHWIzYFRV0+ibW/IH6hmV0CzANaA1cA/QgGNp0ODAfuq0e844FWwHG1PH4ywZh542P374/FNq0e22rUcnI+KElKZ4LRqLc654tIDugELHb3VCZGTFSR+Hwz600w4d9t7v67hPbTzOwGUptUsMrTwOcEA6w+VsPjZwEfArMA3H0VsMrMMjavVq7SHlTj1Ylg+P+KsAMRqa/YVBYHAa9nYPX/A3wW+7mF2PQc42qI6WAzK40dxvvezObF5k6qet56YBJwjJntkPDc3QimZX/ANc5cnVSgGqG4pH417FhEGqgDwbQPi82sScJtq+eetiY24vmRwHOxgpLs87oA8wnmYhoGnESwt/R8wuH08UBTgsOF8YYCBkyob+z5RAWqcapKahUoyXWdYj//BmxIuJU3YL07AM0IDrVtJrEQJjx8S+w5fdz9UXefBpwALCOYgh0Ady8nmOzxrITnnwkscPd3GhB73tA5qMapU+znayHGIJIOnQAnmP01cU9nbdUvZvYTgr2WnQg6VAyLFYmUmNnOwMcJy5q6+0Yza06w13U9UJlQvJ4HzkhY3XjgJjPb293fMbNuwL7A+anGla+0B9U4dYr9fM3MtjOzp8zsHTP7t5k9a2Z7VTU0sytjj1Wa2fGhRCtSu87AMnef7+7lCbf34trdCYx3972BPwITY4e6a/M5sA7YLWH5Z0DX2O3uhMe2BwoJ9pQS9+ZGAtuZWfz/1AeBqlltif2sACYn8XcLKlCNVWfgfXf/kuDb59/dfW93Pxh4Crgnru1zwABgTtajFKlbJ+CNrTUws9bAocA4AHd/juA8T62XWMQmNZwDHGVmRfHLqwog8FHC074kKDj/yw9FbLNbbJr1qnV9RJBfQ2PbOBWY6u5rkaSoQDVOnYgd3nP3L939+bjH5gO7V91x9wXuviybwYkkw8x+DOxM3T34dgM+dvcNccuWs+XeUaKbgR2Bm5KJx92/BV4EDgYW1bBHV9MhxfFAe+CG2LbG19BGaqFzUI1MXFLX1kHiIuBfWQtIpP46x35uNLNDa3j8dXf/rr4rd/cZZnYpcKOZHUTQs+59gs4TexP0wPuW4ChElYsJ9rymm9m9BOerdgS6AIXufmnCZp4EvgZ+B6wGnqlvvPlIBarxqUrqLQqUmf0Z2JPgCnmRqOsU+zk6dovnBOeEviPoVbdLrDND1V7U7tTQQy+Ru99sZvOACwk6P7QmODf1NsG5ojvjR6pw90Vm1hX4M/APYFtgDbCI4DxY4vq/N7MpwLnAQ7FDi5IkFajGp1Ps52vxC83sCuBnwNEN+dYpki3ufiNwYxLt1pjZy8DZwN1mdhTBOahX4tvFet154tBI7j6PYLijZON6iy2vb9pa+2EE10zVKNaZozAWs8TROahGxt1vdHdz9+pBJ2N7TgMJitNX4UUnkjEjgF+a2TsE1yqdkTBSQ3uC3nYzwgiuDpcTxJZ4zVTeM4220biZWUeC8fjeA/4bW7zR3Utij19FcPihNfANweGNQ+MLnEguM7PdCc4TAXzj7m+HGM4WYiOmt4ndXe/umRjWKSepQImISCTpEJ+IiESSCpSIiESSCpSIiESSCpSIiESSCpSIiESSCpSIiESSCpSIiESSCpSIiETS/wee0+NWixPLCQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot all\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2,figsize=(6,6))\n",
    "fntsz = 16\n",
    "\n",
    "f = open('data/fig7a/wpd.json')\n",
    "data = json.load(f)\n",
    "datatmp = data['datasetColl'][0]['data']\n",
    "datatmp2 = np.array([datatmp[i]['value'] for i in range(len(datatmp))])\n",
    "indxs = np.argsort(datatmp2[:,0])\n",
    "axes[0,0].plot(datatmp2[:,0][indxs],datatmp2[:,1][indxs],'k.')\n",
    "#\n",
    "#axes[0,1].plot(lbdlst, lbdlst*0+0.18571428571428572,'k.')\n",
    "#\n",
    "f = open('data/fig7c/wpd.json')\n",
    "data = json.load(f)\n",
    "datatmp = data['datasetColl'][0]['data']\n",
    "datatmp2 = np.array([datatmp[i]['value'] for i in range(len(datatmp))])\n",
    "indxs = np.argsort(datatmp2[:,0])\n",
    "axes[1,0].plot(datatmp2[:,0][indxs],datatmp2[:,1][indxs],'k.')\n",
    "#\n",
    "f = open('data/fig7d/wpd.json')\n",
    "data = json.load(f)\n",
    "datatmp = data['datasetColl'][0]['data']\n",
    "datatmp2 = np.array([datatmp[i]['value'] for i in range(len(datatmp))])\n",
    "indxs = np.argsort(datatmp2[:,0])\n",
    "axes[1,1].plot(datatmp2[:,0][indxs],datatmp2[:,1][indxs],'k.')\n",
    "#\n",
    "\n",
    "axes[0,0].plot(taulst,m*Eplst_tau,'k-')\n",
    "axes[0,0].set_xscale('log')\n",
    "axes[0,0].set_xlim([5,60])\n",
    "axes[0,0].set_ylim([0,0.6])\n",
    "axes[0,0].set_xlabel(r'$\\tau[\\mathrm{fs}]$', size=fntsz)\n",
    "axes[0,0].set_ylabel(r'$E_+[\\mathrm{GeV}]$', size=fntsz)\n",
    "axes[0,0].set_xticks([5,10,20,50])\n",
    "axes[0,0].set_yticks([0,0.1,0.2,0.3,0.4,0.5,0.6])\n",
    "\n",
    "axes[0,1].plot(lbdlst,m*Eplst_lbd,'k-')\n",
    "axes[0,1].set_xlim([0.1,10])\n",
    "axes[0,1].set_ylim([0,0.6])\n",
    "axes[0,1].set_xlabel(r'$\\lambda[\\mu \\mathrm{m}]$', size=fntsz)\n",
    "axes[0,1].set_ylabel(r'$E_+[\\mathrm{GeV}]$', size=fntsz)\n",
    "axes[0,1].set_xticks([0.1,0.5,1,5,10])\n",
    "axes[0,1].set_yticks([0,0.1,0.2,0.3,0.4,0.5,0.6])\n",
    "axes[0,1].set_xscale('log')\n",
    "\n",
    "axes[1,0].plot(I21lst,m*Eplst_I21,'k-')\n",
    "axes[1,0].set_xlim([1,30])\n",
    "axes[1,0].set_ylim([0,0.6])\n",
    "axes[1,0].set_xlabel(r'$I_{21}$', size=fntsz)\n",
    "axes[1,0].set_ylabel(r'$E_+[\\mathrm{GeV}]$', size=fntsz)\n",
    "axes[1,0].set_xticks([1,2,5,10,20])\n",
    "axes[1,0].set_yticks([0,0.1,0.2,0.3,0.4,0.5,0.6])\n",
    "axes[1,0].set_xscale('log')\n",
    "\n",
    "axes[1,1].plot(E0lst,m*Eplst_E0,'k-')\n",
    "axes[1,1].set_xlim([0.4,4])\n",
    "axes[1,1].set_ylim([0,1])\n",
    "axes[1,1].set_xlabel(r'$E_0[\\mathrm{GeV}]$', size=fntsz)\n",
    "axes[1,1].set_ylabel(r'$E_+[\\mathrm{GeV}]$', size=fntsz)\n",
    "axes[1,1].set_xticks([0.5,1,2])\n",
    "axes[1,1].set_yticks([0,0.2,0.4,0.6,0.8,1.0])\n",
    "axes[1,1].set_xscale('log')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}