{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1a342e58",
   "metadata": {},
   "source": [
    "# Astrophysics III: galaxy formation & evolution\n",
    "Moodle: https://go.epfl.ch/PHYS-465\n",
    "\n",
    "## Exercise 5\n",
    "Teaching assistant: Jonathan Petersson <br>\n",
    "Email: jonathan.petersson@epfl.ch"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55040aa8",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20ae81fb",
   "metadata": {},
   "source": [
    "## Load data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3eaf8a3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import binom\n",
    "\n",
    "%matplotlib notebook\n",
    "plt.rcParams['font.size'] = 14\n",
    "\n",
    "\n",
    "# Load data:\n",
    "data  = np.loadtxt('BrinchDave_new.dat')\n",
    "\n",
    "# Only selet galaxies with Mr < -18\n",
    "Mag_r = data[:,20]\n",
    "data  = data[Mag_r < -18]\n",
    "\n",
    "# Assign some data to arrays:\n",
    "ssfr  = data[:,0]       # [log(1/yr)]\n",
    "mass  = data[:,3]       # [log(Msol)]\n",
    "dens_05Mpc = data[:,6]  # [number counts within a projected cylinder]\n",
    "dens_1Mpc  = data[:,7]\n",
    "dens_2Mpc  = data[:,8]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "367d309e",
   "metadata": {},
   "source": [
    "## Exercise 5.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59d96fff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bin your data:\n",
    "bins = np.linspace(-13, -8, 50)\n",
    "h1, bin_edges1 = np.histogram(ssfr[(mass > 9.50) * (mass < 10.0)], bins=bins)\n",
    "h2, bin_edges2 = np.histogram(ssfr[(mass > 10.0) * (mass < 10.5)], bins=bins)\n",
    "h3, bin_edges3 = np.histogram(ssfr[(mass > 10.5) * (mass < 11.0)], bins=bins)\n",
    "h4, bin_edges4 = np.histogram(ssfr[(mass > 11.0) * (mass < 12.0)], bins=bins)\n",
    "\n",
    "# Plot:\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "ax.hist(bin_edges1[:-1], bins=bin_edges1, weights=h1/np.sum(h1), alpha=0.5, label='9.50-10.0')\n",
    "ax.hist(bin_edges2[:-1], bins=bin_edges2, weights=h2/np.sum(h2), alpha=0.5, label='10.0-10.5')\n",
    "ax.hist(bin_edges3[:-1], bins=bin_edges3, weights=h3/np.sum(h3), alpha=0.5, label='10.5-11.0')\n",
    "ax.hist(bin_edges4[:-1], bins=bin_edges4, weights=h4/np.sum(h4), alpha=0.5, label='11.0-12.0')\n",
    "\n",
    "ax.axvline(-11, c='k', ls='--', lw=2, alpha=0.5)\n",
    "\n",
    "ax.set_xlabel('$\\log_{10}(\\mathrm{sSFR} \\ [\\mathrm{yr}^{-1}])$')\n",
    "ax.set_ylabel('$N/N_\\mathrm{Total}$')\n",
    "ax.legend(frameon=False, title='$\\log_{10}(M_\\star \\ [\\mathrm{M}_\\odot])$')\n",
    "fig.savefig('ex1.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52dde892",
   "metadata": {},
   "source": [
    "## Exercise 5.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "367b04b3",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Critical sSFR to distinguish star-forming from quiescent galaxies:\n",
    "ssfr_crit = -11\n",
    "\n",
    "# Bin your data:\n",
    "bins=np.array([9.5, 10.0, 10.5, 11.0, 12.0])\n",
    "hist_ngal, bin_edges = np.histogram(mass, bins=bins)\n",
    "\n",
    "# Fraction of quiescent galaxies:\n",
    "hist_ngal_q, bin_edges = np.histogram(mass[ssfr < -11], bins=bin_edges)\n",
    "frac_q = hist_ngal_q / hist_ngal\n",
    "\n",
    "# Plot:\n",
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "ax.plot((bins[:-1] + bins[1:])/2, frac_q, lw=2)\n",
    "ax.scatter((bins[:-1] + bins[1:]) / 2, frac_q, s=50, marker='o')\n",
    "ax.set_xlabel('$\\log_{10}(M_\\star \\ [\\mathrm{M}_\\odot])$')\n",
    "ax.set_ylabel('$f_\\mathrm{quiescent}$')\n",
    "fig.tight_layout()\n",
    "fig.savefig('ex3.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64b9ba5f",
   "metadata": {},
   "source": [
    "## Exercise 5.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b92c145e",
   "metadata": {},
   "outputs": [],
   "source": [
    "dens_05Mpc = np.nan_to_num(dens_05Mpc)\n",
    "dens_1Mpc = np.nan_to_num(dens_1Mpc)\n",
    "dens_2Mpc = np.nan_to_num(dens_2Mpc)\n",
    "\n",
    "# Mask for quiescent galaxies\n",
    "ssfr_crit = -11\n",
    "mask_q = ssfr < -11\n",
    "\n",
    "# Mass masks (bins):\n",
    "mask_mass1 = (mass > 9.50) * (mass < 10.0)\n",
    "mask_mass2 = (mass > 10.0) * (mass < 10.5)\n",
    "mask_mass3 = (mass > 10.5) * (mass < 11.0)\n",
    "mask_mass4 = (mass > 11.0) * (mass < 12.0)\n",
    "\n",
    "# Density bins:\n",
    "bins_list = [np.linspace(0, 50, 5), np.linspace(0, 20, 5), np.linspace(0, 10, 5)]\n",
    "\n",
    "# Plot:\n",
    "fig, ax = plt.subplots(1, 3, figsize=(12, 3), sharey=True)\n",
    "\n",
    "dens_list = [dens_05Mpc, dens_1Mpc, dens_2Mpc]\n",
    "for i in range(0, 3):\n",
    "    dens = dens_list[i]\n",
    "    bins = bins_list[i]\n",
    "    \n",
    "    hist_ngal, bin_edges = np.histogram(dens[mask_mass1], bins=bins)\n",
    "    hist_ngal_q, bin_edges = np.histogram(dens[mask_mass1 * mask_q], bins=bins)\n",
    "    frac_q = hist_ngal_q / hist_ngal\n",
    "    ax[i].plot((bin_edges[:-1] + bin_edges[1:]) / 2, frac_q, lw=2, marker='o', alpha=1, label='$9.5-10$')\n",
    "    \n",
    "    hist_ngal, bin_edges = np.histogram(dens[mask_mass2], bins=bins)\n",
    "    hist_ngal_q, bin_edges = np.histogram(dens[mask_mass2 * mask_q], bins=bins)\n",
    "    frac_q = hist_ngal_q / hist_ngal\n",
    "    ax[i].plot((bin_edges[:-1] + bin_edges[1:]) / 2, frac_q, lw=2, marker='o', alpha=1, label='10-10.5')\n",
    "\n",
    "    hist_ngal, bin_edges = np.histogram(dens[mask_mass3], bins=bins)\n",
    "    hist_ngal_q, bin_edges = np.histogram(dens[mask_mass3 * mask_q], bins=bins)\n",
    "    frac_q = hist_ngal_q / hist_ngal\n",
    "    ax[i].plot((bin_edges[:-1] + bin_edges[1:]) / 2, frac_q, lw=2, marker='o', alpha=1, label='10.5-11')\n",
    "\n",
    "    hist_ngal, bin_edges = np.histogram(dens[mask_mass4], bins=bins)\n",
    "    hist_ngal_q, bin_edges = np.histogram(dens[mask_mass4 * mask_q], bins=bins)\n",
    "    frac_q = hist_ngal_q / hist_ngal\n",
    "    ax[i].plot((bin_edges[:-1] + bin_edges[1:]) / 2, frac_q, lw=2, marker='o', alpha=1, label='11-12')\n",
    "\n",
    "ax[0].set_ylabel('$f_\\mathrm{quiescent}$')\n",
    "ax[0].set_xlabel('$N_\\mathrm{galaxies}$ within 0.5 Mpc')\n",
    "ax[1].set_xlabel('$N_\\mathrm{galaxies}$ within 1.0 Mpc')\n",
    "ax[2].set_xlabel('$N_\\mathrm{galaxies}$ within 2.0 Mpc')\n",
    "ax[2].legend(title='$\\log_{10}(M_\\star$ [M$_\\odot$])', loc='upper left', bbox_to_anchor=(1, 1), frameon=False)\n",
    "fig.tight_layout()\n",
    "fig.savefig('ex4.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f883a6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "dens_05Mpc = np.nan_to_num(dens_05Mpc)\n",
    "dens_1Mpc = np.nan_to_num(dens_1Mpc)\n",
    "dens_2Mpc = np.nan_to_num(dens_2Mpc)\n",
    "\n",
    "# Mask for quiescent galaxies\n",
    "ssfr_crit = -11\n",
    "mask_q = ssfr < -11\n",
    "\n",
    "# Mass masks (bins):\n",
    "mask_mass1 = (mass > 9.50) * (mass < 10.0)\n",
    "mask_mass2 = (mass > 10.0) * (mass < 10.5)\n",
    "mask_mass3 = (mass > 10.5) * (mass < 11.0)\n",
    "mask_mass4 = (mass > 11.0) * (mass < 12.0)\n",
    "\n",
    "# Density bins:\n",
    "bins_list = [np.linspace(0, 50, 5), np.linspace(0, 20, 5), np.linspace(0, 10, 5)]\n",
    "\n",
    "# Plot:\n",
    "fig, ax = plt.subplots(1, 3, figsize=(12, 3), sharey=True)\n",
    "\n",
    "dens_list = [dens_05Mpc, dens_1Mpc, dens_2Mpc]\n",
    "for i in range(0, 3):\n",
    "    dens = dens_list[i]\n",
    "    bins = bins_list[i]\n",
    "    \n",
    "    hist_ngal, bin_edges = np.histogram(dens[mask_mass1], bins=bins)\n",
    "    hist_ngal_q, bin_edges = np.histogram(dens[mask_mass1 * mask_q], bins=bins)\n",
    "    frac_q = hist_ngal_q / hist_ngal\n",
    "    mean, var, skew, kurt = binom.stats(hist_ngal, frac_q, moments='mvsk')\n",
    "    ax[i].errorbar((bin_edges[:-1] + bin_edges[1:]) / 2, frac_q, yerr=np.sqrt(var)/hist_ngal, lw=1, elinewidth=1, marker='o', ms=5, alpha=1, label='$9.5-10$')\n",
    "    \n",
    "    hist_ngal, bin_edges = np.histogram(dens[mask_mass2], bins=bins)\n",
    "    hist_ngal_q, bin_edges = np.histogram(dens[mask_mass2 * mask_q], bins=bins)\n",
    "    frac_q = hist_ngal_q / hist_ngal\n",
    "    mean, var, skew, kurt = binom.stats(hist_ngal, frac_q, moments='mvsk')\n",
    "    ax[i].errorbar((bin_edges[:-1] + bin_edges[1:]) / 2, frac_q, yerr=np.sqrt(var)/hist_ngal, lw=1, elinewidth=1, marker='o', ms=5, alpha=1, label='10-10.5')\n",
    "\n",
    "    hist_ngal, bin_edges = np.histogram(dens[mask_mass3], bins=bins)\n",
    "    hist_ngal_q, bin_edges = np.histogram(dens[mask_mass3 * mask_q], bins=bins)\n",
    "    frac_q = hist_ngal_q / hist_ngal\n",
    "    mean, var, skew, kurt = binom.stats(hist_ngal, frac_q, moments='mvsk')\n",
    "    ax[i].errorbar((bin_edges[:-1] + bin_edges[1:]) / 2, frac_q, yerr=np.sqrt(var)/hist_ngal, lw=1, elinewidth=1, marker='o', ms=5, alpha=1, label='10.5-11')\n",
    "\n",
    "    hist_ngal, bin_edges = np.histogram(dens[mask_mass4], bins=bins)\n",
    "    hist_ngal_q, bin_edges = np.histogram(dens[mask_mass4 * mask_q], bins=bins)\n",
    "    frac_q = hist_ngal_q / hist_ngal\n",
    "    mean, var, skew, kurt = binom.stats(hist_ngal, frac_q, moments='mvsk')\n",
    "    ax[i].errorbar((bin_edges[:-1] + bin_edges[1:]) / 2, frac_q, yerr=np.sqrt(var)/hist_ngal, lw=1, elinewidth=1, marker='o', ms=5, alpha=1, label='11-12')\n",
    "\n",
    "ax[0].set_ylabel('$f_\\mathrm{quiescent}$')\n",
    "ax[0].set_xlabel('$N_\\mathrm{galaxies}$ within 0.5 Mpc')\n",
    "ax[1].set_xlabel('$N_\\mathrm{galaxies}$ within 1.0 Mpc')\n",
    "ax[2].set_xlabel('$N_\\mathrm{galaxies}$ within 2.0 Mpc')\n",
    "ax[2].legend(title='$\\log_{10}(M_\\star$ [M$_\\odot$])', loc='upper left', bbox_to_anchor=(1, 1), frameon=False)\n",
    "fig.tight_layout()\n",
    "fig.savefig('ex4.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b8d91ba",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}