msfig_mass_rho.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b29751cc-d525-43f5-add2-9f00db1463ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, basic, cosmology, local, tools_multitracer as mass, plottools as pl\n",
    "from matplotlib.pyplot import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7923204d-8db1-428f-bf28-8c78aed84901",
   "metadata": {},
   "outputs": [],
   "source": [
    "lmax = 2048\n",
    "euc_zbin = 5\n",
    "lss_zbin = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "035291a3-9ebe-43e3-85ed-3f50846968f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "frac = mass.galaxy_distribution(np.linspace(0,50,1000),zbn={'euc':euc_zbin,'lss':lss_zbin})[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fe8e577c-4948-4f7f-b71c-4732f8785844",
   "metadata": {},
   "outputs": [],
   "source": [
    "l, camb_cls = mass.read_camb_cls(add_euc=euc_zbin,add_lss=lss_zbin)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "61d159c8-3954-4974-bcad-61f1a6193d73",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(12, 12, 2049)\n"
     ]
    }
   ],
   "source": [
    "# for cross spectrum\n",
    "vecs = {}\n",
    "ocov = {}\n",
    "Scov = mass.get_covariance_signal(lmax,add_cmb=['klb'],add_euc=euc_zbin,add_lss=lss_zbin)\n",
    "print(np.shape(Scov))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d2c76a8b-aabe-4e04-968a-0f84d8e64b46",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CIB covariance\n",
    "scib = mass.get_covariance_signal(lmax,add_cmb=[],add_euc=0,add_lss=0)\n",
    "ncib = mass.get_covariance_noise(lmax,add_cmb=[],add_euc=0,add_lss=0)\n",
    "vecs['cib'] = Scov[0,1:2,:]\n",
    "ocov['cib'] = scib + ncib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "be0c3db9-de9a-4368-98a1-cd87e65138ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Euclid covariance\n",
    "seuc = mass.get_covariance_signal(lmax,add_cmb=[],add_cib=False,add_euc=euc_zbin,add_lss=0)\n",
    "neuc = mass.get_covariance_noise(lmax,frac=frac,add_cmb=[],add_cib=False,add_euc=euc_zbin,add_lss=0)\n",
    "vecs['euc'] = Scov[0,2:2+euc_zbin,:]\n",
    "ocov['euc'] = seuc + neuc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "9578dfd9-83a5-4bc5-9f85-af84b9d3efd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# LSST covariance\n",
    "slss = mass.get_covariance_signal(lmax,add_cmb=[],add_cib=False,add_euc=0,add_lss=lss_zbin)\n",
    "nlss = mass.get_covariance_noise(lmax,frac=frac,add_cmb=[],add_cib=False,add_euc=0,add_lss=lss_zbin)\n",
    "vecs['lss'] = Scov[0,2+euc_zbin:2+euc_zbin+lss_zbin,:]\n",
    "ocov['lss'] = slss + nlss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f745ce88-deb6-4ef2-a69b-9719184eef25",
   "metadata": {},
   "outputs": [],
   "source": [
    "# gal covariance\n",
    "sgal = mass.get_covariance_signal(lmax,add_cmb=[],add_cib=False,add_euc=euc_zbin,add_lss=lss_zbin)\n",
    "ngal = mass.get_covariance_noise(lmax,frac=frac,add_cmb=[],add_cib=False,add_euc=euc_zbin,add_lss=lss_zbin)\n",
    "vecs['gal'] = Scov[0,2:2+euc_zbin+lss_zbin,:]\n",
    "ocov['gal'] = sgal + ngal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "83c6fc26-4605-4b19-abde-4084658ce86a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cib+gal covariance\n",
    "sext = mass.get_covariance_signal(lmax,add_cmb=[],add_euc=euc_zbin,add_lss=lss_zbin)\n",
    "next = mass.get_covariance_noise(lmax,frac=frac,add_cmb=[],add_euc=euc_zbin,add_lss=lss_zbin)\n",
    "vecs['ext'] = Scov[0,1:2+euc_zbin+lss_zbin,:]\n",
    "ocov['ext'] = sext + next"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "bd4e31de-92cb-4ddc-97ba-cc21d8b88b34",
   "metadata": {},
   "outputs": [],
   "source": [
    "icov = {}\n",
    "rho  = {}\n",
    "for t in ['cib','euc','lss','gal','ext']:\n",
    "    icov[t] = ocov[t].copy()\n",
    "    icov[t][:,:,1:] = np.array( [ np.linalg.inv(ocov[t][:,:,l]) for l in range(1,lmax+1) ] ).T\n",
    "    rho[t] = np.array( [ np.dot( vecs[t][:,l], np.dot( icov[t][:,:,l], vecs[t][:,l] ) ) for l in range(lmax+1) ] ) / ( camb_cls['klbklb'][:lmax+1] + 1e-30 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "11e882ae-e6eb-4f56-af8e-0dac2cbf5529",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEPCAYAAABV6CMBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABF2UlEQVR4nO3deVzc133v/9eZfViHHUloA1kbWgHvm2wJ24kTJ3YkK829SW+TSiRNb2rfNCJO+0vTLNeR29vmps1NJDdtHjdNbmzZztKkji3Ju2JbFmixdhnQgiQkEAzDMvv3/P74DggkQIBgmIHP0495DDPzne8cMWbec3altUYIIYQYLstEF0AIIURykeAQQggxIhIcQgghRkSCQwghxIhIcAghhBgRCQ4hhBAjIsEhhBBiRGwTXYAeSikPsBHwaq23DnLMWsALeIB6rXVtvMonhBDClDDBAVQM9WAsWCq11lWx29uByjiUSwghRB8J01Sltd6BWZsYzKNAXZ/bXqVU2bgWSgghxFUSqcZxLR76B0srkD3QgUqpjZjNXqS50ssXlM4f77IJIcSkUlNT06K1zhvosWQKjoF4Broz1keyFeCGGQv0nj174lkmIYRIekqpU4M9ljBNVcPgpX9QZAP1E1ISIYSYwpIpOJ4FSvrc9sioKiGEiL+EaapSSq3BHCXlUUq1aq2fi91fA6zWWnuVUttjxwFsnqiyCiHEVJYwwREbVbVjgPvL+/z8XFwLJYQQ4irJ1FQlhBAiAUhwCCGEGBEJDiGEECMiwSGEEGJEEqZzXAzNMAzCrR0EmtoItXUQ6QgQ7goQ9YeIdIeI+kNEgxEiwSjRYJRo2CAa1hgRiEbAMBSGVmitANBaoQFQgELHLj0/x47CPOqKa2UeaVEGFovGYgGLTWOxKiw2hdWmsNitWB0WrA4rVqcNe5oTe7obhycFR1Ya9qx0nLkZ2DJTsVjk+4sQyUSCY4IYkSiBc5foOn0R/9lWAs0dBFq7CXaEiAQ04TBEIlai2kZUOYla3aCsA5xJAc7YJUZHsUTDKB3GQhgLUfNDXhkoi/mhj9IoFYsNMy9QSqOUQinQaNAKrbUZNhp07yUWI4YiHLFiYEGHrBjY0MqGYbGDZaD/tSKAL3bp+UVEsRhBrDqITYWw2aLYneBIseJIs+P0uHHmpOLKy8RdmIV7dgH2NPcYvQtCiNGQ4BgHRjhCV905Oo6fpfN0K90XOwn6wgQDEA7biOAkYk0DS98gcANuVDSA1QhgJYTNGsVhC2J3BrG7u7Cn2nGku7CnO7GlOLGlObGnurCmu7GnubFnpGDPTMWa4pqof3ovIxwh0ukn0hUg0tFNqK2TkLeLcHs3YV+AcGeAcHeYiD9CJBAlHIRIWBEK2+kOO4n6U6DNBmcADKAtdqnHGunCprtx2MI43eDMsOPOceMu9JA2N5/0hTOxZ6RO6L9fiMlMgmMUDMMgcO4SHcca6Tx5ka5zPvytQQJdimDUSdia0ecbtxu0E2ukGxt+7NYIKa4unGkBnJkuXLkpuAs9pBTlkjIrH0dW+nWVTUfDRLtaMfztGIEOdLADw+/DCHRgBGL3+X0Yfe7XkQA6EoJICB0NoiNhdCSIjprXREPoaOQar6xQNjvK5kTZHCirE3p+jl1bbU5sdiep6emovHQsrtjFmY5yZVy+nZKFSskm4rcRPN+Gv6mNYEsHwdZOAm1+Au1hgn4IhW10d6YQDaTBReBIBF49B5zDGvbhVN243AbuLDuphWmkzskjY8EMUuZOk+YxIa6DBMcQooEg3n31eA+exneyja5LEQIhOyFLOtra863eDuRgiXTh0F2kuPy4M8Kk5KWQOjObjHnTSF9QdF21AB0JEfU1EfGeI9p+jqjvItHOZqKdLRidLUQ7W4h2NmN0tmD42699QosNi9v8oFbODCx2l/kh70rDYstBWc0AwOqIffA7UBZ7rE1rsEIa6GjIDJxI/2siIXSoG6O7DR0JYgQ6Y6HmM9u+BqMUlpQsrKm5ONNySUnPwTItF2taLtaMAmyeGVg9M8CRg/9CiO5TF+lqbKXrfAf+tjABvwVfdyqtoXS4aIEDncAxLJF9OPHhTo2Slu8ioyQXz5I5pC+cicU2UHOgEKIvpYf6w50Ebpi+QJ84d2zIY4xIlI6jp2nb10B73SW6moN0++0ErVmXaw7awBHx4rQFcKUpUnJdpE7PJK24gIxFs3DmZo6qfFprjM5mwi0niVw6SaT1FFHvOSLt54m2nyPiPYfR2Xz1E5UFS2oO1vQ8rGk5WNLyzA/U1Fws7kyUK703HCyuDPPizkA501F2F2qoEIgTrTU61BWrDXWgY7WiaFcrRuely8HY1UK081JvSBrdrVcHjs2JzTMdm6cIq2c6tuzZ2HOLseXOBfc0uhoDdNU30XGyla4LXXR3KAJGOlF7Wu8pVDSE0/CSkhIhY0YKWUumk3PLYpy5GXH+zQgx8ZRSNVrrATfYm1LBYRgG/lNNtNXW0X7iIh3n/XR3WQmozD41CLCF23HbukjNspAx24OntAhP2Q2j7pQ1Qn4irafNYIhdwn1+1qHufsdbUrOxZk7HljnN/BDMnI41c5p5n2ca1oxCLCnZqCna3KKNKNGOZqLes0S8Z82g9TbGAvcckbZGot6zoI3e5yhnGracudjz5prX+TfgKFxExMil/cgF2o830Xmuk+528EfTiNovNxnaw62kuPxkFLrIWlRI/qplo/6iIESymNrBkX+D/tmnv0G3DwI6g6jtcqepNdKFS/lIyYT0GWl4Fk4nq3wervysUb2WEewkfOEY4aZjhJqOxn4+QqT1dL9vyMqRii1nDracOdhz5mDLndN725Y1E4tDRg1dLx0JEWk9TbilnkhLQ+y6nnBLA5FLJyEa7j3WmjULx7SF2AsX4ihchH16KeFuDy176mk/fpGOC2G6QymE7Z7YyQ2ckUukZ4TJmpdF/m0LyFheLP0mYlKZ2sGRt0D/74//Iy7dTkpalPTCFDLn55NVPo/U2YWjOqcR8BFuOk7owlHCTUdiQXGMaNvpywdZHdjzb8BeuABHwQJseSW9QWFJy02IpqKpShtRIpdOEjp/JPb+HY0F/XGIhsyDrA4c00txzFyBs2g5jpkriFJAy9t1tB48R3tThG4jC8Nm1lStkU5SbT6y57gpuHMh2TcvlP4SkdSmdHAsys/Th841jeqPWGtNtP0cwdN7CZ3ZS+jMPkJNR8xmkB42J/b8+TgKF2AvjH1rLViILXcuyipjD5KJjkYIt9QTOnvAfK8b9xE8cwAdiA04sNpxTF+Cc+7NuObejL2oAt+JDprf/QBvQwcdXW7CdrO2ao10kuZsJ2deBtNWlUqNRCSdKR0cC2cW6qNnmoZ1bLTjYm9IBM+Y11HfBfNBixV74SIc00uxFyzAEQsJW84clGVyfbPUWtPVBR0dmo4Ojc93+bqz8/LtYBBCIQiF9BXXYBixeemqZzCWWcOyWMDpBJdL4XJdfe10KjIyFB6PecnMVGRlmffZ7fGvpWmtiVxqIHRmH8Ez+wieqiF0ugYd9gNgy56Fc+4tZpgU34q/LYWmVw5x6Xgbvu50onazY90WbiczvYvC8ulM/8iN1z3sejDdl/yEfGEy56RLrVZclykdHAvmzdTHPjhz1f3RrtZ+ARE8vY+ot9F8UCkzHGauxDlzJY6ZK3HMWJr0fQ/RqKa1VXPxoqa52aC5WXPhgnnd3Gxw8aLm4kWDlhZz5vpQ7HZwucDhUDgcfa/Nn62xLO2Zbd7zs2FAMGiGjt9vXgcCmkBg6JG5AGlpkJtroaBAUVDQ99rC9OmK2bOtZGaO/4eljoYJNR4g0PAOwYZ3CDa82/sFw5KWi/uGu3EtWIVz3l10NgS58MYRLn3QiS+Sjba6UEaYVNVMbkkqRR9aQebSuWNWtt3f30fz4VZSC1IourWQGTcX4s6a+AmhIvlM6eCYv7hYH63dR/DMfkJnantrFJFLJ3uPseWV9AaEc9ZKHDOWYXGNzzfC8aS15tIlzenTBmfOmJfGxsvXFy5ootGrn5eVpcjLU+TlWcjPV+TnW/B4zG/56emK9HRIT+97W+F0ju0HtNZmbSUQAJ9P096uaWszaG/XeL2XLz1hZ16uDjiPRzFrloVZsyzMnm1e5syxcMMNVtLSxidUemolgbq3CZx4g8Dx14j6zFquLXcurvmrcM9fhX3mrVx8/QRNb5+mrdVByJ4NgCPcQk5hlJn3l5Jzx5LratJ68zu7iYYMnBkOWk94QUHeomyKbp1GwfJcrI7JVTsW42dKB0fx7FT9ysOX+xps2bMu1yRmrcRRtBxryuhGUU2UYNAMh7o6g4aGKPX1RuwSpbOz/7H5+YqZMy0UFVmYPt0Mhr4BkZs79iEQL1qbNagLFzTnzhmcOmVeTp82OHUqyrlzGuPyiFyKihTz51tZsMDK/PkWFiywUlxswWod+xAMXzhG4Nhr+I+/RuCDt9DBDrDYcJXcinvxA6SUPkBno+bcywdorgvQZckHZcUebiMnP0TRmoXk3bN8xCHyyl//geySTFb8SSldzd00vt3E2XfO428NYnPbmH5jAUW3FuKZkyFNWWJIUzw4UnTt038bC4sVWNNyJ7pIw6a15tw5zdGjUY4ejXLsmMHRo1FOnTL6fSAWFiqKiy0UF1uZO9f8tl1UZF5crqn74RAKac6eNUP12LEox4+b1w0NRm/NKzUVli61sny5jWXLrCxfbqWgYGw7sXU0QvDUHvyHX6L70EuEzx8CwJY3j5TSB3CX3o/huIHGX71H02EfnboALFZsYS85OQFmfWgxuXcvG1aIvPzlN5heUcCSP1pw+fUNzaXjbTS+fZ7ztc0YYYO0whSKbp3GjJsLcXmcQ5xRTFVTOjhKbsjVdSdaJroY1xSJaI4fNzh4MNonKKL4+iwkO3OmhYULzWaXkhILJSVW5syxkJo6dcNhNIJBTX29wZEjUQ4cMC9HjkSJxJbjKixUrFhh5aabbNxyi4158yxj+u08fOkU/sMv4z/0e/wn3oBoCEtqNilLP0Lqio9D+hLO/HIPTQfb6DDywWLDEW4hf5Zi7tobySidM+B5tda8+MXXKL5vFgs/XjLwa/sjnK+5SOMfztNW3242ZZXmMPPWaeQvy8Vql5FfwjSlg2Pewmn6g6PnJ7oY/WitaWzUHDgQ4cCBKPv3Rzl0KEogYD6emgrz51tZuNDKggUWFi2yMn/++LXRCzNMDh8234v9+6PU1kY4d87828jOVtx8s42bb7Zy66025s4duyAxgp34j75C94Hf0H3wRXSwE0tKFilLHyRl+cchfQmntu2m6XAX3TZz3lFKtIlpi1OZ88k7cE3L6T1XJBjlpb94nYUPl1By/+xrvnbnhW4a3z7P2XeaCHiD2FNiTVm3TSNzlozKmuqmdnCUztYfHDo1oWXweg3efz/Kvn2Xv+G2tpq/d6cTSkutLFt2+TJzpgWLRf5oJ1JPuL/zToTduyO8806EpibzPSsqUqxaZWfVKhs332wbs+ZAIxwgcOwVuvb9ygyRgA+LO5OUpR8lteJRgoEiTv2yhounFSF7DhgRMm0XmXlnEUWfuINwd5SdX93Fkk8tYPZdM4b/bzU0LUdaaXz7PE37WjAiBunTU3ubspwZjjH594nkMqWD44Zl8/SJAx/E9TUvXTLYvTvK7t3mh87x42aHhFJQUmJh+XJrrD3dxvz5lgmZnyBGRmtzQMJbb0V4440Ib78dwe83g/+WW2zcc4+NNWvsY9Y/oiNB/MdeNUPkwG/RwQ6snhmklq0lpWwd7cejnH7xEJfazGV0rGEfWTkRWnzZrPxcKdNvLBjV64a7wpyruUjj2+fxNvhQFkXuwiwKVuRRsDwXV6b0h0wVUzs4lt+gT+w/Ma6v0dxssHt3hHffNcOirs4MipQUKCuzcdNNVlautFFaaiU9XUJiMggGNbt3R3j11QivvRbhzBkDpaC83MoDD9i5/347hYVjEyJGqBv/wRfprHkW/5EdYESwT19KWsWjuBY9xNnfn6DxnSY6mAbKwhz/S8x4+FYy7r8fi3v0c486znXR+M55mmov0t0SAAVZczN7QyQ13y3NWZOYBMcYB0dTkxGrTZhB0dBgBkVqKlRU2LjxRjMsliyxSm1iCtBaU1dn8Pvfh/n978McO2b+/9ATIh/60NjVRKKdLXTtfYHOPc8QOrUHlMJ1w12kVfwRRlo553/zHpZXfkrk9CksqalkPPggnk88gmvZslF/yGut6TjbxYX9zTTta8Z3xhzz7c52krMwm7xF2eQsyJImrUlGguM6g8Pn0+zaFeGtt8K8+645HBbMmcwVFWY79003WVm82IrNJkEx1dXVRXtD5OhRA4sFbr/dxsMP26mstI9Zn0i4uY7OPc/StecXRC6dRDlSSVn+EKkVn0S3p9D+wi/xvfQS2u/HecM8Mh/5BJkPfRRbTs61Tz6E7hY/zYcu0XK0jZZjbUS6zeFoadNTyZqbgac4k6y5GaQVpqKkry5pSXCMMDgMQ3PkiMHrr4d5880Ie/dGiUbNoDBDwrwsWjT2k8fE5FJfH+XXvw7zq1+FOHdOk5YGH/6wnUcecVBWZh2Tph6tNcGGd+jc/XO69v0KHfBh9RSRduN6UhZ9jO49x/A+/zyB/QfAZiP9nlVkPvwwaXfcgXJcXy1BG5r20x20HGmlta4db0M74S4zSGwuKxkz08koSiO9KM28npYqs9eThATHMIKjtdVg164Ib75pXlpazN9LaamFu+6yc9ddNlaskBqFGB3D0OzeHeWFF0K89FKY7m5zoMSnPuXg4YcdY9b3ZYT8dB/8HV3v/QL/0Z3mzpWzK0i78Y+wZayg48VXaP/1r4m2tmL1eMh48EEyP/YQrqVLxyzEui/6aWtop63eh+9MBx3nuogGYzMuFaTmp5ARC5KMmelkzEyTTvcENKWDY97yefqD/VePqopGNe+/H+WNN8xRMgcORNHaXOvojjts3HWXjTvvtJGbKxOixNjq6tL8/vdh/t//C7F/fxS3Gz76UTuf+pST0tKx+zYeaW+iq3Ybnbv/nzlb3eogZckDpK58FOOSm/bf/o7Ona+gQyEcc+aQ+bGHyPjoQziKhj+Udzi0oelu8eM720lHYye+xk58ZzvxtwR6j3FmOHpDJHNmOlnzMiVMJpgERyw4WloM3nzTDIpduyK0tWmUgmXLrNx9txkUS5dapflJxM2hQ1F+/vMg//EfYfx+8//F//pfHTz4oB2HY2z+P9RaEzp7gK73fkFnzbMYnS1Y0nJJLVuHe9FHCb5/Dt9vfkP3nj0AuMvLybj/PtIrK7FPmzYmZRhIuDtshsiZTtrPdOA700Hn+W50bEn+1IIUcuZ7yJ7vIW9RDo40+7iVRVxtSgfHrEXz9GMbDvLGG2EOHTI7tXNyVG+N4o47bGRlSa1CTCyfT/OrX4X4+c9D1NUZ5OcrPv1pB5/8pAOPZ+z+/9TRMP6jO+nc/XO6D74I0RD2aYtJu/GTOKbdRdfr7+L73X8SPGE277qWLSPjvkrS77sPx6xZY1aOwUTDUXyNnbSe8PZeIoEoyqLImpdJ4fI8ClbkkpKT3FscJIMpHRwpuSV6Vm4tK1dauesuG3fdZWfRIpmZLRKT1pq33orwr/8a4q23IrjdsHatgz/+YwezZ49tp3K0u42uvS/Q9d4vCJ7cDcqCa8E9pN34KWzppXS++iYdL79M4JC5KKOjpIS0O24n9Y47SbmxAotr/Pf50Iam/VQHFw40c2F/Cx3nugDIvsFD0a2FTCvLx+aSnTbHw5QOjqyiEn3q8AdkZEhQiORy7FiUf/1XsxkrEoE1a2xs2OBk5cqx/6AMX/yAzj2/oPO9XxBtO4NyppOy5AFSln4EW8ZiOt/4A11vvEn3nj3oUAjldJJSUUHqrbfgLivHtaQUy3WO0BqOruZuzu+5SOM7TXRd6MbqsDCtPJ85q4rInJ0x7q8/lUzp4CgoLtEX6usmuhhCjNrFiwb//u9mM1Z7u+bmm6184QsubrttbIbz9qUNg0DdLrr2/ILugy9idF0CmxP3gntIWfogrpJ7CB49Sddbb9H51i5CdebflnI4cC1dSkrZStwry3CvXIEta/z2udFa423wmYs07r5ANBjFMzeDOauKKCzLl1V+x0DSBIdSai3gBTxAvda6doBj1sQeB/BqrXcMdc6iBSW68ZgEh0h+3d2aZ54J8eMfB7lwQbNsmZUvfMHJvffaxqXpVUcjBBveofv939F14LdE204D4ChajmvBPbgX3IM1Yz6B9w/hr92Lv7YW/+HD9GzL6Jg9G/eKFbhXrsS9cgXOefNQ1rGfwxH2R2h85zynXjtL14VuHOl2Zt0xnVl3zsCdLdvmjlZSBIdSygNs1lpXxW5v11pXDnDMGq31c7Hbm7TWTw113rlL5umGg/Fd5FCI8RQMan75yzBbtwY5c8Zg/nwLVVVOPvxh+7jNMzJHZr2P//BL+I+9SrDhXTAiKLsLZ/FtuEpuwzn3Fuz5iwkdr6d73z78+/bj37eP6KVLAFhSU3EtW4p7xQpSVq7EvXw51szMsSujoWk52sap1xq58L65B0/+khxm3TmD/CU5Mot9hJIlODYCnp4gUEptA57sW+uIBcdOYJ3Wul4ptVlrXT3UeQebxyFEsotENL/7XZgtW4KcOGEwa5aFjRudfPzj9nHfDtgIdhL4YBf+Y68SOP4a4aYj5gMWG44ZS3HOvRnX3JtxzLkJ3WHg37cP/759dO/bR/DYcXq2YHQUF+NeuYK0O+4g9c47saaljUn5ulv8nNl1jjO7zhP0hXBlOZl52zSmlReQPj11TF5jskuW4NiE2fS0NXZ7C7DtyqaoWFPVNmCH1nrdIOfaCGwEKCguKG+qaxrXsgsxkQxDs3NnhB/+MMj770cpKFB87nNOHn3UEbfdIaPdbQRPvkew4V3zcmoPOuwHwOopwjn3JlxzbsI59yZsnnkEDh/tDRP/3r1E29vBbif15ptJv6+SjA99CGt6+nWXy4gaXNjfwuk3z9JytA00pBamULg8j5z5HrJKMmVU1iCSOTi29zRL9TluI7AD2IxZQ6m86mR9SI1DTBVaa/7whyg//GGAd9+N4vEoPvMZB5/+9NjOBRlWWaJhQmffJ9iwm+DJ3QRO7ibadgYAZXfhmFnWGyaOWRWEPmikY+crdOzcQfjUaZTLRfp9lXg+sZaUm24ck0EAgfYgTXubadrbTOsJrznRUEH69FTSClNJK0whJdeNM9OBM9OJK9OBPdU+ZZeOT5bgGE5T1Vr6dJorpTZjhsugHeQSHGIq2rs3wpYtQXbujJCSAp/8pIM/+RPnmO0RMhoR71mCJ3cTbHiPwMl3CTXuh2gYlMIxs4yUxffjWlyJ4bXR/stf4vvt7zA6OnAtXkz2Zz9LxgP3o2xjUzuIBCK01fto/cBrzlhv6qa7xQ9XfBxabAqXx4kry4Ury4k720X69FQyZqaTVpAyqftNkiU4PFy7c3wjZhNVfez2GswgqR/svCXLS3TdfhlVJaam48ejbN0a5Le/DWOxwMMP29mwwcmcORO/Qq0RDhA6s4/AiTfoPvwSodM1oDXWrJmkVTyKe+nD+N89Suu//huhhgbsM2aQ++d/TuZDHx2X0VnRcJRAW5Bge4hA++XrQFsQf1uAQGuQgDfYuySK1WEha56H3IXZ5C3OIn1G2qSqnSRFcEC/4bgA9NQklFI1wGqttbenSSt2SP21huNKcAgBZ84Y/PjHQZ57LkQoBA88YGfjRidLlkx8gPSIdrbgP/wyXXtfuLyy76xy0u/YiG7P5tKPthI4dAjHvBLyH3+ctHvvjfsHtRE16Dzfje9MB96TPi4da6OzqRsAZ6aDvNIc8ktzyF2cjd2d3H0nSRMc40GCQ4jLWloMfvKTED/7WZDOTli50sqnP+3g/vvHblHFsRDxXaCr9jk6//ATwhePY82cRtrtf4olWELLD7YSOnkS98qV5H/lK6SUrZzQsvrbArQcaaP50CWaD7cS8UfMtbVKMslfkkPekhzSp6cmXW1EgkOCQ4h+Ojo0zz8f4t//PcSpUwZ5eYo/+iNzUcW8vMSZda0NA//Rnfhe/wGBY69iSckiY9WXMFryafnhFqLNLaTfdx/5X/4fOGbPnujiYkQNvA0+Lh68RPPBS/gazW12XVlOM0RKc8hdmJUUI7mmdnAsK9F1ByQ4hBiIYWjefDPCT38a4vXXI1itsGqVjUcecbBqlS2haiHBM3vx/ud38B/ZjiU9n4y7Hyd8JMqlf/sJOhQi65OfJPeLfzauS52MVKAtyMXDZoi0HGklEoj2juTyzMnAMyeDzDkZpE9PxWJNnMAGCQ4JDiGG4eTJKNu2mdvcXryoycpSfOxjdj72MQelpZaEaWoJ1L9D2+++SbBuF/Zpi8m4q5rOF/fife45LCkp5FRtJPvTn47L6r0jYUQMWuvaaT3ehvekD+9JX+82uxa7hczZ6b1h4pmTgTvHNaG/cwkOCQ4hhi0SMZd2f/75EDt3RgiHoahIcd99du67z87KldYJ35ZAa033gf+g7Td/TeTSKdxLPkzaso20/uuzdL72GrZp08h/7C/I+OhHUZbE+ibfQ2tzZ0TvSR/eBjNIfKc7MSLmvkGOdHu/WolndkZcN7OS4JDgEGJUvF6DHTsivPRSmF27zBDJy1Pceae5Edrtt0/sRmhGOIDv9f9D+/b/hY4EybizCrtnFS3/tIXAoUM4Fy+i4CtfIfXWWyesjCNhRA06GjvxnurorZV0nu/qnV+SkufuFyaZM9OwOsZnZJwEhwSHENeto0Pz6qthduyI8Ic/RGhvN7deXrrUyh13mLtprlhhxW6Pf20k0t6E9z+/Refun2FxZ5Kx+i+hs4jm//1PRM6dx7VkCZ71j5L54Q9jSU2utaoigQjtfYLEe9JHoC0IgLIo0mdc7i/JmJlu9pfYrj/MJTgkOIQYU9Go5sCBKG+9FeGttyLs2xfFMCA1FW66yayJ3HabjXnz4ts3Ejr7Pq2/+TqBY69gy5lN5r2bCNcZeLdtI3jiA5TbTertt5F+72pSb78de0F+3Mo2lgLeYG+ItJ/04T3VQcRv9pcoq+qd3Z4xM43Mmemkz0gb8bwSCQ4JDiHGlc+neeedCLt2mZdTp8x2+oICxW23mUFy++02cnPj06zlP7qT1t98nfC5g1izZpJx959hcy2n4+WddLzyKpEmc+FT+4wZuMvMjafcy5fjmj8fZY9fP8JY0Yam62I3vsZOfGc6aD9jXoc6wr3HODMcpBakmJf8FNIK3KQWmOtzDVRDmdLBUbysWNcfGHRFEiHEOGhsNHpD5O23I3i95udMaamFD3/YwUc+Ymf69PENEW0Y+A+/RPvO7xFseAdld5Oy9EFSytaiIvkE9h+kO7YBVaS5GQDldOJasgT3smW4ly/HvWI59sLCcS3neNFaE2wP0X6mg46znXRd9NN1oZuui939AgUF7iwn7hw3KbkuUnLduHPdzLxlmgSHEGJiGIbm8GGDXbvCbN8eYf9+cy+OsjIr69aZIeJyjW9zVvDkHjp3/4yufS9gdHvBasc19xacxbdin7EMq2MGoYYmAvsP4N+/n8Dhw+hQCAB7URHpq+8lbfVqUsrKxmyhxYkU7grTddFP58Vuui5047/kp7slQHeLn2C7+e/+yJbVEhxCiMRw+rTBf/5niF//OswHHxh4PIp16+x8+tNOpk0b51pIJEigLrYB1bFXCZ07BNpsVrOk5eKcVYZjZhmOacsgkE7w+Cm6du2i6+230aEQ1qws0tesIf2B+0m9+eZJESJXioai+FsDpE9Lk+AQQiQWrTW7d0f5v/83yI4dEWw2WL/ewec/7yQ/Pz59IUaom9C5Q4TO7CN0Zh/BM7Xmboaxz0V7/nzcpffhnHMXkXNROne8Querr2J0d2P1eEivrJy0ISJ9HBIcQiS0M2cMfvjDAC+8EMZmgz/+Yydf+IKTtLT4D+01gp2EGvcTPLXHrJl8sAuiISypOaSufAT30oeJnO6m46WX6XzllT4hsob0+x8g9eabkrKD/UoSHBIcQiSFU6ei/PM/B/nVr8Lk5ys2bXLx0EMTuwufEezEf+xVumqfx3/oRXQ4gC13Lqlla3Ev/gjhD1rw/f6lfiGSdvfdpNx0Iyk33oh95syEWa5lJCQ4JDiESCr79kX41rcCHDgQpaLCyne+46a4eOL3DjECProP/JbOPc8QOPEGaAP7tMVmTWTxRwgdO4fvxd/TtWsX0bY2AGwFBaTceGPsUoFj7tykCBIJDgkOIZKOYWiefz7M5s0B/H7Nl77k4nOfc2CzJcaHbrTjIl37fk3X3ucJ1r8NgC2vBPfCNbgW3ItFTSdw4BDd771H13vvEW1uAcDq8eBeuRL3ypWklK3EtWRJwi3ICBIcEhxCJLHmZoNvfMPPyy9HKC218OSTKSxaNPG1j74ibY10H/gt/mM7CZx4Ex32g9WBc3YFrnl34Cy5HYsqxL//ffx79+Hfu5dQQ4P5ZJsN1+LFuFcsxzlvHs7iYhwlJRO+PLwEhwSHEEnvxRfD/O3f+mlv12zY4OSLX3TidCZG7aMvIxwgWP82/mOvEPjgLUJn9plDfq0OnLPLcc27A9e8O7Bm3kDw0HH8e/eal4MH0YFA73ksGRnYCwuxTSvEXjgN+7RCbIWF5n05OVhzc7F6POPW7CXBIcEhxKTQ1mbw5JMBfvnLMLNnW/jWt9zcemtiD4M1Aj4C9e8Q+OAtM0ga94ERBasd5ywzSJxzb8YxYxlGR5RQfR3B+nrCp08TPt9EuKmJyPnzRL3eq09ut2PLzsaWm4stNxdrTg6W1FQsLhfK7cLicqOcDvNYjTnMWBugNUYwhA4EMAIB8zoYQPsDGMEgRncXc37yEwkOIcTksWtXhK9/3c/p0wYPP2znq191kZ2dmPtuXMkMkncJ1MWC5MxeM0gAa+Y0HEXLcUxbjC23GHvuXGy5xVgzCtGhEJGmJsIXLhK91EKkpYVIy6XYdTPRlktELl3C8PvRfj86HL5GSUzK5cLidJrXLpd57XYz95lfTN3gmLtsrm440DDRxRBCjLFAQPN//k+Qp58OkpICVVUuPvMZx7gvXzLWjEAHobMHCJ7ZT6hxP6HGfYQvfgBGpPcYZXdj9czAljUDm6cIa1YRNs8MbJ4ZvT9bXOn9zqsjEYxAEB2MNX8pBUqZTVtKoZxOlMMx6EZXU7qpSoJDiORhaIPmYDNn/Gc46z9Lc7CZjkgHnZFODAxsyobT4iTHkUOeM48idxGWC8V8/++svP56hIICxec/7+Thhx2kpiZXgPSloxEibWeIXGog3FxPpKWBiLeRaNtZ89rX1Du7vYfFnYnVU4Q1sxCLKwOLK733WjlTURYbWGxgsVz+GcyAMqJoHQXDQBtR0FE8935JgkMIkXh8YR+HfIc45DvEQd9BjnQcoTva3fu4QpFmSyPNloZVWQkbYYJGkPZwO5rLn11F7iJmhJZz7HfLOfbiUtKt6axb5+CRRxzMn59YI7DGgo6GibafJ+I9S6StkWhbo/mz9yxRXxNGoAMd6DCvQ12jeo25/9snwSGEmFhRHeVk18nekDjkO8Rp/2kArFgpTitmScYSSlJLKHIXUeQuIs+Zh0Vd3ZQSNsK0hFo41XWK453HOdJxhL3evXRFu1BYSL24lPMv30H3ntu4YXoWH/qQnbvvtlFaOvH7pcebjkbQ4W50NALaiF1HzWtAWaxgsZrXyhqrkVixujMlOIQQ8aO15lLoEsc7j3PYd/iq2kSmPZMlGUsozShlScYSFqQvIMWacl2vGdERjviO8E7rO7za/Cpn/GdQ2oK9cTnNr9xOYO+teOwe7rjDxt13m3umJ0uH+kSQPg4JDiHGlKEN/FE/HZEOLgYv0hxs5mLwIucC52joaqChqwFfxAdcrk2UZpSal/RSitxF47rshtaauq46Xm1+lVeaX6HR34jSFtIvLaPtzdu49Oat6I4sli+3smqVjVWr7CxeHN9tbhOdBIcEh5gi3mt9j1eaX6Ej0oEv7MMX8REyQkR0hLARJqIjRIwIYR1Ga41FWS5fuPyz6vtf7MPUgoWQDuGP+PEb/gFfP82axpzUORSnFjM3dS4lqSUsTF+I2+qO56+hn54Qea35NV5tfpXT/tMoFFndCwgcXkLjW0sI1y1kVn4GH/uYnYcesjN79uTrFxkpCQ4JDjFFPLb/MQ60H6DIXUSGPYN0WzpOixObxYZN2bBb7NiVHauyYlEWDG0Q1VE0uvdnAwOtdW/ns0b33rYpGynWFFJsKaRaU0m1pZLvzCffmU+eM480W9oE/waGprXmZPdJXmt+jT1tezjccZiINtv6rZ15dNfPJnx2NkWO2axeNpdH7y1henb6Nc46OU3p4JizbI4+eeDkRBdDiLj40r4vYWDwzyv+eaKLkhQC0QAHfQc52nGU+q56jnnraAyexlCX51DY/TkUuWaxbNocitPmUJpRyvy0+ZO+WWuo4EjsufpCiBHpO0RVXJvL6qIiq4KKrMufjxEd4Vz3Od481sDO/Q0cbT3Jsewz1Ad/j3KZTXQ5jlzuyr2Th6Y9xLy0eRNV/AkjwSHEJKLRKCb3N+HxZlM2ZqXO4r+UzeK/lN1NKKR5880IL/wyyOu1F1Dz9hMuf5dfLf0tvzz3SxY4l/L5BZ/tFz6TnQSHEJOMBMfYcjgUq1fbWb3aTkdHKm+/PYs33vgQr29uxVeyncOVv+bx4ONkti5jrfvPeeimBZN+mK8EhxCTiNYaNcUmuMVTerrivvvs3HefHa2n09DwaV7f9QlebPgPzs7/Of/i/Dzf/8ePUHTsM9xRkcFtt9moqLAl3fpZ1yLBIcQkIn0c8aOUorjYSnFxOn/Cp2gLfITNtVvYtfo/8N36Fj//+Uaefvp2HA5FebmV22+3cfvt5nyRZJ+9nlCjqpRSawEv4AHqtda1gxy3EagHPFrr54Y6p4yqElPJn+39MxwWB99b/r2JLsqUddh3mL8/8fec6DzBAm5hes0XqX0lh2PHDACyshRlZVZKS60sWWJe5+WphBullRSjqpRSHqBSa10Vu70dqBzguG3ABq21N3bMkMEhX8DEVCKd4xNvccZitpZt5bnG5/jxyR9z6sYNfG7d57jb8Qjvvg1/+EOEffui7Nx5echvejrMnWtlzhxL72XGDAsFBRby8xV2e2K9pwkTHMCjQF2f216lVFnfWodSqqzPz8Va66uCRYipLtG+uU5FNmXjkzM/yaq8VfzjiX/kB/U/4OW0l/nKvV/h4x9fBEBHh+bIkSiHD0c5edKgocFgz54Iv/lN/2+7SkFurqKw0EJ2tsLjuXzJyur52UJqKjidCpfLvHY6zecaBkSjl6+jUU0gAH5//+tAQPe7b8h/37j95kbOg9lM1aMVyL7imJ5qUzbgUUpt1lpXj3/RhEgOidT0LKDQVch3l3yX11te53sffI+qvVU8MuMRNszZQHp6KjfdZOOmm/p/DAcCmtOnDZqaDJqaNBcuXP750iVNXV2UtjZN1+hWSx8TiRQcA/EMdFtrXQ+glNocq3n02xs21geyESCvOG/8SylEgpDO8cSjlGJV3ioqsip4uuFpXjj7Am80v8FfzPsL7sq966oaosulmD/fes19REIhjc+naWvTeL2a7m5NMGgGTzAIwaDGMMBqrpSOxaKwWMBuN1/D5QK3u++1wu2+/JjDMfhrJ1JweOkfFNmYHeB91QM5Vzyn+MrjtNZbga0Ac5bOkb8kMWVIH0fiSrOl8fgNj3NfwX38/fG/568P/zW3Zd/GYzc8xjTXtBGfz+FQ5OYqcnPHobDXkEizVJ4FSvrc9gwwqmoHZlD0KAb2DHVS+QYmphoJjsRWmlHK0+VP88XiL7LXu5fPvPcZfnb6Z0T67DGe6BImOLTWXmC7UmqNUmoNsLnnMaVUjVLKEzvmGaXUxlhz1JOx+4QQxGoc0jme8Ho6z39640+pyKrgRw0/4k9r/5T6risbWRJTIjVVMdicDK11+bWOEUJI53iyKXAV8OSSJ3mz5U3+7vjfsaF2A39e/Od8fPrHE/oLQMLUOIQQY0OaqpLPnbl38pOKn1DmKeMfPvgHvnboa3RGOie6WIMaVnAopZ6JXW9QSv2lUmpO7HLv+BZPCDES0jmevLId2Ty15Cn+e8l/5+3Wt6mqreJU96mJLtaAhhUcWuv1sR9bgZ3AVzFHLa0bp3IJIUZB+jiSm1KKR4se5XvLvocv4mNj7UZ2teya6GJdZaRNVbWA1lp/Xmt9n9b6C+NRKCGEmMpWeFbwL2X/wkz3TL566Kv89PRPE6r/akTBobVu0FrvG6eyCCGuk9bSVDVZFLgK+MGKH7A6bzVbG7ay+fhmwkZ4oosFXCM4lFI/vOJ2plLqE0qpFeNaKiHEqMi8pcnFaXXyN4v+hv82+7/xu6bf8Zfv/yUd4Y6JLtY1axxKKZXR57bWWj9P/0l4QogEIZ3jk49Sis/N+Rx/tfCvONB+gM/v/Txn/WcntEzXCo5txBYaVEr9CKhVSp3g6sUHhRAJQjrHJ6cHCh7gH5f9I96wl6raKg60H5iwsgwZHFrrnUDP5Ls6rfU8rfUNWut/Gf+ijQ2puoupRmock9cKzwp+tPJHpNvTeWz/Y2y/sH1CynGtPo4VWuvnlVJ/CkxMCYUQw5ZII2/E+JiZMpMtK7dQmlHKN49+k5+d/lncy3CtpqqvKaVaMdeNWi8T/oRIbNLHMTVk2DP4h2X/wOq81fyo4Uf8oO4Hcf3ScK2mqke11tnAGszJf59XSrUqpV6KS+mEECMiTbNTh91i5+uLvs4j0x/hF42/4LvHv0tEx2eF3WEtcqi13gvs7bmtlMoc6QsppT4B1GitT470uUKI4ZPO8anDoiw8Nu8xMu2Z/Nupf8MX9vGNxd/AaXGO7+uO5kla6/ZRPOd5oEQp9axS6sujeV0hxNBkAuDUo5Tis3M+y1/M+wveuvQWTxx8gkA0MK6vOaYTAK81MTA2SuurwFMjK6YQYjikqWrqWjtjLU8seII9bXuoPliNP+oft9ca6wmA6we5v+8J6oHnh1k+IcQISOf41Pbhwg/zVwv/in3efXzl/a/QHekel9cZ6wmAVUqp/6mUuveKwLnSeyMvqhBiOKSPY2q7v+B+/r9F/x8H2w/y5fe/TFeka8xfY6wnANYD84DngDal1Aml1A+VUo9cESRt111yIcSApMYh1uSv4W8W/w1HOo7wPw78DzoiY7u+1VhPAHyyzxDeeZj9GQqzT6MnSJ5B9vEQYlzIBEDR4568e/jW4m9xvPM4j+9/HF/YN2bnHtMJgLH+j56fG7TWz8f27piH2bz1eeAkUHF9xRZCDEQ6x0Vfd+beybdLv019Vz2PHXhszMIjbhMAtdbtWuudWutq4LujK+7IyR+SmEqkc1xc6fac23lyyZOc6jrFYwceG5Nl2Ye7dexerfXf9QmSR6/zdbde5/OFEIOQznFxpZuzb+Y7pd/hZNdJHj/w+HX3eYx6AqBSakWs4/ul2PXnRvL80byuEGJoMgFQDOaWnFv4dum3qeuq48sHvkxnpHPU5xpVcMSWD6kFSoAcoAp4Wil1SSn18KhLI4S4LtI0K4ZyW85tfHvxtznReYIvHxj9UN1RBQewRmtt0Vrfp7Wu0FpbMIft/gvwY1lSRIiJIX0c4lpuz72dby7+Jsc6j/GX7//lqCYJjjY46q+8I9YPUo05q/x+WYJdiIkhfRziWu7MvZNvLv4mRzqOjCo8RhscOYPNDNdae7XW9wFlozy3EOI6SI1DDMdduXfxjUXf4LDvMJsObqI7OvzwGG1wPAnsVEotH+IY7yjPLYQYJZkAKEZiVd4q/mbx33Cw/eCIVtW9nmXVq4BXlVLvKaU+N0ANJGs0557qsrKyKC8vp7KysvdSW1s7onN4vV7Ky82VYurr6ykpKRnwuMHuF8lLOsfFSN2Tdw9/tfCv2Ovdy9cOfY2gEbzmc4a1kdNAtNa1SqlizFnlTwNblVL1mKOtMjFniU+8JPw72rlzJx6PZ0zOVVxcTF1d3ZicSyQ+6RwXo1FZUElIh/juse/y9UNf59ul3x7y+FEHB5j9GZg1jyql1BrMGeZrMPs36mJBsgN4GdiptR67xVKEEFfRaOkcF6PyYOGDhI0w/+vE/+Jvj/ztkMdeV3D0pbXegRkSPVvL3oi5mGElZrhopVSl1vqVsXrNsfTtb/s5ciQ6rq+xaJGVv/5r96ie6/V6Wb16NTU1NYDZzNRTk9i6dStbtmwBYP369WzcuHHQ51VVVbFnzx4qKipobW29nn+OSERaOsfF6H18+scJG2G+X/f9IY8bs+DoK9YHcmWQVAA14/F6Q5YlCduq1q3rv3jw9u2DL0xcW1vLli1beoPB6/UOeuyOHTuor6/vPfbZZ5+9/sKKhCJNVeJ6rStah0bzFm8Nesy4BMeVYkGyMx6vNVqjrQmMh23btg27j+OZZ55h/frLGy96PJ5Bw2P79u39Qik7e7D9uESySsYvSiLxPFr0KOuH2NB1tMNxhRAJSvo4xHhLqOBQSq1VSq2JXQ85gbDn2HiVLRH0rUnU15uT99evX88zzzxz1f0DqaysZNu2bb3nkj6OyUdqHCIe4tJUNRxKKQ9QqbWuit3ejtmxPtixVZhDgSed8vLyfk1VTzzxBGvXrqW4uJjKykrKysooKzNztaysjKqqqt55G2vWrOGJJ54Y8Lxr1qxh27ZtlJeXU1FRQXFx8bj/W0ScSee4iAOVKDNNlVIbAY/W+qnY7W2YW9FeNfstVtMoA2pjo7kGVbSkSDcebByPIguRcB7c9SBr8tfw+A2PT3RRRJJTStVorQfcrTWRmqo89F+mpBVzu9l+Yk1Ye4Y6kVJqo1Jqj1JqTyAwvCn0QkwG0lQl4iGRgmMgnoHujE08HJTWemtsufcKl8s1HuUSIiHJBEARDwnTx4FZ2/D0uZ3NFcu3x5qzWmNLndyIuUpvvdZ68B5hIaYS6eMQcZBIwfEs/Tu7PVf2b2ite/cqV0rdCGyX0BDiMpkAKOIhYYJDa+1VSm3vM8S2N0SUUjXA6p4mqlg/xxrAc60ah7T5iqlE/n8X8ZAwwQGgtX5ukPvLr7hdi7lVrRDiClLjEOMt0TvHhRAjIDUOEQ8SHEJMJlqWHBHjT4JDiElEOsdFPEhwJJisrKwBV7d97rnnKC8v77303U52sMeeeuqp3u1nS0pKKCkp6b29devWq15DJD9pqhLxkFCd42Jg9fX1VFdX927c1HeBwqEe27RpE5s2bQKgurqanJyc3tticpIah4gHCY4kUF9f37uoIZh7bvQsgjjUY2IKkj4OEQcSHDGXXvgqobPvj+trOGYsJeeR7474eRUVFaxbt47q6mqqqqr6rWo71GNi6pEah4gH6eNIAh6Ph5qaGurr63v7KXr6QYZ6TEw90sch4kFqHDGjqQnEU3Fxce8mTFVVVVRXV7Nly5ZrPiaEEGNNahxJqLq6mj17Bl5ZfqjHxOQnTVUiHiQ4ksCOHTvYsePyflVbtmxhzZo113xMTEHSOS7iYNIHRzK2+a5evbp3vkVPh/e2bdsoKSmhvLwcr9fL5s3mGpBDPSamHqlxiHhImK1jx8v0JdP1uYPnJroYQsTFXa/fxWdmfYY/nfunE10UkeSSZetYIcR1SsYatkg+EhxCTDLSxyHGmwSHEJNET7Oz9HGI8SbBIcQkIc1UIl4kOISYZKTGIcabBIcQk4TUOES8TP7gkL8lMUX0BId0jovxNvmDQ4ipIvYlSZqqxHib9MEh1XcxVfTWOCQ4xDib9MGRrGpra6msrKS8vJzKykqeeuopoP/WsllZWVRWVrJu3TrKy8uprq6ewBKLiSZfkkS8yLLqCai2tpZ169ZRU1PTb6e/gWzbtq33mOrqaqqrq2WtqilK+jhEvEiNIwFt2LCBLVu29NsCdjg7+61fv37QgBFTgFQ4RJxIjSPm+x98nxOdJ8b1NW5Iu4EvzfvSNY+rr68f8dLotbW1bNiwgaeffnq0xROThPRxiPEmwZFg6uvryc7OHvbx69at633e2rVrKSsrG6+iiQQnneMiXiQ4YoZTE4iH4uJiWltbh3183z6OrVu3sm7dut5tZMXUIp3jIl6kjyMBFRcX99vVb7g2btw4queJyUE6x0W8SHAkoKeffpqqqqreYbdg9mFcy44dO0bUzCUmGZkAKOJEmqoSUFlZGdu2bWPdunW9zVbr168fsP9iw4YNAL0hs3379riVUyQW6eMQ8SLBkaDKysoGDIG2trYBfxZC+jhEvEhTlRCThASHiJdJHxzyxySmjJ4+DukcF+MsoZqqlFJrAS/gAeq11lf1CMeOyQbKgW1aaxlGJEQf0schxlvCBIdSygNUaq2rYre3A5VXHFOGGSjPxW63AVlxLqoQCUk6x0W8JFJT1aNAXZ/b3lhQ9JUNVPW53TrAMUJMSRIcIl4SpsaB2Tzl7XO7FTMoesWapfo2TWUP1JwlxFSkL0/kEGJcJVKNYyCewR5QSm0GVg/y2Eal1B6l1J5QKDReZRMioWgtNQ4RH4kUHF76B0U2MOAa4bEO8mcGq21orbdqrSu01hUOh2OsyylEQpPgEOMtkYLjWaCkz23PIKOqyoBarXWtUqpYKXXtjSqEmAKkj0PES8IEh9baC2xXSq1RSq0BerexU0rVKKU8sdDYGTuuDqjRWk+qnYuysrJ6t4vtuQxnnaq+vF4v5eXlgLnceklJyYDHDXb/UGXru35Wj+eee47y8vLeS9/yDvbYU0891fvvKykpoaSkpPf21q1bR1QuYZI5SyJeEqlznJ5htgPcXx77sZYpMPx2586d/Xb/ux7FxcXU1dVd+8BRqq+vp7q6uvc1vF5v7/paQz22adMmNm3aBJhb3ubk5PTeFqPT08chFQ4x3hKmxjFe5FvYyI3kG399fX2/xRc9Hk/vNrdDPSbGjzRVifGWUDWOidT0P/8nwSNHx/U1nIsWUvi1r43quV6vl9WrV1NTUwOYzUw93+S3bt3Kli1bAHMV3Y0bNw76vKqqKvbs2UNFRcWgG0YN1Bw1mIqKCtatW0d1dTVVVVX9gmGox8T4keAQ423S1ziS0bp16/r1cQyltraWLVu2UFNTQ01NTb/QuNKOHTuor6+npqamN2iul8fjoaamprcvpbKysjd4hnpMjD3pHBfxIjWOmNHWBMZD3+1gr+WZZ55h/fr1vbc9Hs+gH87bt2/v3aMcGLNNn4qLi3u3q62qqqK6uro3mIZ6TIwtmQAo4kWCQ+D1evv1a1y5D8jatWuH3cxUXV3dL5yG+5i4fjIBUMSLBEcS6VuTqK83RyGvX7+eDRs29I5Iqq+vH7QmUVlZyebNm9m4cWO/EU4ej+eqEU3DHeHUs8f5mjVrANiyZUvvz0M9JsaPBIcYbxIcCai8vLxfU9UTTzzR+62/srKSsrKy3tFKZWVlVFVV9c7bWLNmDU888cSA512zZg3btm2jvLycioqKUXVWr169ujeYiouLqa6uZvPmzVRVVeHxeKioqOjXTDXYY2LsyQhCES+qd+z3JFVQWqAvHLow0cVIKlu3bh2yk10kpvOB8zz67qN8df5XeXDagxNdHJHklFI1WuuKgR6TUVXiKhIaSU5aqsQ4k+AQYpKQznERLxIcQkwyEhxivE3+4JjcXThC9JIJgCJeJn9wCDFFyARAES+TPjhkiKKYKqSPQ8TLpA8OIaYK+ZIk4kWCQ4hJRmocYrxJcAgxSUjnuIgXCY4EI1vHjnzr2NraWiorK3t/b0899dRV5c3KyqKyspJ169ZRXl5OdXX1iF4jmSglwSHGmdZ6Ul/yFuXpZOLxeHRbW9t1naOtrU2XlZVd87ji4uIRnXegstXV1fU7T1tbm66rq7vmY31t2rRJb968eURl6VFTU6OLi4v7lavnNfqW98qyb9q0SW/atGlUr5moGjob9B2v3aG3X9g+0UURkwCwRw/yuSo1DnGVRNw6drAybdiwgS1btvRbFHI4r7F+/freFYYnC2mqEvEiq+PGHHr2OL4zneP6Ghkz0yh9dP6onjvVt44drEz19fUjXqq9traWDRs28PTTT4+6PIlMgkOMNwmOBHTlZkdXbqzUV9+tY2HoD/2+W8cCPPvss9dd1p7tYaurqykpKeldut3j8Qz52FgYau+RgfT8Xuvr61m7dm2/2tBk0FvjkD4OMc4kOGJGWxMYD7J17PBfd7Ba00D6/l63bt3KunXress2GehJvkWCSBwSHCIht44dbpmKi4vZsWPHiJurNm7cOOlGVskEQBEvEhxJZCptHTvcMj399NOsW7eOmpqa3tpEbW3tNZuhduzYMWY1rkQjfRxivE364EjGb2GydezwlZWVsW3bNtatW9cbhOvXrx8wODZs2ABcDuCh+o6SkYyqEvEy6beOzVucp5sPN090MZJKIm4dm4hlSjQnOk/w2ZrP8p3S73BX7l0TXRyR5GTrWDEiifgBnYhlSjST/UugSBwSHEJMEtJUJeJFgkOISUaCQ4w3CQ4hJolkHAgikpMEhxCTRE8fh9Q4xHiT4BBikpA9x0W8SHAIMUlI57iIFwkOISYZCQ4x3hJq5rhSai3gBTxAvdb6qq3vhnNMX9JhKKYaCQ4x3hImOJRSHqBSa10Vu70dqBzpMUJMVb2d47KsuhhnCRMcwKNAXZ/bXqVU2RU1iuEc00/ICPGlfV8a46IKkXg6o+O7EZkQPRIpODyYTVA9WoErly8dzjEopTYCPWtUBP5p5T8dGqMyAmQC7eP4vOEcd61jBnt8JPfnAi3XKMd4G+3veqzPN5Lnjfa9udbjw77/Fm6R9250z7vev73J8ncHZtlmD/roYJuRx/sCbAI29rm9BVg70mMGOO/WMS7nqM433OcN57hrHTPY4yO5nyE2qo/j/xNJ9d5dz3sj713iv3fXOmayvHfD+V0k0qgqL2aNokc2UD+KY670H9dXrDE733CfN5zjrnXMYI+P9P6Jlmzv3XCOlfcuvucby/fuWsdMlvcOrlG2hFlWPdbxvVn36fjWWg/UOT7kMWJsKKX26EGWVBaJTd675JUs713C9HForb1Kqe1KqZ4t4jb3PKaUqgFWD3WMGHNbr32ISFDy3iWvpHjvEqbGIYQQIjkkTI1DJK5YE+EaoFhr/dQEF0eMQOy9qwCKMSfM7pjYEonhuuK926G1vlZ/btwkUue4SFBaay8w5Ax9kbAexRypsxWonujCiBGpwJxysAeomuCy9CM1jikq9m1mI+CNfaj03D+iJV1E/I3kvet5XClVhoT/hBvhe7dDKVWMGRpb4l/awUmNY+q6auRGnyVddmitn0MGHySq0bx367XWUuOYeCN672LNU9uBtXEr4TBIcExRsbZu7xV3D7ikS9wKJYZlpO9d7Nvsk7Fvr2ICjeS9U0ptUkp5Ys9JqGkHEhyiLw+DL+myBrhRPnwSlocB3rvY0PUq4GkSrJ1c9PIw8N/dDqBYKbWJBKv9Sx+HuBYPQKw9NinGmItenljTh4ykSj6ePiPgEq5vSmocoi8vI1/SRSQGL/LeJSsvSfbeSXCIvp4FSvrc9sioqqQh713ySrr3TpqqpqhY23cl4FFKtWqtn5MlXZKDvHfJa7K8d7LkiBBCiBGRpiohhBAjIsEhhBBiRCQ4hBBCjIgEhxBCiBGR4BBCCDEiEhxCxJlSao1SaptSSiultiilNk50mYQYCZnHIUScxZbLLsPcGEvWjxJJR2ocQkyMSmQNKZGkJDiEmBhrMPdZECLpyMxxIeIs1kxVo7VWE10WIUZDahxCxN8aEnCpbCGGS4JDiPiT/g2R1CQ4hIi/Cgbo35BhuSJZSHAIEUex/o2+u7sJkXQkOISIrwH7N5RS25DmK5EkZAKgEHGglPIAG4EngFal1KbYQyWYYeLVWif0dqFC9JDhuEIIIUZEmqqEEEKMiASHEEKIEZHgEEIIMSISHEIIIUZEgkMIIcSISHAIIYQYEQkOIYQQIyLBIYQQYkQkOIQQQozI/w9YLhaVNLmLCgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pl.plot_1dstyle(usetex=True,xmin=2,xmax=lmax,ymin=0,ymax=1,xlab=r'$L$',ylab=r'$\\rho^{\\kappa\\kappa_i}_L$',ylabsize=22,xlog=True)\n",
    "for t, label in [('euc','Euclid'),('lss','LSST'),('cib','CIB'),('gal','Euclid+LSST'),('ext','Euclid+LSST+CIB')]:\n",
    "    plot(l, np.sqrt(rho[t]),label=label)\n",
    "legend()\n",
    "savefig('fig_rho.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0041e725-e9f2-4db1-9d07-d47c2c388a9d",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}