{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 4 - Model the weak lensing signal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This session assumes that you have your own weak lensing signal in order to model them. We understand coding up the lensing measurement part(Day 3 tutorial) is quite tedious. So in case if you don't have the signal from yesterday's session, we will provide you the output file for modelling and getting the halo mass constraits out of it. Similar to day 3 tutorial, the reader has to fill in the blanks(??) in the code snippets.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In last tutorial/lecture, we have discussed and coded up the R.H.S. of the equation given below, which is nothing but the **observed weak lensing signal** (excess surface mass density).\n",
    "\n",
    "$$\\Delta \\Sigma (R) = \\frac{\\sum_{ls} w_{\\rm ls} e_{t,ls} (\\Sigma^{-1}_{\\rm crit})^{-1}}{2\\mathcal{R} \\sum_{ls} w_{\\rm ls}}$$\n",
    "\n",
    "Now our aim today is to predict $\\Delta\\Sigma(R)$, which can be expressed as -\n",
    "\n",
    "$$\\Delta\\Sigma(R) \\equiv \\bar{\\Sigma}(<R)-\\langle \\Sigma \\rangle (R)$$\n",
    "\n",
    "where - \n",
    "\n",
    "- $\\bar{\\Sigma}(<R)$ : average projected mass density of the matter (dark+baryon) within $R$.\n",
    "- $\\langle \\Sigma \\rangle (R)$ : azimuthally averaged projected mass (dark+baryon) density at $R$.\n",
    "\n",
    "---\n",
    "\n",
    "To model the $\\Delta\\Sigma(R)$, we need to compute the surface mass density, from the 3-dimensional model of halo mass density profile given by **Navarro–Frenk–White (NFW) profile** [arXiv:astro-ph/9508025](https://arxiv.org/abs/astro-ph/9508025) by integrating this along the line-of-sight.\n",
    "  \n",
    "$$\\mathbf{NFW\\ profile:} \\hspace{0.5 cm} \\rho(r) = \\frac{\\rho_0}{\\left(\\frac{r}{r_s}\\right) \\left(1+\\frac{r}{r_s}\\right)^2}$$\n",
    "\n",
    "$$\\rho_0 = \\delta_{c} \\rho_{m,0} \\hspace{0.2 cm} , \\hspace{0.2 cm} \\delta_{c} = \\frac{200}{3} \\frac{c^3}{\\ln{(1+c)}-\\frac{c}{1+c}}$$ \n",
    "\n",
    "\n",
    "If we define $r_{200m}$ to be the radius at which the average density of the halo turns out to be 200 times the present **mean matter density** of the universe, $\\rho_{m,0}$, then -\n",
    "\n",
    "$$M_{200m} = \\frac{4 \\pi r^3_{200m}}{3} (200\\rho_{m,0})$$\n",
    "\n",
    "$r_{200m}$ is related to the concentration parameter by, $$c = r_{200m}/r_s$$\n",
    "\n",
    "\n",
    "Following [Saas-Fee lectures,](https://link.springer.com/book/10.1007/978-3-540-30310-7), this line-of-sight integrated NFW expression i.e. surface mass density, at a given projected radius **R** is given by -\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\Sigma(R) &= {2r_s\\rho_0}\\ f(r/r_s) \\\\\n",
    "f(x) &= \\frac{1}{x^2-1}(1-\\mathcal{F}(x)) \n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\mathcal{F}(x) = \n",
    "\\begin{cases}\n",
    "\\frac{{\\rm acosh}(1/x)}{\\sqrt{1-x^2}}, & \\text{for } x \\lt 1 \\\\\n",
    "\\frac{{\\rm acos}(1/x)}{\\sqrt{x^2-1}}, & \\text{for } x \\gt 1 \\\\\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "$$\n",
    "f(x) = \\frac{1}{3} \\hspace{0.3 cm} , \\hspace{0.3 cm}  \\text{for } x = 1\n",
    "$$\n",
    "\n",
    "\n",
    "**Notice :** NFW profile is a spherically symmetric profile. Hence its line-of-sight (LOS) integration will be symmetric w.r.t an axis passing through the cluster centre and perpendicular to the lens plane. Hence,\n",
    "\n",
    "$$\\Sigma(R) = \\langle \\Sigma \\rangle (R)$$\n",
    "\n",
    "You should use the above fact while implementing the modelling.\n",
    "\n",
    "**Notice:** We are using cluster centre as a proxy for the halo center. But there can always be some offset between them.\n",
    "\n",
    "And the surface mass density within a projected radius, $(<R)$ is given by - \n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\bar{\\Sigma}(<R) &= {2r_s\\rho_0}\\ h(r/r_s) \\\\\n",
    "h(x) &= \\frac{2}{x^2} \\left[ \\mathcal{F}(x) + \\ln\\left( \\frac{x}{2} \\right) \\right] \\\\\n",
    "h(1) &= 2 \\left[ 1 - \\ln (2) \\right]\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "**Thus these expression show that the whole dark matter distribution can be parametrized in terms of only two free parameters, namely - concentration (c) and Halo mass** ($M_{\\rm 200m}$). We suggest the readers to show that these equations are equivalent to the ones described in the lectures.\n",
    "\n",
    "**Note :** Weak gravitational lensing signal probes the total matter distribution around the lens systems (here: galaxy cluster systems). In clusters, the baryonic component only contribute dominantly in the inner regions (small projected radii), where the signal itself is very noisy, resulting in low signal-to-noise ratio. This happens because for smaller projected radii, the radial bin width becomes small due to logarithmic nature of radial binning which reduces the lens-source pair counts in such bins. Plus, the clusters are very bright in the inner regions, due to which a lot of galaxies are not identified by photometric surveys, leading to lack of sources in the inner regions. Thus these signals in the inner regions are not supposed to be very reliable. Thus we altogether ignore the baryonic contribution in the modelling of the weak lensing signals. Moreover, the inner regions don't contribute a lot in the determination of the average halo masses around the lenses.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will first define the class where we will code up all the useful function in it. We will create a class instance then use them to call these functions for model prediction. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have defined two classes one consisting of useful **constants** named constants and the **halo** class contains the functions useful for model predictions. The halo class has couple of functions define in it which are described as \n",
    "\n",
    "- **__init__**: Here we provide initialization parameters for the modelling - for us they are logarithm halo mass **logM200m** and concentration **con_par**. This also sets up other constants required in the other functions.\n",
    "\n",
    "- **nfw** $\\rho(r)$: This will give the standard NFW dark matter distribution for given radial bins.\n",
    "\n",
    "- **sigma_nfw** $\\Sigma (R)$: This will provide the projected NFW dark matter distribution for given radial bins.\n",
    "\n",
    "- **avg_sigma_nfw** $\\bar{\\Sigma}(< R)$: This will provide the mean projected NFW dark matter distribution within the given radial bins.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Please note that the here we are using Halo masses in units of** $h^{-1}M_\\odot$ **and distances in units of** ${\\rm h^{-1}Mpc}.$\n",
    "\n",
    "Useful functions: [arccos](https://numpy.org/doc/stable/reference/generated/numpy.arccos.html), [arccosh](https://numpy.org/doc/stable/reference/generated/numpy.arccosh.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# working with the classes\n",
    "# NFW  --> delta_sigma\n",
    "class constants:\n",
    "    \"\"\"Useful constants\"\"\"\n",
    "    G     = 4.301e-9   #km^2 Mpc M_sun^-1 s^-2 gravitational constant\n",
    "    H0    = 100.       #h kms-1 Mpc-1 hubble constant at present\n",
    "    omg_m = 0.315      #value of omega_matter is consistent with our signal measurement code\n",
    "\n",
    "# Here halo class inherits the constants class \n",
    "class halo(constants):\n",
    "    \"\"\"Useful functions for weak lensing signal modelling\"\"\"\n",
    "    \n",
    "    def __init__(self, log_M200m, con_par):\n",
    "        self.m_tot   = 10**log_M200m # halo Mass h-1 Msun\n",
    "        self.c       = con_par       # concentration parameter\n",
    "        \n",
    "        self.rho_crt = 3 * self.H0**2 /(8.0*np.pi*self.G) # rho critical\n",
    "        \n",
    "        self.r_200   = (3 * self.m_tot /(4*np.pi*200*self.rho_crt*self.omg_m))**(1./3.) # R200m\n",
    "        \n",
    "        # rho_0 is numerator for NFW profile\n",
    "        self.rho_0   = self.c**3 * self.m_tot/(4 * np.pi * self.r_200**3 *(np.log(1 + self.c) - self.c/(1 + self.c)))\n",
    "        #print(\"log_M200m = %s h-1 M_sun, c = %s \\n\"%(np.log10(self.m_tot), self.c))\n",
    "        \n",
    "    def nfw(self,r):\n",
    "        \"\"\"given r, this gives nfw profile as per the instantiated parameters\"\"\"      \n",
    "        r_s    = self.r_200/self.c\n",
    "        value  = self.rho_0/((r/??)*(??+r/??)**??)\n",
    "        return value\n",
    "    \n",
    "    def sigma_nfw(self,r):\n",
    "        \"\"\"analytical projection of NFW\"\"\"\n",
    "        r_s = self.r_200/self.c\n",
    "        k   = 2*r_s*self.rho_0\n",
    "        sig = 0.0*r\n",
    "        c = 0\n",
    "        for i in r:\n",
    "            x = i/r_s\n",
    "            if x < ??:\n",
    "                value = (1 - np.??(1/x)/np.sqrt(1??x**2))/(x**??-1)\n",
    "            elif x > ??:\n",
    "                value = (1 - np.??(1/x)/np.sqrt(x**2??1))/(x**??-1)\n",
    "            else:\n",
    "                value = 1./3.\n",
    "            sig[c] = value*k\n",
    "            c=c+1\n",
    "        return sig\n",
    "    \n",
    "    def avg_sigma_nfw(self,r):           \n",
    "        \"\"\"analytical average projected of NFW\"\"\"\n",
    "        r_s  = self.r_200/self.c\n",
    "        k    = 2*r_s*self.rho_0\n",
    "        sig  = 0.0*r\n",
    "        c=0\n",
    "        for i in r:\n",
    "            x = i/r_s\n",
    "            if x < 1:\n",
    "                value = np.??cosh(1/x)/np.sqrt(1-x**??) + np.??(x/2.0) #numpy default log is base e.\n",
    "                value = value*2.0/x**??\n",
    "            elif x > 1:\n",
    "                value = np.??cos(1/x)/np.sqrt(x**??-1)  + np.??(x/2.0)\n",
    "                value = value*2.0/x**??\n",
    "            else:\n",
    "                value = 2*(1-np.log(2))\n",
    "            sig[c] = value*k\n",
    "            c=c+1\n",
    "        return sig\n",
    "        \n",
    "    def esd(self,r): \n",
    "        \"\"\"ESD profile from analytical predictions\"\"\"\n",
    "        val = self.??nfw(r) - self.??nfw(r) #check with the analytical equation up there\n",
    "        return val\n",
    "\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "Exercise: \n",
    "- Create a instance for the halo class for halo having **logM200m = 14** and **con_par = 5** as: hp = halo(14,5)     \n",
    "- Tell the output of : print(np.log10(hp.m_tot), hp.c, np.log10(hp.rho_crt), hp.r_200, np.log10(hp.rho_0)) \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$\\\\rho(r)$')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEJCAYAAACQZoDoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd1xUV94/8M93qIoIIkUUaWLvNEWNJorpxESJXaMmotmUNdls9smzv2x28zz5JetmsynG2BJLspZI2prdFDF2QARb7CCIEAsggooiCuf5gyEBg0iZmTt35vN+vXxlyp0737y88uGcc885opQCERFRDYPWBRARkXVhMBARUR0MBiIiqoPBQEREdTAYiIioDketC2gpb29vFRwcrHUZRES6kpGRUaSU8qnvPd0HQ3BwMNLT07Uug4hIV0Qk91bvsSuJiIjqYDAQEVEdDAYiIqpDt8EgInEisqS0tFTrUoiIbIpug0EptUEpleDh4aF1KURENkW3wUBERObBYCAiojrsNhgO5Jdg0dYTKLhUrnUpRERWxW6DYdvxQrzxzVHEvP4DnliZju8PncX1yiqtyyIi0pzuZz4319Mju+K+vv5Yn56Pz/bkI+nIOXi3ccbY8ACMjwxAmK+71iUSEWlC9L6DW2RkpGrpkhg3Kquw9XghPk3Pw6YjBbhRpRAe6InxkZ3xQD9/uLs6mahaIiLrICIZSqnIet9jMNRVeOkavtz7E9al5yGr4DJaOTng/r7+mBDVGVHB7SAiJvsuIiKt6CIYRCQBQLZSKqn2cwDFSqk9t/qcqYOhhlIK+/JK8Gl6PjbsP43L124g1NsNE6I6Y2x4AHzcXUz+nURElqKXYAgHEKqUSrw5JBpirmCo7UrFDfznx7NYt/sUdp+8AEeDYFRPX0yMCsTwbj5wMLAVQUT60lAwmH3wWUTiAcxRSo2+6bUSVAfBkno+FgGg2HhcdkMtBkto7eyI+IgAxEcEIKvgMtan5yExIx/fHToHfw9XPBoRgEcjO6OzV2styyQiMgmz366qlEqs/dz4wx61uoxi6/lYCYAkAHsA1Pe+ZsJ82+Cl+3si5aVRWDQ1HN07uOO9zVkY/rfNmLpsFzbsP41rNyq1LpOIqNm0uF01CsA64+NsAOGoDoFYAO1FJAnA6wDGAygGUF+LQnPOjgbc28cf9/bxx+mSq0jMyMe63Xl4Zs1eeLk5Y1x4J0yMDkQXnzZal0pE1CQWGWMQkY01XUkishjAYqXUHmNrYbRS6g9NPF8CgAQACAwMjMjNveVGRBZVVaWwI6sIa9JOYePhc7hRpTAoxAuTogNxb58OcHVy0LpEIiIAGo8x1KMEgFdLTmAcl1gCVA8+m6IoUzAYBMO7+WB4Nx8UXrqGxIx8rN19CvPW7YPnBic8MrATJkUHopsfJ88RkfXSIhh2A/A0Pg4FsLE5JxGROABxYWFhpqrLpHzcXfDknV0wZ3goUrLPY3XaKXySmovlO08iMqgdJkYH4oG+/mjlzFYEEVkXs3clGbuL1gOYXTMQLSIvonpgOVwpNb8l57fE7aqmcv7yNXy2Jx9r0vKQU1SGtq6OGBcRgCmDArkEBxFZlC7mMTSXnoKhhlIKqdnFWJ12Ct8ePIPrldVjEVMGB+Ge3n5wcWQrgojMyyaDoVZX0uzMzEyty2m2osvXsD49H6vTcpFXfBXt3ZzxaGRnTI4ORGB7zosgIvOwyWCooccWQ32qqhS2ZxVh9a5cJB0pQGWVwh1dvTFlUBBie/rC0cFuV0gnIjOwyWCwlRZDfc6WlmPd7jys3X0KZ0rL4dfWBZOiAzE5OhC+bV21Lo+IbIBNBkMNW2kx1OdGZRU2HyvEJ6m52Hq8EI4GwT29O2BaTBAGhXhxpVciajZrm8dAjeToYMDoXn4Y3csPuefL8M9dp/Bpeh7+/eMZdPVtg2kxQXhkYCfuF0FEJqXbFoMtdyU1pPx6JTbsP42PU3NxIL8Ubs4OGBsegKmDg9C9A295JaLGYVeSjdqfV4KPU3Pxr/2nUXGjCtEhXngsJhh39/aDEweriagBDAYbd6GsAusz8vBJ6imcKr6CDm1dMS0mCBOjOqN9G24oRES/xmCwE5VVCluOFWBF8klszyyCs6MBD/XviBlDgtGnk4fW5RGRFbHJYLDXMYbGyiq4hJXJufhsTz6uVFQiIqgdZgwJxr19OrCbiYhsMxhqsMXQsIvl15GYno+VKSeRe/4K/Nq6YMqgIEyKDuS+1UR2jMFAqKpS2HK8ACuSc7HteCGcHQx4aEBHPD4sBD3922pdHhFZGOcxEAwGwcgefhjZww8nCi9jxc6TSMzIR2JGPoZ0aY/Hh4Xgru6+MBg4aY7I3um2xcAxhpYruVKBtbvzsDL5JM6UliPU2w0zhwZjXEQAWjvzdwYiW8auJGrQ9coqfHPwLD7ckYP9eSVo6+qISYMC8VhMMDp6ttK6PCIyAwYDNYpSCntOXcCHO3Lw7cGzEBHc39cfs+8IQb8Az9ufgIh0g2MM1CgigoggL0QEeSGv+ApWpZzE2rQ8bNh/GoNCvDBnRCju7MZxCCJbxxYDNehS+XWs252Hj3bk4HRpObr6tsHs4aEYM6Ajd5oj0jF2JVGLXa+swtcHTmPJthwcOXMRvu4umDE0GFMGBcGjFVd3JdIbmwwG3pWkDaUUdmQVYcm2bGzPLIKbswMmRgdi1rAQdOJANZFu2GQw1GCLQTuHTpdi6bZsbDhwBgAQ188fc+/sgh4dOGGOyNoxGMisfiq5iuU7crA67RSuVFQitqcvnrwzDBFB7bQujYhugcFAFlFypQIrk3OxPDkHJVeuY1CIF566Kwx3dPXmNqREVobBQBZVdu0G1qSdwrLtOTh7sRx9OrXFb+4Mwz29O8CBt7oSWQUGA2ni2o1KfLn3Jyzamo2cojKE+rhh7ogueHhAJzg7culvIi0xGEhTlVUK3x48i/c3Z+HwmYvo6OGKOSO6YEJUZ7g6cS4EkRYYDGQVlFLYerwQCzefQNrJYvi6u2DOiC6YHB2IVs4MCCJLsslg4DwGfUs5cR7vbspESvZ5eLdxRsLwUEwdHMRVXYksxCaDoQZbDPqWllOMdzdlYkdWEbzcnDH7jlBMiwlCGxcGBJE5MRjI6mXkFuPdTVnYerwQnq2d8MSwEEwfEoy2rlxug8gcGAykG/vySvDupkz8cLQAbV0d8fiwUMwaFgx3BgSRSTEYSHd+zC/FO5sykXTkHDxbO2HO8C54bAjHIIhMhcFAuvVjfine2ngMm48VwruNM35zZxgmDwrkba5ELcRgIN3LyC3Gm98dR0r2eXRo64pnRoXh0YjOnChH1EwMBrIZyVlFePP7Y9hzqgSdvVrht6O64eEBHeHowIAgaoqGgoH/mkhXhoR547Mnh2D5zCh4tHLCC+v34+5/bMO/9p9GVZW+f8khshYMBtIdEcFd3X2x4elhWDQ1Ao4OgmfX7MUD7+3A1uOF0HsrmEhrVhMMIpIgIrHGx6Eisl5EErSui6yXiODePh3wzW+H452JA3D52nU89lEapizbhQP5JVqXR6RbVjPGICLhAEKVUokiEgqgWCl123/dHGOgGtduVGL1rlN474csFJdV4IF+/vj93d0R7O2mdWlEVkfTMQYRiReRjfW8FttAi6AYQKjxuFBz10i2wcXRATOHhmDr7+/EsyPD8MORAsS+tRUvf3kQhZeuaV0ekW6YPRiUUom1n4tIvPH1JOPz2Ho+FqmU2gMgG0C8uWsk2+Lu6oTn7+6OrS/eiYnRnbE67RRG/G0z/rHxOC5fu6F1eURWT4sxhihU/8CH8b/hxsexAKJExBNAtrFrKRLAkptPYByPSBeR9MLCQkvUTDrk6+6K/324LzY+Nxx3dffFO5syMWL+ZqxMPonrlVVal0dktSwyxiAiG5VSo42PFwNYrJTaY2wtjFZK/aG55+YYAzXWvrwSvPHNEaRmF6OLjxv++EBP3NXdl/tRk12ytnkMJQC8WnoSEYkTkSWlpaUmKInswYDOnlgzezCWTY+EUsCsFemY9mEajp69qHVpRFZFi2DYDcDT+DgUwMYGjr0lpdQGpVSCh4eHyQoj2yciiO3lh2/nDcefHuyFH38qxf3vbMdLn//IAWoiI0vclRQLILLWoHMiqu84igXgWTMITWRJzo4GzBpWfQfTY0OCsT49D3e9uQULt2Sh/Hql1uURacpq5jE0Fbf2JFM6UXgZr//nCJKOFCCgXSv813098EBff44/kM3iInpEjbQzqwj/8/VhHD17CRFB7fBKXC/0C/C8/QeJdMbaBp+JrNbQMG/8+9k78MbYvsg9X4Yx7+/ES58fQHFZhdalEVmMboOBdyWRuTgYBBOjA/HDC3di1tAQfJqej7ve3IKPU06ikiu4kh1gVxLRbRw/dwmvfHUIKdnn0cu/LV4d0xuRwS2+45pIU+xKImqBbn7uWD17EBZMHogLVyoQvygFz3+6DwWXyrUujcgsdBsM7EoiSxIRPNivIzb9bgSeuqsLvt5/BiPf3Ipl27O5vAbZHHYlETVDTlEZ/rLhELYcK0SYbxu8+lBvDAnz1rosokZjVxKRiYV4u2H5jCgsmx6JazcqMXnZLjy3bh+KLnP2NOmfboOBXUmktZrlNTY+NwLPjgzD1wdOY9Tft2Jt2inuP026xq4kIhPJKriE//7iINJyihEV3A7//5G+6OrnrnVZRPViVxKRBYT5umNdwmDMj++HzILLuP/d7fjbd0e59hLpDoOByIREBOMjO2PT8yPwUP9OeH/zCdz9j23YdpwbSpF+MBiIzKB9Gxf8fXx/rJ49CI4GwfSP0vDsmr2c+0C6oNtg4OAz6cGQLt74Zt4dmBfbFd8ePPvz4LTex/bItnHwmchCThRexh+/+BGp2cUYFuaN18f2RWev1lqXRXaKg89EVqCLTxusfmIwXnukD/blleCet7dhxc4c3tpKVofBQGRBBoNgyqAgfP/ccESHeOHPGw5j/OIUnCi8rHVpRD9jMBBpoKNnKyyfEYW/P9ofmQWXcd872/HBlhO4wXWXyAowGIg0IiIYFxGAjc8Px13dffDXb4/ikYXJOHLmotalkZ3TbTDwriSyFb7urlg0NQLvTw7H6ZKriHtvB/6x8TgqbrD1QNrgXUlEVqS4rAKvbjiEL/edRo8O7vj7+P7o3dFD67LIBvGuJCKd8HJzxtsTB+LDxyJxvqwCD7+/E+9vzuLYA1kUg4HICo3q6Yfv5w3HPb074G/fHcOji1OQU1SmdVlkJxgMRFaqnZszFkwOx7uTBiK7sAz3vbMNq1JOct4DmR2DgcjKPdS/I75/bjgGhbTHn746hMeWp+FM6VWtyyIbxmAg0gG/tq5YMTMKrz3SBxm5F3D3P7bhi735XHOJzILBQKQTItWzpr/57R3o7ueO59btx2/+uQfnuZ0omZhug4HzGMheBbV3w7o5MXjpvh7YdKQA97y9DT8cPad1WWRDdBsMSqkNSqkEDw/e4032x8EgmDOiC/71zFD4uLti1op0/Plfh7hbHJlEk4JBRAaIyFgRGSkiA8xVFBE1To8ObfHlU0Pw+LAQrEg+iTELduLY2Utal0U6d9tgEJEQEXlDRD4AMAGAAGgHYLSILBKR10Uk2LxlEtGtuDg64OUHe2HFzCicL7uGuAU7sCrlJAemqdkaXBJDREYBUEqpHxo8SSOPMwcuiUH0i8JL1/D7xP3YcqwQo3r4Yn58P7Rv46J1WWSFWrIkRnpjftgrpTYByGhOcURkOj7uLlg+IwqvxPXC9swi3PvOdmzPLNS6LNKZBoNBKfXzLT8iMk5E1jXmWCLSjohg5tAQfPX0UHi2csK0D9Pw2r8P49oNDkxT4zRl8DkUwGxzFUJEptXTvy02PDMM0wYHYen2HIxdmIysAu4UR7fXlGBIArBJRJ7gYDORPrg6OeB/Hu6DpdMjf97r4bOMfK3LIivXlGCIBDAewAUA/yUi6cY7lUxCRBJEJPam1/5qqvMT2bPRvfzw7bzh6Bfggd+t348XE/fjagW7lqh+jQ4GpdRSAJ4A9iil5iqlIpVST5qwlnTj+QEAIhJe+zkRtYxfW1f884lBeGZkGNZn5GPM+zuQeY5zHujXmjTBTSm1VymV05TPiEi8iGys57VYEUm4zcdLmvJdRNQwRwcDfnd3d6ycGY3zlyvw0IKd7FqiX2kwGERkbGNPdKtjlVKJNx0Xb3w9yfg89ubP1BwDIFREQhtbAxE1zvBuPvjPb+/4uWvp9+vZtUS/uN3tqp+LyOxbDTiLSLDx/ReUUp838jujAGQbH2cDCDc+jgUQJSKexjDJBruSiMymdtdS4h52LdEvHG93gFJqqYh4AEgQkahab3kA2AjgsyZ2L938w7698Xvm3/S9JQBG13cCYxdUAgAEBgY24auJqLaarqXoEC/MW7sPDy3Yif99uA/GRQRoXRpp6LbBAPw8ee1vJvrOEgBeLTmBUmoJgCVA9ZIYpiiKyJ7d0bW6a+nZNXvxu/X7kZp9Hq+O6YNWzg5al0YaaPTgs4iMEpHvRWS3iHzQgrkMu/FLqyEU1a2OJuN+DESmVdO19Kyxa+mRhTtxsqhM67JIA025KylEKXW3UioKQCKq5zLcdult4+ByZK1B50RUDyrHAvCsGYRuKu7HQGR6jg4GPH93d6yYGY2zF8sRt2AHkg5zEyB70+DqqnUOFBmnlPrsdq9ZiojEAYgLCwubnZmZqUUJRDYtr/gKnvxnBg7+dBFP3xWG50Z3g4NBtC6LTKQlq6vWliQin4rIC7W6kS60tLjmYouByLw6e7VG4twhmBDZGQs2Z2HG8jQUl1VoXRZZQFNmPpcqpcYD2IvqbqQsAH81BgV3cyOyQa5ODvhrfD+8MbYvduUUI+69HTiQz3mntq7Jez4rpTYZl8QIQ/XcgxwAc0SkrcmrawAHn4ksZ2J0IBLnxgAA4j9Iwdq0UxpXRObU6DEGa8Ud3Igsp7isAr9duxfbM4swPjIAr47pA1cn3tKqR6YaYyAiO+fl5owVM6PxzMgwfJqej/hFycgrvqJ1WWRiug0GdiURacPBIPjd3d2xbHokcs9fQdyCHdiZVaR1WWRCug0G3pVEpK3YXn7Y8PQw+Lq7YPpHafhoRw703jVN1XQbDESkvWBvN3z+m6EY1cMXr359GC+sP4Dy61ylVe8YDETUIm1cHLFoagTmxXbFZ3vyMWFJKs5dLNe6LGoB3QYDxxiIrIfBIJgX2w2Lp0Ug69wlPPjeDmTkajb/lVpIt8HAMQYi63NP7w744qmhaO3sgElLUvHp7jytS6Jm0G0wEJF16ubnjq+eGopBoV548bMDeOWrg7heWaV1WdQEDAYiMjnP1s5YPiMKTwwLwcqUXEz7cBfXWdIR3QYDxxiIrJujgwH/78FeeGt8f+w5VYK493bgyJmLWpdFjaDbYOAYA5E+jA0PwPo5MaisUhj3QTI2cn8Hq6fbYCAi/ejf2RNfPT0UYb5tkPBxOhZvPcHJcFaMwUBEFuHX1hXrEmJwf19/vP7NUbyYeAAVNzgobY0ctS6AiOxHK2cHvDdxILr4tMG7mzKRW3wFi6ZGwMvNWevSqBa2GIjIogwGwfOju+GdiQOwL68ED7+/E5nnLmldFtWi22DgXUlE+jZmQCesSxiMKxWVGLswGVuOFWhdEhnpNhh4VxKR/g0MbIevnh6KAK/WmLViN1bs5Aqt1kC3wUBEtqGTZyskzo3ByB5++POGw3iZM6U1x2AgIs25uThi8bQIzBkeik9ST2Hm8t0ovXpd67LsFoOBiKyCg0Hw0v09MT++H1Kzz+PRRcnIv8BtQ7XAYCAiqzI+sjNWzYrGmdJyPLIwGQfyS7Quye4wGIjI6gwJ88bnTw6Bs4MBExan4vtDZ7Uuya4wGIjIKnX1c8cXTw1BN782mPNJBj7akaN1SXZDt8HAeQxEts/X3RVrE2IwuqcfXv36MP78r0OorOLtrOam22DgPAYi+9DK2QEfTI3A48NCsCL5JOZ8nIErFTe0Lsum6TYYiMh+OBgELz/YC6+O6Y0fjp7DhMWpKLhYrnVZNovBQES6MT0mGEunRyKr4DIeWZiMY2e5xpI5MBiISFdG9fTD+rkxuF5ZhfgPkrEzq0jrkmwOg4GIdKdPJw98+dRQdPRshRnL0/Dl3p+0LsmmMBiISJc6erbCp3NjEBHUDvPW7cMi7gpnMgwGItItj1ZOWDkrGg/288cb3xzFXzYc5u2sJsAd3IhI11wcHfDuxIHw93DF0u05OFtajrcnDoCrk4PWpemW1bQYRCRBRGKNj0NFJFxEXhQRT61rIyLrZjAI/vhAL7z8YC98d/gspn24CyVXKrQuS7esJhgApAOoCYFiANnGx17alENEevP4sBAsmBSO/XmliF+UwtVZm8nswSAi8SKysZ7XYkUkob7PKKW4nCIRNcsD/fyx6vFoFFwsx9iFyTh0msvmNJXZg0EplVj7uYjEG19PMj6PvfkzIhJvDIckAPHmrpGIbMvg0PZIfHIIHAyCCYtTsSOTcx2aQouupCj80k2UDSDc+DgWQJRxTGGPiIQbX1ti+RKJSO+6+bnj898MQUC76rkOX+zN17ok3dDirqSbB5PbA4BSan6t12q6kvbUdwJjF1QCAAQGBpq6PiKyEf4e1XMd5qzKwHPr9qPoUgVmDw/Vuiyrp0WLoQQtHFBWSi1RSkUqpSJ9fHxMVBYR2aK2rk5YMSsKD/T1x2v/OYLX/3OEE+FuQ4sWw2780moIBbCxgWNvSUTiAMSFhYWZqi4islEujg54d9JAeLk5Y/G2bBRdrsAb4/rCycGabsy0Hpa4KykWQGStQedEAKHG1z1rBqGbivsxEFFTOBgEr47pjXmxXfHZnnzM+TgDVysqtS7LKonem1SRkZEqPT1d6zKISEc+Sc3Fy18dRHhgO3z4WCQ8WztrXZLFiUiGUiqyvvd0247i1p5E1FxTBwdh4eRw/JhfivGLU3Cm9KrWJVkV3QYDu5KIqCXu6+uPFbOicLqkHOMWJiOr4LLWJVkN3QYDWwxE1FJDunhjbcJgVFRW4dFFydh76oLWJVkF3QYDWwxEZAp9Onkgce4QuLs6YfLSXdhyrEDrkjSn22AgIjKVYG83JD4ZgxBvNzyxMh1f7bPvHeF0GwzsSiIiU/J1d8XaOYN/3hHu49RcrUvSjG6DgV1JRGRqbV2rd4Qb1cMXL395EAt+yLTLWdK6DQYiInNwdXLAB1Mj8MjATnjz++N47d/2t4QGt/YkIrqJk4MBf3+0PzxaOWHZjhyUXr2O18f2haOdLKGh22DgWklEZE4Gg+CVuF7wbO2Et5MycbH8Ot6ZONAu9pLWbfxxjIGIzE1EMC+2G16J64XvDp3D4yt34/K1G1qXZXa6DQYiIkuZOTQEb43vj9TsYkxZtgsXyiq0LsmsGAxERI0wNjwAi6ZG4MiZixi/OAVnS8u1LslsGAxERI00upcfVs6MxpnScsQvSsbJojKtSzIL3QYDJ7gRkRZiurTH6tmDUHbtBuIXpeDo2Ytal2Ryug0GDj4TkVb6BXhi/dwYOBoEE5ekYn9eye0/pCO6DQYiIi2F+bpj/dwYuLs6YsqyXdiVfV7rkkyGwUBE1EydvVpj/Zwh6ODhiukfpdnMyqwMBiKiFujg4Yp1CYMR5tsGs1el45sfz2hdUovpNhg4+ExE1qJ9GxesSRiM/gGeeGr1HiRm5GtdUovoNhg4+ExE1qStqxNWPR6NIV288cL6/ViVclLrkppNt8FARGRtWjs7YtljkRjdyw9/+uoQFm7J0rqkZmEwEBGZkKuTAxZOCcfDAzpi/rfHMP/bo7pbtlu3q6sSEVkrJwcD3ho/AK1dHLFwywmUXbuBV+J6w2AQrUtrFAYDEZEZGAyC1x7uAzdnByzdnoMrFZV4Y1w/OOggHBgMRERmIiL47/t7ws3FEW8nZaL8RhXeGt8fTla+4Q+DgYjIjGr2dHB1csAb3xzFteuVeG/yQLg4Wu+GP9YdWw3gPAYi0pO5I7rgLw/1xveHzyFhVQbKr1dqXdIt6TYYOI+BiPTmsSHB+Ou4vtiWWYiZy3ejzEp3g9NtMBAR6dGEqEC8PWEA0k4WY/pHabhYfl3rkn6FwUBEZGFjBnTCgkkDcSC/BFOWWt9WoQwGIiIN3NfXH4unReDYuUuYtDQVhZeuaV3SzxgMREQaGdnDD8tnRCH3/BVMWGI9+0gzGIiINDQ0zBsrZ0Wj4OI1jF+cgrziK1qXxGAgItJadIgXPnliEEquVGDC4hScLCrTtB4GAxGRFRjQ2RNrEgbj6vVKTFiSghOFlzWrhcFARGQlenf0wNqEGFRWKUxYnIrj5y5pUofVBIOIJIhIrPFxqIiEi8iLIuKpdW1ERJbSvYM71ibEwCDAxCWpOHz6osVrsJpgAJAOoCYEwpVSewAkARivXUlERJYX5tsG6+bEwMXRgMnLUnHwJ8su/WP2YBCReBHZWM9rsSKSUN9nlFKJxoexqA4HIiK7EuLthk/nxMDN2RGTlqZi76kLFvtuswdDrR/yAKpDwfh6kvF5bH2fM76eCKDY3DUSEVmjzl6tsW7OYLRr7YxpH6Yh/aRlfhxq0ZUUBSDb+DgbQLjxcSyAKBHxNIbCH4x/fhUcxvGIdBFJLywstETNRESaCGhXHQ6+7i6Y/lEaUrPPm/07tQiGmweT2wOAUmq+UuoPSqkSpVSSUmq0UmrOzS0O47FLlFKRSqlIHx8fixRNRKQVf49WWJswGB09W2HG8jTszCoy6/dpEQwlALw0+F4iIt3ybeuKtQmDEdzeDbNW7MaWYwVm+y4tgmE3fmk1hALY2MCxt8SNeojI3ni3ccGa2YMR5tsGCasysOnIObN8jyXuSooFEFlr0DkRQKjxdc+aQeim4kY9RGSP2rk5Y/UTgzEo1As+7i5m+Q5RSpnlxOYmInEA4sLCwmZnZmZqXQ4Rka6ISIZSKrK+96xpgluTsMVARGQeug0GIiIyD90GAwefiYjMQ7fBwK4kIiLz0G0wEBGReeg2GNiVRERkHroNBnYlERGZh26DgYiIzEO3E9xqiEghgFwAHgBu1a90q/e8AZh3NaqWaej/yRrO3ZEQPMAAAALrSURBVJxzNPYzjTnudsfc6v2GPmfN14Q5rwdTnF+v10ND71nz9QC07O8sSClV/yqkSimb+ANgSVPfA5Cudd3N/X+yhnM35xyN/UxjjrvdMQ38vTd0rVjtNWHO68EU59fr9dDQe9Z8PZjzmmBXEhER1cFgICKiOmwpGDY08z1rZs66TXHu5pyjsZ9pzHG3O+ZW7/N6MM/59Xo9NKUOa2OWunU/+NwSIpKubrG6INknXhNUm71eD7bUYiAiIhOw92BYonUBZHV4TVBtdnk92HVXEhER/Zq9txjqJSKhIhIuIi+KiOftP0G2TkQSjNvRkp0SEU8RiTf+semfCwyG+oUrpfYASAIwXutiyCqkA7DpHwZ0Wy+p6j3rkwAkaF2MOdlsMBhTfWM9r8WKSIN/qca/fACIRfVFQDrXkuuBbF8jr4/QWm93sVx1lmezwVDrhzuA6r9k4+tJxuexNa/f9Mez1vuJAIotWzmZQ0uvB7Jtjbw+smsdcsJy1Vmeo9YFWFAUgHXGx9kAwgEk3XxBAD9fBH8wHrcR1QFBtqXR14NRLID2IpKklCqxRIGkqfquj8U1gQEbv1vJnoLh5t/82t/qQONvCexCsm2Nvh4AQCk134y1kPX51fWhlMpG3VaDzbLZrqR6lADw0roIshq8Hqghdn192FMw7MYvvwWEorqLiOwXrwdqiF1fHzYbDMZxgshag0iJAEKNr3vWDCqRfeD1QA3h9VEXZz4TEVEdNttiICKi5mEwEBFRHQwGIiKqg8FARER1MBiIiKgOBgMREdXBYCAiojrsaa0kIosxLtWcjurF97IbWJyPyOqwxUBkYiISjupFGJeiehXOCdpWRNQ0nPlMZGIiEorqZZq9lFI2vTwz2Sa2GIhMzLg882hw6XbSKQYDkXlEGgOCSHcYDETmka51AUTNxTEGIiKqgy0GIiKqg8FARER1MBiIiKgOBgMREdXBYCAiojoYDEREVAeDgYiI6mAwEBFRHf8Hi7DOKxRHnWMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# nfw fits only require two parameters logarithm halo mass log10(M200m) and concentration parameter $c$\n",
    "# Here we are using log10(M200m) = 14, c = 5. \n",
    "hp = halo(14,5)  # halo class instance\n",
    "rbin = np.logspace(-2,np.log10(2),30)\n",
    "plt.plot(rbin, hp.nfw(rbin), '-')\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.xlabel(r'$r$')\n",
    "plt.ylabel(r'$\\rho(r)$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "Exercise: \n",
    "    \n",
    "- Make a similar plot as above for projected matter density $\\Sigma(R)$ using **sigma_nfw** and excess surface mass density $\\Delta \\Sigma (R)$ using **esd** function in the halo class. You can use the same radial binning as in the above code.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we want to compare our halo mass estimates with that of the study done on same clusters but with SDSS shape catalog following - [arXiv:1707.01907](https://arxiv.org/abs/1707.01907). So, similar to this work we are also looking cluster masses in units of $10^{14} {\\rm h^{-1} Mpc}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(x, log_M200m, c):\n",
    "    log_M200m = np.log10(log_M200m) + 14 # we are looking in units of 10^14\n",
    "    hp = halo(log_M200m, c);\n",
    "    esd = hp.esd(x)/1e?? # remember we are working with in Mpc but the measurements are in pc, we need a conversion\n",
    "    return esd\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "Exercise: \n",
    "    \n",
    "- print the output of : print(model(np.linspace(0.2,2.0,5), 13, 6))\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will use the above model function to do the fitting to our measurements. At present we are using [curve_fit](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) function from scipy to get the best fit model parameters. More advanced reader can use sophasticated tools like [emcee](https://emcee.readthedocs.io/) to do a MCMC on this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt('/home/idies/workspace/Storage/divyar/IAGRG_2022/DataStore/iagrg_dsigma.dat')\n",
    "x    = data[:,0]\n",
    "y    = data[:,1]\n",
    "yerr = data[:,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(dpi=150)\n",
    "# curve_fit requires initial guess (p0) to get resonable fits, sigma takes in the y errobars \n",
    "popt, pcov = curve_fit(model, x, y, p0=[12,0.5], sigma=yerr)\n",
    "\n",
    "log_M200m, c = popt\n",
    "sig_logMh, sig_c = np.sqrt(np.diag(pcov)) # pcov is the covariance between the parameters\n",
    "\n",
    "# rough estimate of chisq\n",
    "chisq = np.sum((y-model(x, log_M200m, c))**2*1.0/yerr**2)\n",
    "dof = 10 - 2 # 10 datapoints - 2 parameters\n",
    "\n",
    "plt.errorbar(x, y, yerr=yerr, fmt='.', capsize=3, label='Data')\n",
    "plt.errorbar(x, model(x, log_M200m, c) , fmt='-', label=r'Fit, $\\chi^2 /{\\rm dof} = %2.2f/%d$'%(chisq,dof))\n",
    "\n",
    "plt.title(r'$M_{\\rm 200m} = %2.2f \\pm %2.2f \\times 10^{14} \\,h^{-1} M_\\odot,\\, c = %2.2f \\pm %2.2f$'%(log_M200m, sig_logMh, c, sig_c))\n",
    "plt.legend()\n",
    "plt.xlabel(r'$R[{\\rm h^{-1}Mpc}]$')\n",
    "plt.ylabel(r'$\\Delta\\Sigma [{\\rm h M_\\odot pc^{-2}}]$')\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "Exercise: \n",
    "    \n",
    "- Compare your halo mass constraints with that of the weak lensing study by using SDSS shape catalog - [arXiv:1707.01907](https://arxiv.org/abs/1707.01907).\n",
    "- Remember you are studying cluster with redshifts $0.1<z_{\\rm red}<0.33$ and richness $55<\\lambda<100$. So compare the results in the subpanel of fig 4 or fig 5 in [arXiv:1707.01907](https://arxiv.org/abs/1707.01907).    \n",
    "\n",
    " \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given all the machinery to the reader to get the weak lensing signal and model them to constraint halo masses. We encourage reader to try out different cuts on their cluster sample and check the variation in halo masses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**We want to stress on the point that in the current tutorial we did a very simplistic job both in measurements as well as modelling. Please follow the references given in the tutorials as well as lectures to do a more rigorous analysis.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}