F03-Redshift distribution.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "import sys\n",
    "sys.path.insert(0, '/home/yymao/Dropbox/Academia/Collaborations/SAGA-local/saga')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import SAGA\n",
    "from SAGA import ObjectCuts as C\n",
    "from SAGA import HostCuts as H\n",
    "from easyquery import *\n",
    "from SAGA.utils import *\n",
    "from SAGA.utils.display import show_images\n",
    "from astropy.table import vstack, Table, join, unique\n",
    "from astropy.coordinates import Distance, SkyCoord\n",
    "from SAGA.objects.build2 import add_surface_brightness\n",
    "\n",
    "import numexpr as ne\n",
    "from scipy.stats import binned_statistic\n",
    "#from bottleneck import nanmax, nanmin\n",
    "from itertools import product, chain\n",
    "\n",
    "print(SAGA.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saga = SAGA.QuickStart()\n",
    "saga.set_default_base_version(\"paper2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base = SAGA.database.FitsTable(\"/home/yymao/Documents/Research/SAGA/PaperII/data-archive/saga_base_all.fits\").read()\n",
    "base = saga.host_catalog.construct_host_query(\"paper2_complete\").filter(base)\n",
    "\n",
    "assert len(np.unique(base[\"HOSTID\"])) == 36"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "specs_full = Query(C.basic_cut2, C.has_spec, \"RHOST_KPC > 10\", \"SPEC_Z >= 0.002\").filter(base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert not Query(\"SPEC_Z < 0.015\", ~C.paper2_targeting_cut).mask(specs_full).any()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hosts = np.unique(specs_full[\"HOSTID\"])\n",
    "host_dist = saga.host_catalog.load(hosts)[\"DIST\"]\n",
    "\n",
    "full_sq_deg = (np.rad2deg(np.arcsin(0.3 / host_dist))**2 - np.rad2deg(np.arcsin(0.01 / host_dist))**2).sum() * np.pi\n",
    "print(full_sq_deg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "C.paper2_targeting_cut.count(specs_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=2, figsize=(10,3.6), gridspec_kw={\"width_ratios\": (3,2), \"wspace\": 0.27})\n",
    "\n",
    "ax_inset = fig.add_axes([0.34, 0.27, 0.18, 0.34])\n",
    "\n",
    "for ax_this, zbins in zip([ax[0], ax_inset], [np.linspace(0, 1, 41), np.linspace(0, 0.05, 11)]):\n",
    "    zbins[0] = 0.001\n",
    "    for q, kw in (\n",
    "        (Query(), dict(label='Redshifts not in any targeting cuts', color=plt.cm.Oranges(0.25))),\n",
    "        (C.gri_cut, dict(label='Paper I $gri$ cuts', color=plt.cm.Blues(0.35))),\n",
    "        (C.paper2_targeting_cut, dict(label='Primary targets (this work)', color=plt.cm.Greens(0.65))),\n",
    "    ):\n",
    "        ax_this.hist(q.filter(specs_full, 'SPEC_Z'), zbins, weights=np.ones(q.count(specs_full))/full_sq_deg, **kw)\n",
    "\n",
    "        ax_this.set_xlim(0, 1 if zbins[-1] == 1 else 0.052)\n",
    "\n",
    "ax_inset.text(0.003, 5.8, \"zoom-in\")\n",
    "ax[0].set_ylabel(r'$N_{\\rm spec}$ (per sq. deg., per $z$ bin)')\n",
    "ax[0].set_xlabel(r'$z_{\\rm spec}$')\n",
    "ax[0].set_title(\"$N(z)$ split by targeting regions\", fontsize=14)\n",
    "\n",
    "ax_inset.axvline(0.003, ls=':', lw=1, color=\"grey\")\n",
    "ax_inset.axvline(0.013, ls=':', lw=1, color=\"grey\")\n",
    "ax_this = ax[-1]\n",
    "zbins = np.linspace(0, 0.8, 33)\n",
    "zbins[0] = 0.001\n",
    "faint = Query(\"r_mag >= 19\")\n",
    "plt.rcParams['hatch.linewidth'] = 0.25\n",
    "\n",
    "for q, kw in (\n",
    "    (Query(), dict(color=plt.cm.YlGn(0.27), hatch=\"//\", alpha=0.999)),\n",
    "    (Query(C.has_our_specs_only | faint) , dict(label='From literature', color=plt.cm.YlGn(0.2))),\n",
    "    (Query(C.has_our_specs_only) , dict(color=plt.cm.Purples(0.57), hatch=\"//\", alpha=0.999)),\n",
    "    (Query(C.has_our_specs_only, faint) , dict(label='From SAGA', color=plt.cm.Purples(0.5))),\n",
    "):\n",
    "    ax_this.hist(q.filter(specs_full, 'SPEC_Z'), zbins, weights=np.ones(q.count(specs_full))/full_sq_deg, **kw)\n",
    "\n",
    "ax_this.set_xlabel(r'$z_{\\rm spec}$')\n",
    "ax_this.set_ylabel(r'$N_{\\rm spec}$ (per sq. deg., per $z$ bin)')\n",
    "ax_this.set_xlim(0, 0.8)\n",
    "        \n",
    "ax_this.set_title(\"$N(z)$ split by source and magnitude\", fontsize=14)\n",
    "ax_this.text(0.75, 50, \"Hatched: $r_o < 19$\", ha=\"right\")\n",
    "\n",
    "ax_this.legend(loc='upper right', markerfirst=False, handletextpad=0.2, fontsize=13)\n",
    "\n",
    "ax[0].legend(loc='upper right', markerfirst=False, handletextpad=0.2, fontsize=13)\n",
    "plt.savefig('/home/yymao/Downloads/z_distribution.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=2, figsize=(9,4))\n",
    "\n",
    "for ax_this, zbins in zip(ax, [np.linspace(0, 0.8, 41), np.linspace(0, 0.05, 11)]):\n",
    "    zbins[0] = 0.001\n",
    "    for q, kw in (\n",
    "        (Query(), dict(label='All redshifts', color=plt.cm.Oranges(0.25))),\n",
    "        (C.gri_cut, dict(label='Paper I $gri$ cuts', color=plt.cm.Blues(0.35))),\n",
    "        (C.paper2_targeting_cut, dict(label='This work', color=plt.cm.Greens(0.65))),\n",
    "    ):\n",
    "        ax_this.hist(q.filter(specs_full, 'SPEC_Z'), zbins, weights=np.ones(q.count(specs_full))/full_sq_deg, **kw)\n",
    "        ax_this.set_xlim(0, None)\n",
    "\n",
    "ax[1].set_title(\"$N(z)$, zoom-in\", fontsize=14)\n",
    "ax[0].set_ylabel(r'$N_{\\rm spec}$ (per sq. deg., per $z$ bin)')\n",
    "ax[0].set_xlabel(r'$z_{\\rm spec}$')\n",
    "ax[1].set_xlabel(r'$z_{\\rm spec}$')\n",
    "ax[0].set_title(\"$N(z)$ split by targeting regions\", fontsize=14)\n",
    "\n",
    "\n",
    "ax[1].axvline(0.003, ls=':', lw=1, color=\"grey\")\n",
    "ax[1].axvline(0.013, ls=':', lw=1, color=\"grey\")\n",
    "ax[0].legend(loc='upper right', markerfirst=False, handletextpad=0.2, fontsize=13)\n",
    "plt.tight_layout()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## appendix plots "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base1 = saga.host_catalog.construct_host_query(\"COVERAGE_SDSS >= 0.99\").filter(base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=3, figsize=(10, 4))\n",
    "for i, ax_this in enumerate(ax):\n",
    "    \n",
    "    if i == 0:\n",
    "        key = \"r_mag\"\n",
    "        ax_this.set_xlabel(r'$r_o$')\n",
    "        bins = np.linspace(11, 21, 21)\n",
    "        bins[-1] = 20.75\n",
    "    elif i == 1:\n",
    "        key = \"sb_r\"\n",
    "        ax_this.set_xlabel(r'$\\mu_{r_o,{\\rm eff}}$')\n",
    "        bins = np.linspace(18, 26, 17)\n",
    "        cut_other = Query(C.faint_end_limit, C.paper2_targeting_cut)\n",
    "    elif i == 2:\n",
    "        key = \"gr\"\n",
    "        ax_this.set_xlabel(r'$(g-r)_o$')\n",
    "        bins = np.linspace(-0.25, 1.25, 19)\n",
    "        cut_other = Query(C.faint_end_limit, C.paper2_targeting_cut)\n",
    "        \n",
    "    bin_index = np.digitize(base1[key], bins)\n",
    "    cut = Query(C.faint_end_limit, C.paper2_targeting_cut)\n",
    "    \n",
    "    x = bins[:-1]\n",
    "    width = np.ediff1d(bins)\n",
    "    def add_bar(*queries, **kwargs):\n",
    "        y = np.maximum(np.bincount(bin_index[Query(*queries).mask(base1)], minlength=len(bins)+1)[1:-1], 0.01) \n",
    "        return ax_this.bar(x, y, width=width, align='edge', **kwargs)\n",
    "\n",
    "    c_value = 0.35\n",
    "    l0 = add_bar(color=plt.cm.Greys(c_value))\n",
    "    __ = add_bar(C.has_spec, color=plt.cm.Blues(c_value))\n",
    "    __ = add_bar(C.has_spec, C.has_our_specs_only, color=plt.cm.Greens(c_value))\n",
    "    \n",
    "    if False:\n",
    "        c_value = 0.5\n",
    "        __ = add_bar(cut_other, color=plt.cm.Greys(c_value))\n",
    "        __ = add_bar(cut_other, C.has_spec, color=plt.cm.Blues(c_value))\n",
    "        __ = add_bar(cut_other, C.has_spec, C.has_our_specs_only, color=plt.cm.Greens(c_value))\n",
    "    \n",
    "    c_value = 0.65\n",
    "    l1 = add_bar(cut, color=plt.cm.Greys(c_value))\n",
    "    l2 = add_bar(cut, C.has_spec, color=plt.cm.Blues(c_value))\n",
    "    l3 = add_bar(cut, C.has_spec, C.has_our_specs_only, color=plt.cm.Greens(c_value))\n",
    "    l4 = add_bar(cut, C.has_spec, C.is_sat, color=plt.cm.Oranges(c_value))\n",
    "    \n",
    "    ax_this.set_yscale(\"symlog\", linthresh=10, subs=[2,3,4,5,6,7,8,9], linscale=np.log10(2))\n",
    "    ax_this.set_ylim(0, 3e4)\n",
    "    ax_this.grid(\"on\", axis=\"y\", color=\"k\", lw=0.5)\n",
    "    ax_this.set_xlim(bins[0], bins[-1])\n",
    "    if i == 0:\n",
    "        bins[-1] = 21\n",
    "    #ax_this.set_xticks(bins[::5] if i == 2 else bins[::4])\n",
    "    \n",
    "\n",
    "ax[0].set_ylabel(r'$N_{\\rm galaxy}$ (per bin)')\n",
    "\n",
    "x0, y0 = 12, 2e3\n",
    "for mult in (1.1, 2, 5):\n",
    "    ax[0].plot([x0, x0], [y0, y0*mult], c='k', lw=2, marker='')\n",
    "    ax[0].text(x0, y0*0.88, fr'$\\times\\,{mult}$', ha=\"center\", va=\"top\")\n",
    "    x0 += 1.8\n",
    "\n",
    "plt.tight_layout(rect=(0,0,1,0.9), w_pad=0.1)\n",
    "\n",
    "plt.figlegend((l1, l0), (\"Within primary targeting region\", \"Outside of primary targeting region\"), loc=\"upper right\", markerfirst=False, bbox_to_anchor=(0.985,1.02))\n",
    "plt.figlegend((l1, l2, l3, l4), (\"No redshift\", \"Pre-existing redshifts\", \"SAGA redshifts\", \"Satellites\"), loc=\"upper left\", ncol=2, bbox_to_anchor=(0.08,1.02))\n",
    "\n",
    "plt.savefig('/home/yymao/Downloads/target_spec_counts.pdf', bbox_inches='tight')"
   ]
  },
  {
   "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}