04. Detect and Clean Sources -- ciaotools.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from astropy import table\n",
    "from astropy.table import Table, Column\n",
    "from tqdm import tqdm_notebook\n",
    "import tempfile\n",
    "import emcee\n",
    "import numpy.ma as ma\n",
    "import numpy as np\n",
    "import subprocess\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",
    "from astropy.io import fits\n",
    "from astropy import wcs\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "                \n",
    "from regions import read_ds9, write_ds9                \n",
    "from regions import PixCoord, EllipsePixelRegion, CircleAnnulusPixelRegion, CirclePixelRegion\n",
    "\n",
    "from math import sqrt, sin, cos, pow   \n",
    "                \n",
    "# stuff I wrote\n",
    "from utilities import parallel_process, check_exe, system_call, compound_regions\n",
    "from model import beta_model, integ_beta_model, chi2, like, prior, prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Utility functions\n",
    "\n",
    "def pointInEllipse(xo, yo, xp, yp, d, D, angle):\n",
    "    #tests if a point[xp,yp] is within\n",
    "    #boundaries defined by the ellipse\n",
    "    #of center[x,y], diameter d D, and tilted at angle\n",
    "    \n",
    "    #### the angle must be in radians!!!! ####\n",
    "    \n",
    "    ## This is for the detect_vtp function #\n",
    "    \n",
    "    cosa = cos(angle)\n",
    "    sina = sin(angle)\n",
    "    \n",
    "    dd = d**2\n",
    "    DD = D**2\n",
    "\n",
    "    a = pow(cosa * (xp - xo) + sina * (yp - yo), 2)\n",
    "    b = pow(sina * (xp - xo) - cosa * (yp - yo), 2)\n",
    "    ellipse = (a / dd) + (b / DD)\n",
    "    \n",
    "    if ellipse <= 1:\n",
    "        return True\n",
    "    else:\n",
    "        return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def detect_vtp(name, outpath):\n",
    "    # set up all the non-file-specific parameters                              \n",
    "    params = {}\n",
    "    params['limit'] = 1E-6\n",
    "    params['coarse'] = 25\n",
    "    params['maxiter'] = 10\n",
    "    params['clobber'] = 'yes'\n",
    "    params['verbose'] = 1\n",
    "\n",
    "    evts = f'{outpath}/{name}/{name}_events.fits'\n",
    "    expmap = f'{outpath}/{name}/{name}_exp.fits'\n",
    "\n",
    "    # check to make sure the files exist.\n",
    "    if not os.path.isfile(evts) and not os.path.isfile(expmap):\n",
    "        return 0\n",
    "                                                                               \n",
    "    # We are going to run vtpdetect twice... once with a low scale and once with a high(er) scale\n",
    "    params['regfile'] = f'{outpath}/{name}/{name}_vtp_low.reg'\n",
    "    params['log'] = f'{outpath}/{name}/{name}_vtp_low.log'\n",
    "    params['scale'] = 1\n",
    "    outfits = f'{outpath}/{name}/{name}_vtp_low.detect'\n",
    "    with open(f'{outpath}/{name}/{name}_vtp_low.in', 'w') as f:\n",
    "        # build the cmd\n",
    "        cmd = f'vtpdetect {evts}[pi=50:600] {expmap} {outfits} '\n",
    "        for param, value in list(params.items()):\n",
    "            cmd += f'{param}={value} '\n",
    "        f.writelines(f'{cmd}\\n')\n",
    "\n",
    "    os.system(cmd)\n",
    "\n",
    "    # now for the higher level\n",
    "    params['regfile'] = f'{outpath}/{name}/{name}_vtp_high.reg'\n",
    "    params['log'] = f'{outpath}/{name}/{name}_vtp_high.log'\n",
    "    params['scale'] = 1.8\n",
    "    outfits = f'{outpath}/{name}/{name}_vtp_high.detect'\n",
    "    with open(f'{outpath}/{name}/{name}_vtp_high.in', 'w') as f:\n",
    "        # build the cmd\n",
    "        cmd = f'vtpdetect {evts}[pi=50:600] {expmap} {outfits} '\n",
    "        for param, value in list(params.items()):\n",
    "            cmd += f'{param}={value} '\n",
    "        f.writelines(f'{cmd}\\n')\n",
    "\n",
    "    os.system(cmd)\n",
    "\n",
    "    ###\n",
    "    # Now we are building the final catalog by comparing the low and high scale catalogs.\n",
    "    ###\n",
    "    try:\n",
    "        low = Table.read(f'{outpath}/{name}/{name}_vtp_low.detect', hdu=1)\n",
    "        high = Table.read(f'{outpath}/{name}/{name}_vtp_high.detect', hdu=1)\n",
    "    except FileNotFoundError:\n",
    "        return\n",
    "\n",
    "    # create a new table to store our results -- just makes an empty copy of the original\n",
    "    final = Table(dtype=low.dtype)\n",
    "    final.add_column(Column(name='INDEX', dtype='>i4'), index=0)\n",
    "    final.add_column(Column(name='HIGH', dtype='>i4'))\n",
    "                                                                               \n",
    "    # have to add columns to the original tables to make them match\n",
    "    low.add_column(Column(data=np.zeros(len(low)), name='INDEX', dtype='>i4'), index=0)\n",
    "    low.add_column(Column(data=np.zeros(len(low)), name='HIGH', dtype='>i4'))\n",
    "    high.add_column(Column(data=np.zeros(len(high)), name='INDEX', dtype='>i4'), index=0)\n",
    "    high.add_column(Column(data=np.zeros(len(high)), name='HIGH', dtype='>i4'))\n",
    "\n",
    "    index = 1 # vtp is normally 1 indexed\n",
    "    for x_l, y_l, rad_l, rot_l, comp_l in low[['X', 'Y', 'R', 'ROTANG', 'COMPONENT']]:\n",
    "\n",
    "        added_high = False # keep track if we found a high source\n",
    "\n",
    "        for x_h, y_h, comp_h in high[['X', 'Y', 'COMPONENT']]:\n",
    "            if pointInEllipse(x_l, y_l, x_h, y_h, rad_l[0], rad_l[1], np.radians(rot_l)):\n",
    "                final.add_row(high[comp_h - 1]) # add the row\n",
    "                final['INDEX'][index - 1] = index # add the index to the row\n",
    "                final['HIGH'][index - 1] = 1 # say that the source came from high catalog\n",
    "                index += 1\n",
    "                added_high = True\n",
    "\n",
    "        if not added_high:\n",
    "            final.add_row(low[comp_l - 1]) # add the row\n",
    "            final['INDEX'][index - 1] = index # add the index to the row\n",
    "            final['HIGH'][index - 1] = 0 # say that the source came from low catalog\n",
    "            index += 1\n",
    "\n",
    "    # no sources detected\n",
    "    if not len(final):\n",
    "        return final\n",
    "            \n",
    "    # for some reason there are duplicate sources\n",
    "    final = table.unique(final, keys=['X', 'Y'])\n",
    "    final['INDEX'] = [i + 1  for i in range(len(final))]\n",
    "\n",
    "    # remove any source where the center of the region lies inside another region\n",
    "    # gonna do two loops\n",
    "    remove_idx = []\n",
    "    for x1, y1, rad1, rot1, idx1, area1 in final[['X', 'Y', 'R', 'ROTANG', 'INDEX', 'SRC_AREA']]:\n",
    "        for x2, y2, rad2, rot2, idx2, area2 in final[['X', 'Y', 'R', 'ROTANG', 'INDEX', 'SRC_AREA']]:\n",
    "            if idx1 == idx2:\n",
    "                continue\n",
    "            elif any(rad2 < 1e-3): # some sources have radii = 0\n",
    "                remove_idx.append(idx2 - 1) # remember that our source index is 1 indexed\n",
    "            elif pointInEllipse(x1, y1, x2, y2, rad1[0], rad1[1], np.radians(rot1)):\n",
    "                remove_idx.append(idx2 - 1) # remember that our source index is 1 indexed\n",
    "    \n",
    "    # remove sources!\n",
    "    final.remove_rows(remove_idx)\n",
    "    \n",
    "    #reindex the table since we've removed sources\n",
    "    final['INDEX'] = [i + 1  for i in range(len(final))]\n",
    "    \n",
    "    # write detections!\n",
    "    final.write(f'{outpath}/{name}/{name}_vtp.detect', format='fits', overwrite=True)\n",
    "\n",
    "    # write out the regions\n",
    "    with open(f'{outpath}/{name}/{name}_vtp.reg', 'w') as reg:\n",
    "        reg.write(\"# Region file format: DS9 version 4.1\\n\")\n",
    "        reg.write('global color=cyan dashlist=8 3 width=1 font=\"helvetica 10 normal roman\" select=1 '\n",
    "                  'highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1\\n')\n",
    "        reg.write('fk5\\n')\n",
    "        for j, xc, yc, rc, rotc in final[['INDEX', 'RA', 'DEC', 'R', 'ROTANG']]:\n",
    "            reg.write(f'ellipse({xc},{yc},{(rc[0] * 2.36):.3f}\",{(rc[1] * 2.36):.3f}\",{rotc:.3f}) ')\n",
    "            reg.write(f'# text={{{j}}}\\n')\n",
    "\n",
    "    return final"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def centroid_ciao(name, outpath):\n",
    "    \n",
    "    # check to make sure we've started the ciaotools\n",
    "    if not check_exe('dmstat'):\n",
    "        print('dmstat not available. Have you started ciaotools?')\n",
    "        return\n",
    "    \n",
    "    # check for files\n",
    "    # events image\n",
    "    if not os.path.isfile(f'{outpath}/{name}/{name}_img_50-200.fits'):\n",
    "        return\n",
    "    else:\n",
    "        evnts = f'{outpath}/{name}/{name}_img_50-200.fits'\n",
    "    # expo map\n",
    "    if not os.path.isfile(f'{outpath}/{name}/{name}_exp.fits'):\n",
    "        return\n",
    "    else:\n",
    "        expmap = f'{outpath}/{name}/{name}_exp.fits'   \n",
    "    # detections\n",
    "    if not os.path.isfile(f'{outpath}/{name}/{name}_vtp.detect'):\n",
    "        return\n",
    "    else:\n",
    "        srcs = f'{outpath}/{name}/{name}_vtp.detect'\n",
    "    \n",
    "    # WCS info to convert pixels to world coords\n",
    "    hdr = fits.getheader(evnts)\n",
    "    w = wcs.WCS(hdr)\n",
    "    \n",
    "    # now we need to read the individual detections\n",
    "    detects = Table.read(srcs, hdu=1)\n",
    "    \n",
    "    # rename columns\n",
    "    try:\n",
    "        detects['X'].name = 'X_VTP'\n",
    "        detects['Y'].name = 'Y_VTP'      \n",
    "        # add the columns to the detection catalog\n",
    "        detects.add_column(Column(data=np.zeros(len(detects)), name='X', dtype='>f8'))\n",
    "        detects.add_column(Column(data=np.zeros(len(detects)), name='Y', dtype='>f8'))\n",
    "    except KeyError:\n",
    "        pass\n",
    "\n",
    "    try:\n",
    "        detects['X_ERR'].name = 'X_ERR_VTP'\n",
    "        detects['Y_ERR'].name = 'Y_ERR_VTP'      \n",
    "        # add the columns to the detection catalog\n",
    "        detects.add_column(Column(data=np.zeros(len(detects)), name='X_ERR', dtype='>f8'))\n",
    "        detects.add_column(Column(data=np.zeros(len(detects)), name='Y_ERR', dtype='>f8'))\n",
    "    except KeyError:\n",
    "        pass\n",
    "    \n",
    "    try:\n",
    "        detects['RA'].name = 'RA_VTP'\n",
    "        detects['DEC'].name = 'DEC_VTP'\n",
    "        detects.add_column(Column(data=np.zeros(len(detects)), name='RA', dtype='>f8'))\n",
    "        detects.add_column(Column(data=np.zeros(len(detects)), name='DEC', dtype='>f8'))\n",
    "    except KeyError:\n",
    "        pass\n",
    "    \n",
    "    try:\n",
    "        detects['RA_ERR'].name = 'RA_ERR_VTP'\n",
    "        detects['DEC_ERR'].name = 'DEC_ERR_VTP'\n",
    "        detects.add_column(Column(data=np.zeros(len(detects)), name='RA_ERR', dtype='>f8'))\n",
    "        detects.add_column(Column(data=np.zeros(len(detects)), name='DEC_ERR', dtype='>f8'))\n",
    "    except KeyError:\n",
    "        pass\n",
    "    \n",
    "    for i, (xc, yc, rc, rotc) in enumerate(detects[['X_VTP', 'Y_VTP', 'R', 'ROTANG']]):\n",
    "        fd, path = tempfile.mkstemp()\n",
    "        with os.fdopen(fd, 'w') as reg:\n",
    "            reg.write(\"# Region file format: CIAO version 1.0\\n\")\n",
    "            reg.write(f\"+ellipse({xc},{yc},{rc[0]},{rc[1]},{rotc})\\n\")\n",
    "    \n",
    "        cmd = f\"dmstat '{evnts}[(x,y)=region({path})]' centroid=yes clip=yes\"\n",
    "    \n",
    "        result, stderr = system_call(cmd, shlexify=False)\n",
    "#         print(result, stderr)\n",
    "        \n",
    "        center = result.split('cntrd[phys]')[1].split('\\n')[0].split(' ')[1:3]\n",
    "        errors = result.split('sigma_cntrd')[1].split('\\n')[0].split(' ')[1:3]\n",
    "        \n",
    "        # write the new coordinates to the table\n",
    "        detects['X'][i], detects['Y'][i] = float(center[0]), float(center[1])\n",
    "        detects['X_ERR'][i], detects['Y_ERR'][i] = float(errors[0]), float(errors[1])\n",
    "        \n",
    "        ### have to subtract 1 from the region positions, because FITS are 1-index'd and python\n",
    "        ### is 0-index'd.\n",
    "        world = w.wcs_pix2world(float(center[0]) - 1, float(center[1]) - 1, 0)\n",
    "        detects['RA'][i], detects['DEC'][i] = world[0], world[1]\n",
    "        detects['RA_ERR'][i], detects['DEC_ERR'][i] = (detects['X_ERR'][i] * 2.36 / 3600, \n",
    "                                                        detects['Y_ERR'][i] * 2.36 / 3600)\n",
    "        \n",
    "        os.remove(path)\n",
    "        \n",
    "    detects.write(f'{outpath}/{name}/{name}_vtp.detect', format='fits', overwrite=True) \n",
    "    \n",
    "    with open(f'{outpath}/{name}/{name}_vtp.reg', 'w') as reg:\n",
    "        reg.write(\"# Region file format: DS9 version 4.1\\n\")\n",
    "        reg.write('global color=cyan width=1 font=\"helvetica 10 normal roman\" \\n')\n",
    "        reg.write('fk5\\n')\n",
    "        for j, xc, yc, rc, rotc in detects[['INDEX', 'RA', 'DEC', 'R', 'ROTANG']]:\n",
    "            reg.write(f'ellipse({xc},{yc},{(rc[0] * 2.36):.3f}\",{(rc[1] * 2.36):.3f}\",{rotc:.3f}) ')\n",
    "            reg.write(f'# text={{{j}}}\\n')\n",
    "    \n",
    "    return detects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_profile(name, outpath, src_num=0, srcs=None):\n",
    "    # swift pixelscale in degrees\n",
    "    pixsc =  6.548089E-04 * 3600\n",
    "    \n",
    "    # detections\n",
    "    if not srcs:\n",
    "        srcs = f'{outpath}/{name}/{name}_vtp.detect'\n",
    "        # now we need to read the individual detections\n",
    "        detects = Table.read(srcs, hdu=1)\n",
    "    else:\n",
    "        detects = srcs\n",
    "       \n",
    "    # get the source index\n",
    "    i = detects['INDEX'][src_num]\n",
    "    \n",
    "    if os.path.isfile(f'{outpath}/{name}/{name}_vtp_{i}.tmpprof'):\n",
    "        data = Table.read(f'{outpath}/{name}/{name}_vtp_{i}.tmpprof', format='ascii', header_start=2)\n",
    "    else:\n",
    "        return\n",
    "\n",
    "    # x-axis, in arcminutes\n",
    "    x = (data['r1'] + data['r2'])/ 2. / 60. * pixsc\n",
    "\n",
    "    # this is the parameter fitting\n",
    "    ndim = 4  # number of parameters in the model\n",
    "    nwalkers = 100  # number of MCMC walkers\n",
    "    nburn = 100  # \"burn-in\" period to let chains stabilize\n",
    "    nsteps = 500  # number of MCMC steps to take\n",
    "\n",
    "    pos = [np.array([0.001, 1., 1., 0.001]) + 1e-4 * np.random.randn(ndim) for i in range(nwalkers)]\n",
    "    sampler = emcee.EnsembleSampler(nwalkers, ndim, prob, args=(x, data['sb'], data['sb_err']))\n",
    "    sampler.run_mcmc(pos, nsteps)\n",
    "    emcee_trace = sampler.get_chain(discard=nburn, thin=15, flat=True)\n",
    "\n",
    "    #fit -- median (and error) values\n",
    "    So, rc, beta, bg = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]),\n",
    "                         zip(*np.percentile(emcee_trace, [16, 50, 84], axis=0)))\n",
    "\n",
    "     # open a file to write out the results\n",
    "    with open(f'{outpath}/{name}/{name}_{i}_mcmcfits.tmp', 'w') as outfile:\n",
    "        outfile.write('# Field ID So_50 So_84 So_16 rc_50 rc_84 rc_16 '\n",
    "                      'beta_50 beta_84 beta_16 bg_50 bg_84 bg_16\\n')\n",
    "\n",
    "        # write out the fit parameters\n",
    "        outfile.write(f'{name} ')\n",
    "        outfile.write(f'{i} ')\n",
    "        outfile.write(f'{So[0]:.5f} {So[1]:.5f} {So[2]:.5f} ')\n",
    "        outfile.write(f'{rc[0]:.5f} {rc[1]:.5f} {rc[2]:.5f} ')\n",
    "        outfile.write(f'{beta[0]:.5f} {beta[1]:.5f} {beta[2]:.5f} ')\n",
    "        outfile.write(f'{bg[0]:.5f} {bg[1]:.5f} {bg[2]:.5f}\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_radprof(name, outpath, src_num=0, srcs=None):\n",
    "    pixscale = 6.548089E-04 * 3600 # arcsec/pixel\n",
    "\n",
    "    # model parameters\n",
    "    params = {}\n",
    "    params['n'] = 75\n",
    "    params['r0'] = 0 \n",
    "    params['dr'] = 4\n",
    "    \n",
    "    # check for files\n",
    "    # events image\n",
    "    if not os.path.isfile(f'{outpath}/{name}/{name}_img_50-200.fits'):\n",
    "        return\n",
    "    else:\n",
    "        evnts = f'{outpath}/{name}/{name}_img_50-200.fits'\n",
    "\n",
    "    # expo map\n",
    "    if not os.path.isfile(f'{outpath}/{name}/{name}_exp.fits'):\n",
    "        return\n",
    "    else:\n",
    "        expmap = f'{outpath}/{name}/{name}_exp.fits'   \n",
    "        \n",
    "    # detections\n",
    "    if not srcs:\n",
    "        srcs = f'{outpath}/{name}/{name}_vtp.detect'\n",
    "        # now we need to read the individual detections\n",
    "        detects = Table.read(srcs, hdu=1)\n",
    "    else:\n",
    "        detects = srcs\n",
    "        \n",
    "    i = detects['INDEX'][src_num]\n",
    "    xc = detects['X'][src_num]\n",
    "    yc = detects['Y'][src_num]  \n",
    "\n",
    "    # load the data\n",
    "    evnt_data = fits.getdata(evnts)\n",
    "    expmap_data = fits.getdata(expmap)\n",
    "    \n",
    "    # mask out all sources except the source we are working with\n",
    "    if len(detects) < 2:\n",
    "        # here there is only one source in the field. Don't mask anything.\n",
    "        mask_sources = False\n",
    "    else:\n",
    "        mask_sources = True\n",
    "        src_mask = [j for j in range(len(detects)) if j != src_num]\n",
    "\n",
    "    if mask_sources:\n",
    "        regs = []\n",
    "        ### have to subtract 1 from the region positions, because FITS are 1-index'd and python\n",
    "        ### is 0-index'd.\n",
    "        for idx, xc1, yc1, rc, rotc in detects[src_mask][['INDEX', 'X', 'Y', 'R', 'ROTANG']]:\n",
    "            ellipse_region = EllipsePixelRegion(center=PixCoord(x=xc1 -1, y=yc1 -1),\n",
    "                                                width=2 * rc[0], height=2 * rc[1], \n",
    "                                                angle=rotc * u.deg, \n",
    "                                                meta={'ID':idx})\n",
    "            regs.append(ellipse_region)   \n",
    "\n",
    "        # create the compound regions\n",
    "        cmp_reg = compound_regions(regs, evnt_data.shape[0])\n",
    "        # invert the comp region -- we want the places where there aren't regions\n",
    "        cmp_regmask = np.where(cmp_reg == 0, 1, 0)\n",
    "\n",
    "    else:\n",
    "        cmp_regmask = np.ones_like(evnt_data)    \n",
    "    \n",
    "    # output file  \n",
    "    with open(f'{outpath}/{name}/{name}_vtp_{i}.tmpprof', 'w') as radprof:\n",
    "        radprof.write(f\"# source number and position: {name}_{i} {xc:.3f} {yc:.3f}\\n\")\n",
    "        radprof.write(f\"# profile parameters: {params['n']} {params['r0']} {params['dr']}\\n\")\n",
    "        radprof.write(f\"# bin r1 r2 x y w sb sb_err psf area psfsb\\n\")\n",
    "            \n",
    "        for irad in range(params['n']):\n",
    "            r1 = params['r0'] + irad * params['dr']\n",
    "            r2 = params['r0'] + (irad + 1) * params['dr']\n",
    "\n",
    "            center = PixCoord(x=detects['X'][i] - 1, y=detects['Y'][i] - 1)\n",
    "            if r1 == 0:\n",
    "                ann_reg = CirclePixelRegion(center=center, radius=r2)\n",
    "            else:\n",
    "                ann_reg = CircleAnnulusPixelRegion(center=center,\n",
    "                                                   inner_radius=r1,\n",
    "                                                   outer_radius=r2)\n",
    "            # make it into an image\n",
    "            ann_regmask = compound_regions([ann_reg], evnt_data.shape[0])\n",
    "\n",
    "            # combine the regions with the annulus.\n",
    "            final_regmask = (ann_regmask * cmp_regmask)\n",
    "            final_regmask_inv = np.where(\n",
    "                final_regmask == 0, 1, 0)  # this is the final mask for the data\n",
    "\n",
    "            evnt_data_masked = ma.masked_array(evnt_data,\n",
    "                                               mask=final_regmask_inv)\n",
    "            expmap_data_masked = ma.masked_array(expmap_data,\n",
    "                                                 mask=final_regmask_inv)\n",
    "\n",
    "            x = evnt_data_masked.sum()\n",
    "            y = expmap_data_masked.mean()\n",
    "            w = expmap_data_masked.count()\n",
    "\n",
    "            sb = x / (y * w * pixscale**2 / 3600)  #  cts/s/arcmin^2\n",
    "            err_gehrels = sqrt(x + 0.75)\n",
    "            sbe = (1. + err_gehrels) / (y * w * pixscale**2 / 3600\n",
    "                                        )  #  cts/s/arcmin^2\n",
    "\n",
    "            # build the model of the psf\n",
    "            psfr1 = integ_beta_model(r1, rc=(5.6 / 2.36), beta=-1/3)  # cnts\n",
    "            psfr2 = integ_beta_model(r2, rc=(5.6 / 2.36), beta=-1/3)  # cnts\n",
    "            psf = psfr2 - psfr1  # total flux in the annulus\n",
    "\n",
    "            area = np.pi * r2**2 - np.pi * r1**2  # number of square pixels\n",
    "\n",
    "            psfsb = psf / (y * area * pixscale**2 / 3600)\n",
    "\n",
    "            radprof.write(\n",
    "                f'{irad:3d} {r1:7d} {r2:7d} {x:9.1f} {y:9.1f} {w:9.1f}'\n",
    "                f'{sb:12.4e} {sbe:12.4e} '\n",
    "                f'{psf:9.1f} {area:9.1f} {psfsb:12.4e} \\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clean_sources(name, outpath, PSZ_RA, PSZ_DEC):\n",
    "\n",
    "    if not os.path.isfile(f'{outpath}/{name}/{name}_vtp.detect'):\n",
    "        return\n",
    "    else:\n",
    "        srcs = Table.read(f'{outpath}/{name}/{name}_vtp.detect', hdu=1)\n",
    "    \n",
    "    # if there is only a single source, then we don't need to clean anything\n",
    "    if len(srcs) > 1:\n",
    "        clean = True\n",
    "    else:\n",
    "        clean = False\n",
    "    \n",
    "    itter = 1\n",
    "    while clean:\n",
    "        src = np.argmax(srcs['NET_COUNTS'])\n",
    "#         print(name, itter, len(srcs), f'brightest: {src}')\n",
    "        get_radprof(name, outpath, src, srcs)\n",
    "        fit_profile(name, outpath, src, srcs)\n",
    "        \n",
    "        fit = Table.read(f'{outpath}/{name}/{name}_{srcs[\"INDEX\"][src]}_mcmcfits.tmp', \n",
    "                         format='ascii', header_start=0)\n",
    "        \n",
    "        weak_srcs = []\n",
    "        for i in range(len(srcs)):\n",
    "            if i == src:\n",
    "                continue\n",
    "\n",
    "            clust = SkyCoord(srcs['RA'][src], srcs['DEC'][src], unit='deg', frame='icrs')\n",
    "            src1 = SkyCoord(srcs['RA'][i], srcs['DEC'][i], unit='deg', frame='icrs')\n",
    "\n",
    "            sep = clust.separation(src1)\n",
    "            r = sep.arcmin\n",
    "            \n",
    "            flux = beta_model(fit['So_50'], r, fit['rc_50'], fit['beta_50']) + fit['bg_50']\n",
    "            flux = flux.data[0]\n",
    "\n",
    "            cnts = srcs['NET_COUNTS'][i] + srcs['BKG_COUNTS'][i]\n",
    "            exptime = srcs['EXPTIME'][i]\n",
    "            src_area = srcs['SRC_AREA'][i] * 2.36**2/60**2\n",
    "\n",
    "            src_flux = cnts/exptime/src_area\n",
    "            src_flux_err = sqrt(cnts)/exptime/src_area\n",
    "\n",
    "            sigma = (src_flux - flux)/src_flux_err\n",
    "            \n",
    "            if sigma < 5:\n",
    "                weak_srcs.append(i)\n",
    "            #print(srcs['INDEX'][i], sigma)\n",
    "\n",
    "        # print('bad sources', len(weak_srcs))\n",
    "        if not len(weak_srcs):\n",
    "            # remove tmp files.\n",
    "            # print('remove files', srcs[\"INDEX\"][src])\n",
    "            os.remove(f'{outpath}/{name}/{name}_vtp_{srcs[\"INDEX\"][src]}.tmpprof')\n",
    "            os.remove(f'{outpath}/{name}/{name}_{srcs[\"INDEX\"][src]}_mcmcfits.tmp')  \n",
    "            break\n",
    "        else:\n",
    "            srcs.remove_rows(weak_srcs)\n",
    "\n",
    "        itter += 1 \n",
    "        \n",
    "    srcs['INDEX'] = [i + 1  for i in range(len(srcs))]\n",
    "                      \n",
    "    # compute the distance to the PSZ position!\n",
    "    # add some columns to our detect file\n",
    "    try:\n",
    "        srcs.add_column(Column(data=-np.ones(len(srcs)), name='PSZ_dist', dtype='>f8'))\n",
    "    except ValueError:\n",
    "        pass\n",
    "    \n",
    "    for i, (ra, dec) in enumerate(srcs[['RA', 'DEC']]):\n",
    "        PSZ = SkyCoord(PSZ_RA, PSZ_DEC, unit='deg', frame='icrs')\n",
    "        src1 = SkyCoord(ra, dec, unit='deg', frame='icrs')\n",
    "\n",
    "        sep = PSZ.separation(src1)\n",
    "        srcs['PSZ_dist'][i] = sep.arcmin              \n",
    "                      \n",
    "    srcs.write(f'{outpath}/{name}/{name}_vtp.detect', format='fits', overwrite=True)    \n",
    "\n",
    "    with open(f'{outpath}/{name}/{name}_vtp.reg', 'w') as reg:\n",
    "        reg.write(\"# Region file format: DS9 version 4.1\\n\")\n",
    "        reg.write('global color=cyan width=1 font=\"helvetica 10 normal roman\" \\n')\n",
    "        reg.write('fk5\\n')\n",
    "        for j, xc, yc, rc, rotc in srcs[['INDEX', 'RA', 'DEC', 'R', 'ROTANG']]:\n",
    "            reg.write(f'ellipse({xc},{yc},{(rc[0] * 2.36):.3f}\",{(rc[1] * 2.36):.3f}\",{rotc:.3f}) ')\n",
    "            reg.write(f'# text={{{j}}}\\n')\n",
    "            \n",
    "    return srcs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# get file data\n",
    "data = load_PSZcatalog()\n",
    "data = data.sort_values('NAME')\n",
    "\n",
    "outpath = './data_full_new'\n",
    "\n",
    "arr = [{'name':n.replace(' ', '_'), 'outpath':outpath} for n in data['NAME']]\n",
    "\n",
    "###\n",
    "# This is the order you should call the functions. \n",
    "# There are other functions in this notebook, but these are the only ones \n",
    "# that should be called directly.\n",
    "###\n",
    "\n",
    "parallel_process(arr, detect_vtp, use_kwargs=True, n_jobs=1)\n",
    "parallel_process(arr, centroid_ciao, use_kwargs=True, n_jobs=1)\n",
    "\n",
    "arr = [{'name':n.replace(' ', '_'), \n",
    "        'outpath':outpath, \n",
    "        'PSZ_RA': ra,\n",
    "        'PSZ_DEC': dec} for n, ra, dec in zip(data['NAME'], data['RA'], data['DEC'])]\n",
    "parallel_process(arr, clean_sources, use_kwargs=True, n_jobs=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "### For testing ###\n",
    "\n",
    "outpath = './data_full_new'\n",
    "name = 'PSZ1_G031.41+28.75'\n",
    "RA = 256.263152985\n",
    "DEC = 11.4604480691\n",
    "detect_vtp(name, outpath)\n",
    "centroid_ciao(name, outpath)\n",
    "clean_sources(name, outpath, RA, DEC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### For testing ###\n",
    "\n",
    "outpath = './data_full_new'\n",
    "name = 'PSZ2_G313.33+61.13'\n",
    "RA = 197.86033\n",
    "DEC = -1.336351\n",
    "# detect_vtp(name, outpath)\n",
    "# centroid_ciao(name, outpath)\n",
    "clean_sources(name, outpath, RA, DEC)"
   ]
  },
  {
   "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.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
}