F09-11-Sat properties and quenching.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": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from itertools import combinations\n",
    "from collections import Counter, defaultdict\n",
    "from astropy.table import Table, join, vstack\n",
    "from astropy.coordinates import SkyCoord, Distance\n",
    "from SAGA.utils import add_skycoord\n",
    "from easyquery import Query, QueryMaker\n",
    "from astropy.io import ascii\n",
    "\n",
    "import SAGA\n",
    "from SAGA.database import FitsTable, GoogleSheets\n",
    "from SAGA import ObjectCuts as C\n",
    "from SAGA.utils import add_skycoord, fill_values_by_query\n",
    "from SAGA.utils.distance import z2v, d2m, m2d\n",
    "from SAGA.utils.display import show_images\n",
    "from SAGA.objects.object_catalog import calc_fiducial_p_sat, calc_fiducial_p_sat_corrected\n",
    "\n",
    "from matplotlib.transforms import Bbox\n",
    "from scipy.stats import spearmanr\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": [
    "def annotate_base(base):\n",
    "    base[\"Mr_\"] = base[\"r_mag\"] - d2m(base[\"HOST_DIST\"])\n",
    "    base[\"Mr\"] = np.where(np.isfinite(base[\"Mr\"]) & C.has_spec.mask(base), base[\"Mr\"], base[\"Mr_\"])\n",
    "    base[\"log_sm\"] = np.where(np.isfinite(base[\"log_sm\"]), base[\"log_sm\"], 1.254 + 1.0976 * base[\"gr\"] - 0.4 * base[\"Mr\"])\n",
    "    \n",
    "    base[\"human_selected\"] = 0\n",
    "    hs = saga.database[\"human_selected\"].read()\n",
    "    for i in range(1, 4):\n",
    "        base[\"human_selected\"] += np.in1d(base[\"OBJID\"], Query(f\"score == {i}\").filter(hs, \"OBJID\")).astype(np.int) * i\n",
    "\n",
    "    return base"
   ]
  },
  {
   "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",
    "assert len(np.unique(base[\"HOSTID\"])) == 36\n",
    "\n",
    "base = annotate_base(base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsat_missing_all = int(np.ceil(Query(~C.has_spec, C.faint_end_limit).filter(base, \"p_sat_corrected\").sum()))\n",
    "nsat_missing_limit = int(np.ceil(Query(~C.has_spec, C.r_abs_limit).filter(base, \"p_sat_corrected\").sum()))\n",
    "print(nsat_missing_all, nsat_missing_limit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sats = C.is_sat.filter(base)\n",
    "sats_maybe_all = Query(~C.has_spec, \"p_sat_corrected > 0\").filter(base)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 1177 #923 #444 #177 # 114 # 47\n",
    "while True:\n",
    "    #seed += 1\n",
    "    sats_maybe = Query(\n",
    "        (lambda p: p > np.random.RandomState(seed).rand(len(p)), \"p_sat_corrected\"),  # must put in 1st place\n",
    "        ~C.has_spec, \n",
    "    ).filter(base)\n",
    "\n",
    "    if (np.abs((C.faint_end_limit.count(sats_maybe) - nsat_missing_all) < 1.5 and\n",
    "        C.r_abs_limit.count(sats_maybe)) == nsat_missing_limit and\n",
    "        (~Query(C.high_priority_sb_tight, C.gr_cut_tight, \"Mr > -15.5\")).count(sats_maybe) == 0):\n",
    "        print(seed)\n",
    "        break\n",
    "    else:\n",
    "        seed += 1\n",
    "        \n",
    "print(C.r_abs_limit.count(sats_maybe), len(sats_maybe))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert len(sats) == 127\n",
    "assert C.r_abs_limit.count(sats) == 123\n",
    "assert Query(\"EW_Halpha < 2\").count(sats) == 18\n",
    "assert Query(\"EW_Halpha < 2\", C.r_abs_limit).count(sats) == 18\n",
    "assert Query(\"EW_Halpha < 2\", \"Mr < -16\").count(sats) == 2\n",
    "assert len(np.unique(sats[\"HOSTID\"])) == 34"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "is_p1_sats = Query(C.is_sat, saga.host_catalog.construct_host_query(\"paper1_complete\"), QueryMaker.in1d(\"OBJID\", [4443380000000000475, 3115540000000003218], invert=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "\n",
    "EW_Halpha_map = mpl.colors.LinearSegmentedColormap.from_list(\n",
    "    'EW_Halpha_map',\n",
    "    np.vstack((plt.cm.Blues_r(np.linspace(0, 0.75, 256)), plt.cm.Reds(np.linspace(0.25, 1, 256))))[::-1],\n",
    ")\n",
    "\n",
    "def EW_Halpha_normalize(x):\n",
    "    x = np.asarray(x)\n",
    "    with np.errstate(all=\"ignore\"):\n",
    "        return np.where(x > 2, np.log2(x)*2, x)\n",
    "    \n",
    "EW_Halpha_norm = mpl.colors.TwoSlopeNorm(EW_Halpha_normalize(2), EW_Halpha_normalize(-2), EW_Halpha_normalize(200))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_sci(x, n=1):\n",
    "    x = float(x)\n",
    "    s = (f\"{{:.{n}E}}\").format(x)\n",
    "    a, _, b = s.partition(\"E\")\n",
    "    if abs(int(b)) <= 1:\n",
    "        return (f\"{{:.{n}f}}\").format(x)\n",
    "    b = int(b)\n",
    "    return f\"{a} \\\\times 10^{{{b}}}\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(10.5, 6.5), gridspec_kw={'hspace': 0.35, 'wspace':0.3})\n",
    "\n",
    "EW_Halpha_ticks = [-2, 0, 2, 5, 10, 20, 50, 100, 200]\n",
    "EW_Halpha_ticklabels = [f\"$\\\\leq {t}$\" if t == EW_Halpha_ticks[0] else (f\"$\\\\geq {t}$\" if t == EW_Halpha_ticks[-1] else str(t)) for t in EW_Halpha_ticks]\n",
    "\n",
    "####\n",
    "DO_PAPER_I_ONLY = False #or True\n",
    "####\n",
    "\n",
    "sats_this = is_p1_sats.filter(sats) if DO_PAPER_I_ONLY else sats\n",
    "sats_this.sort(\"EW_Halpha\", reverse=True)\n",
    "\n",
    "#####\n",
    "ax_this = ax[0,0]\n",
    "\n",
    "CS = ax_this.scatter(sats_this[\"Mr\"], sats_this[\"gr\"], c=EW_Halpha_normalize(sats_this[\"EW_Halpha\"]), cmap=EW_Halpha_map, norm=EW_Halpha_norm, alpha=0.9, rasterized=True, label=\"Confirmed\")\n",
    "if not DO_PAPER_I_ONLY:\n",
    "    ax_this.scatter(sats_maybe[\"Mr\"], sats_maybe[\"gr\"], marker='o', c='None', edgecolors='k', lw=0.5, label=\"Unconfirmed\")\n",
    "\n",
    "mr = np.linspace(-20, -12, 11)\n",
    "\n",
    "ax_this.set_xlim(-20.5, -11.5)\n",
    "ax_this.set_ylim(-0.05, 0.8)\n",
    "ax_this.set_xlabel(\"$M_{r,o}$\")\n",
    "ax_this.set_ylabel(\"$(g-r)_o$\")\n",
    "ax_this.axvspan(-12.3, -11.5, color=\"grey\", alpha=0.4)\n",
    "ax_this.legend(loc=\"lower left\", handletextpad=0.1, fontsize=12, borderpad=0.2, handlelength=1, frameon=True)\n",
    "\n",
    "s = spearmanr(sats_this[\"Mr\"], sats_this[\"gr\"])\n",
    "#ax_this.text(0.4, 0.15, \"$(\\\\rho_s, p) = ({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue)), fontsize=11, transform=ax_this.transAxes)\n",
    "print(\"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue)))\n",
    "s = spearmanr(np.hstack([sats_this[\"Mr\"], sats_maybe[\"Mr\"]]), np.hstack([sats_this[\"gr\"], sats_maybe[\"gr\"]]))\n",
    "#ax_this.text(0.4, 0.05, \"$(\\\\rho_s, p) = ({:.3f}, {})$\".format(s.correlation, print_sci(s.pvalue)), fontsize=11, transform=ax_this.transAxes)\n",
    "print(\"$({:.3f}, {})$\".format(s.correlation, print_sci(s.pvalue)))\n",
    "\n",
    "#####\n",
    "ax_this = ax[1,0]\n",
    "\n",
    "CS = ax_this.scatter(sats_this[\"Mr\"], sats_this[\"sb_r\"], c=EW_Halpha_normalize(sats_this[\"EW_Halpha\"]), cmap=EW_Halpha_map, norm=EW_Halpha_norm, alpha=0.9, rasterized=True)\n",
    "if not DO_PAPER_I_ONLY:\n",
    "    ax_this.scatter(sats_maybe[\"Mr\"], sats_maybe[\"sb_r\"], marker='o', c='None', edgecolors='k', lw=0.5)\n",
    "\n",
    "ax_this.set_xlim(-20.5, -11.5)\n",
    "ax_this.set_ylim(19, 26.5)\n",
    "ax_this.set_xlabel(\"$M_{r,o}$\")\n",
    "ax_this.set_ylabel(\"$\\\\mu_{r_o, \\\\mathrm{eff}}$ [mag arcsec$^{-2}$]\")\n",
    "ax_this.axvspan(-12.3, -11.5, color=\"grey\", alpha=0.4)\n",
    "\n",
    "s = spearmanr(sats_this[\"Mr\"], sats_this[\"sb_r\"])\n",
    "#ax_this.text(0.57, 0.15, \"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue)), fontsize=11, transform=ax_this.transAxes)\n",
    "print(\"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue)))\n",
    "s = spearmanr(np.hstack([sats_this[\"Mr\"], sats_maybe[\"Mr\"]]), np.hstack([sats_this[\"sb_r\"], sats_maybe[\"sb_r\"]]))\n",
    "#ax_this.text(0.57, 0.05, \"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue)), fontsize=11, transform=ax_this.transAxes)\n",
    "print(\"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue)))\n",
    "\n",
    "\n",
    "####\n",
    "ax_this = ax[1,1]\n",
    "\n",
    "CS = ax_this.scatter(sats_this[\"RHOST_KPC\"], sats_this[\"Mr\"], c=EW_Halpha_normalize(sats_this[\"EW_Halpha\"]), cmap=EW_Halpha_map, norm=EW_Halpha_norm, alpha=0.9, rasterized=True)\n",
    "if not DO_PAPER_I_ONLY:\n",
    "    ax_this.scatter(sats_maybe[\"RHOST_KPC\"], sats_maybe[\"Mr\"], marker='o', c='None', edgecolors='k', lw=0.5)\n",
    "\n",
    "ax_this.set_ylim(-11.5, -20.5)\n",
    "ax_this.set_xlim(-10, 310)\n",
    "ax_this.set_ylabel(\"$M_{r,o}$\")\n",
    "ax_this.set_xlabel(\"Projected distance [kpc]\")\n",
    "ax_this.axhspan(-12.3, -11.5, color=\"grey\", alpha=0.4, rasterized=True)\n",
    "\n",
    "s = spearmanr(sats_this[\"RHOST_KPC\"], sats_this[\"Mr\"])\n",
    "#ax_this.text(0.98, 0.92, \"$({:.3f}, {})$\".format(s.correlation, print_sci(s.pvalue, 2)), fontsize=11, transform=ax_this.transAxes, ha=\"right\")\n",
    "print(\"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue, 2)))\n",
    "\n",
    "s = spearmanr(np.hstack([sats_this[\"RHOST_KPC\"], sats_maybe[\"RHOST_KPC\"]]), np.hstack([sats_this[\"Mr\"], sats_maybe[\"Mr\"]]))\n",
    "#ax_this.text(0.98, 0.82, \"$({:.3f}, {})$\".format(s.correlation, print_sci(s.pvalue, 2)), fontsize=11, transform=ax_this.transAxes, ha=\"right\")\n",
    "print(\"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue, 2)))\n",
    "\n",
    "\n",
    "####\n",
    "\n",
    "####\n",
    "ax_this = ax[0,1]\n",
    "\n",
    "sats_this1 = C.r_abs_limit.filter(sats_this)\n",
    "sats_maybe1 = C.r_abs_limit.filter(sats_maybe)\n",
    "\n",
    "CS = ax_this.scatter(sats_this1[\"RHOST_KPC\"], sats_this1[\"gr\"], c=EW_Halpha_normalize(sats_this1[\"EW_Halpha\"]), cmap=EW_Halpha_map, norm=EW_Halpha_norm, alpha=0.9, rasterized=True)\n",
    "if not DO_PAPER_I_ONLY:\n",
    "    ax_this.scatter(sats_maybe1[\"RHOST_KPC\"], sats_maybe1[\"gr\"], marker='o', c='None', edgecolors='k', lw=0.5)\n",
    "\n",
    "ax_this.set_xlim(-10, 310)\n",
    "ax_this.set_ylim(-0.05, 0.8)\n",
    "ax_this.set_xlabel(\"Projected distance [kpc]\")\n",
    "ax_this.set_ylabel(\"$(g-r)_o$\")\n",
    "ax_this.axhspan(-12.3, -11.5, color=\"grey\", alpha=0.4, rasterized=True)\n",
    "ax_this.text(0, 0, \"($M_{r,o} < -12.3$ only)\", fontsize=12)\n",
    "\n",
    "s = spearmanr(sats_this[\"RHOST_KPC\"], sats_this[\"gr\"])\n",
    "#ax_this.text(0.98, 0.15, \"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue)), fontsize=11, transform=ax_this.transAxes, ha=\"right\")\n",
    "print(\"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue, 2)))\n",
    "\n",
    "s = spearmanr(np.hstack([sats_this[\"RHOST_KPC\"], sats_maybe[\"RHOST_KPC\"]]), np.hstack([sats_this[\"gr\"], sats_maybe[\"gr\"]]))\n",
    "#ax_this.text(0.98, 0.05, \"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue)), fontsize=11, transform=ax_this.transAxes, ha=\"right\")\n",
    "print(\"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue, 2)))\n",
    "\n",
    "\n",
    "####\n",
    "\n",
    "cbar = plt.colorbar(CS, ax=ax, ticks=EW_Halpha_normalize(EW_Halpha_ticks), fraction=0.02, aspect=50, pad=0.03, shrink=0.7)\n",
    "cbar.ax.set_rasterized(True)\n",
    "cbar.ax.set_yticklabels(EW_Halpha_ticklabels)\n",
    "cbar.solids.set_edgecolor(\"face\")\n",
    "cbar.ax.minorticks_off()\n",
    "cbar.ax.set_ylabel(\"EW(H$\\\\alpha$) [$\\\\mathrm{\\\\AA}$]\", fontsize=14, labelpad=-10)\n",
    "\n",
    "if not DO_PAPER_I_ONLY:\n",
    "    plt.savefig(\"/home/yymao/Downloads/halpha.pdf\", dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Query('EW_Halpha < 2', \"Mr < -16\").filter(sats, [\"OBJID\", \"RA\", \"DEC\", \"HOSTID\", \"MASKNAME\", \"r_mag\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# sb--radius-color distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mw = SAGA.database.GoogleSheets(\"1O8tGgnHXRcAT8P78J3V2pWVb3CWE1Hvcowily79tdg4\", 0).read()\n",
    "mw[\"gr\"] = mw[\"Mg_o\"] - mw[\"Mr_o\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 3.6))\n",
    "\n",
    "EW_Halpha_ticks = [-2, 0, 2, 5, 10, 20, 50, 100, 200]\n",
    "EW_Halpha_ticklabels = [f\"$\\\\leq {t}$\" if t == EW_Halpha_ticks[0] else (f\"$\\\\geq {t}$\" if t == EW_Halpha_ticks[-1] else str(t)) for t in EW_Halpha_ticks]\n",
    "\n",
    "sats.sort(\"log_sm\", reverse=True)\n",
    "\n",
    "#####\n",
    "ax_this = ax[0]\n",
    "\n",
    "ax_this.scatter(sats_maybe[\"radius\"]*sats_maybe[\"HOST_DIST\"]*1000/3600/180*np.pi, sats_maybe[\"sb_r\"], marker='o', c='None', edgecolors='k', lw=0.6)\n",
    "CS = ax_this.scatter(sats[\"radius\"]*sats[\"HOST_DIST\"]*1000/3600/180*np.pi, sats[\"sb_r\"], c=sats[\"log_sm\"], vmin=5.5, vmax=10, alpha=0.95, rasterized=True, label=\"SAGA\")\n",
    "ax_this.scatter(mw[\"reff_kpc\"], mw[\"sb_r_eff\"], marker='*', c=mw[\"log_sm\"], vmin=5.5, vmax=10, lw=0.6, edgecolors=plt.cm.gray(0.2), s=50, alpha=1, rasterized=True, label=\"MW\")\n",
    "\n",
    "ax_this.set_xscale(\"log\")\n",
    "ax_this.set_xlim(0.1, 10)\n",
    "ax_this.set_ylim(19, 26.5)\n",
    "ax_this.set_xlabel(\"$R_{r, \\\\mathrm{eff}}$ [kpc]\")\n",
    "ax_this.set_ylabel(\"$\\\\mu_{r_o, \\\\mathrm{eff}}$ [mag arcsec$^{-2}$]\")\n",
    "ax_this.legend(loc=\"lower right\", handletextpad=0.1, fontsize=13, borderpad=0.25, handlelength=1, frameon=True, markerfirst=False)\n",
    "\n",
    "y0, y1 = ax_this.get_ylim()\n",
    "x0, x1 = ax_this.get_xlim()\n",
    "ax_this.axvline(1.5, ymin=(24.7-y0)/(y1-y0), color=plt.cm.gray(0.4), zorder=-99, ls=\":\")\n",
    "ax_this.axhline(24.7, xmin=np.log(1.5/x0)/np.log(x1/x0), color=plt.cm.gray(0.4), zorder=-99, ls=\":\")\n",
    "ax_this.text(1.5*1.05, 26.1, \"UDGs\", fontsize=11)\n",
    "\n",
    "s = spearmanr(sats[\"radius\"]*sats[\"HOST_DIST\"], sats[\"sb_r\"])\n",
    "print(\"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue, 2)))\n",
    "\n",
    "s = spearmanr(np.hstack([sats[\"radius\"]*sats[\"HOST_DIST\"], sats_maybe[\"radius\"]*sats_maybe[\"HOST_DIST\"]]), np.hstack([sats[\"sb_r\"], sats_maybe[\"sb_r\"]]))\n",
    "print(\"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue, 2)))\n",
    "\n",
    "\n",
    "#####\n",
    "ax_this = ax[1]\n",
    "\n",
    "ax_this.scatter(sats_maybe[\"gr\"], sats_maybe[\"sb_r\"], marker='o', c='None', edgecolors='k', lw=0.5)\n",
    "CS = ax_this.scatter(sats[\"gr\"], sats[\"sb_r\"], c=sats[\"log_sm\"], vmin=5.5, vmax=10, alpha=0.95, rasterized=True)\n",
    "ax_this.scatter(mw[\"gr\"], mw[\"sb_r_eff\"], marker='*', c=mw[\"log_sm\"], vmin=5.5, vmax=10, lw=0.5, edgecolors=plt.cm.gray(0.2), s=50, alpha=0.95, rasterized=True, zorder=99)\n",
    "\n",
    "ax_this.set_xlim(0, 0.8)\n",
    "ax_this.set_ylim(19, 26.5)\n",
    "ax_this.set_xlabel(\"$(g-r)_o$\")\n",
    "ax_this.set_ylabel(\"$\\\\mu_{r_o, \\\\mathrm{eff}}$ [mag arcsec$^{-2}$]\")\n",
    "\n",
    "\n",
    "s = spearmanr(sats[\"gr\"], sats[\"sb_r\"])\n",
    "print(\"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue, 2)))\n",
    "\n",
    "s = spearmanr(np.hstack([sats[\"gr\"], sats_maybe[\"gr\"]]), np.hstack([sats[\"sb_r\"], sats_maybe[\"sb_r\"]]))\n",
    "print(\"$({:.2f}, {})$\".format(s.correlation, print_sci(s.pvalue, 2)))\n",
    "\n",
    "\n",
    "fig.tight_layout()\n",
    "cbar = plt.colorbar(CS, ax=ax, pad=0.02)\n",
    "cbar.ax.set_ylabel(\"$\\\\log\\;[M_*/M_\\\\odot]$\", labelpad=0, fontsize=13)\n",
    "cbar.ax.set_rasterized(True)\n",
    "\n",
    "plt.savefig(\"/home/yymao/Downloads/sb-radius-color.pdf\", dpi=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quenched fraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fq_lg = ascii.read(format=\"fast_tab\", names=[\"log_sm\", \"fq\", \"fq_ue\", \"fq_le\"], data_start=0,\n",
    "    table=\"\"\"3.49\t1.00\t1.00\t0.77\n",
    "4.49\t1.00\t1.00\t0.90\n",
    "5.49\t0.94\t0.96\t0.83\n",
    "6.49\t1.00\t1.00\t0.77\n",
    "7.49\t0.80\t0.88\t0.54\n",
    "8.49\t0.66\t0.81\t0.38\n",
    "9.18\t0.00\t0.60\t0.00\n",
    "\"\"\")\n",
    "\n",
    "fq_tinker = ascii.read(format=\"fast_tab\", names=[\"log_sm0\", \"fq\", \"log_sm1\",  \"fq_ue\", \"log_sm2\", \"fq_le\"], data_start=0,\n",
    "    table=\"\"\"9.7632\t0.0625\t9.7655\t0.0455\t9.7633\t0.0788\n",
    "9.8511\t0.1149\t9.8510\t0.0986\t9.8511\t0.1312\n",
    "9.9459\t0.1340\t9.9459\t0.1170\t9.9459\t0.1503\n",
    "10.0479\t0.1835\t10.0479\t0.1637\t10.0480\t0.2019\n",
    "10.1477\t0.2514\t10.1476\t0.2295\t10.1477\t0.2734\n",
    "10.2449\t0.2790\t10.2472\t0.2563\t10.2473\t0.3010\n",
    "10.3494\t0.3618\t10.3493\t0.3363\t10.3495\t0.3859\n",
    "10.4468\t0.4433\t10.4467\t0.4142\t10.4445\t0.4709\n",
    "10.5488\t0.4666\t10.5487\t0.4340\t10.5489\t0.4985\n",
    "10.6485\t0.5381\t10.6484\t0.5033\t10.6486\t0.5728\n",
    "10.7458\t0.5812\t10.7481\t0.5437\t10.7483\t0.6188\n",
    "10.8503\t0.6471\t10.8501\t0.5982\t10.8504\t0.6938\n",
    "10.9523\t0.7079\t10.9522\t0.6526\t10.9525\t0.7632\n",
    "11.0495\t0.7249\t11.0493\t0.6575\t11.0521\t0.7929\n",
    "11.1397\t0.7517\t11.1394\t0.6695\t11.1399\t0.8346\n",
    "11.2347\t0.8211\t11.2319\t0.6943\t11.2327\t0.9472\n",
    "\"\"\")\n",
    "fq_tinker[\"log_sm\"] = (fq_tinker[\"log_sm0\"] + fq_tinker[\"log_sm1\"] + fq_tinker[\"log_sm2\"]) / 3\n",
    "\n",
    "fq_geha = ascii.read(format=\"fast_tab\", names=[\"log_sm\", \"fq\",\"fq_e\"], data_start=0,\n",
    "    table=\"\"\"9.9\t0.084\t0.0052\n",
    "9.7\t0.031\t0.0032\n",
    "9.5\t0.017\t0.0023\n",
    "9.3\t0.007\t0.0016\n",
    "9.1\t0.002\t0.0014\n",
    "8.9\t0.000\t0.0019\n",
    "8.7\t0.000\t0.0028\n",
    "8.5\t0.000\t0.0043\n",
    "8.3\t0.000\t0.0065\n",
    "8.1\t0.000\t0.0091\n",
    "7.8\t0.000\t0.0093\n",
    "7.3\t0.000\t0.0193\n",
    "\"\"\")\n",
    "fq_geha[\"fq_ue\"] = fq_geha[\"fq\"] + fq_geha[\"fq_e\"]\n",
    "fq_geha[\"fq_le\"] = np.maximum(0, fq_geha[\"fq\"] - fq_geha[\"fq_e\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4.8, 3.6))\n",
    "\n",
    "plt.fill_between(fq_tinker[\"log_sm\"], fq_tinker[\"fq_le\"], fq_tinker[\"fq_ue\"], color=\"C4\", alpha=0.7, lw=0, rasterized=True)\n",
    "plt.fill_between(fq_geha[\"log_sm\"], fq_geha[\"fq_le\"], fq_geha[\"fq_ue\"], color=\"C4\", alpha=0.7, rasterized=True)\n",
    "\n",
    "\n",
    "plt.errorbar(fq_lg[\"log_sm\"], fq_lg[\"fq\"], marker='x', ls=\"\", c=\"C1\", ms=6)\n",
    "plt.fill_between(fq_lg[\"log_sm\"], fq_lg[\"fq_le\"], fq_lg[\"fq_ue\"], color=\"C1\", alpha=0.25, lw=0, rasterized=True)\n",
    "\n",
    "pbins = np.linspace(0, 100, 8)\n",
    "log_sm = np.concatenate([sats['log_sm']])\n",
    "bins = np.percentile(log_sm, pbins)\n",
    "x = np.percentile(log_sm, (pbins[1:]+pbins[:-1])*0.5)\n",
    "\n",
    "interloper_corr = lambda d: 0.95 - d*0.001\n",
    "\n",
    "n_q = np.histogram(Query('EW_Halpha < 2').filter(sats, 'log_sm'), bins)[0]\n",
    "sats_q = Query('EW_Halpha < 2').filter(sats, [\"log_sm\", \"RHOST_KPC\"])\n",
    "sats_maybe_q = Query().filter(sats_maybe_all, [\"log_sm\", \"RHOST_KPC\", \"p_sat_corrected\"])\n",
    "n_q_corr = np.histogram(sats_q['log_sm'], bins, weights=interloper_corr(sats_q[\"RHOST_KPC\"]))[0]\n",
    "n_q_corr += np.histogram(sats_maybe_q['log_sm'], bins, weights=sats_maybe_q[\"p_sat_corrected\"]*interloper_corr(sats_maybe_q[\"RHOST_KPC\"]))[0]\n",
    "\n",
    "n_all = np.histogram(sats['log_sm'], bins)[0]\n",
    "n_all_corr = np.histogram(sats['log_sm'], bins, weights=interloper_corr(sats[\"RHOST_KPC\"]))[0]\n",
    "n_all_corr += np.histogram(sats_maybe_all['log_sm'], bins, weights=sats_maybe_all[\"p_sat_corrected\"]*interloper_corr(sats_maybe_all[\"RHOST_KPC\"]))[0]\n",
    "\n",
    "p = n_q/n_all\n",
    "pb = (n_q+1)/(n_all+2)\n",
    "perr = np.sqrt(pb*(1-pb) / n_all)\n",
    "p1 = np.minimum(1, p + perr)\n",
    "p2 = np.maximum(0, p - perr)\n",
    "p_corr = np.maximum(0, (n_q_corr/n_all_corr) - p)\n",
    "\n",
    "plt.errorbar(x, p, yerr=(p-p2, p1-p), ls='', marker='o', lw=2, capsize=2, ms=5, c=\"C2\")\n",
    "plt.errorbar(x, p, yerr=(np.zeros_like(p), p1-p+p_corr), ls='', marker='', alpha=0.6, c=\"C2\", lw=2, rasterized=True)\n",
    "\n",
    "names=[\"m_star\", \"f_quench\", \"e_shot_low\", \"e_shot_high\", \"e_correction\"]\n",
    "t = Table(np.vstack((x, p, p-p2, p1-p, p_corr)).T, names=names)\n",
    "for name in names:\n",
    "    t[name].format = \"%.3f\"\n",
    "t.pprint(-1, -1)\n",
    "\n",
    "this_work_label = \"This work\" #and \"Mao+20\"\n",
    "plt.xlim(5.3, 11)\n",
    "plt.text(6.75, 0.2, f\"SAGA sat.\\n({this_work_label})\", ha=\"right\", va=\"bottom\", color=\"C2\", fontsize=13, fontweight=\"bold\")\n",
    "plt.text(7.2, 0.91, \"MW+M31 sat. (Wetzel+15)\", ha=\"left\", va=\"bottom\", color=\"C1\", fontsize=13, fontweight=\"bold\")\n",
    "plt.text(10.7, 0.6, \"Field gal.\\n(Geha+12)\", ha=\"right\", va=\"bottom\", color=\"C4\", fontsize=13, fontweight=\"bold\")\n",
    "#plt.legend(loc=\"upper right\",  markerfirst=False, fontsize=13,frameon=True, handletextpad=0)\n",
    "plt.xlabel(\"$\\\\log\\;[M_*/M_\\\\odot]$\")\n",
    "plt.ylabel(\"Quenched fraction\")\n",
    "plt.ylim(-0.05, 1.05)\n",
    "plt.axhline(0, lw=0.5, color=\"grey\")\n",
    "plt.axhline(1, lw=0.5, color=\"grey\")\n",
    "plt.tight_layout()\n",
    "\n",
    "if this_work_label == \"This work\":\n",
    "    plt.savefig(\"/home/yymao/Downloads/f-quench.pdf\", bbox_inches='tight', dpi=200)"
   ]
  },
  {
   "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
}