{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gravitational wave events sampling (Analytical formulation)\n",
    "\n",
    "## Initial setup\n",
    "\n",
    "Note: I will interchangeably use terms like observable events and detectable events.\n",
    "\n",
    "Define all the parameters involved.\n",
    "\n",
    "- Source parameters: $\\theta \\in \\{$ $m_1$ (mass of the heavier one), $m_2$ (mass of the lighter one), $\\iota$ (inclination-angle), $\\phi$ (phase-of-coalescence), $\\psi$ (polarization-angle), $ra$ (right-ascension), $dec$ (declination) $\\}$ and $z_s$ : red-shift of the source."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given $d N^U_{obs}(z_s)$ is the number of detectable gravitational wave (GW) events from sources at red-shift $z_s$ in a spherical shell of thickness $d z_s$, then, the rate of observing GWs is given by,\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\mathcal{R}_U &= \\int_{z_{min}}^{z_{max}} \\frac{d N^U_{obs}(z_s)}{d t} \\\\\n",
    "&= \\int_{z_{min}}^{z_{max}} \\frac{d^2 N^U_{obs}(z_s)}{d t \\;d\\mathcal{V}_c}  \\frac{d\\mathcal{V}_c}{d\\Omega dz_s} d\\Omega dz_s \\\\\n",
    "&= \\int_{z_{min}}^{z_{max}} \\frac{d^2 N^U_{obs}(z_s)}{\\left(d t/(1+z_s)\\right) \\;d\\mathcal{V}_c}\\frac{1}{(1+z_s)}\\frac{d\\mathcal{V}_c}{d\\Omega dz_s} d\\Omega dz_s \\\\\n",
    "&= \\int_{z_{min}}^{z_{max}} \\frac{d^2 N^U_{obs}(z_s)}{d\\tau \\;d\\mathcal{V}_c}\\frac{1}{(1+z_s)} \\frac{d\\mathcal{V}_c}{d\\Omega dz_s} d\\Omega dz_s \\\\\n",
    "&= \\int_{z_{min}}^{z_{max}} \\frac{d^2 N^U_{obs}(z_s)}{d\\tau \\;d\\mathcal{V}_c}\\frac{1}{(1+z_s)} \\frac{d\\mathcal{V}_c}{d\\Omega dz_s} \\Bigg\\{ \\int_{all\\; sky} d\\Omega \\Bigg\\}dz_s \\\\\n",
    "&= \\int_{z_{min}}^{z_{max}} \\frac{d^2 N^U_{obs}(z_s)}{d\\tau \\;d\\mathcal{V}_c}\\frac{1}{(1+z_s)} \\Bigg(\\frac{d\\mathcal{V}_c}{d\\Omega dz_s} 4\\pi \\Bigg)dz_s \\\\\n",
    "&= \\int_{z_{min}}^{z_{max}} \\frac{d^2 N^U_{obs}(z_s)}{d\\tau \\;d\\mathcal{V}_c}\\frac{1}{(1+z_s)} \\frac{dV_c}{dz_s} dz_s \\\\\n",
    "\\end{split} \\tag{1}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Note that if $d\\tau=dt/(1+z_s)$ is considered the proper time and it can be converted from the time at detector frame $dt$ using the time-dilation factor $(1+z_s)$. Consequently, $\\frac{d^2 N^U_{obs}(z_s)}{d t \\;d\\mathcal{V}_c}$ and $\\frac{d^2 N^U_{obs}(z_s)}{d \\tau \\;d\\mathcal{V}_c}$ are the observed merger rate density at detector-frame and source-frame respectively. We want to use the $R^U_{obs}(z_s)=\\frac{d^2 N^U_{obs}(z_s)}{d \\tau \\;d\\mathcal{V}_c}$ for our analysis as most observational papers and the output of theoretical predictions are in the source-frame. $\\frac{dV_c}{dz_s}dz_s$ is considered a spherical-shell volume element in co-moving coordinates at red-shift $z_s$. So, the rate expression simplifies to integrating (density) $\\times$ (time dilation effect) over the shell volume element. For more information on the volume element refer to this [page]() of the documentation.\n",
    "\n",
    "Note: $\\frac{dV_c}{dz_s}$ is the differential co-moving volume at red-shift $z_s$ and you can get the value by using `astropy` cosmology package for a given cosmology.\n",
    "\n",
    "Now we have,\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\mathcal{R}_U &= \\int_{z_{min}}^{z_{max}} \\frac{d^2 N^U_{obs}(z_s)}{d \\tau\\;d\\mathcal{V}_c} \\frac{1}{(1+z_s)} \\frac{dV_c}{dz_s} dz_s \\\\\n",
    "&= \\int_{z_{min}}^{z_{max}} R^U_{obs}(z_s) \\frac{1}{(1+z_s)} \\frac{dV_c}{dz_s} dz_s\n",
    "\\end{split} \\tag{2}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Equation 2 is exactly same to the of equation A2 from [WIERDA et al. 2021](https://arxiv.org/pdf/2106.06303.pdf).\n",
    "\n",
    "We want to re-write it in terms intrinsic merger rate distribution $R(z_s)$. $R(z_s)$ is the merger rate density distribution regardless of whether the source is detectable or not.\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\mathcal{R}_U &= \\int_{z_{min}}^{z_{max}} R^U_{obs}(z_s) \\frac{1}{(1+z_s)} \\frac{dV_c}{dz_s} dz_s \\\\\n",
    "&= \\int_{z_{min}}^{z_{max}} R(z_s) P(obs|z_s) \\frac{1}{(1+z_s)} \\frac{dV_c}{dz_s} dz_s \\\\\n",
    "\\end{split} \\tag{3}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\mathcal{R}_U &= \\mathcal{N}^U \\int_{z_{min}}^{z_{max}} P(z_s) P(obs|z_s) dz_s \\\\\n",
    "\\end{split} \\tag{4}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "$P(obs|z_s)$ is the probability of observing a GW event given that the source is at red-shift $z_s$. Here, it is considered marginalized over all the gravitational wave parameters except the red-shift, i.e., $P(obs|z_s) = \\int_{\\theta} d\\theta P(obs|z_s, \\theta) P(\\theta)$.\n",
    "\n",
    "In Eqn(4) the merger rate density becomes the prior distribution $P(z_s)$ for the red-shift, and it is normalized over the red-shift range. So the normalizing factor is given by,\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\mathcal{N}^U = \\int_{z_{min}}^{z_{max}} R(z_s) \\frac{1}{(1+z_s)} \\frac{dV_c}{dz_s} dz_s \n",
    "\\end{split} \\tag{5}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "Therefore,\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "P(z_s) = \\frac{R(z_s)}{\\mathcal{N}^U} \\frac{1}{(1+z_s)} \\frac{dV_c}{dz_s}\n",
    "\\end{split} \\tag{6}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "This is done because we want to sample the red-shift values from this prior distribution, which later will be used in the monte-carlo integration to estimate the rate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rate equation (statistical formulation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now consider all the gravitational waves source parameters $\\theta \\in \n",
    "\\{m_1,m_2,z_s,\\iota,\\phi,\\psi,ra,dec,\\text{spins}\\}$,\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\mathcal{R}_U &= \\mathcal{N}^U \\int_{z_{min}}^{z_{max}} P(z_s) P(obs|z_s) dz_s \\\\\n",
    "& \\text{consider} \\int \\rightarrow \\int_{z_s}\\int_{\\theta} \\rightarrow \\int_{z_s}\\int_{m_1}\\int_{m_2}\\int_{\\iota}\\int_{\\phi}\\int_{\\psi}\\int_{ra}\\int_{dec}\\int_{\\text{spins}} \\\\\n",
    "\\mathcal{R}_U &= \\mathcal{N}^U \\int P(z_s) P(obs|z_s, \\theta) P(\\theta) d\\theta dz_s\n",
    "\\end{split} \\tag{8}\n",
    "\\end{equation}\n",
    "\n",
    "Here, I have assumed $P(\\theta)$ is a normalized multidimensional pdf. All parameters in $\\theta$, other than $m_1$ and $m_2$, are sampled independently. The pdf $P(\\theta)$ is a product of individual pdfs for each parameter.\n",
    "\n",
    "Final expression for the rate of observing GW events is given by,\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\mathcal{R}_U &= \\mathcal{N}^U \\bigg< P(obs|\\theta, z_s) \\bigg>_{z_s\\in P(z_s), \\theta\\in P(\\theta)} \\\\\n",
    "P(\\theta) &= P(m_1,m_2)P(\\iota)P(\\phi)P(\\psi)P(ra)P(dec)P(\\text{spins})\n",
    "\\end{split} \\tag{4}\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The order of sampling and rate calculation steps in LeR are listed below.\n",
    "\n",
    "1. Sample $z_s$ from $P(z_s)$.\n",
    "2. Sample $\\theta$ from $P(\\theta)$.\n",
    "3.  Calculate SNR with [gwsnr](https://gwsnr.readthedocs.io/en/latest/)\n",
    "4.  Apply the SNR threshold and check whether the event is detectable or not.\n",
    "5.  Calculate rate of GW events, $\\mathcal{R}_U$ using equation 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "plaintext"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}