{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Mass Definitions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "``hmf``, as of v3.1, provides a simple way of converting between halo mass definitions, which serves two basic purposes:\n", "\n", "1. the mass of a halo in one definition can be converted to its appropriate mass under another definition\n", "2. the halo mass function can be converted between mass definitions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "By \"mass definition\" we mean the way the extent of a halo is defined. In ``hmf``, we support two main kinds of definition, which themselves can contain different models. In brief, ``hmf`` supports Friends-of-Friends (FoF) halos, which are defined via their linking length, $b$, and spherical-overdensity (SO) halos, which are defined by two criteria: an overdensity $\\Delta_h$, and an indicator of what that overdensity is with respect to (usually mean background density $\\rho_m$, or critical density $\\rho_c$). In addition to being able to provide a precise overdensity for a SO halo, ``hmf`` provides a way to use the so-called \"virial\" overdensity, as defined in Bryan and Norman (1998). \n", "\n", "Converting between SO mass definitions is relatively simple: given a halo profile and concentration for the given halo mass, determine the concentration required to make that profile contain the desired density, and then compute the mass of the halo under such a concentration. \n", "\n", "There is no clear way to perform such a conversion for FoF halos. Nevertheless, if one assumes that all linked particles are the same mass, and that halos are spherical and singular isothermal spheres (cf. White (2001)), one can approximate an FoF halo by an SO halo of density $9 \\rho_m /(2\\pi b^3)$. ``hmf`` will make this approximation if a conversion between mass definitions is desired." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Changing Mass Definitions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Most of the functionality concerning mass definitions is defined in the ``hmf.halos.mass_definitions`` module: " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "init_cell": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using hmf version v3.3.4.dev14+g9a9b63c.d20210610\n" ] } ], "source": [ "%load_ext autoreload\n", "%autoreload 2\n", "import hmf\n", "from hmf.halos import mass_definitions as md\n", "\n", "print(f\"Using hmf version v{hmf.__version__}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While we're at it, import matplotlib and co:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "init_cell": true }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "%matplotlib inline\n", "import inspect" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Different mass definitions exist inside the ``mass_definitions`` module as ``Components``. All definitions are subclassed from the abstract ``MassDefinition`` class:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "[hmf.halos.mass_definitions.FOF,\n", " hmf.halos.mass_definitions.MassDefinition,\n", " hmf.halos.mass_definitions.SOCritical,\n", " hmf.halos.mass_definitions.SOGeneric,\n", " hmf.halos.mass_definitions.SOMean,\n", " hmf.halos.mass_definitions.SOVirial,\n", " hmf.halos.mass_definitions.SphericalOverdensity]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "[x[1] for x in inspect.getmembers(md, inspect.isclass) if issubclass(x[1], md.MassDefinition)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The Mass Definition Component" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "To create an instance of any class, optional ``cosmo`` and ``z`` arguments can be specified. By default, these are the ``Planck18`` cosmology at redshift 0. We'll leave them as default for this example. Let's define two mass definitions, both spherical-overdensity definitions with respect to the mean background density:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "mdef_1 = md.SOMean(overdensity=200)\n", "mdef_2 = md.SOMean(overdensity=500)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each mass definition has its own ``model_parameters``, which define the exact overdensity. For both the ``SOMean`` and ``SOCritical`` definitions, the ``overdensity`` can be provided, as above (default is 200). This must be passed as a named argument. For the ``FOF`` definition, the ``linking_length`` can be passed (default 0.2). For the ``SOVirial``, no parameters are available. Available parameters and their defaults can be checked the same way as *any* ``Component`` within ``hmf``:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'overdensity': 200}" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "md.SOMean._defaults" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The explicit halo density for a given mass definition can be accessed:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(17068502575484.854, 200.0)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "mdef_1.halo_density(), mdef_1.halo_density() / mdef_1.mean_density()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Converting Masses" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To convert a mass in one definition to a mass in another, use the following:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mass in new definition is 0.76 x 10^12 Msun / h\n", "Radius of halo in new definition is 0.16 Mpc/h\n", "Concentration of halo in new definition is 4.81\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/home/steven/miniconda3/envs/hmf/lib/python3.7/site-packages/halomod/halo_exclusion.py:17: UserWarning: Warning: Some Halo-Exclusion models have significant speedup when using Numba\n", " \"Warning: Some Halo-Exclusion models have significant speedup when using Numba\"\n" ] } ], "source": [ "mnew, rnew, cnew = mdef_1.change_definition(m=1e12, mdef=mdef_2)\n", "print(\"Mass in new definition is %.2f x 10^12 Msun / h\" % (mnew / 1e12))\n", "print(f\"Radius of halo in new definition is {rnew:.2f} Mpc/h\")\n", "print(f\"Concentration of halo in new definition is {cnew:.2f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The input mass argument can be a list or array of masses also. To convert between masses, the concentration of the input halos, and their density profile, must be known. By default, an NFW profile is assumed, with a concentration-mass relation from Duffy et al. (2008).\n", "\n", "One can alternatively pass a concentration directly:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mass in new definition is 0.72 x 10^12 Msun / h\n", "Radius of halo in new definition is 0.16 Mpc/h\n", "Concentration of halo in new definition is 3.31\n" ] } ], "source": [ "mnew, rnew, cnew = mdef_1.change_definition(m=1e12, mdef=mdef_2, c=5.0)\n", "print(\"Mass in new definition is %.2f x 10^12 Msun / h\" % (mnew / 1e12))\n", "print(f\"Radius of halo in new definition is {rnew:.2f} Mpc/h\")\n", "print(f\"Concentration of halo in new definition is {cnew:.2f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you have ``halomod`` installed, you can also pass any ``halomod.profiles.Profile`` instance (which itself includes a concentration-mass relation) as the ``profile`` argument." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Converting Mass Functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All halo mass function fits are measured using halos found using some halo definition. While some fits explicitly include a parameterization for the spherical overdensity of the halo (eg. ``Tinker08``, ``Watson`` and ``Bocquet``), others do not.\n", "\n", "By passing a mass definition to the ``MassFunction`` constructor, ``hmf`` can attempt to convert the mass function defined by the chosen fitting function to the appropriate mass definition, using ([Bocquet et al, 2016](http://adsabs.harvard.edu/abs/2016MNRAS.456.2361B)):\n", "\n", "\\begin{equation}\n", " \\frac{dn}{dm'} = \\frac{dn}{dm}\\frac{dm}{dm'} = f(\\nu)\\frac{\\rho_0}{m'}\\frac{d\\nu}{dm'} \\times \\frac{m'}{m}\n", "\\end{equation}\n", "\n", "where the last factor captures all the evolution with mass definition. You should note that this is only ever approximate -- it doesn't account for scatter in the profile of individual halos, nor does it preserve the halo model condition that all mass resides in halos (even if the original mass function did). In effect, the new $m'/m$ term is modifying the $f(\\nu)$.\n", "\n", "By default, the mass definition is ``None``, which uses whatever mass definition was employed by the chosen fit. Nevertheless, care should be taken here: many of the different fits have different mass definitions, and should not be compared without some form of mass conversion (and possibly not even then). To see the mass definition intrinsic to the measurement of each fit, use the following:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "init_cell": true }, "outputs": [], "source": [ "from hmf.mass_function.fitting_functions import SMT, Jenkins, Tinker08" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SO vir\n", "SO *(200m)\n", "FoF 0.2\n" ] } ], "source": [ "print(SMT.sim_definition.halo_finder_type, SMT.sim_definition.halo_overdensity)\n", "print(Tinker08.sim_definition.halo_finder_type, Tinker08.sim_definition.halo_overdensity)\n", "print(Jenkins.sim_definition.halo_finder_type, Jenkins.sim_definition.halo_overdensity)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here \"vir\" corresponds to the virial definition of Bryan and Norman (1998), an asterisk indicates that the fit is itself parameterized for SO mass definitions, and 0.2 is the FoF linking length.\n", "\n", "Let's convert the mass definition of ``Tinker08``, which has an intrinsic parameterization:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEWCAYAAACqitpwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xd8HNW1wPHf0aoXS+5Nrrh3sDGhm9Acug0BQ5JHSWwgQAohjRBK8l56dWiBQJxAggMECMUxmGKKwRUb917lXmVJVtd5f9yRdyVLltaa1a605/v5zGdn7szunmvZOr537twrqooxxhjjp4RoB2CMMab1seRijDHGd5ZcjDHG+M6SizHGGN9ZcjHGGOM7Sy7GGGN8Z8nFGGOM7yy5GGOM8Z0lF2OMMb6z5GKMMcZ3idEOIFo6dOigvXv3PqH3FhUVkZGR4W9AMc7qHB+szvHhROu8aNGifarasTHXxm1y6d27NwsXLjyh986ePZtx48b5G1CMszrHB6tzfDjROovIlsZea91ixhhjfGfJxRhjjO8suRhjjPFdq7jnIiIZwKNAGTBbVf8R5ZCMMSauxWzLRUSeFpE9IrK8Vvl4EVkjIutF5Ade8UTgRVWdDFzR7MEaY4ypIWaTCzANGB9aICIB4BHgC8AQ4HoRGQLkAtu8yyojFVBhaQWHjpRRWqFUVFZF6muMMabFi9luMVX9QER61yoeC6xX1Y0AIjIduBLIwyWYJUQwYf7mzTVM+3izO3j7vyQIpCQGSE5McFsggZQk77W6LDHBXRMIPQ7ZDySQklTzfOjnhH5GSsj52p+THEhARCJVdWOMCUvMJpd6dCfYQgGXVE4DpgIPi8ilwGv1vVlEpgBTADp37szs2bPD+vLNW0trHFcpFJdXUlwescZSWBITICmh+lWO7ieG7CclSMh1oeXB4+r96vKkgFBZVsLSF94mKUFIDnjlCUJSgBplAaHVJLnCwsKw/460dFbn+NAcdW5pyaWu31qqqkXAzQ29WVWfAJ4AGDNmjIb7ENGC/R+ye98mNpa1pUwTUA3r7RFXUeU2JzQ4PwIVoLThqwRSEhNITXItreoWV0pSAqmJAVKSQspCr/Neg+/zrksKfQ2+N/g+V5aa5H/rzR6uiw9W58hoacklD+gRcpwL7GiuL/9ul8V8d9mPqUxNIaHzYLTjICraD6S07QBK2vanOK07ZVVKaUUVZd52dL8ytKzSlVfWuibk2tKKymPeX1p+7OdUnyuvjI1Mpwol5VWUlEfnntTRxJTkEk5oIqtOQqHJKi0pQGrIlpbkEldacoD1uyuQtXtJTUwgLbn6vHtv9fuSArF829KY6GlpyWUB0F9E+gDbgUnADc327XtXAxCoKoWdS5CdS0gGkoEsgKQM6DgAOg6GToPca8eBkN0DEiL7S6iqSr2kVHcyO7pfR5Iqq052tRJeacj7t+/aQ1ZOO0rLKympqKK0vDLkukpKyt1rtJNcddyUVPjzgYvnH/d0YoJ4CSdAWrJrnaUlB0hNDJCaHAgmJq+8OjGlJdWdrILXuiSXnpxIerJLkK2lu9HEh5hNLiLyHDAO6CAiecADqvqUiNwJvAkEgKdVdUWzBZWUBhkdoWhv3efLi2DHYrfVeF+GSzKdBkPHQcHX7FzXj+SDhAQhNcH9gooE14we2+B1lVVKWUUVJeWVRxNPaYVLaCUVlZSWh5RVJ6Wj19Z6X3kdZcd8VvDasiiM4KuoUgpKKygo9SmZ1SNBICM5kbTkABkpiaQlBchICZCWnEhGsktKGV4iqk5I6SmBmsfea/BzXCJLSLCkZfwXs8lFVa+vp3wGMKOZw3Eu/S1c+lvmvPUqZw5oB3tWudbM3jVu/8i+ut9XXgQ7PnVbqORMl3Q6D4XOw9xrpyGQ3i7ydYmQQIKQ5v2ya25HW2/lNVtTtVtXoUmrpLyS4vLgfok3QKOkvIptO3eTmd22RllxWSWlFZUUl7myqmZqqFUpwSRW0PC9r3BUJ57MlABaXkLXtZ+QmZJEm9REMlMTyUxxr1kp1cdJZKYkkpXqturzKYnN/zM3sStmk0ssK09uA73Pcluoon3BhBP6Wnyg7g8qK4Tti9wWKqubl3Cqk84QaN8fEpMjU6FWombrLanJn+daa6fVe15VKa9UissrKQ1NQKFJqqySkopKisuqjpbVvvaYsqPvqeRIWSVHyioi2t14xPuefYXueMvhev6+NiA5kBBMRl7CyUlLIjstiZx095qdnuyOa5VnpSYRsBZUq2LJxU8ZHaDP2W4LVbgX9q6CPatrvhYfrPtzCna4bf2sYFlCUkgrx9s6DYWsLr51rZnwiAjJiUJyYgKkNT2ZHU9ZhUs6R8orKCp1iaeorOJoAqreP/paWkmxd211gqr56vb9HHhRVlnFgaIyDhSVhf1eEWiTWisRpSXRNj2Z9pnJtM9Ipn1mCu0ykumQmUz7jBSy05KsSy+GWXJpDpkd3dbnnGCZqrt3s2cl7F4Ju1fA7uWutVNRcuxnVJW787uX1yxPa1ezhdN5mOtaS0qNbJ1Ms6p+WDbbhxZZqMoq1/I6UlpBYWkF7388jwFDR1JQ4o4LS8op9LrjCksqQsq9stLyo+UVTegjVIX84nLyi8vZ2siGUyBBXPLJcAnIJZ4U2mck07lNKp2zU+ncJoUubVLJTkuyARHNrMnJRUQG4Ob2OgAsBZaq6pGmfm6kiMjlwOX9+vWLdiCQ2cltfccFyysr4MBGL5Gs8JLPcji0te7PKT4Amz9029HPDrgBA12GQ9cR0GWE20/LiWSNTAsUSJCj3VidgK3ZAc7s1yHsz1F1Q/BrJJ8SlywOeUnj0JFyL4GU1Toup+AERvdVVin7CkvZV1gKu49/bXJiwtFE06lNKl3auMTTuU0quw9WMjC/mE5ZqdY15yM/Wi5vAC8CbXFJZriIFKrqQB8+23eq+hrw2pgxYyZHO5Y6BRK94cwDYNjEYHlJvrt/U510qls7ZQXHfoZWwp4Vbls6PVie0ysk2Yxw+1ldrVvNNJmIHH1WqGNWStjvr6h0iSmYiFwCqu5m21dYxv7CUg4UlbG/qIx9haVhJaSyiiq2HShm24HiOs//bN67JAWEbjlpdM9JI7dtGrlt02u8dmmTat1wYfAjuexT1R+GFohIZx8+14RKzYaen3NbNVXXotm9ItittmsZHNhQ92cc2uK2VSEz5KR3qNm66ToS2p0U8edyjAmVGEigbUYybTMaP2iltKKSg0Xl7C8qZX9h2dHXvYWl7Dlcyu7DJew6XMKew6UUNmKoeHmlsmX/Ebbsr7vjJSUxgT4dMmpsfTtm0LdDZlhxxws/kstLInKRqr5VXaCqDTRSjS9EoG0vtw26JFheWgC7lsOupW7budS1eqrKj/2MI/tgw7tuq5aUAV2GuYTTbRR0Oxmpio3504yplpIYoEt2gC7ZDd9fLCytYPfhkpCtlF35bn/1tj0crkhkfwMDEUorqli9q4DVu47tLchJT6JvhwwGdsliYOcsBnZpw6AuWXGddPxILnuAF0XkE+BjYDGwRFXruUlgIi4lC3qd7rZqFWVusEB1stm1zG11dauVF8G2eW7znJWQDBtdojm6degPCfZsg4l9mSmJZHbM5KSOmcecq55n60hZBTsOFbPtYDF5B4vJO3jEey1m24Ejxx0Fd+hIOZ9uPcSnWw/VKO+UlcKgri7RDOqSxYjcHPp2yIiL7jU/ksv/AZfjZkccCVwFPAic4sNnG78kJrvur64j4GSvrKoKDm4KSThe0ik8tuEZqCqDvPluq5aU4T6vOtl0HQXt+1mXmmmR0pMT6dcpi36dsuo8n3+knE37i9i4t5BN+4rYuK+ITXuL2LSvqN6Z0fcUlLKnYC8frA3O6pGVksiw7tmM6JHNyNwcRuRm0z0nrdWNZvMjuaxR1fe9/Q98+DzTXBISoP1Jbhs6IVhesNtLOEtgh7cdzjv2/eVFsPUTt1VLznL3bbqFtHLa9rGEY1q87PQkRqXnMKpHzVGXqsquwyWs213IGq/bbM3uw6zbXejmuauloLSCTzbu55ON+4+WdchMZnSvtozt057T+rRjcNc2LX7kmh/JZZ2I/Ba4T1XrHophWpaszpB1IfS/8GjRnLde4cze6V7C8eZPK9h57HvLCmDLR26rlpIN3Ua6RNN9NHQfA9ndm6EixkSeiNA1O42u2WmcM6Dj0fLKKmXz/qKjCWfF9nw+y8t3Q6dr2VdYxpsrdvPmCtdrkJWSyOjebRnbpx2n9WnPyNxsElvYDNx+PUR5JrBNRNYRvOfyhE+fbWJAeXIODBgHAy4KFh7eGZJslri50+qa1LM0HzZ94LZqWV1doskd45JNt1HuXpExrUQgQTjJu89zyfCugGvl7MgvYem2Q3yWl8/SvEMsy8s/ZuLTgtIKZq/Zy+w17t9Tm9REzurfgXP6d+ScAR3plpPW7PUJ1wklFxFJB9JVdZ+q3uaVCTAQGIW792JauzZd3TbwC+5YFQ7vcMkmtIVzZP+x7y3YCatfdxuAJLgHP0MTTsdB7rkfY1oJEaG79yzNF7yEU1WlrN9byPxNB5i/6QDzNu1n9+GarZvDJRXMWLaLGct2AdC/UybnDujIRUO7MLpX25jsQgvrX66I9AL+CpzjDqUQ+C/wpKq+A6z2tun1f0p0xcwT+q2RiOvuyu4Ogy9zZaqQv80lme2LIG+R2y8vqvlerXKzEexZCYufcWVJGa5FczThjIY23e2hT9OqJCQIAzpnMaBzFl/+XC9UlW0Hipm3aT/zNh1gzvp97MyvOSXUuj2FrNtTyF8+2kSHzGQuHNKFLwzrwuf6tndz3cWAcP9b+A9gCHAfsAHoCowH3hSRJ4Gvq8ba4r81xfwT+q2NCOT0dNuQK11ZVaUbFp23ELYvdAln7yqXYEKVF8GWOW6rltnFJZrcMdDjNHcfJyn2uwiMaSwRoWf7dHq2T+eLY3qgqqzbU8j7a/bywbq9zNt4oMbaRfsKy3hu/laem7+VNqmJXDS0CxNP7s7n+raP6pDncJPLGGCyqj4TUjZVRM4AXsctOfxTv4IzrVRCIDiz8+gbXVlpoetKC004BXWsYF24q2Z3WkKie9izx2nQY6x7tcECphURCbZsJp/Tl+KySuZu3M+sVbt5a8XuGgMEDpdU8OKiPF5clEfX7FSuHNWdiad0Z0Dn5r+fGW5yKQQO1S5U1Y9F5AHg+1hyMSciJfPYNXIO7/C60ha61x2L3Ro4oaoqgguxzXvMlbXJ9RKNt3UZAYHITolvTHNJSw5w3qBOnDeoEz+9chifbj3Im8t3MXPFLvIOBgfs7swv4fH3N/D4+xsY1r0Nk07tyVUndyczpXnuY4b7Lf8B7hKR1+vo/loL2LS7xj9turlt8OXuuKrSrfqZt8Bt2+bDvjXHvu9wHqzIgxUvuePENOh+iks0uV7CyQh/5l9jYk0gQTi1dztO7d2OH106mKV5+by8eDuvfrajxowCy7cf5r7ty/n5jFVcMao7gxIjP51TuMnlHtwUL5+IyE+Ad1S1VETaAHfibu4bExkJAW/NmiHB7rQjB1zLpnq6mu2fHjtYoKL42Hs37U4KdqX1PN0txGYDBUwLJiKM7JHDyB45/OjSwby/Zi8vL97OrFW7KfMe5iwqq+S5+VsJCFxxfllE5z4LK7mo6kEROQv4Na4Vg4gcwrVYdgG3i0i2qub7HqkxdUlv5569qX7+prLCLTWwbX4w4dS1Fs6BDW777J/e57R3SaanNydbFxtNb1qupEACFwzpzAVDOpN/pJyXF+fxj3lbWbfHdSuP7hyI+KSaYXe+qep+4BYRuRu4EBhN8NmWVwEVkS3AYlW92s9gjWlQINFNP9N1JIz1BgQW7ApJNvPdwIHKWpMQHtlfc6BAciYjMvqBXOoSTu4YG5VmWqTs9CRuOrMPN57RmwWbD/LPeVsYlNzI5T6b4ITv7KjqIeAFbwNARLrgEs3JwIgmR2eMH7K6wJAr3AZQXgI7P3PJZutcNzdaca1/bGWFtCtbAu8tcccJSe6+Tc/TodcZrkvNVvY0LYiIMLZPO8b2acfs2bMj/n2+DhtQ1V3ATG8zJjYlpULP09x25jfc7ND71sLWj2HLx7Dlk2Mn6qwqD3azzfkDINB5mLe0wRnQ6yzI7Fjn1xkTjxqdXEQkEzgXGIRb0lhxw5JXA++rauFx3m5M7EpIgE6D3DbmFld2aCurZj7N4PQDrmWzb22tNynsXua2+d40eh0HQ5+zoffZbkh1ertmrYYxsaTB5OLNGfYQcDeQDhwBDgICZAMZwBFvZuQHY/0JfZv+xTRKTk92dxnH4HHj3HHh3uDyAls+dksS1J5RYO8qt81/gqMtm+pk0+sM60YzcaUxLZcHcYnlIWC6qm4LPSkiucAk4AFca+ZBf0P0l03/Yk5IZsea921KDruF07Z84oY45y2stYx0SMtm7qNuYs4uI6DPOW7r+TmbBdq0ao1JLl8D7q5vCn1VzQN+IyKHcQnmQf/CMyZGpbaBfhe4DaCsyA0O2PwhbPrQzSagIQ+qaZUbpbZzCXw8FSTgBgj0OQdO+rx7uDMxftdbN61PY5JLDm6SyoZswJ7QN/EqOQP6ne82cC2brXNh8wcu2dTuRtPK4EwDH/7WzQDd+yw46TyXbDoMsIc6TYvWmOQyF/ieiMxV1aK6LhCRDNy8Yp/Udd6YuJPapubDncWH3L2a6pbN7mU1ry8vgnVvug0gq5tLMiedB33H2XQ1psVpTHK5E3gb2Coib+JGhx3C3V/JwY0euxgoBc6PUJzGtGxpOTDoEreBm7Zm84ewcTZseBcObq55fcEOWPKs28Ddr6lONj1Ph8SU5ozemLA1mFxUdZWIDAVux63dcj5uKDK4UWOrgd8Aj3sPVhpjGpLezq1vU73GzYGNsOE92PgebPzALQ0datdSt835g+tC63su9L/IbbbEgIlBjXrOxUsaP/c2Y4zf2vV126lfdfOj7VjsWjQb33NT1oQODigvgjUz3AZuyHP/i2DAxW55aFsa2sQAPx6iXAV8YA9RGuOTQCL0ONVt477vBgds/sglmw3vuFZOqN3L3fbR7yA1x41gG3AxnHQ+ZLSPTh1M3Iu7hyiNaXFS29S8X7N/A6x9E9a95Z6xCZ2Es+QQLH/RbZLghjgPvgwGXepaRsY0k7h7iNKYFq/9SXD6191WWugGBax7E9bNgoKdweu0CrbNddtb90GnoV6iuQy6DLehziai7CFKY1qylEyXMAZfBqqwa1kw0eQtqPlszZ4Vbnv/l5DT0yWZQZe52QISAtGrg2mV7CFKY1oLEeg6wm3nfBeK9sGa/7o1aja8B5WlwWsPbXXT0sx91C2UNvhyGHZ1zYEDxjRB3D1EaRNXmriR0QFO+YrbSgtg/Tsu0ax9E0oPB687sh8WTYNF0zg9uS2UXAdDJ0LuqW7GaGNOQNw9RGkTV5q4lJIFQ69yW0WZe4Bz9euw+g0o3B28rOwgzHvcbdk9YOgE16LpOtLu0Ziw2EOUxsSbxOTgPGiX/NYtgLb837DyFSjaG7wuf5ubZPPjqdDuJBh1PYyYBDk9ohe7aTHsIUpj4llCgrea5ukw/hcsefURRgXWw8pX3bDmagc2wLv/C+/+n5vJedQN7j5Nckb0YjcxzTpUjTFOIJFDbUfCFX+Ce9bBDS+4lkpy6LozCpveh5dvhd8MgFfugM1z3Eg1Y0L4Nk+EiCQBXVV1q1+faYyJksTk4KzOZUfcvZnP/ulGneElkrLC4OSa7fvDmJth5PW2vLMBGtlyEZE7RGSDiBSIyDwR+Uodl50CbPI3PGNM1CWnw4gvwldehm+vgPMfcOvNhNq/Dt68F343GF6+HbYtsNZMnGswuYjIJOBPuCHJDwE7gGki8qKIpEU4PmNMLMnuDmffDXfMh6+9A2NuqdltVlHiWjhPXQCPnw0Ln3ardJq405iWyz3Ab1T1S6r6G1WdAFwEnAW8JyI2M54x8UYEcsfAZb+H76yGy6e64cqhdi+D178Nvx8Kbz8Eh3fW/VmmVWpMchkIzAgtUNV3gM/hJq78REROikBsxpiWICUTRt8IU96Hye/CyV+GxJBOjeKDbsbmPwyHl6bAjiXRi9U0m8Ykl3zgmDVWVXUzcAawD/gYONXXyIwxLYsIdB8NVz7iWjMX/xxyegXPV5XD0n/BE+fCtMvcjAF2X6bVakxyWQRcVdcJVT2Ie6hyITDVx7iMMS1ZWo6btfkbi+HaZ9zSzKE2fwjPToQnPw+rZ1iSaYUak1yeBfqKSJ3jC1W1GLgC+Atgw5CNMUEJARhyBdwyE772rptKRkJmYN7xKUy/Hh4/C1a8DFU2cWZr0WByUdUXVPUMVT1wnGsqVXWKqvbxNzxjTKuROxqueRq+uQROnQyBlOC53cvhhZvgsTPdMzXWkmnx7Al9Y0zzyukJl/4GvrUUTr8TktKD5/auguk3wFMXuqWdTYtlycUYEx1ZXeDi/4NvLYezvwPJmcFzeQtg2qXwzEQbXdZCxV1yEZHLReSJ/Pz8aIdijAHIaA/n3w/f/Aw+dwcEkoPnNrzjRpe9fDsU7IpejCZscZdcVPU1VZ2SnZ0d7VCMMaEyOsD4n8Fdi2DUl0FCfj199k/402j46PdQUVr/Z5iYEXfJxRgT43J6wlWPwO2fwMBLg+VlhfD2g/DIaW74solpJ5RcRKSbiFwtIpO9125+B2aMiXOdBsH1/4SvvAIdBwfLD25yw5enfwkO74hefOa4wkouIhIQkUeBLcALwJ+91y0i8oiIWEvIGOOvk86D2z6CL/waUnOC5atfh4fHwvwnoaoqevGZOoWbDB4CbgHuBXoDad7rvV75g/6FZowxnkAinDYF7voUTrkxWF5WADPugacvht0roxefOUa4yeV/gPtU9dequlVVS73XXwM/Bm7yPUJjjKmW0R6umAo3zXALlFXLmw9/Psfd8Len/GNCuMmlE7C0nnNLvfPGGBNZvc90XWXnfh8SklxZVbm74T/tUji4OZrRGcJPLmuBSfWcmwSsaVo4xhjTSEmpcN69Lsl0OyVYvvUTN43Mp8/YNDJRFG5y+V/gJhF5W0RuE5EJInKriLwN3OidN8aY5tNpEHz1LRj3w+CkmGWF8Oqd8OItUFoQ3fjiVFjJRVWfB8YDGcAfgX/jptpPB8ar6gu+R2iMMQ0JJMG4H8DXZkH7fsHyFS/BE+Ng1/KohRavwh46rKpvqerpuJFiXYA0b9bkWb5HZ4wx4eg+Gm79EEbfHCzbvx7+cj4s+pt1kzWjcJ9zub/6gUlVrVLVPapa5Z3rKiL3RyJIY4xptOR0uPwPcPVTkJThyipK4LVvwH/utOljmkm4LZcHgNx6znXzzhtjTPQNvwamzIZOQ4JlS56Fv10OhXuiFVXcCDe5CFBfuzIXONi0cIwxxkcdB8DX3oGRNwTLts2DJ86DnZ9FL644kNjQBSJyI24kGLjE8piIHK51WSowHHjL3/CMMaaJktPhqkeh81CY9WPQKjicB0+PhwmPAzZDeiQ0puVyBNjvbQLkhxxXb5uAXwFTIhOmMcY0gQiccSfc8DyktHFl5Ufg+RvpnvdadGNrpRpsuXjDi18AEJG/Aj9R1U2RDswYY3zX/0L42tvw3CQ4sBFQ+q//C8zKhAsedEnI+CLc51xubumJxVaiNCbOdRzo7sPkjg2WzfkDvHI7VJZHL65WJu6myLeVKI0xpLeD//kPDBgfLPvsOfjXl22osk/iLrkYYwzgbvRf9w92drkgWLZ2Jky/AcqLoxdXK2HJxRgTvwKJrBl4J5z5rWDZ+rfhn9dB2ZHoxdUKWHIxxsQ3EXcz/9zvB8s2vQ//+CKUFkYrqhavwdFioUQkBbgZGAgcAJYDS1V1QwRiM8aY5iHipu8PJMG73uTuWz5yXWQ3PO+m9zdhCSu5AP8ErsIllQzcEsciIkXACuAzVb3N1wiNMaa5nPNdCCTDLG+axE3vu2n7r/2bSzym0cLtFrsIuEtVR6pqPyALOB24G1gADPI5PmOMaV5nfhM+f1/weM0b8MrXoaoqejG1QOG2XLbinsYHQFWLgfneZowxrcPZ97hFxub80R0vex5SMuHS39mDlo0UbsvlF8DXIxGIMcbEDBG44CEYc0uwbOHT8NHvohdTCxPuE/rPAJtFZJaIfF5ErBPSGNM6icAlv4XhXwyWvfMTWPZi9GJqQcIdLfYd4A7v8HygXERWA59521JbkdIY02okJMCVj0LBLtj8oSt75XZo0x16nR7d2GJcuN1iPwKexY0SGwr8D/AG0A74JjDTz+CMMSbqEpPhumegwwB3XFkG06+HfeujG1eMCze5lAPTVHWrqq5S1X+p6r2qepmq9gTaRyBGY4yJrrS28KUXIKOjOy4+6GZWLrEJcOsTbnJ5FtcdVidVPdS0cIwxJka17Q3X/wsS09zx/nXw8m02RLke4SaXLcD1InKHiAQiEZAxxsSs3NFw5cPB4zUz4MPfRC+eGBZucvkZ7n7Ln4A9IvIfEXlIRCaKyEm+R2eMMbFm+DVw+p3B4/d+BmvsdnNt4SaXLKA/cDXwR6AMmAQ8D6wTkQJ/wzPGmBh0wUPQ5xzvQOGlybDfplgMFe5zLqqqG1T1ZVX9iap+UVUHApnAWOCuiERpjDGxJJAI1/wVsnu449LDbg6yirLoxhVDGkwuIvJjb2ngXvVdo6olqrpQVaf5Gp0xxsSqjA5w7d8hwXuWfOcSeOeh6MYUQxrTcrkReAXYKCIHRWS2iPxRRG4WkVO8afiNMSb+dD8FLgxJKJ88DOvejl48MaTB5OLNfpwDjAMeALrjur+ews2EXCAiy0TkWRH5bgRjNcaY2PO5r0P/i4LHr9wGBbujF0+MaNQ9F1UtUNUPgRSgGDgb6AaciRtB1gWYiN1zMcbEGxG46jHI7OKOi/bCy7fG/fMv4Y4W+y7wA1Wdo6q7VHWuqj6IW8dlLWAdjsaY+JPRASb+GfCm49/4Hix8KqohRVu4ySUApNcuVNX9wE+B7/gRlDHGtDh9x7mFxqrNuh8ObKrv6lYv3OTyAnC/iOTUca4M94BlTPNGvj2Rn29zAhljfHbevdBxsNsvPwL/uSNuu8dOpFusCPdvgnIlAAAbv0lEQVTA5EMicpaI9BSR83ELia32PUKfqeprqjolOzs72qEYY1qbxBS46lGonh1ryxyY/+foxhQl4T5EWQCcg5v+5evAB7hlj2fhnt6/1e8AjTGmRel+Cpx9d/D47YfgwMboxRMl4bZcUNVyVf0J0Bk4GbgUOA3op6oLfI7PGGNannO+B52Huf2KYnjjHlCNbkzNrDFP6O8Ukb+IyAQRyawuV9UqVf1MVWeq6gJVtXkPjDEG3AJjV0zl6OixDe/A8n9HNaTm1piWyzdxyyE/BuwTkVki8i0R6R/Z0IwxpgXrPhrGTg4ez/whFMfPkleNeUL/eVW9CegKnAt8DHwZWC0ia0Xk9yJygYgkRTZUY4xpYT5/H2R1dftFe+Jq7rFG33PxZkSep6oPqOoYIBf4JdAT+DewX0ReEpFbIhSrMca0LKnZ8IVfBo8X/hW2L4pePM0o7Bv61VR1p6o+papXAx1w079sAb7nV3DGGNPiDb4C+l/sHSj89/txcXM/7OQiIleKyD0icpOInCoiad4IsrdV9duqOigSgRpjTIskAuN/DoFkd5y3AJa9EN2YmkFYyUVEngBeAr4NPALMAw6LyBoReUFE7otAjMYY07K1P8nNnlxt1v1QWhi9eJpBuC2Xa4H7VbW7qmbgljy+DngON+/YjT7HZ4wxrcM590BmZ7dfsBPm/CG68URYuMmlAJhbfeAtefySqj6oqhNV1YYnG2NMXVKy4PwHgsdzpsKhrdGLJ8LCTS7TgPERiMMYY1q/kddDt1PcfmUpvPez6MYTQeEmlzzgShH5hogkRiIgY4xptRIS4KKfBo8/mw67V0QvnggKN7n8FugH/AHYKyKvishPReRqEennf3jGGNPK9D4rZFlkdRNbtkLhJpcs3E38q4HfA6W4m/z/AtaKSIG/4RljTCt0/gMcnXds3ZuweU5Uw4mEsLq2VFWBDd72cnW5iKQCw7zNGGPM8XQZBiOug6XT3fHbD8BXZ7lnYlqJxsyK/GNv9cZe9V2jqiWqulBVp/kanTHGtFbn3Vvzwcq1M6Mbj88a0y12I/AKsFFEDorIbBH5o4jcLCKniEhKhGM0xpjWp20vGPPV4PHsX7SqaWEaMytyPyAHGAc8AHQH7gKeAhYABSKyTESeFZHvRjBWY4xpXc76NiSmuv2dS2Dtm9GNx0eNuqGvqgWq+iGQAhQDZwPdgDOBnwFdcBNX3hWhOI0xpvXJ6gyjbw4ev996Wi/hjhb7LvADVZ2jqrtUda6qPggMAtYCrXNMnTHGRMqZ34SAd3dhx2JYNyu68fgk3OQSANJrF6rqfuCnwHf8CMoYY+JGm64w+qbgcStpvYSbXF4A7heRnDrOlQG9mxyRMcbEm7O+FWy9bF8EG96Jbjw+OJFusSJgnYg8JCJniUhPETkf+AWw2vcIjTGmtWvTDUaHTCo/54/Ri8UnYSUXVS0AzgH+BHwd+ADYBMzCPb1/q98BGmNMXDj9TpCA29/0gbv/0oKFvRKlt+rkT4DOwMnApcBpQD9VXeBzfMYYEx/a9oKhE4LHc6ZGLxYfhJ1cqqlqlap+pqozVXWBqpb5GZgxxsSdM78R3F/5ChzYFL1YmuiEk4sxxhifdR0Jfc9z+1oFcx+NbjxNYMnFGGNiyZnfDO5/+gwU7Y9eLE1gycUYY2JJ33HQZbjbryiGRU9HM5oT1phZkZ8WkT7e/jkikhn5sIwxJk6JwBkh914WPA2V5dGL5wQ1dlbkjt7+e8CQyIVjjDGGIVdBRie3X7ADVr8e3XhOQGOSy05gnNdiESBVRNLr2yIbbtN5a9M8kZ+fH+1QjDGmbonJMCZkQst5T0QvlhPUmOTyBO7p+3xAca2XguNsMU1VX1PVKdnZ2dEOxRhj6jf6ZkjwFgve+jHsWhbdeMLU4DLHqvoTEXkDGAz8Hfhf3DLHxhhjIqVNVxh8Bax4yR3P+zNc+XB0YwpDg8kFQFUXAYu8OcT+qqot98keY4xpKU67NZhclr0AF/4E0ttFN6ZGCndusZtVdZOIdBORq0VksvfaLVIBGmNM3OpxGnQZ4fYrSmDxs9GNJwxhJRcRSRCRR4EtuOn3/+y9bhGRR0TEnpsxxhi/iMDYKcHjT//eYtZ6CTcZ/AS4BbgXt3ZLmvd6r1f+oH+hGWOMYdhESM5y+/vXwdZPohtPI4WbXP4HuE9Vf62qW1W11Hv9NfBj4CbfIzTGmHiWnAHDrwkeL/pb9GIJQ7jJpROwtJ5zS73zxhhj/HTK/wT3V74CxQejF0sjhZtc1gKT6jk3CVjTtHCMMcYco9vJIfONlcCyF6MbTyOEm1z+F7hJRN4WkdtEZIKI3Coib+Omiflf/0M0xpg4JwKnhCyDvOhvMX9jP9yhyM8D44EM4I/Av4GpQDowXlVf8D1CY4wxMPyLkJjm9ncvgx2fRjeeBpzIMsdvqerpuJFiXYA0VT1DVWf5Hp0xxhgnLQeGXhU8XvyP6MXSCE1d5niPqlb5GZAxxph6nPzl4P6Kl6AidleXt4cejTGmpeh5BmT3cPvFB2HdW9GN5zgsuRhjTEuRkAAjrg0eL50evVgaYMnFGGNakhEhT4OsmQlHDkQvluOw5GKMMS1JxwHuuReAqnJY8XJ046nHCSUXmxXZGGOiKLT1svRf0YvjOMKdFTlgsyIbY0yUDbs6uErltnlwYGN046lDuMngIWxWZGOMia7MjtDvguDx0th7ft1mRTbGmJYodNRYDN53sVmRjTGmJRowPjgdzN5VsGdVdOOpxWZFNsaYlig5AwZcHDxe8Ur0YqmDzYpsjDEt1dAJwf0VL8fUTMmJ4Vysqs+LyCHcjf0/AklAObAINyuyTV5pjDHNpf9FkJQO5Udg3xrXNdZ5SLSjAsJMLuBmRQbe8oYddwD2tZbJK8vLy8nLy6OkpOS412VnZ7NqVWz1b0ZadnY2mzZtIjc3l6SkpGiHY4wBSE53XWPVN/RXvNxyk0s1L6Hs8TGWqMvLyyMrK4vevXsjIvVeV1BQQFZWVjNGFn2HDx+mrKyMvLw8+vTpE+1wjDHVhk4IJpeVr8B597rFxaKsweQiIvOBm1R1pYgsAI7bqaeqY/0KrrmVlJQ0mFjilYjQvn179u7dG+1QjDGh+l0ISRlQXgT71sKeldB5aLSjalTLZQVQHLIfO3eMIsASS/3sz8aYGJScDgPHw/J/u+MVL7eM5KKqN4fs3xTRaIwxxoRvyFXB5LLqdfj8fdGNh/DnFru/vkkqRaSriNzvT1jGGGMard/5kJjq9veuiom5xsJ9zuUBILeec92886YJSkpKGDt2LCNHjmTo0KE88EDwj3TmzJkMHDiQfv368Ytf/KLBcmNMnEjOgL7nBY9Xz4heLJ5wk4tQ/z2XXOBg08IxKSkpvPvuu3z22WcsWbKEmTNnMnfuXCorK7njjjv473//y8qVK3nuuedYuXJlveXGmDgz6JLg/uo3oheHp8HkIiI3isi7IvIuLrE8Vn0csn0MPAu8H+mAWzsRITMzE3DP3ZSXlyMizJ8/n379+tG3b1+Sk5OZNGkS//nPf+otr8+jjz7KsGHD6NWrF3/605+aq1rGmEgb8AXc//+BbXOhaF9Uw2lMy+UIsN/bBMgPOa7eNgG/AqZEJsz4UllZyahRo+jUqRMXXnghp512Gtu3b6dHjx5Hr8nNzWX79u31ltfl3//+N7NmzWLx4sXMnTuXhx56iIqKiojXxxjTDDI7Qo/T3L5WwdqZUQ2nMaPFXsAtCIaI/BX4qapG/25RhPX+QeSalZt/celxzwcCAZYsWcKhQ4eYMGECy5cvR+uYM0hE6i2vy9SpU3nyySdJSkqia9euJCUlUVXVKiZXMMaA6xrbNtftr54BJ385aqGEdc9FVW+Oh8QSK3Jychg3bhwzZ84kNzeXbdu2HT2Xl5dHt27d6i2vrby8nKVLlzJgwAAAdu7cSYcOHThw4ADnnnsuv/rVr7jxxhv585//zDXXXMPy5csjX0FjjL8GXRbc3/AulB2JWihhL0ssItd5syJvFZE9tbdIBBlP9u7dy6FDhwAoLi7m7bffZtCgQZx66qmsW7eOTZs2UVZWxvTp07niiivqLa9t5cqV5Ofns3HjRqqqqvjhD3/IXXfdxeLFi5k4cSLf+973yM/PZ/LkyXzxi19ky5YtzV11Y0xTtT8JOgx0+xXFsPG9qIUS1txiInID8DQwDfi8t58AXAEcAv7uc3xRc7yuq0jOLbZz505uvPFGKisrqaqq4tprr+Wyy9z/Rh5++GEuvvhiKisrueWWWxg6dOhxy0MtXryYL33pS1x//fUUFRUxceJEpkyZws9//nMmTJhAeXk57du3JyEhgeXLlzN58uSI1M8YE2GDLoWPvKW1Vs9wx1EQ7sSV3wV+CvwCd/P+UVX9VESygFm4m/+mCUaMGMHixYvrPHfJJZdwySWXNLo81JIlS7jsssu47rrrapSvX7+eAQMGsHTpUgYPHgzA5s2b6dmz5wnWwBgTVQMvgY9+5/bXz4KqKkgIu5OqycL9xv7AHFWtBCqBNgCqWgD8ErjT3/CMX5YsWcKoUaOOKX/qqadISEhg1KhR3HPPPQA888wzzR2eMcYv3U+BtHZuv3A37F4WlTDCTS75QIq3vx0YHHJOgPZ+BGX8N3v2bAYOHBjtMIwxkZYQcNPBVFsXnTUcw00uC4ER3v6rwP0iMllEbgR+DczzMzhjjDEnoN+Fwf0oJZdw77n8HOjl7d/v7T8KBIAFwK3+hWaMMeaE9Dufo7N15c2H4oOQ1rZZQwj3OZe5qvovb/+Qql4JZAI5qnqaqm6IRJDGGGPCkNEBup3s9rUKNjT/kOQmDyFQ1VJVPSwi54nIf/0IyhhjTBP1vyi4v/7tZv/6RiUXEckRkUki8l0RuUZEkkLOfVFEFgLvALa4ujHGxIL+te67NPNUTw3ecxGR4cBbQOeQ4k9F5Grgn8DngJXAl4B/RSJIY4wxYep2shuSXHwAivbArqXQ7djHESKlMS2XnwGHgdOBdNzw4wO4G/jDgBtVdbiqPqeqNguiMcbEgtpDktc376ixxiSXMcCPVXWeqpao6hrgdqAD8B1VfTaiERpjjDkxNe67vNOsX92Y5NIZ2FyrrPr4Mz+DMU7v3r0ZPnw4o0aNYsyYMUfLbZljY0xY+o4L7uctgNLCZvvqxo4Wq29pY1tpKkLee+89lixZwsKFCwFsmWNjTPgyO0HnYW6/qgK2ftJsX93Y5PJmrWn1d3rl79iU+83Dljk2xpyQPucG9zfObravbcwT+g9FPApTg4hw0UUXISLceuutTJkypc7ljOfNm1dveV1Clznet28fw4cP5/bbbycxMdyJGowxLUbfcTD3Ebe/8f1m+9rGLHMc88lFRPoCPwKyVfUaXz70wex6TzV5JZcH8497es6cOXTr1o09e/Zw4YUXMmjQIFvm2BhzYnqdAQmJrlts9zIo3NssX9v8k/zXIiJPe11qy2uVjxeRNSKyXkR+cLzPUNWNqvrVyEbafKqXKe7UqRMTJkxg/vz5EVvmOCEhgbvvvpvvfOc7TJ06NcI1M8Y0u5RMyD01eLypeVovUU8uuFUtx4cWiEgAeAT4AjAEuF5EhojIcBF5vdbWqflDjpyioiIKCgqO7r/11lsMGzYsYsscP/bYY1x55ZX89re/5Rvf+EZzV9cY0xxC77s0U3KJeme7qn4gIr1rFY8F1qvqRgARmQ5cqao/By5rlsCO03UVyWWOd+/ezYQJEwCoqKjghhtuYPx4l3sjsczxzTffzO233x6RuhhjYkTfcfC+95jCxtkw6uqIf6XU1Wff3Lzk8rqqDvOOrwHGq+rXvOOvAKepap0rXYpIe+D/gAuBv3hJqK7rpuCWZ6Zz586jp0+fXuN8dnY2/fr1azDeyspKAoFAo+oWK77//e8zduxYrr665l+qN954gxkzZtC2bVvuvvtu2rVrV+f7q+u8fv168vOPf8+otSgsLCQzMzPaYTQrq3PrJFUVnDnnSyRWlgDw7vDfk9C+b9ifc9555y1S1TENXxkDLZd61HVHut4sqKr7gdsa+lBVfQJ4AmDMmDE6bty4GudXrVrVqBZJJFsukbJy5Uq++c1vHhP3pEmTmDRpUoPvr65zamoqJ598cqTCjCmzZ8+m9t+R1s7q3IrtPAfWvQVA99J1DBx3S0S/LlaTSx7QI+Q4F9gRpVhahdmzZ0c7BGNMNPUddzS5tD0Y+clVYuGGfl0WAP1FpI+IJAOTcMsqG2OMOREhN/XbHlwa8Sn4o95yEZHngHFABxHJAx5Q1adE5E7gTdwSyk+r6ooohmmMMS1bpyHQfTR0Hsq64o4M0Uoi2b6IenJR1evrKZ8BzGjmcIwxpnVKSIDJ7wKwZ/ZshgSSGnhDE78uop/eAsXC6LlYZX82xpjGsuQSIjU1lf3799sv0TqoKvv37yc1NTXaoRhjWoCod4s1NxG5HLi8rudZcnNzycvLY+/e48+9U1JSEne/ZEtKSsjJySE3NzfaoRhjWoC4Sy6q+hrw2pgxYybXPpeUlESfPn0a/IzZs2fHzbMe1eKxzsaYE2fdYsYYY3xnycUYY4zvLLkYY4zxXUxMXBkNIrIXOASEzsKYfZzj0P0OwD4fwqj9fU25tr7zx6tTQ8dWZ6vzibI6N+3aWK1zL1Xt2KgrVTVuN+CJxh7X2l8Yie9vyrX1nQ+njlZnq7PV2ers1xbv3WKvhXFc+1wkvr8p19Z3Ppw61j62OvvD6ty0a63O9ZdHu871ittusaYQkYXayDUNWgurc3ywOseH5qhzvLdcTtQT0Q4gCqzO8cHqHB8iXmdruRhjjPGdtVyMMcb4zpKLMcYY31lyMcYY4ztLLidIRPqKyFMi8mJIWYaI/E1EnhSRL0Uzvkiop87HlLUm9dT5Ku9n/B8RuSia8UVCPXUeLCKPi8iLInJ7NOOLhPr+Hnv/pheJyGXRii1S6vk5jxORD72f9bimfL4llxAi8rSI7BGR5bXKx4vIGhFZLyI/AFDVjar61VofMRF4UVUnA1c0U9hN0tQ61/PnENN8qPMr3s/4JuC6Zgu8CXyo8ypVvQ24FmgRw3Z9+PcM8H3g+eaI1w8+1FmBQiAVyGtKLJZcapoGjA8tEJEA8AjwBWAIcL2IDKnn/bnANm+/MkIx+m0aTatzSzQNf+p8n/eelmAaTayziFwBfAS8E7kwfTWNJtRZRC4AVgK7Ixumr6bRtJ/zh6r6BVxSfagpgVhyCaGqHwAHahWPBdZ7Wb4MmA5cWc9H5OESDLSQP1sf6tziNLXO4vwS+K+qfhrZaP3hx89ZVV9V1TOAFtHl60OdzwM+B9wATBaRmP833dQ6q2qVt3sQSGlKLDH/hxUDuhNsjYBLIN1FpL2IPA6cLCI/9M69BFwtIo/RzFMt+KzRda7nz6ElCufnfBdwAXCNiNzWzHH6KZyf8zgRmSoifwZmRCFWvzS6zqr6I1X9FvBP4MmQX7wtTTg/54nez/gZ4OGmfGncrUR5AqSOMlXV/cBttQqLgJubJarICqfOx5S1UOHUeSowtVmiiqxw6jwbmN0MMUVao+sccnJaRCOKvHB+zi/h/pPcZNZyaVge0CPkOBfYEaVYmovV2ercWlmdm6nOllwatgDoLyJ9RCQZmAS8GuWYIs3qbHVurazOzVRnSy4hROQ54BNgoIjkichXVbUCuBN4E1gFPK+qK6IZp5+szlZnrM5W50jEYhNXGmOM8Zu1XIwxxvjOkosxxhjfWXIxxhjjO0suxhhjfGfJxRhjjO8suRhjjPGdJRdjjDG+s+RijDHGd5ZcTNwQkQdFREVkXT3n13vnH2zm0BpFRH7vxbesjnM5InLAO3+PD991rYjsEhHxjieJSKmIJDX1s018sORi4k0J0EdEaqymKCKnAr2887FqOFAADPAWgAr1PSDZ2z8m+ZyAS4EZGpzCYySwXFXLffhsEwcsuZh4UwS8i5u8L9Qkr7yo2SNqvOG4CQeTgb7VhSLSCfgGwckIlzblS7xFscYDb4QUjwIWN+VzTXyx5GLi0XTg2pAuH8GtDT899CIROV1EXhWRHSJSJCJLROSYVRhFZKiIzPS6pYpEZJWI3NHQuXB4CaQT8Dqu9TIo5PR9uISyBdinqjvD/fxaTgXaArNCykYCq0TkZyKyXUTyReTJlrA6o4kO+4th4tFLQGfgLO/4bKAj8HKt63oBc4CvAZcD/wb+KiLX17ruVaAS+DJwBfAnIKsR58IxwntdilvXfTCAiPQEbgXu9a7xq0vsQ1U97H1HR6ArbgXONOAm4De4P5dWs/y18ZetRGnijqoeEpGZuK6wD73XmV556HVHWzJe6+YD3EJLk4HnvPIOuC6qq1S1+hf7Ow2dOwEjgFJgLbACL7kADwIfqOpsEfk7/qwieCnwj5DjUd7rVFX9nbc/S0RuB/r78H2mFbKWi4lX04FrRCQFuIZaXWIAItLWWzd+C1DubVOAASGXHcCtT/64iFzndV815ly4hgOrvLU5VgCDRGQQ8BXgRyKSjVttsKn3W7oCJ1PzfstI4BCu1VV9nQA5wL6mfJ9pvSy5mHj1KpAJ/B+QAbxWxzXTgOuAXwMX4e5FPA2kVl+gqlXeuV3euV0i8qGInHy8cycQ73CCXV7VLZefAm+o6nzvPNTTLSYiQ0TkaRF5T0QeFpFu9XzPJcBGVV0TUjYSeL/WSLGTcF1krWahLeMvSy4mLqlqEe7m+LeB17zjo0QkFdc99ICqPqyq76rqQur4N6Oqq1X1atz/5C/AJZ83RCTheOcaG6t37RCCiWM5kA1MxN3MB9dtVkUdv+xF5PPAk8ATuPs+rwOvey2f2i6lZqsFXLdY7RbRSO/7lje2Hia+2D0XE88eA1KAx+s4lwIEcPc5ABCRLNwv5zqXb/X+Z/+uiPwO+CcuoRxo6Fwj9Me1EpZ5n7VdRJ4HVqpq9S/34cB6VT1Sx/vvA64JGUU2U0QOeOVfDqlfMi4BXlOrbBDwWa3PHIFr4cTy0G0TRZZcTNxS1dnA7HrO5YvIAuB+ETmM+1/6D4B8oE31dSIyAjdy6l/ARtwQ3u/jfhnnisj0us6p6gHv/eOA94DzvHjqckyXl6peV8c19Y0US1LVnSLSG1itqqmqOt87DnUOrmX2fkjZUNzvidotlxF1lBlzlCUXY+p3A64r6e/AfuBhIB24M+SaXcBu4EdAN9yN7/dwSaTkOOeqpXuve44Tx3DgoKpuP841w4C36jknIpIGbMcb0iwivTi25XQp8LaqloaUjQSOABtqXTsC9+diTJ0kOLuDMaa5ichDwDmqel4Ev+NruO6uKap62LuZ/xzwE1V9J+S6tcCvVfXJSMVi4oe1XIyJrjOA3zV4VROo6l9EpAB4x5t4shD4kaq+X+u6AXV+gDEnwFouxhhjfGdDkY0xxvjOkosxxhjfWXIxxhjjO0suxhhjfGfJxRhjjO8suRhjjPGdJRdjjDG+s+RijDHGd/8P+HY+9bw6q+oAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Default mass function object\n", "mf = hmf.MassFunction()\n", "\n", "dndm0 = mf.dndm\n", "\n", "# Change the mass definition to 300rho_mean\n", "mf.update(mdef_model=\"SOMean\", mdef_params={\"overdensity\": 300})\n", "\n", "plt.plot(mf.m, mf.dndm / dndm0, label=r\"300 $\\rho_m$\", lw=3)\n", "\n", "# Change the mass definition to 500rho_crit\n", "mf.update(mdef_model=\"SOCritical\", mdef_params={\"overdensity\": 500})\n", "\n", "plt.plot(mf.m, mf.dndm / dndm0, label=r\"500 $\\rho_c$\", lw=3)\n", "\n", "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "plt.xlabel(r\"Mass, $M_\\odot/h$\", fontsize=15)\n", "plt.ylabel(r\"Ratio of $dn/dm$ to $200 \\rho_m$\", fontsize=15)\n", "plt.grid(True)\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This did not require any internal conversion, since the `Tinker08` function is able to deal natively with multiple definitions. Let's try converting a fit that has no explicit parameterization for overdensity:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/steven/Documents/Projects/halos/HMF/hmf/src/hmf/mass_function/hmf.py:274: UserWarning: Your input mass definition 'SOMean(300)' does not match the mass definition in which the hmf fit SMT was measured:'SOVirial'. The mass function will be converted to your input definition, but note that some properties do not survive the conversion, eg. the integral of the hmf over mass yielding the total mean density.\n", " f\"Your input mass definition '{mdef}' does not match the mass \"\n", "/home/steven/Documents/Projects/halos/HMF/hmf/src/hmf/mass_function/hmf.py:274: UserWarning: Your input mass definition 'SOCritical(500)' does not match the mass definition in which the hmf fit SMT was measured:'SOVirial'. The mass function will be converted to your input definition, but note that some properties do not survive the conversion, eg. the integral of the hmf over mass yielding the total mean density.\n", " f\"Your input mass definition '{mdef}' does not match the mass \"\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAAEWCAYAAAAtuzN2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XmcFNW5//HPM8O+CAKCILuAiqwKuCvGDXFBjYqY5LoFoz81yb03JjEal2iiMcnNjRrN1bjcJDcSl5i4BcUFF6IBlUUWUWQdFkGQxWGb5fn9cWqYnqaH6Zrpnu5hvu/X67ymq+p01VMzWg+n6tQ55u6IiIjki4JcByAiIpJIiUlERPKKEpOIiOQVJSYREckrSkwiIpJXlJhERCSvKDGJiEheUWISEZG8osQkIiJ5RYlJRETySpNcB9AQderUyXv37l2r7xYXF9O6devMBpTndM6Ng865cajLOb///vufu/t+NdVTYqqF3r17895779Xqu1OnTmX06NGZDSjP6ZwbB51z41CXczazZenU0608ERHJK0pMIiKSV5SYREQkrygxiYhIXlFiEhGRvKJeeSIikpK7U1bulJZX/iwuyf6s50pMIiIplCdckMvcKStzSsvLdy2XljmfFZezaO2XlEfLu+qWl1dZLi2v+L6HuuWVdSqXK0vpbp/Lw3eTkkTycuJxS8tDvNUt71pXVl5lH8nLyfrsU8AZp2T3d6/EJNKIlO+6cIYLYlm5U17Obut2Ffeq30lRt2L7nHWllH/02a6LbVlUtzzFvsoTLnzlCRfusqR/oVdeQMt3u2gnX5hT1U29r5ov+mXueLoNg7feyOrfLN+UZr/BpMQk+cF3Xeyg3D0qUFbuKbeF9ey6QCbWL/c9bysvD/sJF59wAS1PvgBXrE9xIU/cXlZeXvUCHP1cunwHb2yZl/B9qmwvTbpA79q3U+2FvDSpbsU5VV2XGGPl9opzyLr3a/fiueSvAoMmBQU0KTQKC4wWheVZP6YSUz1ZsHozy9YXM2dNKVtmr0q4wFZeFMudXRfL8ugitdtFueJCFV3AKi+wlRe+3S7KFRfOlMch4QK9+8W/6nFCPBVxe8LFsOJ7qRLGjh07KXzz5T0mk73S0qW5jkDqqElBuBhX/AylYNfyzh3badu6VcK25LpGk4KCGraHfRYWUH1dMwoLE+oaFBYWVN1eYDQpDMdrUpBY32haWFBlf8nLFYknebnQjIICq/I7mTp1avZ/71k/ggDwxHsreHTa0rAwa2ZOY8mJkpJcRyCRigtZQQG7LjyV6yovcgXRhbLAwncKrPKCuetzxX4KjE1fbKRTpw40SdheZX9mFBbsvq/EC3PiRb+mi3iVulb1QlxoFv0Lv6DKRTv5Ir6rbkHBru8XWPQz6YKcSmMckqg+KDHVk0Kr+T/yxq7iAmhm4bOFi0RBgVXZVhhtL0i4wJkRPlv4XHHx2/XdXfsi4TuV2yzVRTPhQh0u3ux28d7tQh19XrpkMf37HVjNhbzywlcYHbf6C3nV7bvKbkmkarJJ3lfid7IlXKRHZW3/0ngoMdWTg/Zvy2mHdmH955/TpUvn6CJKygtv1Ytl0kW4ThfluhyH6Du7H6diW+qLP7z7z3c47thjUm5LjG1vMpUVjD6ub67DEGmQGn1iMrPWwP3ATmCqu/9fNo5zwYgeXDCiR/SvysOycYi8tU9zY9/WzXIdhog0EHvlyA9m9oiZrTWzuUnrx5jZQjNbZGY/jFafBzzl7hOBs+s9WBERqWKvTEzAY8CYxBVmVgj8FjgdGAhMMLOBQHdgRVStrB5jFBGRFPbKxOTubwIbklaPAha5+2J33wlMAsYBRYTkBHvp70NEpCExT/v15obFzHoDz7v7oGj5fGCMu38zWv4GcATwA+A+YDvwdnXPmMzsSuBKgC5duhw+adKkWsX15Zdf0qZNm1p9t6HSOTcOOufGoS7nfOKJJ77v7iNqqteYOj+k6vbl7l4MXFbTl939QeBBgBEjRnht311ojO896JwbB51z41Af59yYbl0VAT0SlrsDq3IUi4iIVKMxJaYZQH8z62NmzYCLgGdzHJOIiCTZKxOTmT0OvAMcZGZFZnaFu5cC1wIvAQuAJ9x9Xi7jFBGR3dX4jMnMZgBp95Bw95yPSeLuE6pZ/yLwYj2HIyIiMaTT+WEeMRKTiIhIXdSYmNz90nqIQ0REBNhLnzGJiEjDFfs9JjMbD0wEBgAtkre7e+cMxCUiIo1UrBaTmV0M/C+wiPAe0LPA89F+NhNGUBAREam1uLfyrgduB66Jlu9398uBPsDnwNYMxiYiIo1Q3MTUH5jm7mWEkbj3AXD3LcDPCe8JiYiI1FrcxLQJaB59XgkckrDNgI6ZCEpERBqvuJ0f3gOGEEZPeBa42cxKCbO/3gz8K7PhiYhIYxM3Md0J9Io+3xx9vh8oJIxFd2XmQhMRkcYoVmJy93eBd6PPG4FxZtYcaO7um7MQn4iINDJ1no/J3XcAOzIQi4iIiEZ+EBGR/JLO6OLTgUvdfX46I43nw+jiIiLScKU7uvi26PPcLMYiIiKS1ujilwGYWVPg98BSd1+Z7cBERKRxivOMqQx4DTg4S7HklJn1NbOHzeypXMciItKYpZ2Y3L0c+ATokokDm9l3zGyumc0zs+/WYT+PmNlaM9vtNqOZjTGzhWa2yMx+uKf9uPtid7+itnGIiEhmxO2VdyNhtIfBdTmomQ0iTJ0xChgKnGlm/ZPqdDaztknr+qXY3WPAmBTHKAR+C5wODAQmmNlAMxtsZs8nFU3VISKSJ+K+x3QTYTy8WWa2EviMpF56afbKOwR41923ApjZG8C5wN0JdU4Arjazse6+3cwmRnXGJh3vTTPrneIYo4BF7r44OsYkYJy73wmcmUaMIiKSA3ET01wy0zNvLvBTM+tI6PE3ljAO3y7u/qSZ9QEmmdmTwOXAKTGOcQCwImG5CDiiuspRLD8FhpvZDVECS65zFnBWv36pGm4iIpIJcYckuiwTB3X3BWb2c2AK8CUwGyhNUe/uqKXzAHCgu38Z4zCW6tB7iGk9cFUNcT8HPDdixIiJMeIQEZEYcjbyg7s/7O6HufvxwAZCx4oqzOw4YBDwDHBLzEMUAT0SlrsDq2oZroiI1JOcjfxgZp3dfa2Z9QTOA45K2j4ceAg4A1gC/MnM7nD3m9LZP2G08/7R7cCVwEXAxWl+V0REciTuyA/zqCExxfB09FynBLjG3b9I2t4KuMDdPwUws0uAS5N3YmaPA6OBTmZWBNwStcZKzexawtxRhcAj7j4vQ7GLiEiWpD3yQ/T50kwd2N2Pq2H7tKTlEkILKrnehD3s40XgxdrGKCIi9S/WMyYzO9HMUnUqEBERyYi4nR9eBVaZ2T1mdnQ2AhIRkcYtbmIaTLiddirwtpktN7NfmNnhmQ9NREQao1iJyd3nufvN7n4wcBjwf4TRGGZE49HdkY0gRUSk8aj1e0zuPsvdb3D3fsDZQEvghoxFJiIijVLcIYl2MbMOhPePxhPGtdsG/DlDcYmISCMVKzGZ2T6EW3fjgZMIwwi9QHh59QV335HxCEVEpFGJ22JaS3jB9iXCy67PuntxpoMSEZHGK25iugr4q7tvzkYwIiIicUcXfyxLcYiIiAA5HF1cREQkFSUmERHJK0pMIiKSV5SYREQkr9TqBVsz60aY2K8DYfbZd9xds8OKiEidxX3BthC4F5hImHyvQpmZPQhc5+7lGYxPREQambi38m4DLgd+BPQmjI/XO1q+HLg1c6GJiEhjFPdW3r8BN7n7LxPWLQd+YWYOfBu4OVPBiYhI4xO3xdQZmFPNtjnRdhERkVqLm5g+JgzYmspFwMK6hSMiIo1d3Ft5dwCTzKwn8BTwGaGVdAFwItUnrbxnZn2BG4F27n5+ruMREWms4s5g+wQwBmgN/AZ4GrgHaAWMcfcn092Xmf27mc0zs7lm9riZtYgTS8J+HjGztWY2N8W2MWa2MJpd94d72o+7L3b3K2oTg4iIZE7sF2zd/WV3P4rQI29/oKW7H+3uU9Ldh5kdQOgoMcLdBxG6nl+UVKezmbVNWtcvxe4eIyTL5GMUAr8FTgcGAhPMbKCZDTaz55OKno2JiOSJWInJzG6OXq7F3cvdfW3Fe0tm1tXM4vTIawK0NLMmhBZX8gu6JwB/r2hJmdlEQuusCnd/k/CSb7JRwKKoJbQTmASMc/cP3f3MpLI2RtwiIpJFcVtMtwDdq9nWLdpeI3dfCfyS0NV8NbDJ3V9OqvMkMJnwTOtrhPekLowR6wHAioTlomhdSmbW0cx+Bww3sxuqqXOWmT24adOmGGGIiEgccROTEWawTaU78EVaOzHbFxgH9CEktNZm9vXkeu5+N7AdeAA4292/jBnrbrusrrK7r3f3q9z9QHe/s5o6z7n7le3atYsRhoiIxFFjrzwzuwS4JFp04AEzS57BtgUwGHiZ9JwMLHH3ddEx/gocDfwp6djHAYOAZwitsWvT3D+EFlKPhOXu7H67UERE8kw6LaatwPqoGLApYbmiLAHuBq5M87jLgSPNrJWZGXASsCCxgpkNBx4itKwuAzqY2R1p7h9gBtDfzPqYWTNC54pnY3xfRERyoMYWU/Ss50kAM3sU+Im7L6nLQd39X2b2FPABUArMBB5MqtYKuMDdP42OfQlwafK+zOxxYDTQycyKgFvc/WF3LzWza4GXCL3+HnH3eXWJW0REsi/WC7buflmmDuzut7CHzhLuPi1puYTQgkquN2EP+3gReLEOYYqISD3TRIEiIpJXlJhERCSvKDGJiEheUWISEZG8Endq9eaErtsHEYYBmgvMqeg5JyIiUldxp734M3AOISG1JkyrbmZWDMwDZrv7VRmNUEREGpW4ielU4Dp3vx/AzFoSRnwYklBERERqLW5iWk4Y5QEAd98GTI+KiIhIncXt/HAX8P+yEYiIiAjEn8H2j8BSM5tiZl8xs6ZZiktERBqpuL3y/hO4Jlo8CSgxs4+A2VGZE2cmWxERkWRxnzHdSJia4iZCr7whwNDo53cIE/EVZjJAERFpXOImphLgMXdfHi0vAP5SsdHM2mcqMBERaZzidn74E+EWXkruvrFu4YiISGMXNzEtAyaY2TVmplt2IiKScXFv5f2MMIHfvcBPzOxtYBZR5wcNTSQiInUVNzG1BfoSOjsMjspFhE4RBWZW7O5tMxuiiIg0JnFnsHXg06g8U7HezFoAg6IiIiJSazUmJjP7MeF23Rx3X5aqjrtvB96LioiISK2l02K6BLgVwMw2U/kybcWzpXnuviNbAYqISONSY2Jy935m1hYYBgwHrgOOr9gMlJnZQio7QPwiW8GKiMjeL63u4u6+xd3fApoD24DjgG7AMYSeevsD5xGSloiISK3F7ZV3PXCpu0+LltcA75rZvcCrhG7kIiIitRb3BdtCwntMVbj7euB24D8zEZSIiDRecRPTk8DN1YyJt5Mw1XqDZGZ9zexhM3sq17GIiDRmcRPT9UAx8ImZ3WZmx5pZTzM7iTCJ4Efp7MTMDjKzWQlls5l9N2YsFft6xMzWmtncFNvGmNlCM1tkZj/c037cfbG7X1GbGEREJHPivmC7xcyOB24gdHT4MaFnngFFwKVp7mchoZcf0Zh7K0l4YTda3xnY5u5bEtb1c/dFSbt7DLgP+EPS9wuB3wKnRLHNMLNnCbcj70zax+Xuvjad2EVEJLvidn7A3UsI4+TdQRiSqCuwntBVfGctYjgJ+DTFy7snAFeb2Vh3325mE4FzgbFJ8bxpZr1T7HcUsMjdFwOY2SRgnLvfCZxZizhFRKQe1Hgrz8xWm9nvzexcM2tTsd7dy919trtPdvcZtUxKEMbaezx5pbs/CUwGJpnZ14DLgQtj7PcAYEXCclG0LiUz62hmvwOGm9kN1dQ5y8we3LRpU4wwREQkjnSeMX2H0LJ6APjczKaY2XfNrH9dD25mzYCzCZ0qduPudwPbo2Of7e5fxtl9ql1WV9nd17v7Ve5+YNSqSlXnOXe/sl27djHCEBGROGpMTO7+hLtfSrhldwLwT+DrwEdm9rGZ/drMTjazprU4/unAB+7+WaqNZnYcYWDYZ4BbYu67COiRsNwdWFWLGEVEpB6l3SvPg3+5+y3uPoJwof850BN4GlhvZn81s8tjHH8CKW7jAZjZcOAhYBxwGdAheq6VrhlAfzPrE7XMLgKejfF9ERHJgbjdxXdx99Xu/rC7fxXoRBiSaBnw/XS+b2atCD3m/lpNlVbABe7+qbuXEwaT3W10czN7HHgHOMjMiszsiii+UuBa4CVgAfCEu8+Lc44iIlL/YvfKM7NxQH/gc2AeMNfdtwGvROXf09mPu28FOu5h+7Sk5RJCCyq53oQ97ONF4MV04hERkfwQKzGZ2YPAFYQx8toDLQmjiy8G5hC6jMe53SYiIlJF3Ft5FwI3u/sB7t6a0HIaT3hOVEi43SYiIlJrcW/lbQHerVhw94pp1qt7TiQiIhJL3BbTY8CYLMQhIiICxE9MRcA4M/u2mcXuOCEiIlKTuMnlV4Ru3P8N3GZmbxGmVJ9F6PiQPMCqiIhILHETU1ugLzCEMIDrYEKHiBuAAjMrdve2mQ1xL7FuISx9m302lcLOkdCsda4jEhHJS3GnvXBCZ4dPSZimwsxaEIYOGpTR6PYmn7wML9/EYQAzfwAd+0HXIbD/4KgMgTadcx2liEjO1ZiYzOzHhFt1c1JMTQGAu28H3ouKpLLmw4QFh/WfhDL36crVbfZPSFSDoetQ2LcPFNR6gA4RkQYnnRbTJcCtAGa2mfBMaddzJWCeu+/IVoB7jb6jobyM4sXv0nrbKvDy3et8uQYWrYFFUyrXNW0N+w+qbFXtPxg6D4SmLeorchGRelVjYnL3fmbWljDj7HDCzLXHV2wmjPywkChhufsvshVsgzbsYhh2MTOmTmX00aNg7QJYMzu0pNZ8CJ/Ng5Ktu3+vpBhW/CuUClYInQbsfiuwVYf6Ox8RkSxJ6xlTNL35W2Z2JLANOI7wnKk34b2ma4ADCQlLiakmzVpB98NDqVBeBus/hTVzKpPVmjlQvG7373sZrFsQypy/VK5v3xO6DoNuw6Kfw5WsRKTBidsr73rg0oQBVtcA75rZvcCrwL2ZDK5RKSiE/QaEMvj8yvVb1lQmqTUfwuo5sOHT1PvYuDyUBQmze7TrEZ5VdRsGXYeHn607ZfdcRETqIG5iKiS8x1SFu683s9uB24GHMxGYRNruH0r/UyrX7dgCn82PktWckKzWzoeyFLPbb1oRykfPV67bp3tCqyr62Wa/7J+LiEga4iamJ4GbzewVd9+YtG0n4daeZFvzttDziFAqlO4Mt/ZWzYLVs2D1bFgzF8pS9EvZXBRKlWR1QEKiGho+t+2S/XMREUlSm1t5LwOfmNn9wBRgOWGU8buAjzIbnqStSbMooQxl1yDvZSWw7qPKZLVqFnw2F0q37/79zStDWfhC5bq23eCAw6JyeHhm1aJdvZyOiDRecV+w3WJmxxNGergO+DGhZ54RxtG7NNMBSh0UNq3stcc3wrqykjAKRUWiWj0rtKxKt+3+/S2r4KNVVVtWHfuHJHXA4SFhdRmkrusiklGxB2KNZpL9iZndQRiSqCuwntBVPMVDDskrhU2j96IGwfCvh3VlpfD5x0nJ6sPU3dcrXgyeMyksF0T7O+Bw6Ba1rDr1D505RERqodYjhLt7OZUv20pDVtgEugwMZdjFYV15WWhZrXw/lFUfhHetykurfre8BFbNDKVCs7bhWVV0C7D59p3gDmb1d04i0mBp6gpJraCwMlkdFt0GLNkWWlIr34eVH4Sfqbqu79wCS98KBTgK4MMfVd4C7DEytK5a7FNvpyMiDYcSk6SvaUvoMSqUCls3hNbSyg9Cq6roPSheu/t3i9fCx/8IBQALQyv1GAndR0L3UWFgW40LKNLoKTFJ3bTqAP1OCgXCLbvNK6u0qkpXvEeTsuTOFQ5r54Xy/mNhVYv20H1ESFI9RobWlXoBijQ66Ywu/ghwu7sviXrkfeDuX2Y/NGmQzKBd91AGjgPg7ddfZfShB4RkVTQjtKrWztt9INvtG2HRK6GEncF+B1dtVXUaoFaVyF4u3dHFfwcsAV4nPDKYns2gZC9jhdD54FCGfy2s27EltKiKZlSWreuTvuiVYwJ+8Iewqnm7aJzBUVGyOhxa7luvpyMi2ZVOYloNjDaz+YT3lVqY2W7DElVw9xR9jEWSNG8LfU8IBcItwA2LK5PUiumhF6CXVf3ejk3w6WuhVOg8EHocAT2PCqNhtO+lHoAiDVg6ielBwqgOdxJepn29hvp6gUXiM4OOB4Yy9KKwbmdx6FixYnq4/Vc0PfVo62vnh/L+o2G5bdeERHVkeAm4UI9TRRqKdOZj+omZvQAcAvwBuIMw5YVIdjVrDb2PDQVCq+qLpZVJasX00H09uVW1ZTXM/1soECZb7DESehwZElX3EaHFJiJ5Kd35mN4H3jezk4BH3X1JdsMSScEMOvQJZcgFYd3O4pCoVvwLlr8DK2aE96gSlRTD4qmhAFhBGKapIlH1PBL26VafZyIiexB3rLzLAMysG6ETRAdgA/COu6/KfHj1x8z6AjcC7dz9/JrqS55o1rrqs6rysvBsqiJRLX83dF9P5OVh9PXVs2H6/4R17XuGRNXrKOh1bBhWSc+pRHIiVmIyswLgPmAiVZ8llZnZg8B10VBF6eyrPfB7YBDh2dXl7v5OnHii/TwCnAmsdfdBSdvGAL+JYv29u99V3X7cfTFwhZk9FTcGySMFhWHK+a5DYNTEsG7jipCgVrwbfn42j/CfXIKKSRY/fCIst94Peh0dklSvo0MHC3VTF6kXcZ8I/wS4HPgR8BfgM6ALMD7ath64Oc19/QaY7O7nm1kzkiYgNLPOwLZoWveKdf3cfVHSfh4jJMs/JH2/EPgtcAph5PMZZvYsIUndmbSPy909xXAFsldo3yOUitt/2zeFW34Viarovd1HVy9eB/P/HgqELuk9jw5JqvcxsP8QDVQrkiVxE9O/ATe5+y8T1i0HfmFmDnybNBKTme0DHE80TUY0KnnyyOQnAFeb2Vh3325mE4FzgbGJldz9TTPrneIwo4BFUUsIM5sEjHP3OwktLGmsWrSD/ieHAmEqkNVzYPk/YVlUtifNg7ntizBXVcV8Vc33CT3/eh8DvY4Jc1UVNq3f8xDZS8VNTJ2BOdVsmxNtT0dfYB3wqJkNBd4HvuPuxRUV3P1JM+sDTDKzJwkttVNS7i21A4AVCctFwBHV1MXMOgI/BYab2Q1RAkuucxZwVr9+/WKEIXmvsGn00u7hcPR1UF4eup8vmxbK0mmw9fOq39mxGRZNCQWgaaswhmCvkKgKUk1zLyJpiZuYPgYuIsxim+wiYGGM4x5GeCb1LzP7DfBDwsSDu7j73VFL5wHgwJhDIaV6cu0p1lUcaz1w1Z526O7PAc+NGDFiYow4pKEpKKics+qIb4Vu6p9/AsveDklq2bTQJT1RydYqPf+Otaaw4kjoc3wo3Q4LswyLSI3iJqY7CC2YnsBThGdMnYELgBMJySkdRUCRu/8rWn6KkJiqMLPjCJ0jngFuAa6NEWsR0CNhuTvQoHsOSo6YwX4DQhlxefQ+1ZJwy2/ptJCwNi6v8pUCL6mc+uP1n4YWVc+jokR1HHQdpmdUItWI2138CTPbCNxG6LzQFCgh3Iob4+5T0tzPGjNbYWYHuftC4CRgfmIdMxsOPAScQRin709mdoe735RmuDOA/tHtwJWEpHlxmt8VqZ4ZdOgbSsUswBtXRM+n3g4/1yf10SnZCp++GgqEMf96HwO9jwvJSr3+RHapzdTqLwMvR13HOwGfp9tFPMl1wP9FPfIWA5clbW8FXODunwKY2SVEnSUSmdnjwGigk5kVAbe4+8PuXmpm1wIvEXriPeLu82oRp0jN2veA9uNh6HgA/vnS0xzdtRyWvAFL3oSNy6rW37EJFr4YCkCrjmGEiz7HQ58TwtxUeo9KGqm6Tq1e6y7W7j4LGLGH7dOSlksILajkehP2sI8XgRdrG6NIbe1s3hGGjK7sov7FsnBbb8mboSQ/o9q6vmr39Db7Vz6f6nMc7Nu7PsMXySmNbClSH/btFcrwr4dnVOs/Da2pimSVPOXHl2vCy74VL/zu2xv6nggHnhiSlab6kL2YEpNIfTODTv1CGXlF6J6+bkHUmnoLlr4dbvUl+mJpGD39/UfDWH9dh4Uk1ffE8D6VevzJXkSJSSTXCgqgy6GhHHl1GO9v9ezK1tSyf4bOExW8HFZ9EMpbvwo9/nodU5moOh+i51PSoCkxieSbgkI44LBQjvkOlO4M03x8+josfj3MUZXY36hka9WXfdvsD31HR4lqNLTdv/7PQaQOapWY9sbRxUXyVpNmlfNSnfTjMDzSkjejRDU1vFOV6Ms1MGdSKBC6ovcdHVpTvY8JI7KL5LG4o4sXAveSgdHFRaSWWu4LA8eFAuH5U0VravEbu4/zVzHD77v3Q2GzMBBtv5ND2e9g3faTvBO3xXQbmRtdXEQyYd/eMOKyUMrLYPWsytbU8nehvKSybtnOyqGTXr4J9jkA+p0UklTf0WGAW5Ecy8no4iKSJQWFcMDhoRz/vTDD77J3Qmvq09dCyynR5pXwwR9CscLQw68iUe0/RKNRSE7kanRxEakPzVpXneJj08qQoBa9ElpVid3SvSxM/bH8n/Da7WGyxAOjJHXgV6B1x9ycgzQ6uRpdXERyod0BcNg3QikrhZXvhSS16JXQ2y9R8bqEThQW5pzqdzL0PyW0yDQIrWRJrkYXF5FcK2wCPY8M5Ss3wZfrElpTryaNRuGV7069eTe07BCS1IDTwq0/jUQhGZST0cVFJA+12S8MQjt0fBiNYvUsWPRqSFRF06u+O7VtQ+WQSVYIPY+kR2E/WLs/7HeQevpJneRydHERyVcFBZUv+Z5wfXh3avEb4SXeT14J70pV8DJYNo0DmQb3/y+07wkDxkD/08K7V01b5O48pEHK2ejiItKAtNwXDj0nlPJyWDMHPn4JPnkJVn5AlcmhNy6H6Q+G0rRV6IY+4DTofyrs0y1HJyANSY2JycymA5e6+3wzm8HRFtIsAAAYP0lEQVQepicHcPdRmQpORPJQQQF0GxbK6B/Al2vhkyms/eef6LzpQ9i5pbJuydaq807tPyS0pgacFqabV3d0SSGdFtM8YFvC5z0mJhFpZNp0huFfY/6mA+h87NGw/J3K1lTyTL5r5oTy5t2hO/qA0+CgM0KrqlmrXEQveajGxOTulyV8vjSr0YhIw9akGfQ9IZQxPwvzTn38Enw8OYySnjgKRfE6mPmnUJq0DO9KHTw2tKhad8rdOUjOxR0r72bg96kGbDWzrsBEd/9JpoITkQau44Fw1P8LZfvmMALFxy+H1lTxusp6pdtg4QuhWEEYgeKgsaF06pe7+CUn4nZ+uAWYDKQaSbxbtL1RJqaSkhKKiorYvn37Huu1a9eOBQsW1FNU+aFdu3YsWbKE7t2707Rp01yHI7nSYp/KwWfLy2Hl+yERffQifJ7wbr6Xh9uBy9+BKT+GTgNCgjr4DDhghJ5LNQJxE5NR/TOm7sAXdQun4SoqKqJt27b07t0b28M7HFu2bKFt27b1GFnubd68mZ07d1JUVESfPn1yHY7kg4IC6DEylJNvDbf8PnohdJJY/i5VLjOffxzKtP+G1p3hoDHRc6kToGnLHJ2AZFM6vfIuAS6JFh14wMw2J1VrAQwm9VBFjcL27dtrTEqNlZnRsWNH1q1bV3NlaZw6HgjHfDuU4s/DM6mPXgwjUZRuq6xXvLZy0NmmrcJzqUPODp0oWrbPXfySUem0mLYSprOA0GLaRJgcMNFO4B/A/ZkLreFRUqqefjeSttadYPjXQ9m5NUzRsfAFWDgZtn5eWa9kK3z0fCgFTUML6pCzQmuqzX45C1/qLp1eeU8CTwKY2aPA7e6+ONuBiYjQrFXoqXfw2DDXVNGMylt+iV3Ry0sqB6N9/t+h59EhSR1yVhi4VhqUuGPlXVZzLRGRLCgorBx09tTbYd3H8NFzMP/ZMK5fBS+HZW+HMvkHYST0Q84OSarjgbmLX9IWu3uLmY03s1fMbLmZrU0u2QhS0rN9+3ZGjRrF0KFDOfTQQ7nlllt2bZs8eTIHHXQQ/fr146677qpxvUje228AHPef8K034Lsfwml3Qs+jCE8cEqx8H165Be49DB44BqbeBZ/NB9dYAfkq7ntMFwOPAI8BX4k+FwBnAxuBP2Q4PomhefPmvPbaa7Rp04aSkhKOPfZYTj/9dEaOHMk111zDlClT6N69OyNHjuTss8/moIMOSrl+4MCBuT4VkXja96x8X2rLZ+GZ1PxnYelbUF5aWe+zuaFMvRM6HAgDo5ZUt8M0Inoeidtiuh64HbgmWr7f3S8H+gCfEzpKSI6YGW3atAHCe1UlJSWYGdOnT6dfv3707duXZs2acdFFF/H3v/+92vXVuf/++xk0aBC9evXi3nvvra/TEomnbRcYcTn829/ge5/AOQ+E96AKm1ett+FTePvX8NBX4DdD4OUfhwFp1ZLKubiJqT8wzd3LgDJgHwB33wL8HLg2s+FJXGVlZQwbNozOnTtzyimncMQRR7By5Up69Oixq0737t1ZuXJltetTefrpp5kyZQozZ87k3Xff5bbbbqO0tDRlXZG80aoDDLsYJjwO3/8Uzn8UDj0PmrWpWm/jcvjnPfDQiXDPMHjlVlg9W0kqR+K+YLsJqPhnx0rgEGBqtGxAx8yE1bD1/uELWdv30rvO2OP2wsJCZs2axcaNGzn33HOZO3cunuJ/LjOrdn0q99xzDw899BBNmzala9euNG3alPJyTcMlDUjztjDovFBKtofhkeY/G277bd9UWe+LpaEl9favoUNfOPTcULoM0u2+ehI3Mb0HDAFeAp4FbjazUsJ7TDcD/8pseFJb7du3Z/To0UyePJljjjmGFStW7NpWVFREt27d6N69e8r1yUpKSpgzZw4DBgwAYPXq1XTq1IkNGzYwfvx4zjjjDObNm8fRRx/NlClTuPXWWxk0aFD2T1Kktpq2gINOD6V0Z3hXat4zoSv6joQktWExvPWrUDr2q0xSnQcqSWVR3Ft5dwLLo883A9MJL9U+SnjG9K3MhSZxrVu3jo0bNwKwbds2XnnlFQ4++GBGjhzJJ598wpIlS9i5cyeTJk3i7LPPrnZ9svnz57Np0yYWL15MeXk5N9xwA9dddx0zZ87kvPPO4/vf/z6bNm1i4sSJXHDBBSxbtqy+T12k9po0gwGnwrkPwPWfwIS/wJCLoFnS0GHrF8Gbv4AHjobfjoLXf0ar4uWp9yl1Evc9pneBd6PPG4FxZtYcaO7uycMUNVp7ut2WzbHyVq9ezSWXXEJZWRnl5eVceOGFnHnmmQDcd999nHbaaZSVlXH55Zdz6KGH7nF9opkzZ/K1r32NCRMmUFxczHnnnceVV17JnXfeybnnnktJSQkdO3akoKCAuXPnMnHixKycn0jWNWkejcU3Jtzu+/S10JJa+CLs/LKy3ucfwxs/ZxTA0t+G51aDz9d7UhlS66nVK7j7DmCHmZ0IfN/dT697WFIbQ4YMYebMmSm3jR07lrFjx6a9PtGsWbM488wzGT9+fJX1ixYtYsCAAcyZM4dDDjkEgKVLl9KzZ89anoFIHmnaonLUiZJtsOjVKEn9A0qKK+ut+wim/iyUbsNh8AUhUe3TNXexN3BpJSYzaw+MAXoAS4C/u3tJtO0C4AfAYcDHWYpTcmjWrFlcffXVu61/+OGHARg2bBjDhg0D4I9//GO9xiZSL5q2hEPODGXn1jD00bxnKFvwAoXlOyrrrZoZyks3Qu9jQ5IaeDa03Dd3sTdA6YwuXjFqeJeE1R+Y2VeBPwNHAvOBrwF/yUaQkltTp07NdQgi+aNZq5BsBp7NtFcnc3yXrTD3afjkZSjbGVXy8HLv0rfghf+E/qfAoK+GzhbNWuc0/IYgnRbTz4DNwDnAbKAXcC8wg9B1/BJ3/1PWIhQRyVPlhS1g0JjQBX3bRljwHHz4JCx5k11zSpWXhGdUC1+Epq3DhIeDzw9TdhRq4sxU0klMI4DvuHtFV/CFZnY18Alw5d6SlMysL3Aj0M7dz891PCLSwLRsD4d9I5Qta8LzqA+fDGP1VSgphg+fCKVlBzj0HBh0fhjjTzPz7pLOb6ILsDRpXcXy7Noe2MyWmtmHZjbLzN6rw34eiQaQnZti2xgzW2hmi8zsh3vaj7svdvcrahuHiMgubfeHI6+Gia/BdR/AiTeFKeITbdsA7z0Cj42F/x4EU26BtQtyE2+eSbdXXnXjctR1TJoT3f3zVBvMrDOwLRruqGJdP3dflFT1MeA+kgaQNbNC4LfAKUARMMPMngUKCe9jJbrc3TUyuohkXscD4YTr4fjvwZoPQytq7l9hc1Flnc0rw9Tx0/4bug4N71ENPh/adM5d3DmUbmJ6KRrhIdmryevdPVO/yROAq81srLtvN7OJwLlAlb7N7v6mmfVO8f1RwKKKSQ3NbBIwzt3vBM7MUIwiIukxg65DQjn5NljxbkhS8/4WWk8VVs8O5eWboN9JMGR8eC7VtGXuYq9n6SSm27J0bAdeNjMH/sfdH6yy0f1JM+sDTDKzJ4HLCa2fdB0ArEhYLgKOqK6ymXUEfgoMN7MbogSWXOcs4Kx+/frFCENEJElBAfQ6OpTT7w7dz2dPCu9IlUXdz70s9PT75GVovk/oCTh0Qpiddy9/HpXO1OrZSkzHuPuq6JbdFDP7yN3fTDr23VFL5wHgQHf/MuWeUks1kFW1QwW7+3rgqj3t0N2fA54bMWKEhjYQkcwobFo5bt+2jTD/byFJLX+nss6OzTDzT6G06wlDLoShF0Gn/rmLO4tylnbdfVX0cy3wDOHWWxVmdhwwKNp+S/L2GhQRXgiu0B1YVatgRUTqQ8v2cPilcPlk+M5sOPHGMMJ5ok3L4a1fwn0j4MET4V8PQvH6nISbLTlJTGbW2szaVnwGTgXmJtUZDjwEjAMuAzqY2R0xDjMD6G9mfcysGXARYUT0vVrv3r0ZPHgww4YNY8SIEbvWa2p1kQZm395wwvdDr74rXoGR39x9BIlVH8A/rodfDYDHJ4RpPEp3pNxdQ1LnsfJqqQvwTDT3TxPgz+4+OalOK+ACd/8UwMwuAS5N3pGZPQ6MBjqZWRFwi7s/7O6lZnYtYYqOQuARd5+XpfPJK6+//jqdOnXatVxWVqap1UUaKjPoMTKU0+4Mz5xmPw4fvxRe3oUwfXzFS7wt9w3vRg2b0GCnjM9JYop6yg2toc60pOUSQgsqud6EPezjReDFWoa510icQh3YNYX66NGjU66vLjHdf//93H///WzZsoXvfe97XHfddfV2DiJCmKKjYsy+rRvCS7yzJ0HR9Mo6276AGQ+F0umgkKCGjId9dp9rLV/t3V07GiEz49RTT+Xwww/nwQdDR0dNrS6yF2rVAUZeAd+cEm73HX89tOtRtc7nC8M08b8+FP54Lsx5MgxCm+dydStv73Zru2o31Xkmpls37XHztGnT6NatG2vXruWUU07h4IMP1tTqInu7jgfCV26C0T+CZW/DrMdh/t8rp+fw8jC31KevhQkQDz0Hhl0chkLKw1t9Skx7mYqp0Tt37sy5557L9OnTsza1ekFBAf/xH/+BmdGrVy++/e1vZ/nsRGSPCgqgz/GhjP1FGFR29p9hyVvseltm5xaY+cdQ2vcK70YNvQg69Mlp6Il0K28vUlxczJYtW3Z9fvnllxk0aFDWplZ/4IEHGDduHL/61a+UlETyTfM24fnSJc/Bdz8MLaoOSTPsblwGb9wF9wyDR06HD/4A23M/GblaTNmwh9tt2Zxa/bPPPuPcc88FoLS0lIsvvpgxY8YA2Zla/bLLLks5gaCI5Jn2PcIzqOO+B0UzYNafYd5fYXvCtWr5P0N58fuhc8XQCdB3NBQU1nu4Skx7kb59+zJ7duoB37Mxtfo555zDt771LTp06MANN9xAhw4dah+8iGSfGfQYFcqYu0L38tmTwpBIXhbqlG4LY/h9+CS07RZGmRh2Mex3UL2FqcQkNapuavVx48Yxbty4HEQkInXWtEWY4HDQebDls5CIZj8OnyWMdbBlVeWo590Og2EX06SkS/X7zBAlJqmRplYX2cu17QJHXxvK6jkhQc15ArYmzEq06gNY9QFHWxM4akFWp+RQYhIRkUoVU3Oc8pNwi2/Wn+HjyVC2E4AtbfvTLsvzRCkxiYjI7hJHPd+6AeY+DbMfZ03rI6j+Tc3MUHdxERHZs1YdYNREmPgaq7uemvXDKTFlUKqRFCTQ70ZkL1EPI0UoMWVIixYtWL9+vS7AKbg769evp0WLFrkORUQaAD1jypDu3btTVFTEunXr9lhv+/btje4CvX37dtq3b0/37t1zHYqINABKTBnStGlT+vSpeaypqVOnMnz48HqIKH80xnMWkdrTrTwREckrSkwiIpJXlJhERCSvmHqRxWdm64CNQOIw4u32sJz4uROQMM5HrSUfry51q9u+p3OqaVnnrHOuLZ1z3erm8zn3cvf9aqzl7iq1KMCD6S4nfX4vG8evS93qtsc5R52zzlnnrHPOVNGtvNp7LsZy8rZsHL8udavbHucck5d1zpmhc65bXZ1z9etzfc7V0q28emZm77n7iFzHUZ90zo2DzrlxqI9zVoup/j2Y6wByQOfcOOicG4esn7NaTCIiklfUYhIRkbyixCQiInlFiUlERPKKElMOmFlfM3vYzJ5KWNfazP7XzB4ys6/lMr5sqOacd1u3N6nmnM+J/sZ/N7Psz7hWz6o550PM7Hdm9pSZXZ3L+LKhuv+Oo/+n3zezM3MVW7ZU83cebWZvRX/r0XXZvxJThpjZI2a21szmJq0fY2YLzWyRmf0QwN0Xu/sVSbs4D3jK3ScCZ9dT2HVS13Ou5veQ1zJwzn+L/saXAuPrLfA6yMA5L3D3q4ALgQbRtToD/z8D/AB4oj7izYQMnLMDXwItgKK6xKLElDmPAWMSV5hZIfBb4HRgIDDBzAZW8/3uwIroc1mWYsy0x6jbOTdEj5GZc74p+k5D8Bh1PGczOxt4G3g1e2Fm1GPU4ZzN7GRgPvBZdsPMqMeo29/5LXc/nZCQb6tLIEpMGeLubwIbklaPAhZF/7rYCUwCxlWziyJCcoIG8nfJwDk3OHU9Zwt+DvzD3T/IbrSZkYm/s7s/6+5HAw3iNnUGzvlE4EjgYmCimeX9/9N1PWd3L48+fgE0r0ssef/LauAOoLIVBCH5HGBmHc3sd8BwM7sh2vZX4Ktm9gD1PPxHhqV9ztX8HhqiOH/n64CTgfPN7Kp6jjOT4vydR5vZPWb2P8CLOYg1U9I+Z3e/0d2/C/wZeCjhot3QxPk7nxf9jf8I3FeXg2oG2+yyFOvc3dcDVyWtLAYuq5eosivOOe+2roGKc873APfUS1TZFeecpwJT6yGmbEv7nBM2PpbViLIvzt/5r4R/YNeZWkzZVQT0SFjuDqzKUSz1Reesc95b6Zzr6ZyVmLJrBtDfzPqYWTPgIuDZHMeUbTpnnfPeSudcT+esxJQhZvY48A5wkJkVmdkV7l4KXAu8BCwAnnD3ebmMM5N0zjpndM4652zEokFcRUQkn6jFJCIieUWJSURE8ooSk4iI5BUlJhERyStKTCIikleUmEREJK8oMYmISF5RYhIRkbyixCSSBjO71czczD6pZvuiaPut9RxaWszs11F8H6bY1t7MNkTbv5eBY11oZmvMzKLli8xsh5k1reu+pXFQYhJJ33agj5lVmYXVzEYCvaLt+WowsAUYEE3+luj7QLPo826JqxbOAF70ymFlhgJz3b0kA/uWRkCJSSR9xcBrhIEsE10UrS+u94jSN5gw+GYzoG/FSjPrDHybyoE559TlINGEeGOAFxJWDwNm1mW/0rgoMYnEMwm4MOE2lQEXRut3MbOjzOxZM1tlZsVmNsvMdpu91cwONbPJ0a20YjNbYGbX1LQtjij5dAaeJ7SaDk7YfBMhGS0DPnf31XH3n2QksC8wJWHdUGCBmf3MzFaa2SYze6ghzOoquaH/METi+SvQBTg2Wj4O2A94JqleL2Aa8E3gLOBp4FEzm5BU71mgDPg6cDZwL9A2jW1xDIl+zgHmA4cAmFlP4FvAj6I6mbqN95a7b46OsR/QlTBzb0vgUuCXhN9LtVOxS+OmGWxFYnD3jWY2mXD77q3o5+RofWK9XS2oqFX1JmGStYnA49H6ToTbaue4e0VSeLWmbbUwBNgBfAzMI0pMwK3Am+4+1cz+QGZmHz0D+L+E5WHRz3vc/b+iz1PM7GqgfwaOJ3shtZhE4psEnG9mzYHzSbqNB2Bm+5rZPWa2DCiJypXAgIRqG4AVwO/MbHx0yy2dbXENBhZEc+vMAw42s4OBbwA3mlk7wiyldX2+1BUYTtXnS0OBjYTWXkU9A9oDn9fleLL3UmISie9ZoA3wU6A18FyKOo8B44FfAKcSnr08ArSoqODu5dG2NdG2NWb2lpkN39O2WsQ7mMrbdBUtptuBF9x9erQdqrmVZ2YDzewRM3vdzO4zs27VHGcssNjdFyasGwq8kdQj70DCbb29ZpI9ySwlJpGY3L2Y0JHg34HnouVdzKwF4ZbWLe5+n7u/5u7vkeL/N3f/yN2/SmhBnExIXC+YWcGetqUba1R3IJVJZy7QDjiP0PEBwq2+clIkCjP7CvAQ8CDhOdfzwPNRiyvZGVRtLUG4lZfcEhsaHW9uuuchjYueMYnUzgNAc+B3KbY1BwoJz3UAMLO2hAt7yimjoxbFa2b2X8CfCcloQ03b0tCf0Dr5MNrXSjN7Apjv7hWJYTCwyN23pvj+TcD5Cb31JpvZhmj91xPOrxkheZ6ftO5gYHbSPocQWlb53L1eckiJSaQW3H0qMLWabZvMbAZws5ltJrQOfghsAvapqGdmQwg91P4CLCZ0s/4B4ULe3cwmpdrm7hui748GXgdOjOJJZbfbdO4+PkWd6nrkNXX31WbWG/jI3Vu4+/RoOdHxhBbhGwnrDiVcY5JbTENSrBPZRYlJJDsuJtz++gOwHrgPaAVcm1BnDfAZcCPQjdBJ4HVCAtq+h20VWkU/1+4hjsHAF+6+cg91BgEvV7PNzKwlsJKo27mZ9WL3FtsZwCvuviNh3VBgK/BpUt0hhN+LSEpWOWqIiDQkZnYbcLy7n5jFY3yTcIvuSnffHHV8eBz4ibu/mlDvY+AX7v5QtmKRxkMtJpGG62jgv2qsVQfu/nsz2wK8Gg3C+iVwo7u/kVRvQModiNSCWkwiIpJX1F1cRETyihKTiIjkFSUmERHJK0pMIiKSV5SYREQkrygxiYhIXlFiEhGRvKLEJCIieeX/AwwAGFjXhfLyAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Default mass function object.\n", "mf = hmf.MassFunction(hmf_model=\"SMT\", Mmax=15, disable_mass_conversion=False)\n", "dndm0 = mf.dndm\n", "\n", "# Change the mass definition to 300rho_mean\n", "mf.update(mdef_model=\"SOMean\", mdef_params={\"overdensity\": 300})\n", "plt.plot(mf.m, mf.dndm / dndm0, label=r\"300 $\\rho_m$\", lw=3)\n", "\n", "\n", "# Change the mass definition to 500rho_crit\n", "mf.update(mdef_model=\"SOCritical\", mdef_params={\"overdensity\": 500})\n", "plt.plot(mf.m, mf.dndm / dndm0, label=r\"500 $\\rho_c$\", lw=3)\n", "\n", "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "plt.xlabel(r\"Mass, $M_\\odot/h$\", fontsize=15)\n", "plt.ylabel(r\"Ratio of $dn/dm$ to virial\", fontsize=15)\n", "plt.grid(True)\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, the measured mass function uses \"virial\" halos, which have an overdensity of ~330$\\rho_m$, so that converting to 300$\\rho_m$ is a down-conversion of density. Finally, we convert a FoF mass function:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/steven/Documents/Projects/halos/HMF/hmf/src/hmf/mass_function/hmf.py:291: UserWarning: Your input mass definition 'FoF(l=0.3)' does not match the mass definition in which the hmf fit Jenkins was measured:'FoF(l=0.2)'. The mass function will be converted to your input definition, but note that some properties do not survive the conversion, eg. the integral of the hmf over mass yielding the total mean density.\n", " f\"Your input mass definition '{mdef}' does not match the mass \"\n", "/home/steven/miniconda3/envs/hmf/lib/python3.7/site-packages/halomod/halo_exclusion.py:17: UserWarning: Warning: Some Halo-Exclusion models have significant speedup when using Numba\n", " \"Warning: Some Halo-Exclusion models have significant speedup when using Numba\"\n", "/home/steven/Documents/Projects/halos/HMF/hmf/src/hmf/mass_function/hmf.py:291: UserWarning: Your input mass definition 'SOVirial' does not match the mass definition in which the hmf fit Jenkins was measured:'FoF(l=0.2)'. The mass function will be converted to your input definition, but note that some properties do not survive the conversion, eg. the integral of the hmf over mass yielding the total mean density.\n", " f\"Your input mass definition '{mdef}' does not match the mass \"\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaYAAAEWCAYAAAAtuzN2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl8XVW9///XJ2matBmbpC1t0wlahkKBQimTQFGuIiooFxn0fqUgOFxUvFwHLipSRMUvfPWnOIICTtiLXFFULoNApTIUaG1LSymUTqQttEmbNEOTZvj8/tg77cnJSXJ2cpKe5Lyfj8d+nHP2XuectRvIO2vttdcyd0dERCRdZB3sCoiIiMRSMImISFpRMImISFpRMImISFpRMImISFpRMImISFpRMImISFpRMImISFpRMImISFpRMImISFoZkaoPMrMcYIK7b0nVZ6ar8vJynzZtWp/e29DQQH5+fmorlOZ0zplB55wZ+nPOy5Ytq3L3sb2VSyqYzOwa4DpgHPAK8EN3/3VcsROAZ4HsiHUdcqZNm8ZLL73Up/cuXryY+fPnp7ZCaU7nnBl0zpmhP+dsZpuTKddrV56ZXQrcATwPLAS2Afea2QNmNqpPtRMREelGMteYvgDc7u4fdffb3f1DwLuBdwBPmVnZgNZQREQySjLBdATwcOwOd38COAUoBp4zs8MGoG4iIpKBkrnGVAuUx+90901mdhrwV4JrS99Icd2GlJaWFiorK2lqauqxXHFxMWvXrh2kWqWHVJ9zXl4eFRUV5OTkpOwzRSR9JBNMy4APAg/EH3D33Wb2rvDYD4CMXXWwsrKSwsJCpk2bhpl1W66uro7CwsJBrNnBl8pzdneqq6uprKxk+vTpKflMEUkvyXTl/QY41MxKEx10973A+cDPgWE/VLw7TU1NlJWV9RhK0n9mRllZWa8tUxFJLXfntbfreKuhfcC/q9cWk7v/Hvh9L2XagE+kqlJDlUJpcOjfWWRwVNc384/1VSx5vYolr+/k7T3NnFkxgksH+HtTdoOtiIgMbc2tbSzbtJunwyBas21PlzJrqtpw9wH9AzGVMz9MBiwTZn4QERkO3J3Xd9Tz9Gs7+cf6Kp7fUE1TS/dddUV5I5he7DTuayM/d+DaNan85A0E16yG/cwP6Sw7O5vZs2fvf/3HP/6R7qZPSqbs4sWL+cUvfsGvfx0/0Ud0jzzyCNdeey1tbW1cddVVXH/99Z2ONzU1ceaZZ9Lc3ExraysXXXQRCxcu7Pf3isgBO+uaefaNzt1z3cnOMk6YUsIZM8dyxsxyZk8q5h9Lnh7QUILUBtPHAXX+H2SjRo1ixYoVKSu7YsUK5syZ0+96tbW1cc011/D4449TUVHBSSedxPnnn8+sWbP2l8nNzeXJJ5+koKCAlpYW3vGOd/De976XU045pd/fL5Kp6ppaWLphF8+8UcWz66tZ93Zdj+Wnl+fzjhnlnDGznFMPK6Mwb/Bvy0hZMLn7r1L1WZJa3/3ud7n77rsBuOqqq/j85z+f9HtXrlxJaWkpJ598Mjt37uTuu+/u0zxZL7zwAjNmzODQQw8F4NJLL+VPf/pTp2AyMwoKCoDgvrCWlhYNdBCJqLm1jeWba3j2jSqeWV/Fyspa2tq7v5OnKG8Ep88o398qmlw6ehBrm5gGPwyAadf/dcA+e9Ot7+vx+N69ezn++OMBmD59Ol/96le55557WLp0Ke7OySefzFlnncWcOXO6lH3wwQe7fN6KFSu44IILWLp0KY899hhf+9rXWLJkSacyZ5xxBnV1Xf8Ku/322znnnHMA2Lp1K5MnT95/rKKigqVLl3Z5T1tbGyeeeCLr16/nmmuu4eSTT+7lX0Qks7W1O69s28MzYRC9uGlXj9eJcrKNOVPGcPph5ZxxeDnHTipmRHZ6rYCUdDCZWQFwFnAkMIbgZtoa4FXg7+5ePyA1lEjiu+e+//3v86EPfWj/NPUXXnghS5YsYc6cOb125bW2tlJdXc0NN9wAwPHHH09VVVWXcvFBlYh717/YErWGsrOzWbFiBTU1NXzoQx9i9erVHHPMMb1+vkimcHc2VjXwzBvVPPN6Fc9tqKZ2b0u35c1g1oQiTp9RzmmHlTFveimjR6Z3m6TX2lnw22MhwbIXo4FGYDfB9aRiIB9oNLP/B9zkiX4DyUHTnx/HK6+8wowZMxg5ciQAy5cv57jjjutSLpkWU0VFBW+++eb+Y5WVlUycOLHb7y4pKWH+/Pk88sgjCibJeG/VNvHchir+8Xo1z75Rxfbanm8wn1Y2mtNmlPOOGeWcemgZY/JHDlJNUyOZ2LyJIJQWAovc/c3Yg2ZWAVwKfJ2gFXVTaqs49PTU3TbYUxKdeeaZLFiwgOuvvx5358EHH0x6hN3KlSvZuHEjzc3NtLS0sHDhQr73ve91KZdMi+mkk07i9ddfZ+PGjUyaNIlFixZx3333dSqzc+dOcnJyKCkpYe/evfztb3/jy1/+cnInKjKMvL2niec3VIfbLjZWNfRYvrwgl9NnlHH6jHJOn1HOpJKhvSJRMsF0FXCdu9+Z6KC7VwK3m9kegnC6KXXVk/464YQTWLBgAfPmzQOCwQ/JjrJbuXIlH/3oRznttNPYu3cvX/va1/o8Qm7EiBH88Ic/5D3veQ9tbW1ceeWVHH300QCcd955/PznP6eqqorLL7+ctrY22tvbufjii3n/+9/fp+8TGUoOBNEunt9Q3WsQFeaO4ORDy/aH0cxxBcNqoFAywVQCvJFEuTfCsnIQ1dd3vdR33XXXcd111yVVNtbtt98OwC233JKSup133nmcd955XfY//HCwqsrEiRP55z//mZLvEklnO/Y08fzGXTz3RjVLN1SzoZcgyh2RxQlTxnD6jDJOm5GeAxZSKZlgeh74kpk97+4J//XMLB/4MvBcKisnIjIc7Khr2t8aen5DNRt29hxEI0dkceKUMZxyaBmnHFrK8VNKyB2ROXMXJBNMnwH+Bmwxs0cJRuHVEFxPKiEYpfceoBl41wDVU0RkyNhR18TSmCB6I4kgOmFKCaceWs4ph5Zy3OQS8nIyJ4jiJTO7+FozOxr4NHAuQfiMCQ/vJgiq24GfunvNQFVURCQduTtbdjWydOMuXty4ixc37WJTdWOP7+kIoqBFVMbxGR5E8ZIazB4GzrfDTUQkY7W3O6++VceLm3bx1xVNfOmZJ9hR1/18cwAjs7OYExNEc6YoiHqS3ndZiYgcZPta23l5aw0vbNzNi5t28dKmXexpao0p0dblPSNHZHH85JL914hOmDJGQRSBgklEJEZDcyvLt+zmxY27eGHTLla8WdPjFD8QDN8+cdoYTppWysnTS5ldUZxRgxVSTcE0jHzzm9/kvvvuIzs7m6ysLH72s59x8skns2/fPr70pS/x5z//maysLGbNmsWPfvQjKioqOr1/wYIFnHrqqXzyk5/cv++Pf/wjd955Jw8//DCnnXYazz77bMLv7ulYh4KCgl6HqIsMtp11zSzbvJuXNgXXh1Zv29PjpKcQ3NA6b/oYxrTu4iP/Mo8jDykiO2v43Ed0sCmYhonnnnuOv/zlLyxfvpzc3FyqqqrYt28fADfccAN1dXW89tprZGdnc88993DhhReydOnSTjflXXbZZdx6662dgmnRokVcdtllAAmDp62tjezs7F5DSSQdtLc7r+2oY9nm3fu3zb0MVACYUjqaedNLmTetlJOmlzKtbDRmxuLFizl6YvEg1DyzRA4mM8sG9gFzgVXh85PcfXmK6yYRbN++nfLycnJzcwEoLy8HoLGxkXvuuYeNGzeSnR10LVxxxRXcfffdPPnkk7zrXQdG+J9zzjksWLCA7du3M2HCBBobG/nb3/7GXXfdBRxo8SxevJiFCxcyYcIEVqxYwSuvvLL/WH19PRdccAG7d++mpaWFW265hXe+852D/K8hEqhvbmXlmzW8tGk3y7bs5p9bdlPX6fpQV2ZwxPhC5k0v5aRppcybXsr4orxBqrFA31tMxoFFAdV+jXdT939B9XuWvJtqE+5+97vfzc0338zhhx/OOeecwyWXXMJZZ53F+vXrmTJlCkVFRZ3Kz507lzVr1nQKpuzsbC688ELuv/9+rr32Wh566CHOPvvshHP7vfDCC6xevZrp06d32p+Xl8eDDz5IUVERVVVVnHLKKSxfrr9ZZOC5O1tr9nZqDa3dvodeeuUYOSKLYycVc+K0McybVsrcqaUUjx78xfHkAHXlDRMFBQUsW7aMJUuW8NRTT3HJJZdw6623MmfOnIRzaLl7wv2XXXYZX/ziF7n22mtZtGgRH/vYxxJ+37x587qEUsfn3nDDDTz99NNkZWWxdetWduzY0SUYRfqrpa2dNdv2sGzzbpZv3s1Lm3f1uEx4h/KCXOZOHcPcaWM4YeoYjp5YpIEKaUbBNIxkZ2czf/585s+fz+zZs/nlL3/Jhz/8YTZv3txlVvPly5fzgQ98oMtnnH766Wzfvp2VK1fy7LPPsmjRooTf1bG+U7zf/va37Ny5k2XLlpGTk8O0adNoaup5in6RZOzY08TyLTWseLOG5Vt2s6qy99FyZnDkIUWcOLWEE6eOYe7UUirGjBpWE54ORwqmgdBNdxsM3LIX69atIysri5kzZwLByrNTp04lPz+fyy+/nOuuu46f/vSnZGdn86tf/YrGxsaE137MjIsvvpjLL7+c8847j7y8aH3rtbW1jBs3jpycHJ566ik2b96ckvOTzNLU0sbqrbX8Mwyif27ZzbZe1iACKMgdwZwpJZwwJWgRHT+5hMI8dcsNNQqmYaK+vp7Pfvaz1NTUMGLECGbMmMGddwYrlXz729/mC1/4AocffjhZWVkceeSRPPjgg93+1XjZZZdx2223ceutt0aux0c/+lE+8IEPMHfuXI4//niOPPLIfp2XDH8dK7IGARQE0drte2jt7eIQMLl0FHOnlnLC1DGcOGUMRxxSqGHbw4CCaZg48cQTux2ynZubyx133MEdd9yR1GfNmTMn4cq3HfcgdXQXJjpWXl7Oc891nmS+Y3Vb3cMkALWNLayoDFpB/9xSw8rKGmoau18avMOonGxmVxQzZ3IJx08OuubGabTcsKRgEpEB09LWzrq36vhn2B234s2aXpd86HDY2HzmTAm64+ZMKeGI8YXDeg0iOUDBJCIp0e7O+h31rKqsYVVlLSsra3hl2x6aW3seoABQMjonbAmNYc6UEo6rKNGQ7QwWOZjcvc3MrgA2xj5PfdVEJF25O5W797KqspZVlUF33IrNjTQ9+vde3zsiy5g1sSgIoiklzJk8hqnhTAoi0McWk7v/MtHzTNfdvUGSWomuf8nA2lHXxKo3O0Kolpe31rKrYV9S751UMioMoKBL7uiJxZppW3qkrrwUycvLo7q6mrKyMoXTAHJ3qqurIw9jl+TVNO7j5a21QXfcm0G33Ft7krsXrbxgJMdWlDB7UjHHTS5m9qQSxhbmDnCNZbhRMKVIRUUFlZWV7Ny5s8dyTU1NGfdLNdXnnJeX12VmdOmb2sYW1myrZfW2Wl7euodVlTVJTWoKUJg3gmMrgvA5rqKYxsq1XHju2frDTPpNwZQiOTk5Cafoibd48WLmzJkzCDVKH5l4zumoqr6ZNdv2sHprbbBtq+XNXXuTem9eThbHTCzm2IoSjq0o5tiKYqaV5ZMVc8/Q4up1CiVJCQWTyDDj7ry9p3l/+Kzeuoc122rZnsTMCQA52caRhxRxbEUxx1WUMLuimJnjCjRUWwZNX5a9uAL4IDAW2Ar8DfiNuyd3c4KIpEzH6Lg1YQB1BFFVfe+TmUIQQkccUsgxE4s5emIRx1aUcOSEQk1qKgdVpGAys28B1wPPAq8Ak4HvAN8wsyvd/S+pr6KIALS1B1P3vLJ9TxhEQQjV7u191gSA3BFZzJpYxDETizlmUhFHTyzm8PGFjByhlpCkl6gtpiuBb7v7Vzp2mFkR8B/AA2Z2gbs/msoKimSi+uZWXt2+h7Xb9/DK9j28sr2OdW/t6XU27Q75I7M5emIxR0/qCKJiDhubr+44GRKiBpMBj8fucPc9wEIzGwV8E1AwiSTJ3dlW28TabUEAdQRRsiPjAIpH5XBMGEBHTyrmmIlFXQYmiAwlvQaTmeW4e0dfwa+AfwEWJyj6CPC51FVNZHjZ19rO6zvqeGXbHtZur+OV7bWs3V6XdFccwPiiXI6aUMSsCUXMnhS0hLS+kAw3ybSYGsxsDfBP4FXg382sAfiBu8dOF30WsGoA6igy5NTtc55ZXxW0gMLW0Pod9Ukt5QCQnWXMGFvArIlFHDWhkFkTijlqQiFlBbpZVYa/ZILpI8CxwPHANQQDHm4BvmhmzwGbgcPC7dwBqqdIWtq7r43Xd9Tx6lt1rAu3V9+qC0bFPbk0qc8ozBvBrAlFQUtoYtAamjGuQNP2SMbqNZjc/QHggY7XZlZCEFLHhdvJwCxgJLAcSP3yrCIHWVu7s7m6YX/wrHurjnVv17G5uoEkG0EATCkd3akFNGtiEZNK1BUnEqsvs4vXEFxjWtyxz8yygaMIgkpkyHJ3dtY3d2r9rHurjtd31CU9Ig5gZBYcObG4U0voyEMKtcy3SBJSMvODu7cBq8NNZEhoaG7ltbfrurSCkp01GyDLYFpZPkccUsgRhxRy5CGFHHFIERtffoF3nv2OAay9yPClKYlk2GtobuX1HfW8/nbd/sfX3q5na01y88R1GFuYGwTP+I4QKmLm+MTXgjara06kzxRMMmzUN7d2Cp/gMXoAjR6ZzeHjO1o/B0KoNH/kANVcRGIpmGTI2dPUwvod9ax/u57XYoJoW5KTlHYYkWVMKw+64Y6MaQVVjBmlm1NFDiIFk6St2r1BAHW0fl57u471O+qTniW7w4gsY3p5PjPHFzBzXCEzxxdw+PhCppXla544kTSUzMwPjwGfdfd1MfveCSzVjOLSX+3tzrbavbyxs4E3dtTzxs6OrYGddcnNkN1hRJZx6Nj8/eEzc1whh48vYKoCSGRISabFdA5Q3PEiHBr+OHASwX1LIr1qamljY1VDEDo7GvYH0IadDextaYv0WTnZHS2gQmaOC1o/HQGUo0lKRYa8vnblqQNeEqqubw5aPzvr97eAVm9ppOrRR/AIN6ICjMzO6tQFd/j4AmYqgESGPV1jksha2tp5c1djwhbQ7sbkJyTtUJo/ksPG5nPY2IJgGxc8rxgzmmwNQhDJOMkGU6K/dSP+/StDSXu7s31PExt3NrCxqp6NVY3hYwNv7t5LW5R5eAhuRJ1cOjoMnzCExgVBpGHYIhIr2WB61Mxa4/Y9kWAf7j6u/9WSweDuVDfsY1NVAxuqGthY1RAGUQObqhtobk1+Cp4Oo3Ky97d4ZoThs2vTWi469yxNSioiSUkmmBYOeC1kQNU1tbCpqpENYYtnUxhCG6oaqGvq8rdFUiYW5zF9bD7Ty/P3B9BhYws4pCivyz1Ai6vXKZREJGnJzC6uYBoC6ptb2VzdwObqRjZXN3YKn6r6aMOuO5Tmj2R6ef7+7dDyfKaV5zOtLJ9RIxU0IjIw+jT4wcwmAqcCpcAu4Dl335bKikln7s7uxpb94bOpuoEtHY+7GqmqT37i0Vj5I7OZPjYIm0PL88NWUAHTy/IpHq2ZsEVk8EUKpvAepjuAq4HYP5nbzOxOghtxo1+YECAYcLCjrrlT6Gze1bg/jPra7ZaTbUwt69zq6Xg+tjBXawGJSFqJ2mJaCFwJ3AD8N/A2MB64BLgZqAZuTGUFh5u2du8cOlWdw6cvAw4gCJ/JpaOZWjqaqWX5TC0bHYZPARNL8hih+35EZIiIGkwfA77q7rfH7NsC3GZmDnwOBVNCtz+6jr+s2sabuxppe+ypPn3G6JHZTCkdzbQweKbufxzNhOJRuudHRIaFqME0DljVzbFV4XFJoGbvPjZVN/ZarmR0ThA4paOZVjaaKWX54eNoxhao201Ehr+owfQacCnwWIJjlwLrEuwXglVOO4wvymVq6YHWzv6WT6kGHIiIRA2mW4BFZjYFeIDgGtM44MPA2QThJAmcf9xEzpg5lk1rXuI97zr7YFdHRCRtRQomd7/fzGoIBkF8H8gBWoBlwLnu/njqqzg8jCvKY1xRHttfVVeciEhPIt/H5O6PAY+ZWRZQDlRpiLiIiKRKn2cXD8NoRwrrIiIiQq83t5jZY2Z2RNy+d5pZfnfvERER6atk7rrsbgXbI7p9h4iISB/1dToAXcEXEZEBoXlqREQkrSQbTFrBVkREBoVWsBURkbSiFWxFRCStaAVbERFJKxr8ICIiaUXBJCIiaUXBJCIiaUXBJCIiaUXBJCIiaaVPs4ub2UTgVKAU2AU85+7bUlkxERHJTJGCKZzA9Q7gaiA75lCbmd0JfFZrM4mISH9E7cpbCFwJ3ABMA0aFjzeE+29KXdVERCQTRe3K+xjwVXe/PWbfFuA2M3Pgc8CNqaqciIhknqgtpnHAqm6OrQqPi4iI9FnUYHoNuLSbY5cC6/pXHRERyXRRu/JuARaZ2RTgAeBtglbSh4Gz6T60REREkhIpmNz9fjOrIRgE8X0gB2gBlgHnuvvjqa+iiIhkksj3Mbn7Y8BjZpYFlANVGiIuIiKpEukak5ndGN5ci7u3u/uOjlAyswlmphF5IiLSL1EHP3wdqOjm2MTwuIiISJ9FDSYDvJtjFcDu/lVHREQyXa/XmMzscuDy8KUDPzGzPXHF8oDZwGOprZ6IiGSaZAY/NALV4XMDagkmbo21D/hf4Mepq5qIiGSiXoPJ3X8P/B7AzO4Bbnb3jQNdMRERyUxR72O6YqAqIiIiAlooUERE0oyCSURE0oqCSURE0oqCSURE0krUpdVzgSuAIwiGjK8GVrn7GwNQNxERyUBRJ3G9D/ggQSDlEyyrbmbWAKwBVrr7p1JaQxERyShRg+ndwGfd/ccAZjaKYMaHY2M2ERGRPosaTFuA/TfXuvte4IVwExER6beogx9uBf59ICoiIiICEYPJ3X8NbDKzx83snWaWM0D1EhGRDBV1VN5/AteEL98FtJjZq8DKcFul5dVFRKQ/ol5j+grwG+CrBKPyjgWOCx+vBSYB2amsoIiIZJaowdQC3OvuW8LXa4H/7jhoZiWpqpiIiGSmqIMffkPQhZeQu9f0rzoiIpLpogbTZuAyM7vGzNRlJyIiKRe1K+9bwGjgDuBmM/sHsIJw8IOmJhIRkf6KGkyFwKEEgx1mh9ulBIMissyswd0LU1tFERHJJFFXsHXgjXB7sGO/meUBx4SbiIhIn/UaTGb2NYLuulXuvjlRGXdvAl4KNxERkT5LpsV0OXATgJnt4cDNtB3Xlta4e/NAVVBERDJLr8Hk7jPMrBA4HpgDfBY4s+Mw0GZm6zgwAOK2gaqsiIgMf0kNF3f3OndfAuQCe4EzgInA6QQj9Q4BLiQILRERkT6LOirvi8ACd38mfP0W8LyZ3QE8QTCMXEREpM+i3mCbTXAfUyfuXg18A/jPVFRKREQyV9Rg+j1wYzdz4u0jWGpdRESkz6IG0xeBBuB1M1toZu8wsylm9i6CRQRfTXkNRUQko0S9wbbOzM4E/otgoMPXCEbmGVAJLEh1BUVEJLNEHfyAu7cQzJN3C8GURBOAaoKh4vtSXD8REckwycz8sB34a7g97u71AO7ezoGbbUVERFIimWtM1xIE2E+AKjN73Mw+b2YzB7ZqIiKSiXoNJne/390XEHTZnQU8C/wb8KqZvWZm3zOzc8wsZ2CrKiIimSDpUXkeWOruX3f3uUAF8B1gCvA/QLWZ/cHMrhyguoqISAaIOlx8P3ff7u6/cPd/BcoJpiTaDHwpVZUTEZHME3lUnpldAMwEqoA1wGp33wv8Ldz+I6U1FBGRjBIpmMzsTuDjBHPklQCjCGYX3wCsIhgyfkvKaykiIhkjalfexcCN7j7J3fMJWk6XAL8jmEfv8hTXT0REMkzUrrw64PmOF+7escz6H1JZKRERyVxRW0z3AucOQD1ERESA6MFUCVxgZp8zs8gDJ0RERHoTNVz+H8F6TP8fsNDMlhBMSbSCYODD+hTXT0REMkzUYCoEDgWOJZjAdTbBgIj/ArLMrMHdC1NbRRERySRRl71wgsEObwAPduw3szzgmHATERHps2RmF/8aQVfdKnffnKiMuzcBL4XbkGJm+cCPCVbgXezuvz3IVRIRyWjJDH64HPgjsMHMdpvZYjP7vpldYWYnmFnuANcxMjO728x2mNnquP3nmtk6M1tvZteHuy8EHnD3q4HzB72yIiLSSa8tJnefYWaFwPHAHIKVa8/sOEww88M6wrWZ3P22gapsBPcCPwR+1bHDzLKBHwH/QjC68EUze4hgMtqXw2Jtg1tNERGJl9RwcXevc/clQC6wFzgDmAicDnwLOISg5fHZAapnJO7+NLArbvc8YL27bwhX2l0EXEAQUhVhmT5PaisiIqlhwXiGJAub7QAWuPvDcfvLgCeAO9z9F6mtYt+Y2TTgL+5+TPj6IuBcd78qfP1/gJOBLxO0rpqAf3R3jcnMPgF8AmD8+PEnLlq0qE/1qq+vp6CgoE/vHap0zplB55wZ+nPOZ5999rJw2aQeRR0unk1wH1Mn7l5tZt8AvgGkRTAlYAn2ubs3AFf09mZ3vxO4E2Du3Lk+f/78PlVi8eLF9PW9Q5XOOTPonDPDYJxz1K6r3wM3mllJgmP7gGn9rtHAqQQmx7yuALYdpLqIiEg3ogbTF4EG4HUzW2hm7zCzKWb2LuBW4NWU1zB1XgRmmtl0MxsJXAo8dJDrJCIicSIFk7vXEYzIuwP4d+BpYCPwOMGsEJ9MdQX7wsx+BzwHHGFmlWb2cXdvBT4DPAqsBe539zUHs54iItJV5IlY3b0FuNnMbiGYkmgCUE0wVHxfiuvXJ+5+WTf7HwYeTnRMRETSQ59nCHf3dsJ7l1JXHRERyXS6b0dERNKKgklERNKKgklERNJKr8EUTog6PXx+ppll1m3OIiIyqJKdXXxs+PwpYNbAVUdERDJdMqPytgPzzewVgml98sysy7REHdy9MVWVExGRzJNMi+lOglkdagmWuXgKqOthExER6bNk1mO62cz+ChxFsL7RLQRLq4uIiKRcUjfYuvsyYFk4J9497r470mKiAAAVVElEQVRxYKslIiKZKtLMD+5+BYCZTQROBUoJFuR7zt01U7eIiPRbpGAysyyCRfWuJlibqUObmd0JfDacqkhERKRPot5gezNwJXADwdpLo8LHG8L9N6WuaiIikomiTuL6MeCr7n57zL4twG1m5sDngBtTVTkREck8UVtM44BV3RxbFR4XERHps6jB9BrByq+JXAqs6191REQk00XtyrsFWGRmU4AHgLcJWkkfBs6m+9ASERFJStTh4vebWQ2wEPg+kAO0AMuAc9398dRXUUREMklfllZ/DHgsHDpeDlRpiLiIiKRKn9djcvd2d98xXELJzA41s1+Y2QMHuy4iIpnsoC0UaGbXmtlqM1tjZp/vx+fcbWY7zGx1gmPnmtk6M1tvZtf39DnuvsHdP97XeoiISGoclGAys2MIZo+YBxwHvN/MZsaVGWdmhXH7ZiT4uHuBcxN8RzbwI+C9BGtIXWZms8xstpn9JW7TMHcRkTRxsFpMRwHPu3uju7cCfwc+FFfmLOBPZpYHYGZXAz+I/yB3f5pgvr5484D1YUtoH7AIuMDdX3b398dtO1J4biIi0g8HK5hWA2eaWVm46OB5wOTYAu7+e+ARguHpHyWY8ujiCN8xCXgz5nVluC+hsC4/BeaY2X91U+YDZnZnbW1thGqIiEgUkUflQf9nF3f3tWb2HeBxoB5YCbQmKPd/zWwR8BPgMHevj1LNRF/dQ52qgU/1Uu8/A3+eO3fu1RHqISIiEURqMZlZtpn9GNgM/B74Wfi42cx+FA4hT4q7/8LdT3D3MwnC7fUE33cGcAzwIPD1KHUlaCHFtsIqAC3NISKS5qJ25S0kRbOLdww4CGeRuBD4XdzxOcBdwAXAFUCpmd0Soa4vAjPNbLqZjSSYleKhCO8XEZGD4GDOLv4/ZlZGMHPENe6+O+74aODD7v4GgJldDiyI/xAz+x0wHyg3s0rg62FrrNXMPgM8SrB21N3uvibJuqXesl/Cqv/mmPp9UP0bGFkAuYWQWxQ+9rDl5EPWQRvZLyIyqKIGU8pmF3f3M3o5/kzc6xaCFlR8uct6+IyHgYeTrdOAqn4dNj9DOUD1ixHfbEFA7Q+z2K0j2AoS7CvsGoAjRg7AyYmIpE7UYOqYXfyxBMc0u3hPmuv68WaH5j3B1p+PAcjOjQmwgu5bbCNjnucVQ15RUDavOHhUwInIANHs4oPl9M/DMRfx8rJnmT1zWhg0dcG2r/7A89j9zXXQXA8tDamrR1szNDZDY1X/PmfEqJiwigms2AAL95VVbYGN2Z0DLrcIsvs0KFREhjnNLj5YSqdD6XSqN7fBcfOjvbetNS686mBfXVyAxQZbfTf768DbUnM+rXuhfi/Uv91r0dkQ3LkWLyc/rjUWG2pxLbRE5XKLdO1NZBjS7OJDQfYIGFUSbP3hDq1N3bTO6hO01sLXTXuguRaaasPneyAVP/KWhmCr6+sofovpaiyGvJLgcVRJ3PPwWPzzEXlgiW53E5GDqc99KWEYaSqfocQMckYFW0E/pgd0h30NQUDFhlVTbYJ9e6jeuoGygpxO+2jeQw/3OydbkQPX3mrf7L14vOyRiQOrx5ArOdCKU2tNZED0Gkxm9gKwwN1fMbMX6eW3ibvPS1XlJE2ZhQMnCqBoYq/FX168mPnz53fe2d4edEfGhlovAddl374oE4Ek0LYPGnYEW1SWdaCLMUELbcr23fDi+gNhN2rMgS23WKEm0oNkWkxrgL0xz/v7Z65I8Iu5owuur9paY8KqY6uBvTUHnjfVhq/jnu+tgfaWvn+3t4efWQM1m7scPhRg428Sv9fCcx9V2jmwEm2jY8rkFUNWdt/rLDJE9BpM7n5FzPMFA1obkSiyRwS/uEeXRn+vO7TsTRxmyQRbf1pr3g57dwdbVHnFCQKst4ArgeycvtdXZJBFusZkZjcCP080YauZTQCudvebU1U5kQFjBiNHB1vRhOjvb2uNCa+uwbbltZeZMq4opoUWBtHemvD6Wh91tAx3b4r2vtyirl2Ko0phdFkY7uHj/n1lMDJfg0PkoIg6+OHrBEtRJBpGNTE8rmCS4S97BOSXBVsCG9oWMyX+ulqHtpaw9bUbGnfFhFZ3W1imqR/LrXQMEqnZkvx7snO7D60E+7Nb9wYtUYWZ9FPUYDK6v8ZUAfShb0Ikw2TnQH55sEXR3nYg0GK33sKtqaZvw/vbmqFue7Al4QyA55IMs9FlB46pZSZxkhmVdzlwefjSgZ+YWXxfRB7BfZSJpioSkVTIyu7bNbX29uA+tE5hFrbEGndBY3Ww7e14Hj62NkWvY8QwA+JaZmVhaI+F0eVBi3R0+Dq/PDieV6JRjcNcMi2mRqA6fG5ALV2XMt8H/C/w49RVTURSIivrwHWlKPY1JgisREFWDY27aavfQXb7vuj1ixpmWSPCEAuDa3+IhcG1P8TCfQqyISeZUXm/J1gMEDO7B/iGu28Y6IqJyEHWMTikZHLvZYElixcz/7R5cYHVc5jRWBW9ZdbeGkyFlcR0WABYdkxLLCaw8sce2B/bKlOQHXRR58q7ovdSIpKxOsKsuCL59+xrDAKroSoIqobq8HFnuK+68/Oooxq9LdqN1PuDbCwUjIX8ccFMKfljw8dx+/dbe2u0ukhSIk9JZGaXAFcDhxNcW+rE3fsx142IZJyoYdbaHBNiVZ2fJ9rXnyDrJcvOAnixtNvg6rx/LIzIjVaXDBX1PqaPAHcD9wLvDJ9nAecDNcCvUlw/EZHORuRC8aRgS0Zrc9dWV0PYIutoocU+b444LH/vrmDb+WrvZfOKe2yBddqfMypaPYaRqC2mLwLfAG4FPgH82N2Xm1kh8DjBQAkRkfQxIjeY0zGJeR2BA0FWvyMIrPqw9VS/M3w8sN8bq7Eos7R13CBd/XrvZUcWBoFVMD4IqoJDoHB8+Lrj+SFBt+MwuyYWNZhmAs+4e5uZtQFFAO5eZ2bfAb4H3J7iOoqIDJ4IQfb0k09w1knHhF1/OxOG14FQ2xltPbR9dbCrDnb1MtbMssPgCkOrcHw3ITZ+yHQlRg2mWqDjzLYCRwGLw9cGJL4NXkRkGPKs7OCXfuH43gu3h3MkdhtccfuTnWTY25Ifbp9XAoWHdA2x+EDLLTqoNz1HDaaXgGOBR4GHgBvNrJXgPqYbgaWprZ6IyDCRlXVgGqtxR/Vc1j2YsaN+54Gh8XVvJX4eZTLgjrkde7seNiJcs63wkPBxQhheEyjZXQXMT/47+yBqMH0bmBo+vzF8/mMgG3gR+GTqqja4zOxQ4CtAsbtfdLDrIyIZzOzATdFjD++5bGtzGFI7wsB6C+reThBiO5LvSmzdGyznkmBJlyNzxwLXRj+nCKLex/Q88Hz4vAa4wMxygVx3jzQm08z+A7iKYJqjl4Er3D3yHChmdjfwfmCHux8Td+xc4PsEwflzd7+1u88Jbxr+uJk9ELUOIiIHzYhcKJkSbD1pbw8HdYRBVfd29yHW0v04tubc0q73CaVYn5dW7+DuzUCzmZ0NfMnd39vbe8xsEvA5YJa77zWz+4FLCYahd5QZB+x197qYfTPcfX3cx90L/JC4oepmlg38CPgXoBJ40cweIgipb8d9xpXurmXiRWT4ysoKR/mNJZjatBvuwXpj+4PrrQOhVfcWu/cY/VjeMylJBZOZlQDnApOBjcCf3L0lPPZh4MvACcBrEb97lJm1AKPpupTGWcCnzew8d28ys6uBDwHnxRZy96fNbFqCz58HrO+YPsnMFgEXuPu3CVpYIiISzwxyC4OtfEaXw5sWL2baAFeh18HvZjYbWAvcB3wHuB94zsymmtkzwCKCkXofBWYl86XuvpVgWPkWYDtQ6+6PxZX5PcHaT4vM7KPAlcDFSZ4XwCTgzZjXleG+hMyszMx+Cswxs//qpswHzOzO2tp+rIsjIiI9SuaurG8Be4BTCVo2RxHMLv4icAxwubvPdvffuSe36IuZjQEuAKYTLDCYb2b/Fl/O3f8v0AT8BDjf3aOsZ51orGO3d8K5e7W7f8rdDwtbVYnK/NndP1FcPNANWRGRzJVMMM0FvubuS929yd3XAZ8GyoH/dPff9OF7zwE2uvvOsEvwD8Bp8YXM7AyC8HuQYHXcKCoJuh47VJB45V0REUkjyQTTeGBT3L6O1yv7+L1bgFPMbLSZGfAugu7C/cxsDnAXQcvqCqDUzG6J8B0vAjPNbLqZjSQYXPFQH+srIiKDJNkJlrrrAuvTnO/uvhR4AFhOMFQ8C7gzrtho4MPu/kbYRXg50GVQvZn9DngOOMLMKs3s4+F3tAKfIbgZeC1wv7uv6Ut9RURk8CQ7XPzRcIaHeE/E70922Qt3/zo9dM+5+zNxr1sIWlDx5S7r4TMeBh5Opj4iIpIekgmmhQNeCxERkZC5R5iyXQAws50E60/Fjhsv7uF17PNyoCoF1Yj/vv6U7e54T+fU22uds865r3TO/Subzuc81d3H9lrK3bX1YQPuTPZ13POXBuL7+1O2u+NRzlHnrHPWOeucU7UNr9WlBtefI7yOPzYQ39+fst0dj3KO8a91zqmhc+5fWZ1z9/sP9jl3S115g8zMXnL3uQe7HoNJ55wZdM6ZYTDOWS2mwRc/LD4T6Jwzg845Mwz4OavFJCIiaUUtJhERSSsKJhERSSsKJhERSSsKpoPAzA41s1/ELuNuZvlm9kszuytcf2pY6eacu+wbTro55w+GP+M/mdm7D2b9BkI353yUmf3UzB4ws08fzPoNhO7+Ow7/n15mZsNuYdJufs7zzWxJ+LOe35/PVzCliJndbWY7zGx13P5zzWydma03s+sB3H2Du3887iMuBB5w96uB8wep2v3S33Pu5t8hraXgnP8Y/owXAJcMWsX7IQXnvNbdP0Ww0OeQGFqdgv+fIVjZ+/7BqG8qpOCcHagH8giWHeozBVPq3Euw/Px+ZpYN/Ah4L8HqvpeZWXer/FZwYMXdtgGqY6rdS//OeSi6l9Sc81fD9wwF99LPczaz84F/AE8MXDVT6l76cc5mdg7wCvD2wFYzpe6lfz/nJe7+XoJA7tccqwqmFHH3pwlW9o01D1gf/nWxj2AZ+gu6+YhKgnCCIfJzScE5Dzn9PWcLfAf4X3dfPrC1TY1U/Jzd/SF3Pw0YEt3UKTjns4FTgI8AV5tZ2v8/3d9z9gMrmO8GcvtTl7T/xxriJnGgFQRB+EwyszIz+ykwx8z+Kzz2B+BfzewnDPL0HymW9Dl38+8wFEX5OX+WYAXni8zsU4Ncz1SK8nOeb2Y/MLOfMbSXoUn6nN39K+7+eeA+4K6YX9pDTZSf84Xhz/jXwA/786XJrsckfWMJ9rm7VwOfitvZQLBS71AX5Zy77BuiopzzD4AfDEqtBlaUc14MLB6EOg20pM855uC9A1qjgRfl5/wHgj+w+00tpoFVCUyOeV0BbDtIdRksOmed83Clcx6kc1YwDawXgZlmNt3MRgKXAg8d5DoNNJ2zznm40jkP0jkrmFLEzH4HPAccYWaVZvZxd28FPgM8CqwF7nf3NQeznqmkc9Y5o3PWOQ9EXTSJq4iIpBO1mEREJK0omEREJK0omEREJK0omEREJK0omEREJK0omEREJK0omEREJK0omEREJK0omESSYGY3mZmb2evdHF8fHr9pkKuWFDP7Xli/lxMcKzGzXeHxL6Tguy42s7fMzMLXl5pZs5nl9PezJTMomESS1wRMN7NOq7Ca2UnA1PB4upoN1AGHh4u/xfoSMDJ83iW4+uB9wMN+YFqZ44DV7t6Sgs+WDKBgEkleA/AkwUSWsS4N9zcMeo2SN5tg8s2RwKEdO81sHPA5DkzMuao/XxIuiHcu8NeY3ccD/+zP50pmUTCJRLMIuDimm8qAi8P9+5nZqWb2kJltM7MGM1thZl1WbzWzo83skbArrcHM1prZNb0diyIMn3HAXwhaTUfGHP4qQRhtBqrcfXvUz49zEjAGeDxm33HAWjP7lpltNbNaM7trKKzqKgeH/sMQieYPwHjgHeHrM4CxwINx5aYCzwBXAR8A/ge4x8wuiyv3ENAG/BtwPnAHUJjEsSiODR9XAa8ARwGY2RTgk8ANYZlUdeMtcfc94XeMBSYQrNw7ClgA3E7w79LtUuyS2bSCrUgE7l5jZo8QdN8tCR8fCffHltvfggpbVU8TLLJ2NfC7cH85QbfaB929IxSe6O1YHxwLNAOvAWsIgwm4CXja3Reb2a9Izeqj7wN+G/P6+PDxB+7+3fD542b2aWBmCr5PhiG1mESiWwRcZGa5wEXEdeMBmNkYM/uBmW0GWsLtE8DhMcV2AW8CPzWzS8Iut2SORTUbWBuurbMGONLMjgT+D/AVMysmWKW0v9eXJgBz6Hx96TighqC111HOgBKgqj/fJ8OXgkkkuoeAAuCbQD7w5wRl7gUuAW4D3k1w7eVuIK+jgLu3h8feCo+9ZWZLzGxOT8f6UN/ZHOim62gxfQP4q7u/EB6HbrryzGyWmd1tZk+Z2Q/NbGI333MesMHd18XsOw74e9yIvMMIuvWGzSJ7kloKJpGI3L2BYCDBfwB/Dl/vZ2Z5BF1aX3f3H7r7k+7+Egn+f3P3V939XwlaEOcQBNdfzSyrp2PJ1jUsO4sDobMaKAYuJBj4AEFXXzsJgsLM3gncBdxJcJ3rL8BfwhZXvPfRubUEQVdefEvsuPD7Vid7HpJZdI1JpG9+AuQCP01wLBfIJriuA4CZFRL8Yk+4ZHTYonjSzL4L3EcQRrt6O5aEmQStk5fDz9pqZvcDr7h7RzDMBta7e2OC938VuChmtN4jZrYr3P9vMec3kiA8L4rbdySwMu4zjyVoWaXz8Ho5iBRMIn3g7ouBxd0cqzWzF4EbzWwPQevgeqAWKOooZ2bHEoxQ+29gA8Ew6y8T/CKvMLNFiY65+67w/fOBp4Czw/ok0qWbzt0vSVCmuxF5Oe6+3cymAa+6e567vxC+jnUmQYvw7zH7jib4HRPfYjo2wT6R/RRMIgPjIwTdX78CqoEfAqOBz8SUeQt4G/gKMJFgkMBTBAHU1MOxDqPDxx091GM2sNvdt/ZQ5hjgsW6OmZmNArYSDjs3s6l0bbG9D/ibuzfH7DsOaATeiCt7LMG/i0hCdmDWEBEZSsxsIXCmu589gN9xFUEX3SfcfU848OF3wM3u/kRMudeA29z9roGqi2QOtZhEhq7TgO/2Wqof3P3nZlYHPBFOwloPfMXd/x5X7vCEHyDSB2oxiYhIWtFwcRERSSsKJhERSSsKJhERSSsKJhERSSsKJhERSSsKJhERSSsKJhERSSsKJhERSSv/P8Cm7uEaujoHAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Default mass function object.\n", "mf = hmf.MassFunction(hmf_model=\"Jenkins\", disable_mass_conversion=False)\n", "dndm0 = mf.dndm\n", "\n", "# Change the mass definition to 300rho_mean\n", "mf.update(mdef_model=\"FOF\", mdef_params={\"linking_length\": 0.3})\n", "plt.plot(mf.m, mf.dndm / dndm0, label=r\"FoF $b=0.3$\", lw=3)\n", "\n", "\n", "# Change the mass definition to 500rho_crit\n", "mf.update(\n", " mdef_model=\"SOVirial\",\n", " mdef_params={}, # NOTE: we need to pass an empty dict here\n", " # so that the \"linking_length\" argument is removed.\n", ")\n", "plt.plot(mf.m, mf.dndm / dndm0, label=r\"SO Virial\", lw=3)\n", "\n", "plt.xscale(\"log\")\n", "plt.yscale(\"log\")\n", "plt.xlabel(r\"Mass, $M_\\odot/h$\", fontsize=15)\n", "plt.ylabel(r\"Ratio of $dn/dm$ to FoF $b=0.2$\", fontsize=15)\n", "plt.grid(True)\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To see the difference between an assumed mass conversion and a properly re-fit mass function, we can compute an automatically-converted mass function for `Bocquet2016`, and compare to the refit:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/steven/Documents/Projects/halos/HMF/hmf/src/hmf/mass_function/hmf.py:291: UserWarning: Your input mass definition 'SOCritical(200)' does not match the mass definition in which the hmf fit Bocquet200mDMOnly was measured:'SOMean(200.0)'. The mass function will be converted to your input definition, but note that some properties do not survive the conversion, eg. the integral of the hmf over mass yielding the total mean density.\n", " f\"Your input mass definition '{mdef}' does not match the mass \"\n", "/home/steven/Documents/Projects/halos/HMF/hmf/src/hmf/mass_function/hmf.py:291: UserWarning: Your input mass definition 'SOCritical(500)' does not match the mass definition in which the hmf fit Bocquet200mDMOnly was measured:'SOMean(200.0)'. The mass function will be converted to your input definition, but note that some properties do not survive the conversion, eg. the integral of the hmf over mass yielding the total mean density.\n", " f\"Your input mass definition '{mdef}' does not match the mass \"\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEiCAYAAADTSFSPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzsnXl4Tcf7wD8j9iWWhBBCEPsWEqpF0VLaWlpb7dReqopW0Zb2W6W11VLaquVXe9UaeymKqiXUUkvFLrKIJQgh2/z+mJO4iSxXcm/uTTKf5znPPcvcOe+cc+95zzsz7/sKKSUajUaj0aSWbLYWQKPRaDQZG61INBqNRpMmtCLRaDQaTZrQikSj0Wg0aUIrEo1Go9GkCa1INBqNRpMmtCLRWBUhRDchxO+2liMWIUQeIcRGIcQ9IcRvtpYnsyOECBNClLO1HJZGCDFBCHFLCBEkhChttNPB1nLZCq1IMghCiK5CCF/jBxsohNgqhGhoa7lSQkq5TEr5mq3lMKED4AI4SSk7JjwohPhCCCGFEB8k2P+hsf+LdJITIUROQx4/IcRDIcQVIcRCIYR7esmQVqSU+aWUlyxdr3EtmqXh+9K4pmFCiBtCiOnmKgIhhBswEqgqpSwupbxmtDPaOL5HCNEvtbJlRLQiyQAIIUYAM4CJqIdgaWAu0NaWcqWEECK7rWVIhDLAeSllVDJlzgO9EuzraexPT1YDbYCuQEGgFnAUeDWd5UgSO73H5lJLSpkfaAy8A/Qx83tlgNtSyptWkyyjIaXUix0vqAdIGNAxmTK5UIomwFhmALmMY00Af2AUcBMIBN4C3kA9GO8AY03q+gL1APsVeAAcQ/3hYo+PBi4ax84Ab5sc6w38BXxn1DvB2LffOC6MYzeBe8BJoLpJOxcDIcBV4DMgm0m9+4GpwF3gMvB6MtejCrAHCAVOA22M/V8CEUCkcU37JvLdL4ClwFmgmrGvmrG9FPjC2FcY2GTIe9dYL5XgWlwyrtNloJux3wP402j/LeDXJNrQDAgH3JJppyvgY1zrC0D/BO1YZVzTB8Z18Da5h6sT1DUTmGVyLxYYv5Ubxn10SOYeJ9kmQAIelrzHwBIgxrg+YcAoY38bo52hxv2vksy1i5PL2F4FzEnwv3vmGpjclxjj3P8HuBv1ZQe+BqKBx8bx7239DEmX55StBdBLCjcIWgJRQPZkyvwPOAgUA4oCB4CvjGNNjO+PA3IA/Y0/8nKgAOoh+RgoZ5T/AvWg7WCU/8j4U+cwjndEPcCyod7iHgIljGO9jXMNNf5UeYivSFqg3qgLoZRKFZPvLgY2GDK5o5RcX5N6Iw3ZHYD3UApTJHItcqAeqmOBnMArqAdpJZP2LU3mWn6BUhhjgW+NfZOBMcRXJE5AeyCvIfNvwHrjWD7gvsk5S/BUKa0APjWuX26gYRJyfAP8mcJv40+UZZob8DTu66sm7XiMemFwACYBB41jZYBHgKOx7YB6YNY3ttcDPxntKAYcBgYmc4+TbBPxFYlF7rFR/grQzGS7Iuq32Nz4DYwyfgc5k/i+qVyVjfYPNzme3DVoAviblHU36stubO8B+tn62ZGuzylbC6CXFG4QdAOCUihzEXjDZLsFcMVYb4J6g4p9oyxg/OhfMCl/FHjLWP8i9oFjbGcz/mSNkjj3caCtsd4buJbgeG+eKpJXjIdHfYw3UWO/A/AE1eccu28gsMekjgsmx/IabSieiDyNgKAE9a/gqQL4AvMUSWngmvFQuga4YaJIEvmeJ3DXWM+HeituD+RJUG4xMA8T6yWJ+n4GViZz3A315lvAZN8k4P9M2rHT5FhVINxkez/Q01hvDlw01l2Me5HHpGwXYHcy9zjJNhn3ycOS99g4foX4iuRzYFWC3+0NoEkS35coZf/QWF/BUys+pWvQBK1I4i16jMT+uQ04p9AX7YrqKojlqrEvrg5pDASilApAsMnxcCC/yfb12BUpZQyqa8wVQAjRUwhxXAgRKoQIBaoDzol9NyFSyl3A98AcIFgIMU8I4Wh8P2cibShpsh1kUs8jY9VU5lhcgeuG3EnVlSJSymuoN9qJgJ+UMl67hBB5hRA/CSGuCiHuA3uBQkIIBynlQ5S1NggIFEJsFkJUNr46CmWNHRZCnBZCJNUvfxtlySSFK3BHSvkgmXYGmaw/AnKb/I6Wox6OoMZglhvrZVDKM9DkHv+EeiuPJeE9NqdNlrzHiRHvP2Dc/+skf9/rGPW/A7yAegEA866BxgStSOyfv1FdFG8lUyYA9eOPpbSxL7W4xa4IIbIBpYAAIUQZ1Jvy+6hZT4WAf1EPkVhkchVLKWdJKb1QXWoVgY9R/eqRibThRipkDwDcDLnTWtdi1OycxYkcGwlUQll2jsDLxn4BIKXcLqVsjlIG51DXDSllkJSyv5TSFfVGPlcI4ZFI/TuBekKIUknIFgAUEUIUMNn3PO38DWhi1P82TxXJddTbuLOUspCxOEopq5l8N949NrNNlrzHz8hAgv+AEEKgfsfJ1i8Vq1D/s3HGbnOuwfPIlunRisTOkVLeQ/3A5wgh3jLehHMIIV4XQkw2iq0APhNCFBVCOBvll6bhtF5CiHbG2+uHqD/VQdQbm0T1xSOEeBdlkZiFEKKuEOIFIUQOVJfCYyDasJZWAV8LIQoYCmtEKttwyKh7lHGdmgCtgZWpqOtX4DVDtoQUQFlyoUKIIsD42ANCCBchRBshRD7UtQtDdUMhhOhoohzuoq5nNAmQUu4EdgDrhBBeQojsxrUZJIToY1hIB4BJQojcQoiaQF9gmTkNk1KGoLpgFgGXpZRnjf2BwO/ANCGEoxAimxCivBCicVJ1mdMmC99jUBa1qX/KKuBNIcSrxu9rJOraHzCzvm+AAUKI4qm5BinIlunRiiQDIKWcjvrTfYZ6iF9HWQXrjSITAF/ULKhTqJlWE9Jwyg0oc/8u0ANoJ6WMlFKeAaah3t6CgRqoGTzm4oh6M7+L6oa4jZqlA2rw9iFqptN+1BvywucVXEoZgZq98zrqLXguaizgXCrqCpdS7pRShidyeAZqoPkWSsluMzmWDfUgC0DNbGoMDDaO1QUOCSHCUDOuhkkpLychQgdgC0qh3UNZf94oawVU15S7cZ51wHgp5Y7naOJy1Cyk5Qn290R1Q51B3avVJN/NZm6bLHKPDSahXp5ChRAfSSn/A7oDs1H3pDXQ2vg9pIiU8hRq8sLHxq7nvQamzAQ6CCHuCiFmmdugjIwwBoc0GkA55KFms3S3tSwajSZjoC0SjUaj0aQJrUg0Go1GkyZ015ZGo9Fo0oS2SDQajUaTJrQi0Wg0Gk2ayMiRO83G2dlZuru7p/r7Dx8+JF++fCkXzEToNmcNdJvNJCoKQkKgSBHIlcsqcl26dIkHDx5Qy80NLl+GypXBQvcmtff56NGjt6SURVMsaOsYLemxeHl5ybSwe/fuNH0/I6LbnDXIKm1evXq1vHjxopQylW3etElKkHLvXssKZhARESELFiwo+/TpI2WvXlI6OUkZFWWx+lN7nwFfqWNtaTSarM6DBw/o1q0bs2fPTn0lx46pT09PywiVgH379nHv3j1at2oFv/8OzZqBQ8ZJuKgViUajydRs3ryZJ0+e0L59+9RXcuwYVKwIBQqkXDYVbNy4kZw5c/KaqysEBkKLFlY5j7XQikSj0WRq1q5di4uLCy+++GLqKzl2DLy8LCeUCVJKNm7cyCuvvELe/fvVzubNrXIua5ElBtsTIzIyEn9/fx4/fpxi2YIFC3L27Nl0kMp65M6dm1KlSpEjRw5bi6LRpBvh4eFs2bKFHj164JDarqL79+HWLahTx7LCGZw7d46LFy8yYsQIWL8eqlaFUkkFfbZPsqwi8ff3p0CBAri7u6MiTifNgwcPKGAlkzY9kFJy+/Zt/P39KVu2rK3F0WjSjSNHjvD48eO0dWs5OkJoKESYFf/xudm0aRMArZs0geHD4f33rXIea5JlFcnjx4/NUiKZASEETk5OhISE2FoUjSZdefnllwkICMDJySltFeXIoRYrsHHjRmrVqoXbf/8pZdW6tVXOY02y9BhJVlAisWSltmo0phQrViz13VqgLIRp0ywnkAnBwcH89ddftGnTBjZuhIIFoUEDq5zLmmRpRWJrrl+/TtOmTalSpQrVqlVj5syZANy5c4fmzZtToUIFmjdvzt27dwHVRfXBBx/g4eFBzZo1ORY7JVGj0TzD/PnzadKkCaGhoamvJDoafvkFLl2ynGAmrF+/npiYGDq0awebN8Prr1vN8rEmWpHYkOzZszNt2jTOnj3LwYMHmTNnDmfOnOGbb77h1Vdfxc/Pj1dffZVvvvkGgK1bt+Ln54efnx/z5s3jvffes3ELNBr7ZfHixQQHB1OwYMHUV3LmDISFQf36lhPMhDVr1uDh4UGN8HC4eTNDdmuBViQ2pUSJEtQxZoIUKFCAKlWqcOPGDTZs2ECvXr0A6NWrF+vXq0SIGzZsoGfPngghqF+/PqGhoQQGBgIwefJkatSoQa1atRg9erRtGqTR2AlXr15l3759dOvWLW3dugcPqk8rKJLbt2+za9cuOnTogNi0STkgtmxp8fOkB1l2sN2UDz/8kOPHjyd5PDo6+rn7WD09PZkxY4bZ5a9cucI///zDCy+8QHBwMCVKqKyeJUqU4ObNmwDcuHEDNze3uO+UKlWKGzducPz4cdavX8+hQ4fImzcvd+7ceS5ZNZrMxooVKwDo2rVr2io6eBCcnMDDwwJSxWfDhg1ER0fToUMH6NNHjY0UKWLx86QH2iKxA8LCwmjfvj0zZszA0dExyXIykdwxQgh27tzJu+++S968eQEokkF/jBqNJZBSsnDhQho2bEi5cuXSVlmhQtC2LVhhssqaNWtwd3enjrMznDyZYbu1QFskAClaDtb0I4mMjKR9+/Z069aNdu3aAeDi4kJgYCAlSpQgMDCQYsWKAcoCuX79etx3/f39cXV1RUqpZ2VpNAZRUVEMHTrUMj5TVpqtFRoayo4dO/jggw8QGzaonRlYkWiLxIZIKenbty9VqlRRXq0Gbdq04ZdffgHgl19+oW3btnH7Fy9ejJSSgwcPUrBgQUqUKMFrr73GwoULefToEYDu2tJkaXLkyMHQoUNp1apV2iqKjraMQIng4+MT9xLJ6tVQvTpUqmS181kbrUhsyF9//cWSJUvYtWsXnp6eeHp6smXLFkaPHs2OHTuoUKECO3bsiBs8f+ONNyhXrhweHh7079+fuXPnAtCyZUvatGmDt7c3np6eTJ061ZbN0mhsRmhoKAsWLCAsLCztlU2YAOXLw5Mnaa8rAcuWLcPd3Z367u6wfz+kxfPeDtBdWzakYcOGiY57APzxxx/P7BNCMGfOnETLjx49Ws/W0mR5li1bxvvvv0/t2rXjZkSmmn37VHgUCyeyCgoKYufOnYwePRqxfj1ICR06WPQc6Y3dWSRCiJZCiP+EEBeEEEk+GYUQHYQQUgjhnZ7yaTQa+0RKyfz58y2jRCIi4MABePllywhnwqpVq4iJiaFbt26wZo0KT1+tmsXPk57YlSIRQjgAc4DXgapAFyFE1UTKFQA+AA6lr4QajcZeOXToEMePH6d///5pr+zoUQgPh8aN015XApYtW0atWrWoWqwY7NmjrJEMPlnGrhQJUA+4IKW8JKWMAFYCbRMp9xUwGUg5BrxGo8kSzJo1i4IFC9KjR4+0V/bnn+qzUaO012WCn58fhw8fVtbIhg1qQD+Dj4+A/Y2RlASum2z7Ay+YFhBC1AbcpJSbhBAfJVWREGIAMADUdNo9e/bEO16wYEEePHhgllDR0dFml7VnHj9+/Mx1SIqwsDCzy2YWdJszLtHR0Zw6dYrXXnsNX1/fZMua0+ZCuXJRuHt3Lp8+bUEp1SxMIQTu7u7cmTKFPCVKcOjePWWZWBGr32dzErun1wJ0BOabbPcAZptsZwP2AO7G9h7AO6V6vby8nklqf+bMmRQT38dy//59s8vaM8/T5t27d1tPEDtFtznj8+TJkxTL2KrNMTEx0sPDQzZp0kTK4GApHRykHD06Xc6d2jYDvtKMZ7e9dW35A24m26WAAJPtAkB1YI8Q4gpQH/DRA+4aTdblyZMncb5TOXPmTHuFQUFw6hTExKS9LhP27dvHhQsXePfdd+HXX1W3VrduFj2HrbA3RXIEqCCEKCuEyAl0BnxiD0op70kpnaWU7lJKd+Ag0EZKmbwta8e4u7tTo0YNPD098fZW+lCHkddozGfFihWULFmSc+fOWabC5cuhZk2lUCzIggULKFCggHJCjD1H9eoWPYetsCtFIqWMAt4HtgNngVVSytNCiP8JIdrYVjrrsXv3bo4fPx7Xt6vDyGs05hETE8PkyZOpWLEilSzlGf7nnypIo6urZeoD7t+/z2+//UaXLl3IFxSkgkGmNaCkHWFXigRASrlFSllRSlleSvm1sW+clNInkbJNMrI1khQ6jLxGYx4+Pj6cPXtWOfdZYgptdLRyRLSw/8ivv/5KeHg4ffv2VdYIQJcuFj2HLbG3WVs2o0mTJs/s69SpE4MHD+bRo0e0TiSgWu/evenduze3bt1SoaBNMHeGhBCC1157DSEEAwcOZMCAATqMvEZjBlJKJk2aRLly5ejYsaNlKj12DO7ehWbNLFOfwYIFC6hWrRp1vb2hZ0+lqEqXtug5bIlWJDbmr7/+wtXVlZs3b9K8eXMqV66cZFmpw8hrNHH8+++/HD58mB9++IHs2S30KPv9d/VpQUVy+vRpDh06xPTp0xHHjsF//8Hw4Rar3x7QisQgOQsib968yR53dnZO9RxtV6MftlixYrz99tscPnxYh5HXaMygRo0anD59Ou05R0wZOlRlQyxa1GJVzp07l1y5cilHyXHjIHdueOcdi9VvD9jdGElW4uHDh3GOjg8fPuT333+nevXqOoy8RpMCERERAFStWpXcuXNbrmJHR3j1VYtVd//+fRYvXsw777yDc758anykfXuVMCsToS0SGxIcHMzbb78NqGQ8Xbt2pWXLltStW5dOnTqxYMECSpcuzW+//QaoMPJbtmzBw8ODvHnzsmjRIkCFkT9+/Dje3t7kzJmTN954g4kTJ9qsXRqNNZFS0rJlS6pXr86sWbMsV/GBA7B7t7JKkslU+jwsXbqUsLAwhgwZAmvXwr170LevReq2J7QisSHlypXjxIkTz+x3cnLSYeQ1miTYtWsXu3fv5q233rJsxStXwvz5MHKkRaqTUjJnzhy8vb2pV68ejBkDZctaJRCkrdFdWxqNJsMgpeTTTz/Fzc2NgQMHWrbyHTvUQ95CXWV//vknZ86cYfDgwXD5MuzaBX36QLbM99jNfC3SaDSZlo0bN3Lo0CE+//xzclky4dS1a3DuHLz2msWqnD17NoULF6Zz586waJEKFW/4h2U2tCLRaDQZAiklX331FZUqVaJ3796WrTx22m/z5hap7uLFi6xbt45BgwaRx8EBfv4ZXn8d3NxS/nIGJEuPkWSlabOJ+aBoNBkJIQSbNm0iICCAHDlyWLbygACLZiqcPn06OXLkYOjQobB+vYrbNWSIReq2R7KsRZI7d25u376dJR6wUkpu375t2WmSGk068uTJE6SUuLi4ULt2bcufYNw4OHPGIpkKb926xaJFi+jevbuKUDFnjhpkb9HCAoLaJ1nWIilVqhT+/v6EhISkWPbx48cZ/iGcO3duSpUqZWsxNJpUMXr0aHx9fdm1a5flrZGYGDUA7uBgkep++OEHwsPDGTlypApHv3cvTJ5ssfrtkSyrSHLkyEHZsmXNKrtnzx7rvAVpNJoUOXnyJLNnz6Zv376WVyKgupwuXoTt29NskYSHhzN79mzefPNNqlatCoMHQ65c8O67FhLWPsmyXVsajcb+iYmJYdCgQRQuXJhJkyZZ4wTg4wMFC1qkW2vJkiWEhITw0UcfKefDJUugc2dwdraAsPZLlrVINBqN/bNgwQL+/vtvfvnlF+sEIz12TA20JxLd+3mJiopi8uTJeHt707hxY5g6FcLClKd8JkcrEo1GY5dIKZk3bx6NGzdWAQ+twcaNanzkjTfSXNXSpUu5ePEiPj4+iKgomDULmjQBL6+0y2nnaEWi0WjsEiEEe/fu5e7du9abpr9hA7z4Ypq7nqKiopgwYQJ16tShVatWKjijvz/88IOFBLVvtCLRaDR2x+HDh6levTp58+YlT5481jlJTAwMHAhGmoa0EM8aAZg2DSpXtoilkxHQg+0ajcauuHHjBi1atLB8LK2EZMsG772nwrqngWeskT174J9/YMSITBlXKzG0RaLRaOwGKSX9+vUjIiKC8ePHW/dkv/6qco+ksVsrnjUiBEycqKyc7t0tJKj9kzXUpUajyRDMmTOHbdu2MXnyZDw8PKx3ov/+U9Nyly1LUzVPnjzhyy+/fGqNHDgAO3fCxx+Dtbrk7BBtkWg0GrvA19eXESNG8Oabb/Lee+9Z92RGsjg6dEhTNXPmzOHKlSv8/PPPyhr56itl4QwaZAEhMw7aItFoNHaBk5MTrVq14pdffiGbtccWVq2Chg2hZMlUV3H37l0mTJhAixYtaNasGRw+DNu2qcRY+fNbUFj7R1skGo3GpsQGTi1btixr1661+vnyXrumYmClMU3vxIkTCQ0NZfLkyWrHhAlQpEimjvKbFNoi0Wg0NmXy5Mm0a9eO8PDwdDlfoX/+UQEU09CtdeXKFWbNmkXv3r2pWbOmmqW1cSMMHw4FClhQ2oyBViQajcZmbNiwgTFjxpArV650i7Ad0LatyohYokSq6/j0009xcHDgf//7n9rx+ecqXlcWCIeSGFqRaDQam3DixAm6deuGt7c3ixYtSp8kc7H5h1xdU13F3r17Wb58OR999JFKzfDnn7B5M4wZo5RJFkSPkWg0mnQnKCiI1q1bU6hQITZs2GA97/WEDBhApRs3VAysVBAZGcmQIUNwd3dn9OjRSjF98okatP/gA8vKmoFIsyIRQtQFhJTysAXk0Wg0WYCgoCCyZcvG+vXrVRbB9ODBA1i+HPnKK6mu4vvvv+fff/9lw4YN5M2bF9auhUOHYP78LOU3khBLWCQHUF1kmTf9l0ajsQjR0dE4ODjg6enJ+fPnyZkzZ/qdfM0aePSIoBYtSE3HVkBAAOPHj+fNN9+kdevWEBWlurOqVoVevSwubkbCEorkDfRYi0ajSYHIyEjat29P7dq1+fLLL9NXiQAsWgQVKnC/WrVUff2jjz4iIiKCmTNnqvGcBQvg/HkVQTh71h4lSFIBCCG6CiFSzCQjpdwhpdxuWbE0Gk1mIjIyku7du7Nx40aKFy+e/gKcPq1yp/fvn6pMiBs3bmTFihWMGTOG8uXLw9278Nln0KiRRZJiZXSSsySWAB4AQogIYyxEo9FonovIyEi6dOnCqlWrmDJlivXDnySGs7N68Kcid3poaCiDBg2iRo0ajBkzRu0cNw7u3IHZsy2Sojejk5w9FgoUN6OcRqPRJIqUkq5du7JmzRqmT5/O8OHDbSOIi4uKg5UKRowYQXBwMD4+Pqo77sQJmDtXhaCvVcvCgmZMklMQfwBLhBBnje0FQoiwpApLKV+yhEBCiJbATNTg/Xwp5TcJjo8A+gFRQAjQR0p51RLn1mg0lkUIQZs2bXj55ZcZaitnvW3bICJCdUE9p/Wwbds2Fi1axJgxY/Dy8lLTfYcOVaFQYp0RNckqkj7A+0AloB5wA7hlTWGEEA7AHKA54A8cEUL4SCnPmBT7B/CWUj4SQrwHTAbesaZcGo3m+QgJCeH48eM0b97cevnWzUFKNbMqJua5xzLu3r1L//79qVKlCuPGjVM7ly+Hfftg3jylTDRAMopEShkGfAMghGgGjJZSnrCyPPWAC1LKS8Z5VwJtgThFIqXcbVL+IJB1ssdoNBmAS5cu0aJFC27fvs3ly5cpaEtv7wMH4PhxlTv9OawRKSX9+/cnKCiItWvXqvAtt2+rWFp160KfPlYUOuOR3KytCCGEt7G5DbiXDvKUBK6bbPsb+5KiL7DVqhJpNBqz2b17N3Xr1uXOnTts2rTJtkoEYOpUcHKCnj2f62sLFixgzZo1fP3119Sta8wzGj5czdaaP18FfdTEkVzXVjSQy1jvA/wMXLGyPIm9MshECwrRHfAGGidxfAAwAMDFxYU9e/akWqiwsLA0fT8jotucNbBkm9evX8/s2bNxc3Pj66+/JiIiwqbXM8/169TbsIGr3btz5fDTwBsptfnatWsMHTqUOnXq4O3tzZ49eyhy6BA1lyzhSo8eXLlzR+Vlz0BY/bctpUx0AY4B24FBQAwwCfVgTmzpn1Q9z7MALwLbTbbHAGMSKdcMOAsUM6deLy8vmRZ2796dpu9nRHSbswaWbPOIESNk69at5b179yxWZ5r44w8pK1WSMigo3u7k2vzo0SPp6ekpnZyc5I0bN9TO+/eldHOTskoVKR8/tqLA1iO19xnwlWY8Y5OzSD5AWSHfo6yCT5LTR0bZtHIEqCCEKIsa3O8MdDUtIISoDfwEtJRS3rTAOTUaTSq5cOECoaGheHt7M3nyZIQQ1s9uaC6vvAJnz5o9NiKlZMiQIRw/fpyNGzfiGhsh+JNPwN8f/voLcuVKvpIsSpJ3XEq5X0pZBciJ6nJqAORIYrFIrAMpZRRqpth2lMWxSkp5WgjxPyFEG6PYFCA/8JsQ4rgQwscS59ZoNM/HypUrqVOnDv3790dKiYODg/0oEV9fePLkuQbY582bx6JFi/j8889p1aqV2rl5sxqoHz4cXnzRSsJmfFJ0NJRSxgghmgOnpJTR1hZISrkF2JJg3ziT9WbWlkGj0STNgwcPGD58OAsWLKBBgwYsX748fXKJmEtoKDRvDu3aqXhYZnDo0CGGDh1Ky5YtGT9+vNoZHKxmZ9WsCRMnWlHgjI9ZHutSyj+sLYhGo7F/rl+/TqNGjbh+/Tpjx47lyy+/JLu9BSycNUspk/ffN6t4QEAA7du3p2TJkixbtgwHBwflf9KnD9y/D7t26S6tFEjyFyCECADekFIeF0IEksTsqViklKlPOabRaOwaKSVCCEqWLEnz5s3p06cPL9pjV09oKHz3Hbz1FtSunWLxhw8f0rp1a0IlIglUAAAgAElEQVRDQ9m/fz9FYp0M586FLVtULK1URgvOSiT3KrEAuGmynqwi0Wg0mZOtW7fy+eef4+Pjg6urKz//bIl5NVYi1hoZNy7FotHR0XTr1o3jx4+zYcMGPD091QFfXxgxAt54A4YMsbLAmYPkPNs/N1n/LH3E0Wg09sL169f58MMPWbt2LZUqVeLmzZtPZzLZK3//bbY1MmrUKDZs2MDMmTOfDq7fvg0dOkDx4rB4sY7sayZmd24KIVwBF2MzWEoZYB2RNBqNLZFSMmXKFL788kuklEycOJERI0aQKyOME2zZAvdSDsIxdepUpk+fzvvvv88HsbnWo6Ohe3cIDIT9+5VHvMYsklUkRhDFUcB7JAhVIoTwB+YCU6SUMVaTUKPRpCtCCP7991+aNWvGzJkzcXd3t7VIKRMYqMKWFCsGhQolW3Tz5s1MnTqVTp06MWPGjKcHvvpKRQr+8UcVT0tjNsnF2hLABuAr4BBKmbQG2hjrvsDXwHrri6nRaKzJsWPHaN68Of/88w8A8+fPZ8OGDRlDiYAa06hZEx4/TrZYbF6Uli1bsmTJEjVDC2DVKvjyS5V7fcCAdBA4c5GcRdIFeA14UyaeSneekTtkgxCis5RypVUk1Gg0ViMoKIgePXqwdOlSnJycuHbtGrVr107/fOpp4cgRWLlSZUDMnTvJYhs3bqRr165UrVqV1atXP23j4cNKgTRoAD/9pMdFUkFybqhdgV+SUCIASCm3AYvRodw1mgzH+PHj6dmzJ6tXr2b06NFcvHiRtm3b2lqs5yMmRvmLuLjAqFFJFlu3bh3t27enVq1aTJo0iXz58qkD165BmzZQogSsW6f9RVJJcoqkNip8fEpsNcpqNBo759atW7GBTxFC8Morr/Dff/8xadIk24d8Tw0LFiiLYupUKFAg0SKrV6+mU6dOeHt7s2PHDvLnz68OhIZCq1YQHg4bN0LRoukoeOYiOUXiBASaUUeQUVaj0dgpt27dYvTo0bi7u+Pjo8LTjR8/ntGjR1O6dGkbS5cGTp6Exo2hW7dED69YsYLOnTtTv359tm/f/lRZPnqkMiaeOwerV2unwzSS3BhJTlRe9JSIRgVu1Gg0dkZgYCAzZsxg7ty5PHz4kM6dO1O1alUA+4qPlVpmz1YD7Im0ZcaMGQwfPpzGjRuzadOmOEtEREXBO++oaL6//qricmnSREp+JIOFEG+lUMbOPZQ0mqyJlJKmTZvi5+dHx44dGTduXJwSyfD8/Tfky6dmaiUYYI+JiWHUqFFMmzaNdu3asWzZMpUqFyA6mkqTJ8OOHSqqb8eONhA+85GcIgkAXjWzHu2cqNHYAUePHuWHH37g+++/J3fu3Pzwww+4ubnh4eFha9Esx4MH0KUL5M+vurZMQtc/efKEd999lxUrVvD+++8zY8aMp1N8o6OhTx+K79gBEybAoEE2akDmI7kQKaXSUxCNRpM6pJTs2rWLb7/9lh07duDo6MiAAQOoV68eTZs2tbV4lueTT9Rsq3374imRoKAg2rdvz4EDB/j222/5+OOPn3bfRUfDu+/CkiVc7t2bsp9+aiPhMyd2Fv9Zo9E8D7dv345zJCxevDjffvstAwcOzJgzsMwhtktqxAjl92Hg6+vLW2+9xd27d/n111/p1KnT0+9ER0Pv3rB0KXz1FVcbNqRs+kueqbGTdGYajcZcgoKC2Lx5MwBFihShcuXKzJs3j8uXLzNq1KjMq0RCQqBnT6hSRYUzMViyZAkNGzYke/bsHDhwIL4SefJEdYMtXaq6sz7T8WetgbZINJoMwrFjx5g5cyYrV64kZ86cBAUFkS9fPpYvX25r0dKHQoWgb1/o1Any5uXRo0d8+OGH/PzzzzRt2pRVq1bh7Oz8tPy9eyoS8J49MG2asmI0VkFbJBqNnXP06FEaN26Ml5cXa9asoX///vj6+j71zs4KPHkCOXIoq6JmTU6fPk29evX4+eefGT16NNu3b4+vRAIDlX/J/v3KGtFKxKpoRaLR2CGhoaH4+/sDkCtXLq5du8bUqVPx9/fn+++/p1KlSjaWMB3ZvRsqVoR//0VKybx586hbty4hISFs376dSZMmkSOHiSvb6dPw0ktw4QJs3pyks6LGcuiuLY3Gjjh//jyzZ89m0aJFtGnThuXLl1O9enUuXrxItmxZ8L3vyhXl61GsGNeFoF/Llvz+++80a9aMJUuWULx48fjlfXyU4sifX3VpeXvbQuosR3I5258nlrKUUtpx/k2Nxr75888/mTJlCps3byZnzpx07tyZYcOGxR3Pkkrk4UNo2xYZHc1v3brR78UXiYmJYc6cOQwaNCj+NZESJk6Ezz8HLy9Yvx5Klky6bo1FSc4i+fE56pGAViQazXMQHh5Orly5yJYtG1u3buXIkSN88cUXDBo0CBcXl5QryMwYU3blv//yqacnkz77jCZNmrBw4ULKlk0weff+fejXD377TVkjP/8MefLYRu4sSnKvOTmeY8lAyQs0Gtvi7+/P2LFjKVWqFNu2qQDbY8eO5dq1a4wfP14rEeBRaCh+//zDJ0LwvZ8fc+bM4Y8//nhWiRw7BnXqwNq1MHkyLFmilYgNSFKRSCmjn2dJT6E1mozIwYMH6dKlC+7u7nz77bc0btyYEiVKAODo6JgxcqJbGSklG9eupZq3N5UvXiSgc2f+++8/Bg8e/GxX1uzZ8OKLakbXnj3w8cc6KZWNeK7BdiFECaAC8EwaMinl75YSSqPJbERFRdGpUyfu37/PsGHDGDp0aMZJY5tOHD16lD3duvHyf/9RolIl/m/PHho3bvxswYAAGDgQNm2CN9+EX34BJ53JwpaYpUiEEPmBFcAbprtRYyOxOFhQLo0mQ3Pr1i1++ukn1q1bx4EDB8iZMyfr16+nYsWKTxMraQC4dOkSn376KblXrmQBcK1KFf709SVH3rzxC0oJy5fD0KEqGdV338EHH8SLt6WxDebegUmAB9AUpUA6Ac2AX4DLQIOkv6rRZB1OnTpFv379cHNz47PPPsPJyYmQkBAA6tSpo5WICVevXmXw4MFUrlyZwmvWsACIadoU96NHn1UiwcHQvj107w6VK8OJE/Dhh1qJ2Anm3oU3gQnAX8b2VSnlLillH2AT8KE1hNNoMhK+vr7UrFmT5cuX06tXL06fPs327dspqaehxuPixYv069cPDw8P5s+fz4+NGjE3MpJsLVuSffPm+IPl0dEqSGPlyrBlC0yZoqL+VqxouwZonsHcMRIX4JqUMloI8ZD4qXU3AastLplGY+c8ePCARYsW8eTJEz7++GO8vLz48ccf6dixI0WKFLG1eHbHv//+y+TJk1m+fDnZs2dn0KBBjBo1CrdYH5AZM+InqTpyBAYPBl9faNIE5s5VARs1doe5Fsl1niqPC8QfK/EGHltSKI3Gnrl06RLDhw+nVKlSDBs2jF27diGlRAjBwIEDtRIxISYmhk2bNtGsWTNq1KjBmjVrGDZsGJf9/JhdowZurq5QujT8+ONTJRIcrJJOvfAC3LihxkV27dJKxI4x1yLZiRoTWQ/MABYJIWoDT1DjJjOtI55GY1/MmjWLDz/8EAcHBzp16sSwYcOoV6+ercWyO+7du8fixYuZNWsWFy5coFSpUnzzzTf079+fIg4O0LWr6qoqWVLNvAIIC4Pp01X3VXg4DBsGX34Jjo62bYwmRcxVJKOBfABSyl+EEI+ADkAeYDgw1zriaTS2JSIigpUrV1KnTh2qV6/Oyy+/zNixYxk8eDCurq62Fs+ukFKyf/9+5s+fz2+//UZ4eDj169dnwoQJtGvXTgVWPHsW3n4bLl5UVsibb0JkJCxcCOPHPx1UnzhRj4NkIMxSJFLKMCDMZPs34DdrCaXR2JrQ0FB++uknZs2aRUBAAKNGjeLbb7/F09MTT09PW4tnVwQFBbFkyRLmz5/P+fPnKVCgAD179qRfv354mwZN3LwZOneGvHnhjz9U19WPP8I338DVq9CwIaxbp5wMNRmK547+K1QS5BwJ90spIywikUZjY8aNG8d3331HWFgYzZo1Y+HChbz22mu2FsuuuHv3LmvXrmXFihXs3r2bmJgYGjZsyNixY+nQoUPiuVKcnKB2bViwALZtU91bN24ohTJnDrzxhvZMz6A8j0PiBKAdUILEB+kt4pAohGiJGnNxAOZLKb9JcDwXsBjwAm4D70gpr1ji3Jqsy/Hjx5FS+ddGRETw9ttvM3LkSGrVqmVjyeyH+/fvs3nzZlasWMG2bduIjIykfPnyjB07lm7dulG5cuVnv3T8uMqz/vHHUK4cNGsGjRqpLqxGjWDRIrVPK5AMjbkWyY/AW8Ai4AxgFetDCOEAzAGaA/7AESGEj5TyjEmxvsBdKaWHEKIz8C3wjjXk0WRuYmcUTZ06lX379jFlyhSaNm3KpEmTEPrBBsD169fx8fHBx8eH3bt3ExkZScmSJRk6dChdunTBy8sr8WsVGakGzseNg4IFlQPh6tUqLlaLFjBmjMpgqMkUmKtIXgeGp0POkXrABSnlJQAhxEqgLUp5xdIW+MJYXw18L4QQMvZ1UqNJgYiICBYtWsT06dM5f/48pUuX5rvvvqOiMbiblZVIdHQ0x44dY+vWrWzYsIFjx44BUKFCBYYNG0bbtm156aWXks+PcuSICut+8iQ4O0NIiBr76NNHhTfR03gzHcKc568Qwh/oJ6XcZlVhhOgAtJRS9jO2ewAvSCnfNynzr1HG39i+aJS5laCuAcAAABcXF6+VK1emWq6wsLAsF9oiM7Y5OjoaBwcHIiMj6dKlC05OTrzzzjs0btwYBweHTNnmlAgLCyMsLAxfX198fX05duwYDx48QAhB1apVadCgAQ0aNKB06dJm1ed46hSeI0YgoqMRUhLu4kJA27YEvvkmUXYyjTer3ufUtLlp06ZHpZQpp5mUUqa4ACOAdRiKx1oL0BE1LhK73QOYnaDMaaCUyfZFwCm5er28vGRa2L17d5q+nxHJTG0+d+6cHDhwoKxYsaKMiIiQUkrp7+8vY2Ji4pXLTG1Ojhs3bsiVK1fKwYMHSzc3N4kKvipLliwpe/fuLZcvXy6Dg4PNr/DaNSnffVfK2rWlBClz5JCyfXspd+6UMjraeg1JJVnlPpuS2jYDvtKMZ7e5XVtFgTrAWSHELiD0WX0kPzWzruTwB9xMtksBAUmU8RdCZAcKAncscG5NJkIaPg1Tp07Fx8eHXLly0bNnT8LCwihcuHCWiX8lpcTPz499+/bFLZcuXQIgX758VK9enZEjR9K8eXOqVKlifrferVtqzOOHH1QXFkD58ipHSNeuoL37sxTmKpLuxmc+oHUixyVgCUVyBKgghCgL3AA6A10TlPEBegF/o5widxmaU6OJY+/evTRp0gQnJyfGjRvHkCFDKFasmK3Fsjp37tzhyJEjccvBgwe5efMmAM7OzjRq1Ij333+fRo0a4enpyf79+2nSpIl5lfv7w8aNKh/6zp0QE6P2Fy4MX3yhxj+y8PhSVsZch0S3lEulHSlllBDifWA7avrvQinlaSHE/1Amlg+wAFgihLiAskQ6p4dsGvsmLCyMhQsXEh0dzfDhw2nUqBGLFy+mffv25E0YkjyT8PDhQ/755x+OHDnC4cOHOXLkCBcvXow7XqlSJVq0aEHDhg1p1KgRlStXfr6JBFKqVLY+PkqB/POP2l++POTPr2JjjRunBtZ1dscszXM7JFobKeUWYEuCfeNM1h+jxlI0GgICApg9ezY//vgjoaGhtGrViuHDh5MtWzZ69Ohha/EsgpSSq1evcvLkSU6cOMHJkyc5efIkfn5+cb4vbm5u1K1bl379+lGvXj28vLwoWLDg858sMFB5ne/cqfw/AgKUlVGxolpWrYKaNVWok/LltQLRAM+hSIQQ7sBIoCFQBGUN7AOmS+0QqLEBP/30E0OHDiU6Opp27doxcuRI6tevb2uxUo2UkpCQEM6dO8eZM2c4depUnNK4f/9+XDkPDw9q1qwZ58dRt25dihcvnqpzOjx8qFLW7typltOn1QEnJ+UwmCuXyv/x339QqZKyQoSAqlUt0WRNJsFcz/bawB4gCtgMBKNylHQFegohGkspT1hLSI0G1IN2586duLu7U6FCBby9vRk0aBAffvgh5cqVs7V4ZhMdHc2VK1c4e/Ys586di/d59+7duHKOjo7UrFmTHj16ULNmTWrWrEn16tXTNnX12jX466+4peHJk2qsI3duePll6NVLeZo7OkKtWvDwITRtCvPmweuv64yEmkQx1yKZCpxE+W88jN0phMgHbAOmocLMazQWJzYC77Rp0zh58iRDhw5l1qxZeHl54eXlZWvxEiUqKorr169z8eLFuOXChQtcuHCB8+fP8+TJk7iyLi4uVKlShXfeeYfKlStTpUoVKleujJubW9qcI8PDlUe5r69SHPv3qwFzgHz5oH59rnbvjnuvXspJcPVqpThq11bjI8OGQYcOalujSQZzFUl9VEyrh6Y7pZQPhRCTgdR7+2k0yTB79my++eYbAgICqFq1KgsXLqRr14QT+WxDeHg4ly5diqcsYhXGlStXiIqKiiubK1cuypUrh4eHBy1btoxTFpUrV6Zw4cKWEEZNw/X1haNH1XL6tEpVC1CqFDRo8HSpWROk5MHUqTB/vvI8f/xYWR+ffKK6r77+Ou1yabIE5iqSx0ChJI4VQmdI1FiQGzdu4OrqihCC8+fPU6VKFRYsWECLFi3SLXyJlJLQ0FCuXr3K1atXuXbt2jPrwcHB8b5TsGBBypcvT+3atenYsSPly5ePW0qWLJl8WBHzBVPdU6dOwb//quXUqfhKo2hR8PKCNm3Up5cXuBkTL2NilJIQAj7+mBpTpyqfj3ffhf79tfWhSRXmKpItwDdCiAtSyoOxO4UQ9YGJqHETjSZNHDlyhGnTprF69Wr++OMPGjduzHfffUf27JafXBgZGUlgYCD+/v5xCuLvv/9mypQpccriwYMH8b6TO3duSpcuTZkyZWjdujWlS5eOUxQeHh4UKVLEcopOSjVj6tw5pSRiFcfp02AqV+nSUK0atG4N3t5KaZQqFd+fIyYGDh6EtWvh119h6VI1kN63L6cKF6bGRx9BzpyWkVuTJTH3HzoC2Aj8JYQIRA22FwNcgcPGcY3muUkYgdfR0ZHhw4dToUIFgOdWIlJK7t27x40bN5Jdbt68SUI/VkdHR8qXL0+FChV49dVXKVOmDGXKlIlTHkWLFrW8RXT7Npw/D35+6jN28fODR4+elnNygho1oHdvqF5dLdWqqci6SREaCp9+qhwIAwIgRw41kB6rNCpX5vZLL2klokkz5jokhgD1hRCtgLqonCSBwCHD70OjeS6klAghiIiIoH///uTJk4fp06fTt29fHJMI7vfkyROCg4MJDAxMVkk8fPjwme86OTlRsmRJSpYsSZ06deLWS5UqFacsfH19zffyNpeICLh+Ha5cgcuX1eeVKyrV7PnzcMckuo+Dg8rZUaGCGquI9d2oXh2KFUvZazw8XPl+PH4MnTqpAXUfH6hfH9q1U2ltCyXVQ63RpJ7net2TUm4CNllJFk0W4ObNm8yZM4ft27fz119/kTNnTtauXUu+fPm4desWPj4+BAUFERgYSFBQUNwSGBgYb2psLDly5MDV1ZVSpUrh6enJm2++GackYhdXV1dy585tnQY9eqSy/F2/rtLFmiqLy5fVMVPLx8FBjVeULQsdOz5VFhUrqn05nkk+mjxXr6oUtps3w65dSol4eSlFkiOHksPBIjnnNJokSVKRCCFySiN9rhAiRdtX6lS7GhMePXoUTxEcP36cTZs2cfLkSaKjoylYsCBubm7cunWLyMjIZ76fJ08eSpQoQfHixalcuTJNmzalePHicUusknB2drbMIHZCpIS7d5Ui8Pd/9jN2PaFyE0KNUbi7wyuvqE93d6Uk3N2hZMnnVxamREaqfB8vvqjONXYsLF+uvMwHDIBWrcDUqtJKRJMOJGeRhAshXpRSHkbNykopMKL+xWZiYmJiCA0N5ebNm8kuwcHBBAUFxfPENiVPnjyUKVMGd3f3OKUQqzBMlwIFClhnhlZkpEq0FBSk0r2afFY5dUoFH7xxQy3h4fG/KwS4uChlUK6ccuArWVIpjpIllaJwc7P8mENICGzdqqyO7dvh3j3laV6xInz2mYp3VbGiDpiosRnJKZIBqFwfses6wm4m49GjR4SEhMQpAFOFcOrUKSZOnBi3HRISEs8vwhQnJyeKFStGsWLFqFWrFs2bNyckJIQ8efLQuXNnihUrxsaNGxkwYAAlSpSwfEMeP1ZhzW/dUkohgYKI93nrVuJ1FChAAUdHZTl4eUHbtvGVRKlSUKJE2qwJc4mJUQovVy6lOF5/XVlIxYtD+/ZqrCM2DL7ONqixA5JUJFLKBSbr89NHHE1aiIqK4vbt2ylaDbFLWFhYovXky5cPR0dHSpcuTenSpfH29o5TFAkXZ2fnuJlV9+7d4+eff2bmzJn4+/vToEGDON+POnXqmNsINZMpVjEktoSExN9OZHAdgLx51cPXxUW9sTdq9HQ74WfevBzes8fyg+3mcu+einW1ebOyPkaOhI8+gnr14MsvlfLw9NQhSjR2ibmxts4DHaSUJxM5Vg1YJ6WsaGnhsjrh4eGEhIRw69atJD9j10NCQrh9+/YzU1oBHBwcKFq0aNzDv1y5cri4uCSqGIoWLUq+fPnY85wP1SVLljBkyBAePHhAkyZN+OH773njpZcQFy6ocYQ7d+J/xq4nVBqJDKjHUaCAcrZzdlYP/2rV1Lqz89P9RYs+VQ4ZIZ1qTAy0aAF79iglWrAgtGyp4lyByvXx+ec2FVGjSQlzZ215AElNe8kHuFtEmkxMdHQ0d+/eTVIRJPb5yNSPwIRs2bLh7Owct1StWhVnZ+cklUPhwoVTNyAtpXJ+S0QRyDt32HviBG5SUi4qivKXL9Mqe3ZGliiB17Fj8NZbydedL596SMYqgjJlnq6bKobYxckp44csDwtTIdq3bIH792HFCmVhlCunnAnfeEMNolvBAVOjsSbJzdrKD5hO6HcWQrgmKJYb6MSz6XAzPbHjC+ZYC7du3eLOnTvExGaUS0D+/PlxdnaOsxqqVasWtx2rLGLXixYtSqFChcxTDI8fqy6TCxfUZ2io+jRdktjXILbLKDbshsETVGC1GcBxYFjevMwoXZqXihThpZdeUuE2Chd++mm6brovKznBrV2rUtLu3av8SgoUUOMeMTFKkfz0k60l1GjSRHKvPiOB8ahBdonybE8MAYyysFx2QWhoKF988QWnT59m0qRJ8RRDeMIZPQYODg7xrIXq1as/owgSKolnfBykVP4JDx6oN9cHD9QD/syZ51IGRJgxI9vRUXWnFCqkPkuUgMqVuRkWRsnq1eM9/Cdv3860337j5p07VKtShZ+HD6db9+6QJ48FrnYm4dEj1U21ZYsa23ByUr4cAQHwwQfK6mjQIGspUk2mJzlFshL10imAtcAnwPkEZSKAc1LKy9YRz/YsXLiQ/Pnz4+bmRvHixalRo0bi1kKRIhTNm5eC2bKRLSwsvhIwXf/338T3m+578OBpPuzkyJcvvhJwdlb+BLHbsUtS2wUKJDl467dnDyWbNOHMmTNUqVIFIQSXf/+dei+9xLBhw3j11VfTLYCi3XP7tuqm2rIFdu9WlmDevMrhsHFj+PBDGKGjCGkyL8nN2voP+A9ACNEcFQ4l8Wk+mZRCDg7cHz2aq6dPU6ZQIfWADw1VXswJH/5hYfE9mJPCwUE9wB0d1WeBAuqhXqrUs/sTljFVAo6OVutLj46OZu/evYwfP569e/dy4MABXnzxRebMmWMd57+MRlQU/P23UuR16qhxo6FD1cywgQOV1fHyyypZFOiZVppMj7lPomigHbA44QEhRA/gmpTyT0sKZhdERMCnn+KWPfvTN/jYh7uTk/I5SOyhn9x6bKpSO+TRo0f8+OOPzJ49mytXrlCmTBmmTZtGFcNXIUsrkfv3lU/Hxo3K8rh9Gzp3VpZIhQpw6ZL6PWg0WRBzFclEwCeJY8WBQUADi0hkTxQpAo8fs/fvv23nX5AO3Lt3j4IFCyKE4Ouvv6Z69eq8++67jB071ioh3DMMd++q8SGAhg1VKPciRZTF0aaNmrYbi1YimiyMuU+J6sC4JI4dAz61jDh2hhAZf8ppEkRERLB27VrmzJlDcHAw586dI0+ePJw5cwYXFxf27NmT9ZRIbN6OjRvVcv063LypvNknTVIWpZ6eq9E8g7n/iBggqXygTkAW7vPIWAQGBvLjjz8yb948goKCKFeuHIMHDyYqKoqcOXPi4uJiaxFtw6pVvDRokLJCHByUBdK7t+rezJFDeZZrNJpEMVcB/AWMFELECzRkbA8H9ltaMI3lkFISYUwF/vvvv/nqq6+oU6cOW7Zswc/Pj5EjR5IzK01HvX0b/u//VDyt/cZP182NUE9PWLZMWSF79qiZVvny2VJSjSZDYK5FMhalLPyEECtQSa1KAJ2BIkAj64inSQthYWEsWbKEuXPn0qFDB8aPH0+bNm3w8/OjfPnythYvfQkPh/nzYd065RgYHa0i9d68qY6/+CJnxo2jWCYeC9NorIVZFomU8gRQHzgC9Ae+Mz4PAy8kFoNLYztOnDjB4MGDcXV1ZfDgweTMmZOqVasCKnVtllEiZ8+qZE+gxjXGj1cRgEePBl9flRSqXTvbyqjRZALMHjWUUp4GOlpRFk0aePLkCbmMiQFffPEFW7dupVOnTgwePJgXXnghazgPxsSopE/r1qk85f/9Bx4eKv95jhxw7pxKWavRaCyKHiTP4Jw8eZIhQ4bg4uLChQsXAJg+fToBAQEsXryY+vXrZ24lYpojZcgQlZ986lTl4DlnjhrriEUrEY3GKphtkQgh6gJ9gYokEglYSvmSBeXSJEN4eDgrV67kp59+4tChQzlgYXYAABN6SURBVOTKlYuOHZ8ai2Uzu0/Do0fKOXDdOti0SXmZV6oE3bur6bmtWil/D41Gky6Ym4/kVWAb8CfQGNgJ5EGNm1xHzerSWJmHDx+SL18+Hj58yKBBgyhXrhzfffcdPXv2pEhWeHBevqxmUm3frgbPCxeG1q2fHm/QQC0ajSZdMdci+QqYDXwMRAJjpJTHhBDlUApmu5Xky/KEhoaycuVKFi1ahIODAwcOHMDZ2ZmTJ09SsWLFzN1tdeOGGuuITTFbuDCcOAF9+sDbb6t4VumR+laj0SSLuYqkGvA5yjFRopJZIaW8JIQYjwo3v8wqEmZRDh8+zMyZM1m7di2PHz+mRo0a9OnTh5iYGLJly0alSpVsLaJ1uHBBdVmtXau8zAG6dlWKpFAhuHjRbmOVaTRZFXMVyRNASCmlECIQKAvsM46FAm7WEC6rceXKFYoUKYKjoyPHjh1jy5Yt9OnThz59+lCnTp3MaX1IqZSDh4fa7t9fDZB7ecHXXyvLwwgaCWglotHYIeYqkhNAJdTYyG5gjBDiOiofyZfAaeuIl/kJDw9n3bp1LFy4kD/++IO5c+fy3nvv0atXL3r37v1s0qvMQGxMq7Vr1XLtmvLvcHKC775TXVhlythaSo1GYybmKpKZKCsEYAywGfjD2A4A3k6rIEKIIsCvqPzvV4BOUsq7Ccp4Aj+gUgBHA19LKX9N67ltQVRUFB988AHLly/n3r17lC1blv/973+0atUKgDyZNevgzp3QowcEBanxjWbNYOzYp8ExPT1tK59Go3luzFIkUspNJuv+QojaKAslD3BaSvnEArKMBv6QUn4jhBhtbH+SoMwjoKeU0s/IH39UCLFdShlqgfNbHT8/P3x9fenSpQvZs2fn/PnztG7dmj59+tC4cePMl+/j0SP4/XdldbRrB2+9BeXKqYCI7dqpcOwFC9paSo1Gk0ZSVCRCiNyoUPHDpZTbAaSUMcBZC8vSFmhirP8C7CGBIpFSnjdZDxBC3ASKosZp7JKQkBB+/fVXli5dyqFDh8iTJw+tW7cmf/787NixI/ONe8TEqGRP69bB1q1KmRQuDC+8oI6XKwe//WZbGTUajUVJ8RVYSvkYcEbN1rImLlLKQOOcgUCybshCiHpATuCileVKNUuXLsXV1ZWhQ4fy+PFjpkyZgp+fH/nz5wfIPEokKAj+MHo6s2WDCROUk2Dv3qorKzhYeZ1rNJpMiZBm5BkXQswEnKSU3dN0MiF2ojIqJuRT4BcpZSGTsnellInmQBFClEBZLL2klAeTKDMAGADg4uLitXLlylTLHRYWFvfwT4ro6GhOnDjBzp07afL/7d17kJX1fcfx9ycCErLeigRJvIWqxSs4EMc4YkDRGLHqEGXNzdjxglqTqZNOteK0WqcZTbxkLCbWW1HTiI5KRbBEwkW3TlvFVnGNosZLEdcrKoFMiLDf/vF7Vo6Hc3bP7nMuu3s+r5kz5zy/32+f5/vdA+e7z+U8vylTOOyww1i7di0LFy7k2GOPZezYsX3efiP0lPPwN99k18cfZ9e2NnZqb2fLiBE8Pn8+MXQow959lz+OHDng5iqv5H0ebJxzc+hrzlOnTn0qIib1ODAienwAPwDWAv9FmilxFulDuutxTiXr6WEbq4Ex2esxwOoy43YkHWo7rdJ1T5w4MfJYvnx5yfbOzs5oa2uLCy+8MEaPHh1AtLS0xI033phre/3BNjlv2ZIeERE//nFEunA3Yvz4iCuuiFi1KqKzs+5xVlO593kwc87Noa85Ayujgs/YSq/a+mn2PAY4rFQ9Am6pcF3lLAC+B1yVPT9YPEDSMGA+cGdENORAe0SwZs0a9txzTwDOOOMMOjo6OPHEE2ltbeWEE05gxIgRjQit+jZtSrdhX7AgPe66C44+Go47Lt2W/eST0zkPM2tqlRaSetyH4irgXklnAf9Hdst6SZOA8yLibGAmcBQwUtKZ2c+dGRFP1zKwiOCZZ55h3rx53HPPPaxfv56Ojg6GDh3K/PnzGTt2LDvssEMtQ6iv997jgMsvh6eegg0b0iyBX/sadO0ajx+fHmZmVH7575ZaBxIR7wPHlGhfCZydvf4F8Itax1Jo0aJFnH/++axZs4btttuOadOm0draSmdnJwDjB8MH6quvpj2O4cNh1izYeWdGrFmTbk1y8slpL2QwfjHSzKqibCGR9Ajw/YhYXdB2NPDfEbGxHsH1BzvuuCMjR45k9uzZzJgxg1GjRjU6pOp4+mm4/3548EF49tnUNn16KiRDhrDyttuY4mlnzawC3V1WMw345NtikrYDlpC+iNg0Jk+ezPXXX8+sWbMGdhFZvz7N3dF1ld6118KPfpS+43HttelmiQsXdr8OM7MSKp7YKjNIvvjQBCLgxRdh0aL0aGuDjz9O85iPGwdXXpnua7Xrro2O1MwGuN4WEuvPNm1KxaKlBR56KJ3fADjwQLjoonToqusuu3vv3bAwzWxw6ekbY6W+rVjrb7hbb7z5Jtx6a7rd+siR8LOfpfavfjXNWf7aa9DeDldfnSaCGuK/Hcysunr6VPmVpM1FbUtLtBER3d7SxKokIs3J0dkJhx8OTz6Z2vfYI91V98gj0/JOO8EFFzQuTjNrGt0VkivqFoWV13Wu45FH0lzlmzfD4sXp9iNdd9GdPh0OOsiTPplZQ5QtJBHhQtJoc+bANdfA66+n5X32SUWja6/kuusaG5+ZGT7Z3j9s3py+Rb5kSdrzuP9+GDUqTfw0YQJcfHH6ZrlvR2Jm/ZALSSO1t8Ps2WmO8vXrU9vEidDRkQrJrFnpYWbWj7mQ1Murr6Y5O5YtS1dYnXZauu1Iezu0tsIxx6RbkQzkLz2aWVNyIamlzZvhvPNSAXnttdS2224weXJ6vc8+8Nt+Oy+XmVlFXEiqpaMDHnssPYYNS98aHzIkfZN8wgT44Q/TXse4cb66yswGFReSvK6+Gm67DV56KS23tMAJJ2ztf/zxxsRlZlYnLiSViIDVq+HRR9MexxNPwKpV6RzHli2w//7ppPhRR8Ghh/rb42bWVPyJ15P58znirLPggw/S8ujR6fYjH36Yzndcemlj4zMzazAXkp7suSfrJk1it5kzUwHZZx+f4zAzK+BC0pOJE3nh0kvZzZM8mZmV1NPdf83MzLrlQmJmZrm4kJiZWS4uJGZmlosLiZmZ5eJCYmZmubiQmJlZLi4kZmaWiwuJmZnl4kJiZma5uJCYmVkuLiRmZpaLC4mZmeXiQmJmZrm4kJiZWS4uJGZmlosLiZmZ5dJvComkP5G0RNJL2fMu3YzdUdJaSXPqGaOZmW2r3xQS4BJgaUTsCyzNlsu5Eni0LlGZmVm3+lMhORm4I3t9B3BKqUGSJgKjgUfqFJeZmXVDEdHoGACQ9GFE7Fyw/EFE7FI05jPAMuC7wDHApIi4sMz6zgXOBRg9evTEefPm9Tm2DRs20NLS0uefH4icc3Nwzs2hrzlPnTr1qYiY1NO4IX2Kqo8k/RrYrUTX7ApXcQHwcESskdTtwIi4GbgZYNKkSTFlypReRPppK1asIM/PD0TOuTk45+ZQ65zrWkgiYlq5PklvSxoTER2SxgDvlBj2FWCypAuAFmCYpA0R0d35FDMzq6G6FpIeLAC+B1yVPT9YPCAivt31WtKZpENbLiJmZg3Un062XwUcK+kl4NhsGUmTJN3a0MjMzKysfrNHEhHvk06gF7evBM4u0T4XmFvzwMzMrFv9aY/EzMwGIBcSMzPLxYXEzMxycSExM7NcXEjMzCwXFxIzM8vFhcTMzHJxITEzs1xcSMzMLBcXEjMzy8WFxMzMcnEhMTOzXFxIzMwsFxcSMzPLxYXEzMxycSExM7NcFBGNjqHmJL0LvA7sBHxU0NXdcuHrXYH3qhBK8fb6OrZcX6l25+yci5eds3Ou1F4RMarHURHRNA/g5kqXi16vrMX2+zq2XF+pdufsnJ2zc65WzuUezXZo66FeLBf31WL7fR1brq9Uu3N2zsXLzrk6mjHnkpri0FZeklZGxKRGx1FPzrk5OOfmUOucm22PpK9ubnQADeCcm4Nzbg41zdl7JGZmlov3SMzMLBcXEjMzy8WFxMzMcnEhqZCksZJuk3RfQdvnJN0h6RZJ325kfLVQJudt2gaTMjmfkr3HD0o6rpHx1UKZnPeXdJOk+ySd38j4aqHcv+Ps//RTkk5sVGy1UuZ9niKpLXuvp/R13U1dSCTdLukdSe1F7cdLWi3pZUmXAETEKxFxVtEqZgD3RcQ5wEl1CjuXvDmX+T30a1XI+d+y9/hMoLVugedQhZyfj4jzgJnAgLhUtgr/nwEuBu6tR7zVUIWcA9gADAfe6GscTV1IgLnA8YUNkrYDbgS+DhwAfFPSAWV+fndgTfZ6S41irLa55Mt5IJpLdXK+LPuZgWAuOXOWdBLwH8DS2oVZVXPJkbOkacBvgLdrG2ZVzSXf+9wWEV8nFdAr+hpEUxeSiHgMWFfUfBjwcla9/wjMA04us4o3SMUEBsjvsgo5Dzh5c1ZyNfDvEfE/tY22OqrxPkfEgog4AhgQh22rkPNU4HDgW8A5kvr9/+m8OUdEZ/byA2D7vsbR739RDfBFtu5lQCoWX5Q0UtJNwKGS/jbrewD4hqSfU8fbEdRAxTmX+T0MRL15n78PTANOlXReneOspt68z1Mk3SDpn4GHGxBrtVScc0TMjoi/An4J3FLwITvQ9OZ9npG9x3cBc/q6wSF5oh2kVKItIuJ94Lyixo3AX9QlqtrqTc7btA1Qvcn5BuCGukRVW73JeQWwog4x1VrFORd0zq1pRLXXm/f5AdIfxLl4j2RbbwB7FCzvDrzZoFjqxTk758HKOdchZxeSbT0J7CvpS5KGAacDCxocU605Z+c8WDnnOuTc1IVE0t3AfwJ/JukNSWdFxGbgQuBXwPPAvRHxXCPjrCbn7Jxxzs652nH4po1mZpZHU++RmJlZfi4kZmaWiwuJmZnl4kJiZma5uJCYmVkuLiRmZpaLC4mZmeXiQmJmZrm4kNigJOlySSHppTL9L2f9l9c5tIpIuj6L79kSfTtLWpf1/3UVtjVT0luSlC2fLmmTpKF5123NwYXEBrM/AF+S9KkZ/iR9Gdgr6++vDgZ+B+yXTVRU6G+AYdnrbQpNH0wHHo6tt7kYD7RHxMdVWLc1ARcSG8w2AstIN60rdHrWvrHuEVXuYNKN9oYBY7saJX0e+AFbb8K3Ks9GssmbjgcWFTRPAP43z3qtubiQ2GA3D5hZcNhGpHnI5xUOkvQVSQskvSlpo6SnJW0zM6CkAyUtzg4tbZT0vKS/7KmvN7Ji8XlgIWmvZFxB92Wk4vE68F5EdPR2/UW+DOwCLCloGw88L+lHktZK+kjSLQNhxkBrDP/DsMHuAWA0cGS2PBkYBcwvGrcX8DhwNvDnwP3Av0j6ZtG4BcAW4DvAScA/ATtU0Ncbh2TPq0hziO8PIGlPYBZwaTamWoe12iJifbaNUcAY0qyQnwXOBK4h/V4GzfTLVl2eIdEGtYj4UNJi0uGstux5cdZeOO6TPZRsr+Ux0oRA5wB3Z+27kg4znRIRXR/iS3vq64NDgE3Ai8BzZIUEuBx4LCJWSLqTKsxsRyok/1qwPCF7viEirsteL5F0PrBvFbZng5D3SKwZzCPNt749cCpFh7UAJO2SzVH+OvBx9jgX2K9g2DrSXNg3SWrNDkFV0tdbBwPPZ/NKPAeMkzQO+C4wW9JOpBnw8p4fGQMcyqfPj4wHPiTtTXWNE7Az8F6e7dng5UJizWAB0AL8I/A54KESY+YCrcBPgONI5w5uB4Z3DYiIzqzvrazvLUltkg7trq8P8R7M1sNWXXskVwKLIuKJrB/KHNqSdICk2yUtlzRH0hfKbOcE4JWIWF3QNh54tOiKrT8lHeYaNBNCWXW5kNigFxEbSSeuLwIeypY/IWk46RDP30fEnIhYFhErKfH/IyJeiIhvkP5Cn0YqNIskfaa7vkpjzcYewNYi0Q7sBMwgnWiHdOirkxIf7JKOBm4Bbiadp1kILMz2aIpN59N7I5AObRXv6YzPttdeaR7WXHyOxJrFz4HtgZtK9G0PbEc6LwGApB1IH8QlpxDN/mJfJuk64Jek4rGup74K7Ev66//ZbF1rJd0L/CYiuj7IDwZejojfl/j5y4BTC67mWixpXdb+nYL8hpGK3alFbeOAZ4rWeQhpz6U/Xy5tDeRCYk0hIlYAK8r0fSTpSeDvJK0n/fV9CfARsGPXOEmHkK5gugd4hXTZ7MWkD97dJc0r1RcR67KfnwIsB6Zm8ZSyzWGriGgtMabcFVtDI6JD0t7ACxExPCKeyJYLHUXa43q0oO1A0mdC8R7JISXazD7hQmKWfIt0OOhO4H1gDjACuLBgzFvA28Bs4Aukk9LLSQXjD930dRmRPb/TTRwHAx9ExNpuxhwEPFKmT5I+C6wlu4xY0l5su0c0Hfh1RGwqaBsP/B74bdHYQ0i/F7OStPWuCGZWS5KuAI6KiKk13MbZpENW50bE+uxE+93AP0TE0oJxLwI/iYhbahWLNQ/vkZjVzxHAdT2OyiEibpX0O2BpdtPFDcDsiHi0aNx+JVdg1gfeIzEzs1x8+a+ZmeXiQmJmZrm4kJiZWS4uJGZmlosLiZmZ5eJCYmZmubiQmJlZLi4kZmaWy/8DS1mXL5xJ5NMAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Default mass function object.\n", "mf200c_conv = hmf.MassFunction(\n", " hmf_model=\"Bocquet200mDMOnly\",\n", " mdef_model=\"SOCritical\",\n", " mdef_params={\"overdensity\": 200},\n", " disable_mass_conversion=False,\n", ")\n", "mf500c_conv = hmf.MassFunction(\n", " hmf_model=\"Bocquet200mDMOnly\",\n", " mdef_model=\"SOCritical\",\n", " mdef_params={\"overdensity\": 500},\n", " disable_mass_conversion=False,\n", ")\n", "mf200c_refit = hmf.MassFunction(hmf_model=\"Bocquet200cDMOnly\")\n", "mf500c_refit = hmf.MassFunction(hmf_model=\"Bocquet500cDMOnly\")\n", "\n", "m = mf200c_conv.m\n", "\n", "fig, ax = plt.subplots(1, 1, sharex=True)\n", "ax = [ax]\n", "ax[0].plot(m, mf200c_conv.dndm / mf200c_refit.dndm - 1, color=\"k\", label=\"200c\")\n", "ax[0].plot(\n", " m,\n", " mf200c_conv.fsigma / (mf200c_refit.fsigma / mf200c_refit.hmf.convert_mass()) - 1,\n", " color=\"r\",\n", ")\n", "\n", "ax[0].plot(m, mf500c_conv.dndm / mf500c_refit.dndm - 1, ls=\"--\", color=\"k\", label=\"500c\")\n", "ax[0].plot(\n", " m,\n", " mf500c_conv.fsigma / (mf500c_refit.fsigma / mf500c_refit.hmf.convert_mass()) - 1,\n", " color=\"r\",\n", " ls=\"--\",\n", ")\n", "\n", "ax[0].set_xscale(\"log\")\n", "ax[0].set_xlabel(r\"Mass, $M_\\odot/h$\", fontsize=15)\n", "ax[0].set_ylabel(r\"dn/dm\", fontsize=15)\n", "ax[0].grid(True)\n", "ax[0].legend()\n", "ax[0].set_ylim(-0.5, 0.5)\n", "ax[0].set_ylabel(\"Fractional Diff.\")\n", "ax[0].set_title(\"Comparison of Mass Conversion to Refit\")\n", "plt.savefig(\"/home/steven/Desktop/comparison_of_massconv_to_refit.pdf\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the conversion is only accurate to ~20% at lower masses, and blows up at high masses. The red lines here are directly using the fits for $m'/m$ given in Bocquet et al (2016) (instead of actually converting mass definitions using a halo profile). There is some discrepancy between that fit and the conversion, resulting in the discrepancy here. " ] } ], "metadata": { "celltoolbar": "Initialization Cell", "kernelspec": { "display_name": "Python [conda env:hmf]", "language": "python", "name": "conda-env-hmf-py" }, "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.0" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 1 }