{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AKARI-SEP - Merging HELP data products\n",
    "\n",
    "This notebook merges the various HELP data products on AKARI-SEP.\n",
    "\n",
    "It is first used to create a catalogue that will be used for SED fitting by CIGALE by merging the optical master list, the photo-z and the XID+ far infrared fluxes.  Then, this notebook is used to incorporate the CIGALE physical parameter estimations and generate the final HELP data product on the field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set this to true to produce only the catalogue for CIGALE and to false \n",
    "# to continue and merge the CIGALE results too.\n",
    "MAKE_CIGALE_INPUT = True\n",
    "MAKE_FINAL_CAT = True\n",
    "\n",
    "SUFFIX = '20180221'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from astropy.table import Column, MaskedColumn, Table, join, vstack\n",
    "\n",
    "from herschelhelp.filters import get_filter_meta_table\n",
    "\n",
    "from herschelhelp_internal.utils import add_column_meta\n",
    "\n",
    "filter_mean_lambda = {\n",
    "    item['filter_id']: item['mean_wavelength'] for item in\n",
    "    get_filter_meta_table()\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reading the masterlist, XID+, and photo-z catalogues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Master list\n",
    "\n",
    "ml = Table.read(\n",
    "    \"../../dmu1/dmu1_ml_AKARI-SEP/data/master_catalogue_akari-sep_{}.fits\".format(SUFFIX))\n",
    "ml.meta = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# # XID+ MIPS24\n",
    "\n",
    "xid_mips24 = Table.read(\"../../dmu26/dmu26_XID+MIPS_AKARI-SEP/data/\"\n",
    "                          \"dmu26_XID+MIPS_AKARI-SEP_cat_20190909.fits\")\n",
    "xid_mips24.meta = None\n",
    "\n",
    "# Adding the error column\n",
    "xid_mips24.add_column(Column(\n",
    "    data=np.max([xid_mips24['FErr_MIPS_24_u'] - xid_mips24['F_MIPS_24'],\n",
    "                 xid_mips24['F_MIPS_24'] - xid_mips24['FErr_MIPS_24_l']],\n",
    "                axis=0),\n",
    "    name=\"ferr_mips_24\"\n",
    "))\n",
    "xid_mips24['F_MIPS_24'].name = \"f_mips_24\"\n",
    "xid_mips24 = xid_mips24['help_id', 'f_mips_24', 'ferr_mips_24', 'flag_mips_24']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# XID+ PACS\n",
    "\n",
    "xid_pacs = Table.read(\"../../dmu26/dmu26_XID+PACS_AKARI-SEP/data/\"\n",
    "                        \"dmu26_XID+PACS_AKARI-SEP_cat_20190924.fits\")\n",
    "xid_pacs.meta = None\n",
    "\n",
    "# Convert from mJy to μJy\n",
    "for col in [\"F_PACS_100\", \"FErr_PACS_100_u\", \"FErr_PACS_100_l\",\n",
    "            \"F_PACS_160\", \"FErr_PACS_160_u\", \"FErr_PACS_160_l\"]:\n",
    "    xid_pacs[col] *= 1000\n",
    "\n",
    "xid_pacs.add_column(Column(\n",
    "    data=np.max([xid_pacs['FErr_PACS_100_u'] - xid_pacs['F_PACS_100'],\n",
    "                 xid_pacs['F_PACS_100'] - xid_pacs['FErr_PACS_100_l']],\n",
    "                axis=0),\n",
    "   name=\"ferr_pacs_green\"\n",
    "))\n",
    "xid_pacs['F_PACS_100'].name = \"f_pacs_green\"\n",
    "xid_pacs['flag_PACS_100'].name = \"flag_pacs_green\"\n",
    "\n",
    "xid_pacs.add_column(Column(\n",
    "    data=np.max([xid_pacs['FErr_PACS_160_u'] - xid_pacs['F_PACS_160'],\n",
    "                 xid_pacs['F_PACS_160'] - xid_pacs['FErr_PACS_160_l']],\n",
    "                axis=0),\n",
    "    name=\"ferr_pacs_red\"\n",
    "))\n",
    "xid_pacs['F_PACS_160'].name = \"f_pacs_red\"\n",
    "xid_pacs['flag_PACS_160'].name = \"flag_pacs_red\"\n",
    "\n",
    "xid_pacs = xid_pacs['help_id', 'f_pacs_green', 'ferr_pacs_green',\n",
    "                    'flag_pacs_green', 'f_pacs_red', 'ferr_pacs_red',\n",
    "                    'flag_pacs_red']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# # XID+ SPIRE\n",
    "\n",
    "xid_spire = Table.read(\"../../dmu26/dmu26_XID+SPIRE_AKARI-SEP/data/\"\n",
    "                         \"dmu26_XID+SPIRE_AKARI-SEP_cat_20190916.fits\")\n",
    "xid_spire.meta = None\n",
    "\n",
    "xid_spire['HELP_ID'].name = \"help_id\"\n",
    "\n",
    "# Convert from mJy to μJy\n",
    "for col in [\"F_SPIRE_250\", \"FErr_SPIRE_250_u\", \"FErr_SPIRE_250_l\",\n",
    "            \"F_SPIRE_350\", \"FErr_SPIRE_350_u\", \"FErr_SPIRE_350_l\",\n",
    "            \"F_SPIRE_500\", \"FErr_SPIRE_500_u\", \"FErr_SPIRE_500_l\"]:\n",
    "    xid_spire[col] *= 1000\n",
    "\n",
    "xid_spire.add_column(Column(\n",
    "    data=np.max([xid_spire['FErr_SPIRE_250_u'] - xid_spire['F_SPIRE_250'],\n",
    "                 xid_spire['F_SPIRE_250'] - xid_spire['FErr_SPIRE_250_l']],\n",
    "                axis=0),\n",
    "    name=\"ferr_spire_250\"\n",
    "))\n",
    "xid_spire['F_SPIRE_250'].name = \"f_spire_250\"\n",
    "xid_spire.add_column(Column(\n",
    "    data=np.max([xid_spire['FErr_SPIRE_350_u'] - xid_spire['F_SPIRE_350'],\n",
    "                 xid_spire['F_SPIRE_350'] - xid_spire['FErr_SPIRE_350_l']],\n",
    "                axis=0),\n",
    "    name=\"ferr_spire_350\"\n",
    "))\n",
    "xid_spire['F_SPIRE_350'].name = \"f_spire_350\"\n",
    "xid_spire.add_column(Column(\n",
    "    data=np.max([xid_spire['FErr_SPIRE_500_u'] - xid_spire['F_SPIRE_500'],\n",
    "                 xid_spire['F_SPIRE_500'] - xid_spire['FErr_SPIRE_500_l']],\n",
    "                axis=0),\n",
    "    name=\"ferr_spire_500\"\n",
    "))\n",
    "xid_spire['F_SPIRE_500'].name = \"f_spire_500\"\n",
    "\n",
    "xid_spire = xid_spire['help_id',\n",
    "                      'f_spire_250', 'ferr_spire_250', 'flag_spire_250',\n",
    "                      'f_spire_350', 'ferr_spire_350', 'flag_spire_350',\n",
    "                      'f_spire_500', 'ferr_spire_500', 'flag_spire_500']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Photo-z\n",
    "\n",
    "photoz = Table.read(\"../../dmu24/dmu24_AKARI-SEP/data/AKARI-SEP_DESI-DR8_Zou_et_al_2019_photo_z_withhelp_id.fits\")\n",
    "photoz.meta = None\n",
    "\n",
    "photoz = photoz['help_id', 'photo_z']\n",
    "photoz['photo_z'].name = 'redshift'\n",
    "\n",
    "photoz['redshift'][photoz['redshift'] < 0] = np.nan  # -99 used for missing values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/pyenv/versions/3.7.2/lib/python3.7/site-packages/astropy/table/column.py:965: RuntimeWarning: invalid value encountered in less\n",
      "  return getattr(self.data, op)(other)\n"
     ]
    }
   ],
   "source": [
    "# Temp spec-z\n",
    "\n",
    "\n",
    "\n",
    "ml['zspec'][ml['zspec'] < 0] = np.nan  # -99 used for missing values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Flags\n",
    "flags = Table.read(\"../../dmu6/dmu6_v_AKARI-SEP/data/akari-sep_20180221_flags.fits\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merging"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "merged_table = join(ml, xid_mips24, join_type='left')\n",
    "\n",
    "# Fill values\n",
    "for col in xid_mips24.colnames:\n",
    "    if col.startswith(\"f_\") or col.startswith(\"ferr_\"):\n",
    "        merged_table[col].fill_value = np.nan\n",
    "    elif col.startswith(\"flag_\"):\n",
    "        merged_table[col].fill_value = False\n",
    "merged_table = merged_table.filled()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "merged_table = join(merged_table, xid_pacs, join_type='left')\n",
    "        \n",
    "# Fill values\n",
    "for col in xid_pacs.colnames:\n",
    "    if col.startswith(\"f_\") or col.startswith(\"ferr_\"):\n",
    "        merged_table[col].fill_value = np.nan\n",
    "    elif col.startswith(\"flag_\"):\n",
    "        merged_table[col].fill_value = False\n",
    "merged_table = merged_table.filled()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "merged_table = join(merged_table, xid_spire, join_type='left')\n",
    "        \n",
    "# Fill values\n",
    "for col in xid_spire.colnames:\n",
    "    if col.startswith(\"f_\") or col.startswith(\"ferr_\"):\n",
    "        merged_table[col].fill_value = np.nan\n",
    "    elif col.startswith(\"flag_\"):\n",
    "        merged_table[col].fill_value = False\n",
    "merged_table = merged_table.filled()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "merged_table = join(merged_table, photoz, join_type='left')\n",
    "\n",
    "# Fill values\n",
    "merged_table['redshift'].fill_value = np.nan\n",
    "merged_table = merged_table.filled()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for col in flags.colnames:\n",
    "    if 'flag' in col:\n",
    "        try:\n",
    "            merged_table.remove_column(col)\n",
    "        except KeyError:\n",
    "            print(\"Column: {} not in masterlist.\".format(col))\n",
    "        \n",
    "merged_table = join(merged_table, flags, join_type='left')\n",
    "# Fill values\n",
    "for col in merged_table.colnames:\n",
    "    if 'flag' in col:\n",
    "        merged_table[col].fill_value = False\n",
    "merged_table = merged_table.filled()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Saving the catalogue for CIGALE (first run)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "set()\n",
      "set()\n"
     ]
    }
   ],
   "source": [
    "if MAKE_CIGALE_INPUT:\n",
    "    \n",
    "    # Sorting the columns\n",
    "    bands_tot = [col[2:] for col in merged_table.colnames\n",
    "             if col.startswith('f_') and not col.startswith('f_ap')]\n",
    "    bands_ap = [col[5:] for col in merged_table.colnames\n",
    "             if col.startswith('f_ap_') ]\n",
    "    bands = list(set(bands_tot) | set(bands_ap))\n",
    "    bands.sort(key=lambda x: filter_mean_lambda[x])\n",
    "    \n",
    "    columns = ['help_id', 'field', 'ra', 'dec', 'hp_idx', 'ebv', 'redshift', \n",
    "               'zspec']\n",
    "    for band in bands:\n",
    "        for col_tpl in ['f_{}', 'ferr_{}', 'f_ap_{}', 'ferr_ap_{}',\n",
    "                        'm_{}', 'merr_{}', 'm_ap_{}', 'merr_ap_{}',\n",
    "                        'flag_{}']:\n",
    "            colname = col_tpl.format(band)\n",
    "            if colname in merged_table.colnames:\n",
    "                columns.append(colname)\n",
    "    columns += ['stellarity', 'stellarity_origin', 'flag_cleaned',\n",
    "                'flag_merged', 'flag_gaia', 'flag_optnir_obs',\n",
    "                'flag_optnir_det', 'zspec_qual', 'zspec_association_flag']\n",
    "    \n",
    "\n",
    "    # Check that we did not forget any column\n",
    "    print(set(columns) - set(merged_table.colnames))\n",
    "    print( set(merged_table.colnames) - set(columns) )\n",
    "    \n",
    "    merged_table = add_column_meta(merged_table, '../columns.yml')\n",
    "    merged_table[columns].write(\"data/AKARI-SEP_{}_cigale.fits\".format(SUFFIX), overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Merging CIGALE outputs\n",
    "\n",
    "We merge the CIGALE outputs to the main catalogue. The CIGALE products provides several χ² with associated thresholds. For simplicity, we convert these two values to flags."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if MAKE_FINAL_CAT:\n",
    "\n",
    "    # Cigale outputs\n",
    "    cigale = Table.read(\"../../dmu28/dmu28_AKARI-SEP/data/zphot/HELP_final_results.fits\")\n",
    "    cigale['id'].name = \"help_id\"\n",
    "\n",
    "    # We convert the various Chi2 and threshold to flags\n",
    "    # UPDATE - we no longer do this and instead provide the chi2 values \n",
    "    #flag_cigale_opt = cigale[\"UVoptIR_OPTchi2\"] <= cigale[\"UVoptIR_OPTchi2_threshold\"]\n",
    "    #flag_cigale_ir = cigale[\"UVoptIR_IRchi2\"] <= cigale[\"UVoptIR_IRchi2_threshold\"]\n",
    "    #flag_cigale = (\n",
    "    #    (cigale[\"UVoptIR_best.reduced_chi_square\"] \n",
    "    #         <=  cigale[\"UVoptIR_best.reduced_chi_square_threshold\"]) &\n",
    "    #    flag_cigale_opt & flag_cigale_ir)\n",
    "    #flag_cigale_ironly = cigale[\"IRonly_IRchi2\"] <= cigale[\"IRonly_IRchi2_threshold\"]\n",
    "\n",
    "    #cigale.add_columns([\n",
    "    #    MaskedColumn(flag_cigale, \"flag_cigale\", \n",
    "    #                 dtype=int, fill_value=-1),\n",
    "    #    MaskedColumn(flag_cigale_opt, \"flag_cigale_opt\", \n",
    "    #                 dtype=int, fill_value=-1),\n",
    "    #    MaskedColumn(flag_cigale_ir, \"flag_cigale_ir\", \n",
    "    #                 dtype=int, fill_value=-1),\n",
    "    #    MaskedColumn(flag_cigale_ironly, \"flag_cigale_ironly\", \n",
    "    #                 dtype=int, fill_value=-1)\n",
    "    #])\n",
    "    \n",
    "    #UPDATE Instead of flags above we include all cols used to make \n",
    "    #deprecated flags\n",
    "    cigale[\"UVoptIR_OPTchi2\"].name                 =  \"cigale_chi2_opt\"\n",
    "    cigale[\"UVoptIR_IRchi2\"].name                  =  \"cigale_chi2_ir\"\n",
    "    cigale[\"UVoptIR_best.reduced_chi_square\"].name =  \"cigale_chi2_red\"\n",
    "    cigale[\"IRonly_IRchi2\"].name                   =  \"cigale_chi2_ironly\"\n",
    "    \n",
    "    cigale['UVoptIR_bayes.stellar.m_star'].name =  \"cigale_mstar\"\n",
    "    cigale['UVoptIR_bayes.stellar.m_star_err'].name = \"cigale_mstar_err\"\n",
    "    cigale['UVoptIR_bayes.sfh.sfr10Myrs'].name = \"cigale_sfr\"\n",
    "    cigale['UVoptIR_bayes.sfh.sfr10Myrs_err'].name = \"cigale_sfr_err\"\n",
    "    cigale['UVoptIR_bayes.dust.luminosity'].name = \"cigale_dustlumin\"\n",
    "    cigale['UVoptIR_bayes.dust.luminosity_err'].name = \"cigale_dustlumin_err\"\n",
    "    cigale['IR_bayes.dust.luminosity'].name = \"cigale_dustlumin_ironly\"\n",
    "    cigale['IR_bayes.dust.luminosity_err'].name = \"cigale_dustlumin_ironly_err\"\n",
    "\n",
    "    cigale = cigale['help_id', 'cigale_mstar', 'cigale_mstar_err', 'cigale_sfr',\n",
    "                    'cigale_sfr_err', 'cigale_dustlumin', 'cigale_dustlumin_err', \n",
    "                    'cigale_dustlumin_ironly', 'cigale_dustlumin_ironly_err', \n",
    "                    #Old flags\n",
    "                    #'flag_cigale', 'flag_cigale_opt', 'flag_cigale_ir', 'flag_cigale_ironly'\n",
    "                    #New chi2 values\n",
    "                    \"cigale_chi2_opt\", \"cigale_chi2_ir\", \"cigale_chi2_red\", \"cigale_chi2_ironly\"\n",
    "                   ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if MAKE_FINAL_CAT:\n",
    "\n",
    "    merged_table = join(merged_table, cigale, join_type='left')\n",
    "\n",
    "    # Fill values\n",
    "    for col in cigale.colnames:\n",
    "        if col.startswith(\"cigale_\"):\n",
    "            merged_table[col].fill_value = np.nan\n",
    "        elif col.startswith(\"flag_\"):\n",
    "            merged_table[col].fill_value = -1\n",
    "    merged_table = merged_table.filled()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sorting columns\n",
    "\n",
    "We sort the columns by increasing band wavelength."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if MAKE_FINAL_CAT:\n",
    "\n",
    "    bands = [col[2:] for col in merged_table.colnames\n",
    "             if col.startswith('f_') and not col.startswith('f_ap')]\n",
    "    bands.sort(key=lambda x: filter_mean_lambda[x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if MAKE_FINAL_CAT:\n",
    "\n",
    "    columns = ['help_id', 'field', 'ra', 'dec', 'hp_idx', 'ebv', 'redshift', 'zspec']\n",
    "    for band in bands:\n",
    "        for col_tpl in ['f_{}', 'ferr_{}', 'f_ap_{}', 'ferr_ap_{}',\n",
    "                        'm_{}', 'merr_{}', 'm_ap_{}', 'merr_ap_{}',\n",
    "                        'flag_{}']:\n",
    "            colname = col_tpl.format(band)\n",
    "            if colname in merged_table.colnames:\n",
    "                columns.append(colname)\n",
    "    columns += ['cigale_mstar', 'cigale_mstar_err', 'cigale_sfr', 'cigale_sfr_err',\n",
    "                'cigale_dustlumin', 'cigale_dustlumin_err', 'cigale_dustlumin_ironly', \n",
    "                'cigale_dustlumin_ironly_err', \n",
    "                #Old flags\n",
    "                #'flag_cigale', 'flag_cigale_opt', 'flag_cigale_ir', 'flag_cigale_ironly'\n",
    "                #New chi2 values\n",
    "                'cigale_chi2_opt', 'cigale_chi2_ir', 'cigale_chi2_red', 'cigale_chi2_ironly',\n",
    "                'stellarity', \n",
    "                'stellarity_origin', 'flag_cleaned', 'flag_merged', 'flag_gaia', \n",
    "                'flag_optnir_obs', 'flag_optnir_det', 'zspec_qual', \n",
    "                'zspec_association_flag']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if MAKE_FINAL_CAT:\n",
    "\n",
    "    # Check that we did not forget any column\n",
    "    assert set(columns) == set(merged_table.colnames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if MAKE_FINAL_CAT:\n",
    "\n",
    "    merged_table = add_column_meta(merged_table, '../columns.yml')\n",
    "    merged_table[columns].write(\"data/AKARI-SEP.fits\", overwrite=True)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python (herschelhelp_internal)",
   "language": "python",
   "name": "helpint"
  },
  "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}