{ "cells": [ { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy.integrate as integrate\n", "from scipy.integrate import quad\n", "from scipy.special import kv\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "#Plot particle density as a function of m for T=.1565 GeV using quantum statistics\n", "g=1 \n", "mu=0 #assuming no chemical potential\n", "mass=np.linspace(0.1,1,1000) # mass in GeV \n", "T=.1565\n", "i=[]\n", "for m in mass:\n", " def n(p):\n", " return (g*(4*np.pi)/((2*np.pi)**3))*(p**2)/(np.exp((np.sqrt(p**2+m**2) - mu)/T)-1)\n", "\n", " solution1, abserr = quad(n,0,20.)\n", " \n", " i=np.append(i, solution1)\n", "\n", "# Plot particle density using Boltzmann approximation\n", "k=[]\n", "\n", "for m in mass:\n", " def b(j):\n", " return (g/(2*np.pi**2)*(m**2)*T*np.exp(mu/T)*kv(2,m/T))\n", " k=np.append(k, b(m))\n", " \n", "#Plot particle density using Boltzmann & large argument approximation\n", "l=[]\n", "\n", "for m in mass:\n", " def b(j):\n", " return (g/(2*np.pi**2)*(m**2)*T*np.exp(mu/T))*np.sqrt(np.pi*T/(2*m))*np.exp(-m/T)\n", " l=np.append(l, b(m))\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0,0.5,'Density')" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEKCAYAAAAiizNaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmcTfX/wPHXe3YME2MnZmQfsxhjkC3Zd7KnMhVaUFL90KKiQgllK6WikqWQ+kplXxLGvjMYDGLsy5j98/vjXtMYs2Hu3Bnzfj4e9+Hecz7nc97n1sx7zjmf83mLMQallFIquzjYOwCllFJ5iyYepZRS2UoTj1JKqWyliUcppVS20sSjlFIqW2niUUopla008SillMpWmniUUkplK008SimlspWTvQPIiYoWLWq8vLzsHYZSSuUqW7ZsOWeMKZZRO008qfDy8iI0NNTeYSilVK4iIscy004vtSmllMpWmniUUkplK008SimlspXe41Eqk+Li4oiIiCA6OtreoShlV25ubpQtWxZnZ+e72l4Tj1KZFBERQcGCBfHy8kJE7B2OUnZhjOH8+fNERETg7e19V33Y9FKbiLQSkQMiEiYiw1JZ7yoic63rN4qIV7J1w63LD4hIyzvoc5KIXMvMPpS6E9HR0Xh6emrSUXmaiODp6XlPZ/42Szwi4ghMAVoD1YFeIlI9RbNngYvGmIrABGCsddvqQE/AB2gFTBURx4z6FJEg4IHM7EOpu6FJR6l7/zmw5RlPMBBmjDlijIkF5gAdU7TpCMy0vv8JaCqWI+oIzDHGxBhjjgJh1v7S7NOalD4G/i+T+8hy/16O5r1f9xCXkGiL7pVS6r5gy8RTBjiR7HOEdVmqbYwx8cBlwDOdbdPrcyCw2BhzOpP7yHLbT1zim/XhfLb8kC26V3lcREQEHTt2pFKlSlSoUIGBAwcSExOT5ftZtWoVf//9d5b3C3DmzBnatWuHv78/1atXp02bNgCEh4cze/bsDLdP2S40NJSXXnopy9qr7GHLxJPaWYXJZJs7Wi4ipYFuwKS7jAMR6S8ioSISGhkZmcomGWtVoyRda5Vlysowthy7cFd9KJUaYwyPPfYYnTp14tChQxw6dIgbN27wf/+X8gT/3tky8YwYMYLmzZuzY8cO9u7dy5gxY4C7TzxBQUF89tlnWdZeZQ9bJp4I4MFkn8sCp9JqIyJOgAdwIZ1t01peE6gIhIlIOJBfRMIy2MctjDHTjTFBxpigYsUynGooTe+0r07pB/LxytwdXIuJv+t+lEpuxYoVuLm58fTTTwPg6OjIhAkTmDVrFteuXePbb79l4MCBSe3btWvHqlWrAHjhhRcICgrCx8eHd955J6mNl5cX77zzDoGBgfj6+rJ//37Cw8P5/PPPmTBhAgEBAaxdu5aQkBB++umnpO3c3d0BS4Jq3Lgx3bt3p3LlygwbNowffviB4OBgfH19OXz48G3Hcfr0acqWLZv02c/PD4Bhw4axdu1aAgICmDBhAuHh4TRs2JDAwEACAwOTEmHKdqtWraJdu3YArF69moCAAAICAqhZsyZXr15Nt/21a9d4+umn8fX1xc/Pj59//pmEhARCQkKoUaMGvr6+TJgw4Z7/26nb2XI49Wagkoh4AyexDBZ4PEWbxUAfYAPQFVhhjDEishiYLSLjgdJAJWATlrOX2/o0xuwBSt7sVESuWQcTpLkPWxwwQEE3Zyb0CKDHFxsY9etexnb1s9WulB299+se9p66kqV9Vi9diHfa+6S6bs+ePdSqVeuWZYUKFcLLy4uwsLBUt7npgw8+oEiRIiQkJNC0aVN27tyZ9Au/aNGibN26lalTpzJu3Di++uornn/+edzd3XnttdcAmDFjRpp979ixg3379lGkSBEqVKhA37592bRpE59++imTJk1i4sSJt7QfMGAAPXr0YPLkyTRr1oynn36a0qVLM2bMGMaNG8dvv/0GQFRUFH/99Rdubm4cOnSIXr16ERoaelu7m8kVYNy4cUyZMoX69etz7do13Nzc0m0/atQoPDw82LVrFwAXL15k+/btnDx5kt27dwNw6dKldL9bdXdsdsZjvZ8yEPgD2AfMM8bsEZGRItLB2mwG4Gk9OxkCDLNuuweYB+wFlgIDjDEJafWZQSip7sOWansV4fnGDzE39AR/7PnX1rtTeYAxJtWRRJn5G2revHkEBgZSs2ZN9uzZw969e5PWPfbYYwDUqlWL8PDwO46rdu3alCpVCldXVx566CFatGgBgK+vb6r9tWzZkiNHjtCvXz/2799PzZo1Se3SdlxcHP369cPX15du3brdEnNa6tevz5AhQ/jss8+4dOkSTk7p/129bNkyBgwYkPS5cOHCVKhQgSNHjjBo0CCWLl1KoUKFMtyvunM2fYDUGLMEWJJi2Yhk76Ox3JtJbdsPgA8y02cqbdwzsw9bGtysMmsORTJ8wS5qlnuA4gXdsjsEZUNpnZnYio+PDz///PMty65cucKZM2eoUqUKu3fvJjHxv9GUN5+xOHr0KOPGjWPz5s0ULlyYkJCQW56/cHV1BSyX7uLjU7807OTklNS3MYbY2NjbtgdwcHBI+uzg4JBmf0WKFOHxxx/n8ccfp127dqxZswZPz1vH+0yYMIESJUqwY8cOEhMTcXPL+Odn2LBhtG3bliVLllC3bl2WLVuWbvvUknnhwoXZsWMHf/zxB1OmTGHevHl8/fXXGe5b3Rmdq81GXJwcmNgjgOsx8Qz9aWem/jJVKi1NmzYlKiqKWbNmAZCQkMCrr77KwIEDyZcvH15eXmzfvp3ExEROnDjBpk2bAEtyKlCgAB4eHpw5c4bff/89w30VLFiQq1evJn328vJiy5YtAPzyyy/ExcXd9XGsWLGCqKgoAK5evcrhw4cpV67cbfu8fPkypUqVwsHBge+++46EhIRUY0vu8OHD+Pr6MnToUIKCgti/f3+67Vu0aMHkyZOTPl+8eJFz586RmJhIly5dGDVqFFu3br3rY1Vp08RjQxWLF2R466qsPBDJDxuP2zsclYuJCAsXLuSnn36iUqVKeHp64uDgwJtvvglYLjN5e3vj6+vLa6+9RmBgIAD+/v7UrFkTHx8fnnnmGerXr5/hvtq3b8/ChQuTBhf069eP1atXExwczMaNGylQoMBdH8eWLVsICgrCz8+PevXq0bdvX2rXro2fnx9OTk74+/szYcIEXnzxRWbOnEndunU5ePBg0j5Ttktu4sSJ1KhRA39/f/Lly0fr1q3Tbf/WW29x8eLFpG1WrlzJyZMneeSRRwgICCAkJITRo0ff9bGqtIn+JX67oKAgk1WF4BITDX2+2URo+EX+91IDKhRzz3gjlSPt27ePatWq2TsMAP7++2969erFggULbht0oFR2SO3nQUS2GGOCMtpWz3hszMFBGNfNH1dnB16Zu11nNVBZ4uGHH+bYsWOadFSupIknG5Qo5MaHnX3ZEXGZqStvf7ZBKaXyEk082aSNbyk6BpRm0opD7D552d7hKKWU3WjiyUYjO9TA092FIfO2ExOfYO9wlFLKLjTxZCOP/M6M6eLHwTPXmPCXTiSqlMqbNPFksyZVitMr+EGmrzmsE4kqpfIkTTx28GZby0Sir87bQVSsTiSqMufm5Jz2Eh8fT9GiRRk+fLhd47gTEydOTHpgNTWRkZE4OzvzxRdfZGNUt1u8eHHSTN33KuUxt2nTJsfNOaeJxw7cXZ34uKs/4eejGPv7fnuHo+5zaU1dc6f+/PNPqlSpwrx589KciePmDAM5RUaJZ/78+dStW5cff/wxy/Z5N993hw4dGDYsa6aRTHnMS5Ys4YEHUhZmti9NPHZS7yFPnqnvzcwNx1gfds7e4ahc6tdff6VOnTrUrFmTZs2acebMGQDeffdd+vfvT4sWLXjqqaeIioqie/fu+Pn50aNHD+rUqcPNh6T//PNP6tWrR2BgIN26dePatWup7uvHH3/k5Zdfply5cvzzzz9Jy728vBg5ciQNGjRg/vz5bN68OWlmgtdff50aNWoApFu6wd3dnaFDh1KrVi2aNWvGpk2beOSRR6hQoQKLFy8GLEnt9ddfT5rp4OZZyqpVq3jkkUfo2rUrVatWpXfv3hhj+Oyzzzh16hRNmjShSZMmaR7TJ598QkREBCdPnkxa7u7uzquvvkpgYCBNmzZNmsj0kUceYfDgwTz88MPUqFEjaWqilN93dHR0UsmFmjVrsnLlSgDGjx/PM888A8CuXbuoUaMGUVFRt3w3ISEhvPDCCzRp0oQKFSqwevVqnnnmGapVq0ZISEhSjKmVu0jtmL28vDh37lzS/mvUqEGNGjWSZg4PDw+nWrVq9OvXDx8fH1q0aMGNGzdS/b6yjDFGXyletWrVMtnhRmy8aTJupan34TJz+UZstuxT3b29e/f+92HJUGO+bpO1ryVD091/gQIFblt24cIFk5iYaIwx5ssvvzRDhgwxxhjzzjvvmMDAQBMVFWWMMebjjz82/fv3N8YYs2vXLuPo6Gg2b95sIiMjTcOGDc21a9eMMcaMGTPGvPfee7ftJyoqypQqVcpcv37dfPHFF2bQoEFJ68qXL2/Gjh2b9NnHx8esX7/eGGPM0KFDjY+PjzHGmG+++cYMGDAgqV3btm3NypUrjTHGAGbJkiXGGGM6depkmjdvbmJjY8327duNv7+/McaYL774wowaNcoYY0x0dLSpVauWOXLkiFm5cqUpVKiQOXHihElISDB169Y1a9euTYotMjIy1e/z+PHjpmLFisYYY4YPH24++eSTpHWA+f77740xxrz33ntJcTdu3Nj07dvXGGPM6tWrk44t5fc9btw4ExISYowxZt++febBBx80N27cMAkJCaZhw4ZmwYIFplatWmbdunW3fTd9+vQxPXr0MImJiWbRokWmYMGCZufOnSYhIcEEBgaabdu2GWOMOX/+vDHGmPj4eNO4cWOzY8eOVI/55ufQ0FBTo0YNc+3aNXP16lVTvXp1s3XrVnP06FHj6OiY1G+3bt3Md999l+p3ltwtPw//fW+hJhO/Y/WMx47cnB35pJs//16JZtSvGU/7rlRKERERtGzZEl9fXz7++GP27PmvSkiHDh3Ily8fAOvWraNnz54A1KhRI6kezz///MPevXupX78+AQEBzJw5k2PHjt22n99++40mTZqQP39+unTpwsKFC2+5rNajRw/AUr/m6tWrPPzwwwA8/njKElypc3FxoVWrVoClpELjxo1xdna+pbzCn3/+yaxZswgICKBOnTqcP3+eQ4cso0ODg4MpW7YsDg4OBAQEZKrEw5w5c+jevTsAPXv2vOVym4ODQ9IxPfHEE6xbty5pXa9evQBo1KgRV65cSbp/kvL7fvLJJwGoWrUq5cuX5+DBgzg4OPDtt9/y5JNP0rhx4zTnzmvfvj0igq+vLyVKlMDX1xcHBwd8fHySji29chepWbduHZ07d6ZAgQK4u7vz2GOPsXbtWgC8vb0JCAgA7r5Exp2waVkElbGa5QrzwiMPMWXlYVr6lKRZ9RL2DkllRuusuRF8rwYNGsSQIUPo0KEDq1at4t13301al3wyT5PGPRljDM2bN8/wHsePP/7I+vXr8fLyAuD8+fOsXLmSZs2a3bKvtPYDt5ZXAG4pz+Ds7JxUoiCt8grGGCZNmkTLli1v6XfVqlW3lGdIr8RDymM6c+YMP/zwAwCnTp3i0KFDVKpU6ba2ycsnpCylcPNzZr5vgEOHDuHu7s6pUykLMv8n+fGnLD0RHx+fYbmL1KQXU8rvz9aX2vSMJwd4uWllqpYsyLAFu7hwPTbjDZSyunz5MmXKlAFg5syZabZr0KAB8+bNA2Dv3r1JVTfr1q3L+vXrk6qYRkVFcfDgwVu2vXLlCuvWreP48eOEh4cTHh7OlClTUk1WhQsXpmDBgkn3gObMmZO0Lq3SDZnVsmVLpk2bllSW4eDBg1y/fj3dbdIqi3DgwAGuX7/OyZMnk45p+PDhSfEmJiYmlfuePXs2DRo0SNp27ty5gOUMwsPDAw8Pj9v6b9SoUVJCO3jwIMePH6dKlSpcvnyZl19+mTVr1nD+/PlbSorfifTKXaR1zI0aNWLRokVERUVx/fp1Fi5cSMOGDe9q//fKpolHRFqJyAERCROR24ZsiIiriMy1rt8oIl7J1g23Lj8gIi0z6lNEZojIDhHZKSI/iYi7dXmIiESKyHbrq6+tjjchMYEVx1fc8XYuTg6M7x7A5RuxvP3LbhtEpu4HUVFRlC1bNuk1fvx43n33Xbp160bDhg0pWrRomtu++OKLREZG4ufnx9ixY/Hz88PDw4NixYrx7bff0qtXL/z8/Khbty7799860nLBggU8+uijt/xV3LFjRxYvXkxMTMxt+5oxYwb9+/enXr16GGOSfjGnVbohs/r27Uv16tUJDAykRo0aPPfccxme2fTv35/WrVvfNrjgxx9/pHPnzrcs69KlS1IyLVCgQFK58RUrVjBiRFL9SgoXLszDDz/M888/n2ZZ8BdffJGEhAR8fX3p0aMH3377La6urrzyyiu8+OKLVK5cmRkzZjBs2DDOnj17R98DpF/uIq1jDgwMJCQkhODgYOrUqUPfvn2pWbPmHe87S2TmRtDdvABH4DBQAXABdgDVU7R5Efjc+r4nMNf6vrq1vSvgbe3HMb0+gULJ+h0PDLO+DwEm30nsdzu4YP6B+abGtzXMxC0Tk2743onJKw6Z8kN/M4u3n7yr/SvbSu1mam4RHx9vbty4YYwxJiwszJQvX97ExMTYZF9Xr15Nej969Gjz0ksv2WQ/tpTaQA5jLIMLNm/enM3R5Ez3MrjAlvd4goEwY8wRABGZA3QEkt8B6wi8a33/EzBZLBdMOwJzjDExwFERCbP2R1p9GmOuWJcJkA/I9kJDnSt2Zs/5PXy16yuuxl7ljTpv4CCZP6l8rlEF/tp7hrd/2U0d7yIUL6TlslXWiIqKokmTJsTFxWGMYdq0abi4uNhkX//73/8YPXo08fHxlC9fnm+//dYm+1G5ly0TTxngRLLPEUCdtNoYY+JF5DLgaV3+T4pty1jfp9mniHwDtMGS3F5N1q6LiDQCDgKvGGOS95FlHB0cGVF3BIVcCvH17q+5EnuFDxp8gLODc6a2d3J04JPu/rT5dC3DFuxiRp+g225kKnU3ChYsSFYVN8xIjx49kkaE5VZpPct087kjdW9seY8ntd+YKc9C0mpzp8stb4x5GigN7ANu/p//K+BljPEDlgGp3oEVkf4iEioioTcfFrsbIsIrtV7h5cCX+f3o7wxeOZjo+PRHmyT3UDF3hraqyor9Z5kfGnHXcSilVE5ly8QTATyY7HNZIOX4waQ2IuIEeAAX0tk2wz6NMQnAXKCL9fN56yU7gC+BVEs2GmOmG2OCjDFBxYoVy+Qhpq2vb1/ervs2ayPW8uLyF4mKS3vajpRCHvaijncRRv62l4iLmd9OKaVyA1smns1AJRHxFhEXLIMHFqdosxjoY33fFVhhvUG1GOhpHfXmDVQCNqXVp1hUhKR7PO2B/dbPpZLtrwOWs6Fs0b1Kd0Y3HM2WM1sYsHxAppPPzXLZxhhen7+TxMRsv12llFI2Y7PEY4yJBwYCf2D5ZT/PGLNHREaKSAdrsxmAp3XwwBBgmHXbPcA8LPdqlgIDjDEJafWJ5RLcTBHZBewCSgEjrft4SUT2iMgO4CUso9yyTdsKbRndYDRbz269ozOfB4vk5+121dlw5Dzf/h1u2yCVUiobiUnnada8KigoyGT1jdjfj/7OsLXDCCgWwLRm08jvnD/DbYwxPDszlPVh5/jfSw2pWNy+0+Lndfv27aNatWp227+joyO+vr4YY3B0dGTy5MlJU9OkJjw8nHbt2rF79262b9/OqVOnaNOmTTZGfOfi4+MpWbIk/fr1Y/To0fYOR6UjtZ8HEdlijAnKaFuduSCbtPZuzZiGY9geuZ0Xlr2QqTMfEWFMF1/yuzgyZN524hISM9xG3b/y5cvH9u3b2bFjB6NHj76jujjbt29nyZIlNowua2Sm9MKdyqqyECrraOLJRq29WzO24Vi2R25nwPIB3IjPeD6k4gXdeL+TLzsjLjN15eFsiFLlBleuXKFw4cKA5cz4ZvkBX1/fpCldboqNjWXEiBHMnTuXgIAA5s6dS5s2bQgICCAgIAAPDw9mzpzJt99+S6dOnWjfvj3e3t5MnjyZ8ePHU7NmTerWrcuFC5aKuV9++SW1a9fG39+fLl26JNV+CQkJ4aWXXuLhhx+mQoUKSdPBpFW2IDXplV4YOnQowcHBBAcHJ03xExISwvPPP0/Dhg2pXLkyv/32G2ApwdCtWzfat29PixYt0vyOFi5cSLNmzTDGcPr0aSpXrsy///6bVf+ZVBp0ktBs1sq7FQkmgeFrhzNk1RA+a/IZzo7pP+fT1q8Uf+4tzaQVh3i0anF8y94+N5TKXmM3jWX/hawt4le1SFWGBg9Nc/2NGzcICAggOjqa06dPs2KFZXqmBQsWJJ0JnTt3jtq1a9OoUaOk7VxcXBg5ciShoaFMnjwZ+G826S1btvD000/TqVMnFi5cyO7du9m2bRvR0dFUrFiRsWPHsm3bNl555RVmzZrF4MGDeeyxx+jXrx8Ab731FjNmzGDQoEEAnD59mnXr1rF//346dOhA165dAdi2bRt79uyhdOnS1K9fn/Xr198y/9nN41u+fDlffPEFly5d4scff6RevXpJ6wsVKsSmTZuS4riZZMLDw1m9ejWHDx+mSZMmSUlpw4YN7Ny5kyJFivDzzz+n+h117tyZn3/+mSlTprB06VLee+89SpYseff/EVWm6BmPHbSt0Ja3673NupPrGLp2KPGJGV8KGNmhBp7uLrwybzvRcTmryqPKHjcvte3fv5+lS5fy1FNPYYxh3bp19OrVC0dHR0qUKEHjxo3ZvHlzhv2dO3eOJ598ktmzZyfNp9akSRMKFixIsWLF8PDwoH379gC3lCfYvXs3DRs2xNfXlx9++OGWUgydOnXCwcGB6tWrJxWlg8yVLcio9MLNcgS9evViw4YNScu7d++Og4MDlSpVokKFCklzzTVv3pwiRYoApPsdTZo0idGjR+Pq6pq0D2VbesZjJ90qdyMqLopxoeN45+93GFV/VLrT63jkd+bjrv489fUmxv1xgLfaVc/GaFVK6Z2ZZId69epx7tw5IiMj7+peSEJCAj179mTEiBFJFUKB26bgT608QUhICIsWLcLf359vv/32lqf5k2+fPK7MlC3IqPRCWqUJ7rVMwcmTJ3FwcODMmTMkJibi4KB/j9uafsN21MenDy/6v8jiw4sZvXF0hr9AGlUuxhN1yzFj/VH+OXI+m6JUOdH+/ftJSEjA09OTRo0aMXfuXBISEoiMjGTNmjUEBwff0j7lVPnDhg3Dz88vqTjcnbh69SqlSpUiLi4uaer/e5WZ0gs378vMnTv3lktw8+fPJzExkcOHD3PkyBGqVKlyW/9pfUfx8fE8/fTTzJ49m2rVqjF+/PgsOR6VPj3jsbPn/Z/netx1Zu6dSQHnAgyuNTjd9m+0qcbaQ+d4bf4Ofn+5IQXdMjcPnMr9bt7jActf8DNnzsTR0ZHOnTuzYcMG/P39ERE++ugjSpYsecvlrCZNmjBmzBgCAgIYPnw448aNw8fHJ6m/kSNHprbLVI0aNYo6depQvnx5fH19U639cqfSKr3wf//3f0mlF2JiYqhTpw6JiYm3JKQqVarQuHFjzpw5w+eff46b2+2T66b1HY0cOZKGDRvSsGFDAgICqF27Nm3btsXDw4O+ffvmipGAuZE+x5MKWzzHkx5jDKP+GcX8g/MZHDiYZ32fTbf9lmMX6Pb5BrrVepCxXf2yKUpl7+d48jIvLy9CQ0NvqzkUEhJCu3btkgYxqOyjz/HkciLCW3XforVXayZunciisEXptq9VvgjPNX6IuaEnWL7vTLptlVIqp9FLbTmEgzjwQYMPuBhzkXf/fpcibkVoVLZRmu0HN6vEyv1nGfrzLv58pTBFCtimtopSOUFqo+AArfWTS+kZTw7i7OjMxCYTqVKkCq+uepXtZ7en2dbVyZEJPSzlst9atCvLnvJW6dPvWal7/znQxJPDFHAuwNSmUymevzgDVwzkyKUjabatVqoQrzSvzJJd/7J4R8qKEyqrubm5cf78eU0+Kk8zxnD+/PlUB3Fklg4uSEV2Dy5IzYmrJ3hyyZM4OzrzXevvKFkg9aepExIN3b/YwKEzV/njlUaU8siXzZHmHXFxcURERBAdnfnCfkrdj9zc3ChbtizOzreOqs3s4AJNPKnICYkHYP+F/YQsDaFk/pLMbD0TD9fUp8oJP3ed1p+uJcirMLOeCdZy2Uopu9BRbfeBqkWq8lmTzzh+9TiDVgxKc1JRr6IFeLOt5fme7zcez+YolVLqzmjiyeGCSwVbyimc3c7wtcNJSEx9nrbedcrRqHIxPvzfPo5EXsvmKJVSKvNsmnhEpJWIHBCRMBEZlsp6VxGZa12/UUS8kq0bbl1+QERaZtSniMwQkR0islNEfhIR94z2kVu08GrB0OChLD++nHGh41JtIyJ81MUPV2cHBs/V2j1KqZzLZolHRByBKUBroDrQS0RSzmz5LHDRGFMRmACMtW5bHegJ+ACtgKki4phBn68YY/yNMX7AcSwlstPcR27Tu1pvnqz+JN/v+57v9n6XapuSHm6M7myp3TNx2cFsjlAppTLHlmc8wUCYMeaIMSYWmAN0TNGmIzDT+v4noKlY7ox3BOYYY2KMMUeBMGt/afZpjLkCYN0+H2Ay2Eeu81rQazQv35yPN3/MX8f+SrVNa99SdA8qy9RVh9moE4kqpXIgWyaeMsCJZJ8jrMtSbWOMiQcuA57pbJtunyLyDfAvUBWYlME+ch0HceDDBh/iV8yP4WuHp/mA6TvtfShfJD9D5u3g8o24bI5SKaXSZ8vEk9pZRcqx22m1udPlljfGPA2UBvYBPe4gDkSkv4iEikhoZGRkKpvkDG5Obkx6dBIl8pdg0IpBHLty7LY2BVydmNAjgH+vRDPil912iFIppdJmy8QTATyY7HNZIOXj9UltRMQJ8AAupLNthn0aYxKAuUCXDPZBiu2mG2PI0ZvQAAAgAElEQVSCjDFBxYoVy/RB2kNht8JMazYNQXhh2QtciL7tcKhZrjCDm1bil+2nWLTtpB2iVEqp1Nky8WwGKomIt4i4YBkssDhFm8VAH+v7rsAKY3midTHQ0zoizRuoBGxKq0+xqAhJ93jaA/sz2EeuVq5QOSY1ncTZqLMMWp76Mz4vNqlIUPnCvL1oNycuRNkhSqWUup3NEo/1fspA4A8sl77mGWP2iMhIEelgbTYD8BSRMGAIMMy67R5gHrAXWAoMMMYkpNUnlstpM0VkF7ALKAWMTG8f9wP/Yv6MbTiWXed2pfqMj6ODMKGHpdDXkHnbSUjM9flWKXUf0ClzUpFTpszJrB/2/cCYTWN4otoTDA0eetv6hdsieGXuDl5rUZmBj1ayQ4RKqbwgs1PmaD2e+0Dvar05ee0k3+39jlIFSvGUz1O3rO8UUIYV+yOZuOwQDSoVI+DBB+wUqVJK6ZQ5943Xgl6jWblmjAsdd9szPiLC+51qUKKQG4PnbON6TLydolRKKU089w0HcWB0w9H4FvNN9Rkfj3zOfNLdn2MXonjv1z12ilIppTTx3FduPuNTPH9xXlrxEieunLhlfd0Kngx4pCLzQiP4ZbsOsVZK2YcmnvtMEbciTG06lUQSeWH5C1yKvnTL+sHNKhFUvjBvLtzNsfPX7RSlUiov08RzH/Ly8GLSo5M4fe00L618iZiEmKR1To4OfNqrJg4Cg37cRmy8zmKtlMpemnjuUzWL1+TDhh+y7ew23lz3JonmvwRT5oF8fNTVn50Rl/n4j/3p9KKUUllPE899rKVXS4bUGsIf4X8wcevEW9a1qlGSp+qV58u1R1m5/6ydIlRK5UWaeO5zIT4h9KjSg292f8O8A/NuWfdGm2pULVmQV+fv4MyVaDtFqJTKazTx3OdEhGHBw2hUthEfbPyANRFrkta5OTsy+fFAbsQmMHiOTqmjlMoemnjyACcHJz5u9DFVClfhtdWvsef8f8/xVCzuzsiOPmw4cp4pK8PsGKVSKq/QxJNH5HfOz5SmU3jA9QEGLh/IqWv/VZPoWqssnQJKM3HZQTYdvb3EglJKZSVNPHlIsfzFmNp0KjHxMby47EWuxF4BrFPqdPalXJH8vDxnGxeux9o5UqXU/UwTTx5TsXBFJjaZyLGrx3hl5SvEJVhKY7u7OjH58UDOX49l8NztJOr9HqWUjWjiyYOCSwUz8uGRbPp3E+/8/Q43S2PUKOPBu+19WHMwksl6v0cpZSNaFiGPav9Qe05eO8mU7VMoU7AMAwIGANAr+EFCwy8wYdlBAssVpkGlonaOVCl1v7HpGY+ItBKRAyISJiK3Vf60lraea12/UUS8kq0bbl1+QERaZtSniPxgXb5bRL4WEWfr8kdE5LKIbLe+RtjymHOT5/yeo3PFzny+43MWHloI3LzfU4NKxd15ec42/r2sz/copbKWzRKPiDgCU4DWQHWgl4hUT9HsWeCiMaYiMAEYa922OtAT8AFaAVNFxDGDPn8AqgK+QD6gb7L9rDXGBFhfI1GAJcm8Xe9t6pWqx8gNI9lwagMA+V2cmNq7FtFxCQycvZW4BJ3PTSmVdWx5xhMMhBljjhhjYoE5QMcUbToCM63vfwKaiohYl88xxsQYY44CYdb+0uzTGLPEWAGbgLI2PLb7hrODM+MfGY/3A94MWTWEgxcPApbne0Z38SP02EU+WqrzuSmlso4tE08ZIHlBmAjrslTbGGPigcuAZzrbZtin9RLbk8DSZIvricgOEfldRHzu9oDuV+4u7kxtOpX8Tvl5cdmLnLl+BoAO/qXpY53Pbenuf+0cpVLqfmHLxCOpLEs5RjetNne6PLmpwBpjzFrr561AeWOMPzAJWJRqsCL9RSRUREIjIyNTa3JfK1mgJFOaTeFq7FUGLB/A9ThLrZ432lbD/8EHeH3+Dq3fo5TKErZMPBHAg8k+lwVOpdVGRJwAD+BCOtum26eIvAMUA4bcXGaMuWKMuWZ9vwRwFpHbhmoZY6YbY4KMMUHFihW7syO9T1QtUpVPHvmEsEthvLr6VeIS43B1cmTK4zVxcBCe/34rN2IT7B2mUiqXs2Xi2QxUEhFvEXHBMlhgcYo2i4E+1vddgRXWezSLgZ7WUW/eQCUs923S7FNE+gItgV7G/Fd8RkRKWu8bISLBWI75vE2O+D7QoEwD3qr7FutPrueDfz7AGEPZwvmZ2DOA/f9eYfiCnUnP/Sil1N2w2XM8xph4ERkI/AE4Al8bY/aIyEgg1BizGJgBfCciYVjOdHpat90jIvOAvUA8MMAYkwCQWp/WXX4OHAM2WPPMAusItq7ACyISD9wAehr9zZmurpW7cvLaSb7a9RVlC5alr29fmlQpzqvNKzPuz4P4ln2AZxt42ztMpVQuJfo7+HZBQUEmNDTU3mHYVaJJZNjaYfx+9HfGNhxLmwptSEw0vPDDFpbtO8t3zwbz8EP6cKlS6j8issUYE5RRO50yR6XKQRx4v/771CpRi7fWv0Xov6E4OAifdA/Au2gBBs7eRsTFKHuHqZTKhTKVeETkZxFpKyKaqPIQF0cXPm3yKWXcy/Dyypc5cvkI7q5OTH+yFnEJiTz//Rai43SwgVLqzmQ2kUwDHgcOicgYEalqw5hUDuLh6sG0ZtNwcnDihb9e4GzUWSoUc+fTngHsOXWF4Qt26WADpdQdyVTiMcYsM8b0BgKBcOAvEflbRJ6+OSeaun+VLViWqc2mcinmEs8ve54rsVd4tGoJXmlWmYXbTvLN+nB7h6iUykUyfelMRDyBECxzoG0DPsWSiP6ySWQqR/Hx9GFik4kcvXyUQcsHER0fzcAmFWlRvQQfLNnH34fP2TtEpVQukdl7PAuAtUB+oL0xpoMxZq4xZhDgbssAVc5Rr3Q9RjcYzbaz2xi6ZiiJJPBJd3+8PPMz4IetOrOBUipTMnvG85UxproxZrQx5jRYShoAZGbonLp/tPJuxdDgoaw4sYL3/3kfd1cnvupTm0QDz84M5Up0nL1DVErlcJlNPO+nsmxDVgaico/e1XrTz7cfPx/6mcnbJ+NdtADTnggk/Nx1Bs3eRryWUVBKpSPdmQtEpCSW2Z/ziUhN/puksxCWy24qjxpUcxAXoi8wfed0irgVoXe13rzX0Yc3F+7mwyX7GdE+ZeklpZSyyGjKnJZYBhSUBcYnW34VeMNGMalcQER4q+5bXIi+wNhNY/F086R3nVaEnb3G1+uPUrG4O4/XKWfvMJVSOVC6iccYMxOYKSJdjDE/Z1NMKpdwcnDio0Yf8dxfzzF83XA8XD14s00djkReZ8Qvu/Eqml+n1VFK3Sbdezwi8oT1rZeIDEn5yob4VA7n5uTGpKaT8PbwZvDKwey/uJdJj9fEq2gBXvh+K0fP6Ug3pdStMhpcUMD6rztQMJWXUhRyKcTnzT6nsFthnl/2PGduhDOjTxAOAs/O3MzlKB3pppT6j85OnQqdnfrunLh6gj6/98FgmNVqFqfO5efJGZsIKPcA3z0bjKuTo71DVErZUJbOTi0iH4lIIRFxFpHlInIu2WU4pQB4sOCDTG8+nfjEePr91Y/yJeL4uJsfm45e4LX5O0lM1D9ylFKZf46nhTHmCtAOS/npysDrNotK5VoVC1fk8+afcynmEv3/6k/Dqvn4v1ZV+HXHKT7+84C9w1NK5QCZTTw3JwJtA/xojLmQmY1EpJWIHBCRMBEZlsp6VxGZa12/UUS8kq0bbl1+QERaZtSniPxgXb5bRL6+OXmpWHxmbb9TRAIzeczqLvl4+jCl6RROXzvN8389T+96xXi8TjmmrTrM9/8cs3d4Sik7y2zi+VVE9gNBwHIRKQZEp7eBiDgCU4DWQHWgl4ikfKrwWeCiMaYiMAEYa922OpYy2D5AK2CqiDhm0OcPQFXAF8iHZTJTrG0rWV/9sZR4UDZWq0QtJjSZwKFLhxi4fCDD23jzaNXijPhlN8v3nbF3eEopO8psWYRhQD0gyBgTB1wHOmawWTAQZow5YoyJBeaksk1HYKb1/U9AUxER6/I5xpgYY8xRIMzaX5p9GmOWGCtgE5aHXm/uY5Z11T/AAyJSKjPHre5NgzINGNtwLDvP7eTV1UP4pLsPPqU9GDh7GzsjLtk7PKWUndxJRdFqQA8ReQroCrTIoH0Z4ESyzxHWZam2McbEA5cBz3S2zbBP6yW2J4GldxCHspEWXi147+H32HB6A+/8M4wvngrA092FZ77dzPHzWjpbqbwos6PavgPGAQ2A2tZXRkPmJJVlKYc1pdXmTpcnNxVYY4xZewdxICL9RSRUREIjIyNT2UTdrU4VOzE8eDgrT6xk3LYRfNUnkPhEwxMzNnL2arpXbJVS96GM5mq7KQiobu7soZ8I4MFkn8sCp9JoEyEiToAHcCGDbdPsU0TeAYoBz91hHBhjpgPTwfIcT/qHpu7U49UeJz4xno9DP8ZRHPnqqf/jqa+30OfrzczpXxePfFrIVqm8IrOX2nYDJe+w781AJRHxFhEXLIMFFqdosxjoY33fFVhhTW6LgZ7WUW/eWAYGbEqvTxHpi2VS017GmMQU+3jKOrqtLnD5Zk0hlb2e8nmKV2u9ytLwpcw/Po6pvQMIO3uVfjNDiY5LsHd4SqlsktkznqLAXhHZBMTcXGiM6ZDWBsaYeBEZCPwBOAJfG2P2iMhIINQYsxiYAXwnImFYznR6WrfdIyLzgL1APDDAGJMAkFqf1l1+DhwDNljGJ7DAGDMSWIJlGHgYEAU8ncljVjYQUiOEeBPPp1s/xVEcGdftBQbP3cnA2Vv5/IlaODneyW1HpVRulKkpc0SkcWrLjTGrszyiHECnzLG9L3Z8weTtk+lUsROVHJ5hxC976RJYlo+7+uHgkNptOaVUTpfZKXMydcZjjFktIuWBSsaYZSKSH8sZh1J35Tn/50gwCUzbMY0ulRwZ3LQXE5eHUaSAM2+0qYb1rFUpdR/KVOIRkX5YHr4sAjyEZTjy50BT24Wm7ncv+L9Agklg+s7pdKvswFN1O/Hl2qMUdHPmpaaV7B2eUspGMnuPZwCWhzc3AhhjDolIcZtFpfIEEWFgwEASEhOYsXsGj1VMoHPNjoz/6yCuTg481/ghe4eolLKBzCaeGGNM7M3LH9ahzzrkWN0zEeHlwJdxEAe+3PUlbb1jaRv/GKN/34+rkwMh9b3tHaJSKotlNvGsFpE3gHwi0hx4EfjVdmGpvEREeCnwJVwdXZm8fTItysfRLL4r7/66F1dnR3oFl7N3iEqpLJTZxDMMy4Seu7A8nLkE+MpWQam86Tn/53BxdGH8lvE8UjaWRok9eGPhLlwcHehSq2zGHSilcoXMjmpLFJFFwCJjjM4no2zm6RpP4+LowphNY6hfNo66CY/z+k87cHFyoL1/aXuHp5TKAukmHutM0e8AA7HMeSYikgBMsj6cqVSW612tNy6OLozaMIraZeMIjH+CwXO34+QgtPbVicWVyu0yekx8MFAfqG2M8TTGFAHqAPVF5BWbR6fyrG6VuzGq/ihCz24mX7lvqFHWhYE/buPXHbdNs6eUymUySjxPYZn77OjNBcaYI8AT1nVK2UzHih0Z03AMu8/txKnsF/iXc+TlOdtYtO2kvUNTSt2DjBKPszHmXMqF1vs8Op2wsrnW3q2Z1HQSx6+EE138MwK8Da/M285PWyLsHZpS6i5llHhi73KdUlmmQZkGfNniSy7FXOTyAxOo9VAsr/+0gzmbjts7NKXUXcgo8fiLyJVUXlcB3+wIUCmAgOIBfNPyGxJJ5N8Cn1Cr0lWGLdjFd/8cs3doSqk7lG7iMcY4GmMKpfIqaIzRS20qW1UpUoVZrWdR0KUgx13HE1Q1krcX7eartUfsHZpS6g5o8ROVqzxY8EFmtZ5F2YJlOez4KUHVj/P+//bxyZ8HuLMCuUope9HEo3KdYvmL8U3Lb/Ar6scBM5Xa/ruYtOIQI37ZQ2KiJh+lcjqbJh4RaSUiB0QkTESGpbLeVUTmWtdvFBGvZOuGW5cfEJGWGfUpIgOty4yIFE22/BERuSwi262vEbY7YpVdPFw9mN5iOs3LN2d/7A/UClzDd/8cZfDc7cTGJ2bcgVLKbmyWeETEEZgCtAaqA71EpHqKZs8CF40xFYEJwFjrttWxlMH2AVoBU0XEMYM+1wPNsJS/TmmtMSbA+tIZF+4Tro6ujGs8jqeqP8XBG7/jF7iYxTvD6f9dKDdiE+wdnlIqDbY84wkGwowxR4wxscAcoGOKNh2Bmdb3PwFNrdP0dATmGGNirA+vhln7S7NPY8w2Y0y4DY9H5UAO4sDrtV9nWPAwwm9spGrN2aw5fJQnZ2zk8o04e4enlEqFLRNPGeBEss8R1mWptjHGxAOXAc90ts1Mn6mpJyI7ROR3EfG5k4NQuUPvar2Z8MgEzsUepVyNr9l5Joxun//NyUs37B2aUioFWyYeSWVZyju/abW50+Xp2QqUN8b4A5OARak1EpH+IhIqIqGRkToBd27UtHxTvmrxFQlEUaTSF5yO3kPnKevZffKyvUNTSiVjy8QTATyY7HNZIOUMj0ltrFVNPYAL6WybmT5vYYy5Yoy5Zn2/BHBOPvggWbvpxpggY0xQsWLFMj46lSMFFA/ghzY/UCx/EZzKfEmi+0Z6fLGBVQfO2js0pZSVLRPPZqCSiHiLiAuWwQKLU7RZDPSxvu8KrDCWhzEWAz2to968gUrApkz2eQsRKWm9b4SIBGM55vNZcoQqRypXqBw/tP2B4FLBxDwwh0JlfufZmZt0ih2lcgibJR7rPZuBwB/APmCeMWaPiIwUkQ7WZjMATxEJA4ZgqXSKMWYPMA/YCywFBhhjEtLqE0BEXhKRCCxnQTtF5GaF1K7AbhHZAXwG9DT6pOF9r5BLIaY0nULvar255raCkpVnM2zRZsb9oQ+aKmVvoj+EtwsKCjKhoaH2DkNlkfkH5/PhPx/iRnH+PdSb9tV9+airH27OjvYOTan7iohsMcYEZdROZy5Q971ulbvxRfMvcHS+jmflafzv0Bp6fLGBM1ei7R2aUnmSJh6VJwSXCmZ229mULViCAuW/Jiz2V9pPXsuOE5fsHZpSeY4mHpVnlCtUjtltZ9O8fDMcPJcQ5zmTbtNX8st2rWiqVHbSxKPylPzO+RnXeByvBb1GotsuClaYxisL/uSjpft1glGlsokmHpXniAh9fPowvcV08ueLxuOhaUzf8gv9ZoXqNDtKZQNNPCrPqlOqDvPbz6Nq0YfIV/Z71l+YSbtJq9lzSmc6UMqWNPGoPK1kgZLMbDWTbpW74ey5misPfMZj0//H/NATGW+slLormnhUnufi6MKIeiP4qNFHuBU4Sz6vTxm2dA7DF+wkOk7LKyiV1TTxKGXV2rs189vPo7JnefI/OIsF4VPp+vlaTlyIsndoSt1XNPEolUy5QuX4vs339K7WGxfP9RxzHUvbaYtYuvu0vUNT6r6hiUepFFwcXRgWPIyJTSZSwP0ylJnAwMVf8YZeelMqS2jiUSoNTcs1ZUGHnwgoXp18Zeax4OQY2k1ZyoF/r9o7NKVyNU08SqWjtHtpvmn1NYMDB+PmsY8zhT6k49cz+O6fYzrLtVJ3SROPUhlwdHDkWd9n+bHtbMo/UBTnMjN4f8MHPDvrb85fi7F3eErlOpp4lMqkap7V+LnjPJ6o+gQuRTawMfZtmk35nj/2/Gvv0JTKVTTxKHUHXB1dGVpnKF+2+JJihYT4Ep/y0h+jeHnOJp1uR6lMsmniEZFWInJARMJEZFgq611FZK51/UYR8Uq2brh1+QERaZlRnyIy0LrMiEjRZMtFRD6zrtspIoG2O2KVV9QtVZdfOy+ic6XHcPFcw7KrQ2k+9RvWHoq0d2hK5Xg2Szwi4ghMAVoD1YFeIlI9RbNngYvGmIrABGCsddvqQE/AB2gFTBURxwz6XA80A46l2EdroJL11R+YlpXHqfKugi4FGVn/Xb5s8SUlCjlzo+gk+v72JkMXbuJqtJ79KJUWW57xBANhxpgjxphYYA7QMUWbjsBM6/ufgKYiItblc4wxMcaYo0CYtb80+zTGbDPGhKcSR0dglrH4B3hAREpl6ZGqPK1uqbr89tgielTuhUuRf/jt/Ks8OnU6y/aesXdoSuVITjbsuwyQfKbFCKBOWm2MMfEichnwtC7/J8W2ZazvM+ozM3GUAfLGo+iJiRB1Hq6egiunIfoSxFyFmCsQcw1MAhgDWIcGO7qASwFwcbf861oQChQD9xLgXtzyWd0mv3N+3qr3Bm0fas3/rXqTf50/Z+CyDTTY9gwfdniY4gXd7B2iUjmGLROPpLIs5YMPabVJa3lqZ2gZPUyRmTgQkf5YLsVRrly5DLrMgRIT4dxBOL0DIvdB5AGI3A+XIyAhNvVtxBEcnW9+ABGIj7Eko7Q4F4CCJaCwNxTxhiIVrO8rgOdDyfrLm2oWr8lvXRYyfcdXfLXrK/6JGUrTL1vxRoNn6FnbC8sJvVJ5my0TTwTwYLLPZYFTabSJEBEnwAO4kMG2GfV5N3FgjJkOTAcICgrK+U8GJibCqW1weAWc2AgRmyDaWkfGwQk8K0JJX6jWHgqVgUKloWBpyPcAuBaynLk4uVqSTXLGWJJP7HWIvWY5M7oeCdfOwtV/Lf9eOQkXwyEiFGKS1a5xcIZiVaCEj+VV3AdK1oCCJbPta8kJXB1dGRQ4gA4V2/Hm2pHscFzEqK2bmbOjD590bEfF4nrWqPI2sdXT19ZEchBoCpwENgOPG2P2JGszAPA1xjwvIj2Bx4wx3UXEB5iN5Z5OaWA5lsEBkok+w4EgY8w56+e2wECgDZbLcp8ZY4LTiz0oKMiEhobe+5eQ1WKj4NAfcGAphC2DqHOAQPFqULY2PFgHygRakk52nHkYAzcuwoWjcD4Mzu61vM7ssSSnmwqVgbJBUCbI8m+pAHDJb/v4cgBjDEuP/sHIDaO5FneR+Et16F6xP683r4m7qy3/7lMq+4nIFmNMUIbtbDnth4i0ASYCjsDXxpgPRGQkEGqMWSwibsB3QE0sZzo9jTFHrNu+CTwDxAODjTG/p9WndflLwP8BJYGzwBJjTF/rYIXJWEbHRQFPG2PSzSo5KvEkJsDR1bBzPuxbbDkLye8JFZtBpRbw0KOQv4i9o7xd1AVLEjq9E06GWs6OLlkHHIqj5YyoXD3wqg/l60OBoun3l8tdi73GuM2f8XPYHEx8PlyuteKtRs/SOeBBvfym7hs5IvHkVjki8URdgK0zYfMMuHzCcnmsekfw6275Re3gaN/47sa1yP+SUMRmyyvOWuumWDVLEvJqYDk+9+L2jdVGDlw4wIi1H7L30lYSYorjLb0Y374HVUrq5TeV+2niuQd2TTwXjsC6CbBzHsRHg1dDqP0sVG4NzvfZyKj4WDi9HcLXQvh6OP4PxF23rCtWzXI299CjUP7h++rSnDGG5cdX8t76MVyKO03Ctao0L9mXt1s+QlF3V3uHp9Rd08RzD+ySeM4fhjXjYOdcy/0Z/14Q3B9KpHzm9j6WEGcZlRe+Fo6sgmMbICEGHF2hfD14qKklEZXwuX1QRC4UmxDLVztm8eWu6cSZGLjyMCHV+zKwcQBuzrnwjFbleZp47kG2Jp5rkbDyfdg6y/ILNugZqP9SnhsJlqrYKDj+N4StsIzei9xnWe5ewpKAKjaDik0hX2H7xnmPzt84z4cbJvDnicWYBGdcrz/K6/X60T2wEg4OuT/BqrxDE889yJbEExcNG6fBmk8g/gYEPQuNXrtv721kicsn4chKSxI6vBJuXLAMVChXDyq3hMqtoGilXHs2dOTyEd5b+wlbz68hMb4ARePbMLJJXxpXLm3v0JTKFE0898DmiefoWvj1Zbhw2HLvpsUoyy9MlXmJCXByKxxcCgf/gDO7LMsLe1uTUEvLIAWn3HfPZMfZnby95iOOXt9BYmxhyjl0ZmTTJ6jtVczeoSmVLk0898BmiefGRfjzbdj2HRT2grbjLZeK1L27HGFJQAf/sAw/j4+2TPvzUBPLmVClFrnqbNIYw+oT6xj19zjOxhwhMaYoDzl3YlSzJ/B/0NPe4SmVKk0898AmiSdsOSx6Aa6fg4cHQuNh99VIrRwlNsoyQOHm2dDNh1lLB0KV1pZEVNI3V1ySSzSJLDnyJ+M2TuF8XDiJMUWp4taZ95s9SfXSufvelrr/aOK5B1maeOJjYPlI2DAZilWFzl9A6YCs6VtlzBg4s9sy28PBpXByC2AssylUbmm51OndKMcPVU80ifzv8F+M2zSFC3FHSYz1pIprZ0Y88gQB5fQMSOUMmnjuQZYlnnOH4Ken4d9dULsvtHgfnPPde7/q7l07C4f+hAO/WwYoxF0H5/xQocl/AxQKlrB3lGkyxvBr2F98smkyF+KPkhhbhHJOrRne4EkaVSqTcQdK2ZAmnnuQJYln36+w8HnLze2OUyyXeFTOEhcN4eusl+SWWmaIgGSX5FpCSb8ceUnOGMPSI8v5ZPPnnIk5QGJ8foqZRxlS92na+1TSaXiUXWjiuQf3lHgSE2DVaFjzMZSpBd2/Aw/9SzTHM8YyuenNJBQRyu2X5BrmyDPWjadCGf335xy+vhGT6Ix7XD2eqdGHPsG1cHXSB1FV9tHEcw/uOvHcuAQL+lku5dR8Atp8kuPvHag0pHlJ7hHL5bjKLXPcQ74Hzofx/rppbL+4AkMCjjd8aV2uC681akNRLUSnsoEmnntw14ln5zxY9CK0HmuZgUAvd9wf4mMso+QOpLwkV9NyJlSlVY66JHf2+lnGbviK5RGLSZDrJMaUxK9QG4Y37I1fmdwzpFzlPpp47sE9XWo7f9hSiVPdn4yxlHs48Putl+QKlracBVW5OUrO/pfkouOjmblzIbP2zOZKYjgmwY2ipgEhfo/zeM1auDilVtBXqbuniece5IiyCCp3uBZpuSR30HpJLvZajrskZ4xh7fFQJv0VaaoAABK/SURBVGz6mrDrf2MwOERXoVGpNgyp/xgVinrYNT51/9DEcw808ai7kgsuyZ2+dobx/8xixcnfiP3/9u49OK6zvOP499nV1ZKsuyVblizZlh1fYnzDcUKaGKdpTEpioBAbCm1opgwplCl0GMgwQxloCwQYCC2FJiHDbZIQaIodSHBJYkggsWPHjuNLfJctyZas+8W6rLS7T/94jyxZluS1pb1Ifj4zZ3b3nLPnvPuO4l/Oe97zvrSgwQyK/e/go9dvYtPy1ST57SrIXL2ECB4R2QA8hJst9FFV/dqw7anAT4BVQDOwSVVPedseAO4DQsCnVHXbWMcUkQrgSSAP2AN8RFX7RORe4Bu4qbIB/lNVHx2r3BY8ZtzGapKb703vMHdd3GaPDYVD/PrY73nkjSc43bMLJIwEKlhT8C4+ueZ9LJ+duM8ymcQV9+ARET9wFLgdqAV2AR9U1UND9vkHYJmqflxENgPvVdVNIrIYeAJYA8wCngcWeF8b8Zgi8hTwtKo+KSI/APap6ve94Fmtqp+MtOwWPGbCXWiS+y2c/AME2kF87pmheetdGJWsBn9SzItWf76R7+x8ghdqn6GXetclO7iCDeV38om176Iwy4Z2MpFJhOC5EfiSqt7hfX4AQFW/OmSfbd4+r4pIElAPFAKfH7rvwH7e1y45JvA1oBEoVtXg0HNb8JiEEwrC2T1u/L4TL7hhfDQMqdkw9xZv5tXbIHdOTIulqmw/vZNH9/6SA20vob4eNJjJrKS1vH/hRj6y8mbSU2IfjGbyiDR4ovlXVALUDPlcC9ww2j5eYLQD+d76HcO+O/AU5kjHzAfaVDU4wv4AfyUit+Culj6tqkOPYUxs+ZOgdI1b3vmAG7X85B9cCB1/0Y16AZA3z10JzV3npnhIz4lqsUSE9eVrWV++lr5QH4/v38bP39pCbWA7/3Hkeb67v4C56TfzgUXv5gPLVpGWbCFkrk40/3JGuoM6/PJqtH1GWz/Snc+x9gd4BnhCVQMi8nHgx8D6Swor8jHgYwBlZWUjHM6YKEnPhSXvcYuqG+PvxAtuwru9P4PXHnbNcsXL3OgJFbdC2VpIzYpakVL8Kdy7/C7uXX4XrT1tPPz6Fp499Ruqglt48MCveHBPIXOn3cQHFt/JPdffQIqNkGCuwJRvahtWJj/Qoqpj9h+1pjaTMIIB1zHh1MtQ9RLU7oJQn5t5tWSVF0S3QOkNMXl26EznOX645xleqP4dzaG3EFHoz2dO2g3cNX8Df738ZjLTkqNeDpOYEuEeTxKuaes2XI+yXcCHVPXgkH0+AVw/pHPB+1T1HhFZAjzOYOeCF4BK3JXNiMcUkV8A/zOkc8GbqvpfIjJTVeu8870X+Jyqrh2r7BY8JmH1dUPNzsEgOrMHNAT+FJi9BspvdldDs98OqZlRLUpdZyOPvP4ML9Y8T3PoIEgYDU5nhn85t5bewkdX3kFZbnSbB01iiXvweIW4E/gOruvzY6r6byLyZWC3qm4VkTTgp8AKoAXYrKonve9+Afg7IAj8k6o+N9oxvfVzGexOvRf4sNe89lXgbu84LcD9qnp4rHJb8JhJI9AJp1+FUy+5IKrf7zoqiB9mLoOyG10Qld0Y1RlYW7rb+NEbz/K709up7d0LvgAaTiIjvJAVBTdxz+I7WDdvAT5fYgwrZKIjIYJnsrLgMZNWb7trjqve4QLpzG43DTi4zgoDQTTnJsibG5WHWfuCfWw5/CeePvw7DnfsJOhvAED6iilNX8G6snfwoWXrKMmxEROmGguecbDgMVNGsA/q9kH1q4NLT6vbNi3fPTs0e7W7X1SyKio95/bWHeXx/dvYWf9HWsNHQEJo2E96eB6LclZxx9w/4z2LbyAjNWXCz21iy4JnHCx4zJQVDkPzMTj9irsyqt0NTUcGt+dXXhxERUshaeICoauvm/996488d/xljnS8TsDnnmzQUDrTuY6leSu5fe6NvHvRStKTrZPCZGPBMw4WPOaa0tvuOimc2Q21r7vXrka3zZ8KM9/mQmjmMve+YOGEjbBQ3dbAUwdf5OWaVzjds4+QrwXwgkjmsyh3OevL13LXwjVMT7c5hRKdBc84WPCYa5oqtFVfHER1b0Kwx21PSoMZi10IzXybC6QZSyZk0sMjTad5+tDL7Di7m+ruAwT951yRwsmkhyuoyLyet89awbsXruG6wmKb4jvBWPCMgwWPMcOEQ+7B1vo33T2jun0ujALtbrv4YcYi95Br8fXufdESyCgcVweGqpZ6thz+E6+c2UXV+f30SI17dgiQYAEFyZUsylvCzaUruaNyJXnTMibi15qrZMEzDhY8xkRAFVpPDQujfYPNdOA6MMxY7JYi73XGoqsedaG9t5Ntx3bzh+rXeavlIE39x1B/u1ccP8mhEopTK7kubxE3lS7jtnnLLIxiyIJnHCx4jBmH843QcBAa3oJz3mvDW9DfNbhPdpkLosLroGCBt8x3wwddoYMN1fzf8dfYXbePqvOH6QxXgS8AgKqP5FAxhSkVzMtZyMriJbyzYjnz8gutmS4KLHjGwYLHmAkWDkN7NZw75OYpajjkwqjpGIT7B/fLKHQ96woGlgXuNWcO+CIbDy6sYfbVneDFqjfYU3+A053HaQ+fBu/KCIBgLpkym1nTylmYV8mqWddxa8VSCjKiN/7dtcCCZxwseIyJkVAQ2k5D01EXQk1Hofm4e+1uHtzPn+IeeM2bB3kVkFsOuRXufXZpRF2+TzTXsb1qH6/XH+Bk+1GaAtUEfPWIhABQFfyhPLL8JRSnlzM/dx7LixaypnQhFbn5doUUAQuecbDgMSYBdLe4MGo+NhhMLVXuvtJADztwI3dPnw155S6McssHwylnjmu+GyU0+oL97DpzjFdrDnKg8SjVnVW0BWvo8527EEgAhDJIo4jclFmUZs5hQV4FS4vmsXpWJUVZNgLDAAuecbDgMSaBqUJnvQug1qrBMBp439108f4pmZA9e9hSOvg+a9YlV0x9wX721J3gtdpDHG4+SXVHNU2BM3TruQudGS4ITieNInJSiimaNpOy6SUsyC9jWfFclswoJTXp2nkQ1oJnHCx4jJnEAp0uiFqqoL0G2msHX9tqLg0mBLKKLw6irGLImglZRd5r8YWeeI1dnbxWc4QDDSc43lrFma4aWvrO0h1uvCSUVH34wzmkSyHZyTMoTJ/J7KwSKrJLWFAwm8VFpRRlTo9NvcSABc84WPAYM4X190D7mSGhVDsknGrc1VR/96XfS84YOZAyi93njEK6kzPZ33Ge/Y3VHGuuprqzlobuOjqCDQRoIuzruPAc0gANpZGsOaT788hOLiA/rZCiaUWUTC+iIqeE+fmzqMyfOSlmfLXgGQcLHmOuYaruqqmzHjrr4Pw599pZP2TxPg+913SBuOeXMgpcL71p+e41o4Ce1GyOhX0cC/RxItDHqd7znO3toLWvifOhFvq0FfV3IhIeViQfEsogiemk+bLJTMolJzWX/PR8ijIKmJlZSElWIeW5RZTnFJGVFp/hhSINnsSPUGOMiSURSJvulsIFo++nCoEO6DwH5+uhq8kt3U3uIdquRuhqds8ydTVCbxvpwDJvGXJCSMt2nSDSc+lPqqAuOZ0aSaEaHzXhMGfDQRqln+ZQgI5wG419Z6kLdiE9/W6WseFFC6Xj1yySJZM0XxbT/FlkpWSTk5pDXnoOBdNyKcrIozgzj5Lp+ZRmF5CdlhGznnsWPMYYczXEC4y07LEDakCo33UR7xoIJi+kelovWpJ7WilrP01ZTyvv6GkDRm6VUqA7ZRoNqVk0JE/jXFIqDb5kGnx+GgWaCdIuLXSGG2gPBTkX7CfcG4L2EQ+HhpPwaQY35G/kkY2fvepqiYQFjzHGxII/2btHVBz5d8JhNx7ehWBqu/BeetrI6G2jItBBRW+Hax4MdLqrsEAn9HZAX+dFhwsItPv8tPl8tPt9tPt8tPr9tPpTaUlKpsXXz8KeExP8wy8V1eARkQ3AQ7hpqh9V1a8N254K/ARYBTQDm1T1lLftAeA+IAR8SlW3jXVMEalgcOrrPcBHVLVvrHMYY0xC8/kuNMFdlXAY+s5fCKPU3g5mBDqZEegYDKhAJwTOu5Dq64IF6yf2N4wgasEjIn7ge8DtQC2wS0S2quqhIbvdB7Sq6nwR2Qx8HdgkIouBzcASYBbwvIgMXMuOdsyvA99W1SdF5Afesb8/2jmi9buNMSZh+HyD96sSiC+Kx14DHFfVk6rah7sa2Thsn43Aj733vwRuE3d3ayPwpKoGVLUKOO4db8Rjet9Z7x0D75jvucw5jDHGxEE0g6cEqBnyudZbN+I+qhrE3fbKH+O7o63PB9q8Yww/12jnuIiIfExEdovI7sbGxuGbjTHGTJBoBs9IVxXDu2eMts9ErY+0HKjqw6q6WlVXFxYWjvAVY4wxEyGawVMLlA75PBs4O9o+IpIEZON6pY/23dHWNwE53jGGn2u0cxhjjImDaAbPLqBSRCpEJAXXWWDrsH22An/rvX8/8KK6oRS2AptFJNXrrVYJvDbaMb3vbPeOgXfMLZc5hzHGmDiIWq82VQ2KyCeBbbiuz4+p6kER+TKwW1W3Aj8Efioix3FXIZu97x4UkaeAQ0AQ+ISqhgBGOqZ3ys8BT4rIvwJ7vWMz2jmMMcbEh43VNgIbq80YY65cpGO1RbOpzRhjjLmEXfGMQEQagdPxLsc4FeA6XRjH6uNiVh+DrC4uNp76mKOql+0WbMEzRYnI7kguea8VVh8Xs/oYZHVxsVjUhzW1GWOMiSkLHmOMMTFlwTN1PRzvAiQYq4+LWX0Msrq4WNTrw+7xGGOMiSm74jHGGBNTFjyTnIhsEJEjInJcRD4/wvbPiMghEXlTRF4QkTnxKGesXK4+huz3fhFREZmyvZkiqQsRucf7+zgoIo/HuoyxFMF/K2Uisl1E9nr/vdwZj3LGgog8JiINInJglO0iIt/16upNEVk5oQVQVVsm6YIbNugEMBdIAfYBi4ft805gmvf+fuDn8S53POvD2y8LeAnYAayOd7nj+LdRiRteKtf7PCPe5Y5zfTwM3O+9Xwycine5o1gftwArgQOjbL8TeA43uv9aYOdEnt+ueCa3y062p6rbVbXb+7gDN3L3VBXJ5IMAXwEeBHpjWbgYi6Qu/h74nqq2AqhqQ4zLGEuR1IcCA1N1ZnPpaPpThqq+xNij9G8EfqLODtzo/zMn6vwWPJNbJJPtDXUf7v9ipqrL1oeIrABKVfXXsSxYHETyt7EAWCAifxKRHSKyIWali71I6uNLwIdFpBZ4FvjH2BQtIV3pvy1XJGqjU5uYiGiSOwAR+TCwGrg1qiWKrzHrQ0R8wLeBe2NVoDiK5G8jCdfctg53JfyyiCxV1bYoly0eIqmPDwI/UtVviciNuFHtl6pqOPrFSzgR/9tyNeyKZ3KLZLI9ROTPgS8Ad6tqIEZli4fL1UcWsBT4vYicwrVdb52iHQwinYhxi6r2q2oVcAQXRFNRJPVxH/AUgKq+CqThxi27FkX0b8vVsuCZ3C472Z7XtPTfuNCZym34cJn6UNV2VS1Q1XJVLcfd87pbVafiHBiRTMT4K1znE0SkANf0djKmpYydSOqjGrgNQEQW4YKnMaalTBxbgb/xeretBdpVtW6iDm5NbZOYRjbZ3jeATOAXIgJQrap3x63QURRhfVwTIqyLbcBfiMghIAR8VlWb41fq6ImwPv4ZeEREPo1rVrpXvS5eU42IPIFrYi3w7mn9C5AMoKo/wN3juhM4DnQDH53Q80/RejXGGJOgrKnNGGNMTFnwGGOMiSkLHmOMMTFlwWOMMSamLHiMMcbElAWPMXHmjZL90yGfk0SkUUSm+rA+5hplwWNM/HUBS0Uk3ft8O3AmjuUxJqoseIxJDM8Bf+m9/yDwxMAGEVkjIq9488S8IiILvfVLROQ1EXnDmzOlUkQyROQ3IrJPRA6IyKY4/BZjxmTBY0xieBLYLCJpwDJg55Bth4FbVHUF8EXg3731HwceUtXluAFga4ENwFlVfZuqLgV+G6sfYEykbMgcYxKAqr4pIuW4q51nh23OBn4sIpW4oVySvfWvAl8QkdnA06p6TET2A98Uka8Dv1bVl2PyA4y5AnbFY0zi2Ap8kyHNbJ6vANu9K5i7cINXoqqPA3cDPcA2EVmvqkeBVcB+4Ksi8sVYFd6YSNkVjzGJ4zHcKMD7RWTdkPXZDHY2uHdgpYjMBU6q6ne998tE5DDQoqo/E5HzXBtzD5lJxoLHmAShqrXAQyNsehDX1PYZ4MUh6zfhZszsB+qBLwNvB74hImGgH7g/uqU25srZ6NTGGGNiyu7xGGOMiSkLHmOMMTFlwWOMMSamLHiMMcbElAWPMcaYmLLgMcYYE1MWPMYYY2LKgscYY0xM/T+hc1ub+7hzRQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(mass,i, label=\"Quantum Statistics\")\n", "plt.plot(mass,l, label=\"Large Argument Approximation \")\n", "plt.plot(mass,k, label=\"Boltzmann Approx.\")\n", "plt.legend(loc=\"upper right\")\n", "plt.ylim()\n", "plt.xlabel(\"Mass\")\n", "plt.ylabel(\"Density\")" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Abundance Calculated from Statistical Model: [0.9970017500294499, 0.0029919632572737082, 6.2867132764680795e-06]\n", "Abundance of Nuclei from ALICE : [0.9868356129235992, 0.013157808172314654, 6.578904086157327e-06]\n", "0.003000960888168281 : 0.013333333333333332 6.305619098745192e-06 : 6.666666666666666e-06\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:12: RuntimeWarning: overflow encountered in exp\n", " if sys.path[0] == '':\n" ] } ], "source": [ "#masses in GeV/c^2\n", "md = 1.8756\n", "m3He = 2.8084\n", "m4He = 3.7274\n", "\n", "#Abundance of species (given on Page 8 of lecture) for fixed T\n", "T=.1565 #GeV \n", "mu=0\n", "\n", "\n", "def n(p,m,g,s):\n", " return (g*(4*np.pi)/((2*np.pi)**3))*(p**2)/(np.exp((np.sqrt(p**2+m**2) - mu)/T)-s*1)\n", "\n", "n_md,err=quad(n,0,np.inf,args=(md,3,1))\n", "n_m3He, err=quad(n,0,np.inf,args=(m3He,2,-1))\n", "n_m4He, err=quad(n,0,np.inf,args=(m4He,1,1)) \n", "\n", "raw=[n_md,n_m3He,n_m4He]\n", "normed = [i/sum(raw) for i in raw]\n", " \n", "\n", "# Estimate of abundances from graph on page 25 \n", "a_d=6*10**-1\n", "a_3He=8*10**-3\n", "a_4He=4*10**-6\n", "estimate=[a_d, a_3He, a_4He]\n", "normed_e = [i/sum(estimate) for i in estimate]\n", "\n", "\n", "print('Abundance Calculated from Statistical Model: ', normed)\n", "print('Abundance of Nuclei from ALICE : ', normed_e)\n", "print(n_m3He/n_md, ':' ,a_3He/a_d, n_m4He/n_md, ':', a_4He/a_d,)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }