{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting Galaxy Cluster Member Colors in Extragalactic Catalogs\n",
    "\n",
    "Owners: **Dan Korytov [@dkorytov](https://github.com/LSSTDESC/DC2-analysis/issues/new?body=@dkorytov)**\n",
    "\n",
    "Last verified run: Mar 8, 2024 (by @patricialarsen)\n",
    "\n",
    "This notebook demonstrates how to access the extra galactic catalog through the Generic Catalog Reader (GCR, https://github.com/yymao/generic-catalog-reader) as well as how filter on galaxy features and cluster membership.\n",
    "\n",
    "__Objectives__:\n",
    "\n",
    "After working through and studying this Notebook you should be able to\n",
    "\n",
    "1. Access extragalactic catalogs (cosmoDC2, skysim, roman_rubin) through the GCR\n",
    "2. Filter on galaxy properties\n",
    "3. Select and plot cluster members\n",
    "\n",
    "\n",
    "__Logistics__: This notebook is intended to be run through the JupyterHub NERSC interface available here: https://jupyter.nersc.gov. To setup your NERSC environment, please follow the instructions available here: https://confluence.slac.stanford.edu/display/LSSTDESC/Using+Jupyter+at+NERSC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import GCRCatalogs\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as clr\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gc = GCRCatalogs.load_catalog('cosmoDC2_v1.1.4_small')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other catalog options:\n",
    "Some alternate catalogs that you might want to look at include skysim, which is a 5000 square degree extragalactic catalog, and the Roman Rubin catalog which is a smaller extragalactic catalog and image simulation. To load these instead execute one of the two cells below.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gc = GCRCatalogs.load_catalog('skysim5000_v1.1.2_small')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gc = GCRCatalogs.load_catalog('roman_rubin_2023_v1.1.3_elais')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "data = gc.get_quantities(['halo_mass', 'redshift',\n",
    "                        'mag_u', 'mag_g', 'mag_r', \n",
    "                        'mag_i', 'mag_z'], filters=['halo_mass > 3e13'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading catalog\n",
    "We load in the catalog with the \"load_catalog\" command, and then the values with the \"get_quantities\" command using filters to select sub-samples of the catalog. For this case we only need the magnitudes in several filters and the redshift. Galaxies are filtered on host halo mass to be at least 3e13 h$^{-1}$M$_\\odot$.\n",
    "\n",
    "### Help for error messages:\n",
    "If this fails to find the appropriate quantities, check that the desc-python kernel is being used and if this is not available source the kernels by running the following command on a terminal at nersc: \"source \n",
    "/global/common/software/lsst/common/miniconda/kernels/setup.sh\"\n",
    "\n",
    "We are loading in a smaller version of the full cosmoDC2 catalog - this contains the same information as the full catalog but with a smaller sky area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "h,xbins = np.histogram(np.log10(data['halo_mass']),bins=40)\n",
    "xbins_avg = (xbins[1:]+xbins[:-1])/2.0\n",
    "plt.semilogy(xbins_avg, h)\n",
    "plt.ylabel(r'Galaxy Count')\n",
    "plt.xlabel(r'log10( M$_{\\rm{halo}}$ / M$_\\odot)$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a sanity check, we made sure no galaxies have a host halo below 3e13 h$^{-1}$ M$_\\odot$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "gal_clr = data['mag_g']-data['mag_r']\n",
    "plt.hist2d(data['redshift'], gal_clr, bins=100, cmap='PuBu', norm=clr.LogNorm())\n",
    "plt.colorbar(label='population density')\n",
    "plt.ylabel('Observed g-r')\n",
    "plt.xlabel('redshift')\n",
    "plt.title('Galaxy Colors in Clusters')\n",
    "plt.tight_layout()\n",
    "\n",
    "plt.figure()\n",
    "gal_clr = data['mag_r']-data['mag_i']\n",
    "plt.hist2d(data['redshift'], gal_clr, bins=100, cmap='PuBu',norm=clr.LogNorm())\n",
    "plt.colorbar(label='population density')\n",
    "plt.ylabel('r-i')\n",
    "plt.xlabel('redshift')\n",
    "plt.title('Galaxy Colors in Clusters')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After loading the catalog and data, we are able to plot the color distribution as a function of redshift for cluster galaxies with very little effort!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "desc-python",
   "language": "python",
   "name": "desc-python"
  },
  "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}