03. Combine Reduced Products.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from glob import glob\n",
    "\n",
    "import os\n",
    "import sys\n",
    "sys.path.append(f'{os.environ[\"HOME\"]}/Projects/planckClusters/catalogs')\n",
    "from load_catalogs import load_PSZcatalog\n",
    "\n",
    "# parallel processor\n",
    "from utilities import parallel_process, system_call_env, system_call"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def combineXRT(name, outpath):\n",
    "    ''' This combines the individual observations\n",
    "    \n",
    "    IMPORTANT! The scripts that this (and other functions) creat are designed\n",
    "    to be run from the same directory as this notebook. They WILL NOT work \n",
    "    if you try to run them from the individual data directories.\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    if not os.path.isdir(f'{outpath}/{name}'):\n",
    "        return\n",
    "    \n",
    "    # find the x-ray files\n",
    "    files = glob(f'{outpath}/{name}/reduced/**/*xpcw[1-4]po_cl.evt', recursive=True)\n",
    "\n",
    "    # sort the files so we can control their input order\n",
    "    # this is used to make sure all of the data products are the same\n",
    "    # when we are combining observations\n",
    "    \n",
    "    files = np.sort(files)\n",
    "    \n",
    "    if len(files) < 1:\n",
    "        return\n",
    "\n",
    "    # write xsel.in\n",
    "    with open(f'{outpath}/{name}/{name}_xsel.in', 'w') as f:\n",
    "        for i, f_in in enumerate(files):\n",
    "            f_parts = f_in.split('/')\n",
    "            if i == 0:\n",
    "                # this is the session name... specify random letters to run many at once.\n",
    "                f.writelines(f'{name}\\n')\n",
    "                f.writelines('read events\\n')\n",
    "                # set the data directory\n",
    "                f.writelines('/'.join(f_parts[:3]) + '\\n')\n",
    "                # first entry\n",
    "                f.writelines('/'.join(f_parts[3:]) + '\\n')\n",
    "                f.writelines('yes\\n')\n",
    "                continue\n",
    "\n",
    "            f.writelines('read events\\n')\n",
    "            f.writelines('/'.join(f_parts[3:]) + '\\n')\n",
    "            # if you try to read more than 20 exposures, it says \"more?\"\n",
    "            if i >= 19:\n",
    "                f.writelines('\\n')\n",
    "            if i >= 42:\n",
    "                f.writelines('\\n')\n",
    "            if i >= 65:\n",
    "                f.writelines('\\n')\n",
    "            if i >= 88:\n",
    "                f.writelines('\\n')\n",
    "            if i >= 111:\n",
    "                f.writelines('\\n')\n",
    "            if i >= 134:\n",
    "                f.writelines('\\n')\n",
    "            if i >= 157:\n",
    "                f.writelines('\\n')\n",
    "            if i >= 180:\n",
    "                f.writelines('\\n')\n",
    "            if i >= 203:\n",
    "                f.writelines('\\n')\n",
    "            if i >= 226:\n",
    "                f.writelines('\\n')\n",
    "            if i >= 249:\n",
    "                f.writelines('\\n')\n",
    "        f.writelines('extract events\\n')\n",
    "        f.writelines(f'save events {outpath}/{name}/{name}_events.fits\\n')\n",
    "        if os.path.isfile(f'{outpath}/{name}/{name}_events.fits'):\n",
    "            f.writelines('yes\\n')\n",
    "        f.writelines('yes\\n')\n",
    "\n",
    "                     \n",
    "        f.writelines('set phaname PI\\n')\n",
    "        # here we are going to make a few binned images for a few different energy ranges\n",
    "        # energies in loop\n",
    "        for eng in [200, 300, 400, 500, 600]:\n",
    "            f.writelines(f'filter pha_cutoff 50 {eng}\\n')\n",
    "            \n",
    "            # save non-binned image -- the yes's are to overwrite if file is already there\n",
    "            f.writelines('set xybinsize 1\\n')\n",
    "            f.writelines('extract image\\n')\n",
    "            f.writelines(f'save image {\"/\".join(f_parts[:3])}/{name}_img_50-{eng}.fits\\n')\n",
    "            if os.path.isfile(f'{outpath}/{name}/{name}_img_50-{eng}.fits'):\n",
    "                f.writelines('yes\\n')\n",
    "\n",
    "            # save binned image -- see above\n",
    "            f.writelines('set xybinsize 8\\n')\n",
    "            f.writelines('extract image\\n')\n",
    "            f.writelines(f'save image {\"/\".join(f_parts[:3])}/{name}_img_50-{eng}_bl8.fits\\n')\n",
    "            if os.path.isfile(f'{outpath}/{name}/{name}_img_50-{eng}_bl8.fits'):\n",
    "                f.writelines('yes\\n')\n",
    "\n",
    "            f.writelines('set xybinsize 4\\n')\n",
    "            f.writelines('extract image\\n')\n",
    "            f.writelines(f'save image {\"/\".join(f_parts[:3])}/{name}_img_50-{eng}_bl4.fits\\n')\n",
    "            if os.path.isfile(f'{outpath}/{name}/{name}_img_50-{eng}_bl4.fits'):\n",
    "                f.writelines('yes\\n')\n",
    "\n",
    "        f.writelines('exit\\n')\n",
    "        f.writelines('no\\n')\n",
    "\n",
    "                         \n",
    "    # call\n",
    "    cmd = f'xselect @{outpath}/{name}/{name}_xsel.in'\n",
    "#     stdout, stderr = system_call_env(cmd)\n",
    "\n",
    "    # log the output\n",
    "    log_file = f'{outpath}/{name}/{name}_xsel.log'\n",
    "#     with open(log_file, 'w') as f:\n",
    "#         f.writelines(stdout)\n",
    "\n",
    "    os.system(f'{cmd} > {log_file}')                  \n",
    "                         \n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def combineXRT_exp(name, outpath):\n",
    "    ''' This combines the exposure maps'''\n",
    "    \n",
    "    if not os.path.isdir(f'{outpath}/{name}'):\n",
    "        return\n",
    "    \n",
    "    # find the x-ray files\n",
    "    files = glob(f'{outpath}/{name}/reduced/**/*xpcw[1-4]po_ex.img', recursive=True)\n",
    "\n",
    "    # sort the observations -- see above\n",
    "    files = np.sort(files)\n",
    "    \n",
    "    if len(files) < 1:\n",
    "        return name\n",
    "\n",
    "    # remove the old file if it is there\n",
    "    if os.path.isfile(f'{outpath}/{name}/{name}_exp.fits'):\n",
    "        os.remove(f'{outpath}/{name}/{name}_exp.fits')\n",
    "\n",
    "    # write xsel.in\n",
    "    with open(f'{outpath}/{name}/{name}_ximg_exp.in', 'w') as f:\n",
    "        for i, f_in in enumerate(files):\n",
    "            f_parts = f_in.split('/')\n",
    "            f.writelines(f'read {f_in}\\n')\n",
    "            if i == 0:\n",
    "                continue\n",
    "            f.writelines('sum\\n')\n",
    "            f.writelines('save\\n')\n",
    "\n",
    "\n",
    "        f.writelines(f'write/fits {\"/\".join(f_parts[:3])}/{name}_exp.fits\\n')\n",
    "\n",
    "        f.writelines('exit\\n')\n",
    "\n",
    "#     # call\n",
    "    cmd = f'ximage @{outpath}/{name}/{name}_ximg_exp.in'\n",
    "    stdout, stderr = system_call_env(cmd)\n",
    "    print(stdout)\n",
    "    # log the output\n",
    "    log_file = f'{outpath}/{name}/{name}_ximg_exp.log'\n",
    "    with open(log_file, 'w') as f:\n",
    "        f.writelines(stdout)                 \n",
    "    \n",
    "#     os.system(f'{cmd} > {log_file}')                 \n",
    "                     \n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# get file data\n",
    "data = load_PSZcatalog()\n",
    "data = data.sort_index(axis=1)\n",
    "\n",
    "outpath = './data_full_new'\n",
    "\n",
    "arr = [{'name':n.replace(' ', '_'), 'outpath':outpath} for n in data['NAME']]\n",
    "parallel_process(arr, combineXRT, use_kwargs=True, n_jobs=7)\n",
    "parallel_process(arr, combineXRT_exp, use_kwargs=True, n_jobs=7) #4.5 hours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# outpath = './data_full'\n",
    "# name = 'PSZ2_G313.33+61.13'\n",
    "# combineXRT(name, outpath)\n",
    "# combineXRT_exp(name, outpath)"
   ]
  }
 ],
 "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.7.1"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}