segment_info.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Checking and playing with NICER data\n",
    "Looking through and plotting segment information from quicklook_segments.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy.table import Table, Column, unique, vstack\n",
    "from astropy.time import Time\n",
    "import os\n",
    "# import subprocess\n",
    "from astropy.io import fits\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.font_manager as font_manager\n",
    "from matplotlib.ticker import MultipleLocator, FixedLocator\n",
    "from matplotlib.ticker import ScalarFormatter, NullFormatter\n",
    "import matplotlib.gridspec as gridspec\n",
    "from matplotlib.colors import cnames as mcolors\n",
    "from matplotlib._color_data import XKCD_COLORS as xkcdcolor\n",
    "from scipy.stats import binned_statistic\n",
    "# import matplotlib.patches as patches\n",
    "# from matplotlib.collections import PatchCollection\n",
    "from xcor_tools_nicer import find_nearest as find_nearest\n",
    "from xcor_tools_nicer import clock_to_mjd as clock_to_mjd\n",
    "font_prop = font_manager.FontProperties(size=20)\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "homedir = os.path.expanduser(\"~\")\n",
    "obj_name = \"Swift_J1728.9-3613\"\n",
    "obj_prefix = \"SwiftJ1728\"\n",
    "obj_dir = \"%s/Documents/Research/%s\" % (homedir, obj_prefix)\n",
    "out_dir = \"%s/out\" % (obj_dir)\n",
    "# obj_dir = \"%s/Documents/Research/NICER_exploration\" % (homedir)\n",
    "# out_dir = \"%s/out/%s\" % (obj_dir, obj_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## This seg_info_file was made in analyze_all_segments.py\n",
    "seg_list_file = \"%s/%s_seg-info-list.txt\" % (out_dir, obj_prefix)\n",
    "seg_list = [line.strip() for line in open(seg_list_file)]\n",
    "seg_list = seg_list[1:]\n",
    "seg_tab = Table()\n",
    "for seg_info_file in seg_list:\n",
    "    temp_tab = Table.read(seg_info_file, format='ascii', \n",
    "                     names=('obsID', 'start_time', 'end_time', \n",
    "                            'total_rate', 'broad_rate', 'rms', \n",
    "                            'hard_rate', 'soft_rate', 'hr', 'fe62', \n",
    "                            'fe63', 'fe64', 'fe65', 'fe66'))\n",
    "    seg_tab = vstack([seg_tab,temp_tab])\n",
    "#     print(seg_tab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Unique obsIDs in this table\n",
    "obsIDs = unique(seg_tab, keys='obsID')['obsID']\n",
    "print(\"%d obsIDs\" % len(obsIDs))\n",
    "print(obsIDs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "times = seg_tab['start_time'] - seg_tab['start_time'][0]\n",
    "times /= (60*60*24)  # now it's in days? elapsed\n",
    "nanrms_mask = seg_tab['rms'] < 90.\n",
    "seg_tab = seg_tab[nanrms_mask]\n",
    "times = times[nanrms_mask]\n",
    "mjdstart = Time(clock_to_mjd(seg_tab['start_time'][0]), format='mjd')\n",
    "isostart = mjdstart.to_value('iso', 'date_hm')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "binning_bins = np.arange(0,np.ceil(times[-1]),step=0.1)\n",
    "bin_rate, bin_days, temp = binned_statistic(times, seg_tab['broad_rate'], statistic='mean', bins=binning_bins)\n",
    "bin_hr, bin_days, temp = binned_statistic(times, seg_tab['hr'], statistic='mean', bins=binning_bins)\n",
    "bin_rms, bin_days, temp = binned_statistic(times, seg_tab['rms'], statistic='mean', bins=binning_bins)\n",
    "bin_fe62, bin_days, temp = binned_statistic(times, seg_tab['fe62'], statistic='mean', bins=binning_bins)\n",
    "bin_fe63, bin_days, temp = binned_statistic(times, seg_tab['fe63'], statistic='mean', bins=binning_bins)\n",
    "bin_fe64, bin_days, temp = binned_statistic(times, seg_tab['fe64'], statistic='mean', bins=binning_bins)\n",
    "bin_fe65, bin_days, temp = binned_statistic(times, seg_tab['fe65'], statistic='mean', bins=binning_bins)\n",
    "bin_fe66, bin_days, temp = binned_statistic(times, seg_tab['fe66'], statistic='mean', bins=binning_bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_fe = np.average([bin_fe62, bin_fe63, bin_fe64, bin_fe65, bin_fe66], axis=0)\n",
    "print(np.shape(bin_fe))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Rate: \", np.nanmin(bin_rate), np.nanmax(bin_rate))\n",
    "print(\"HR: \", np.nanmin(bin_hr), np.nanmax(bin_hr))\n",
    "print(\"rms: \", np.nanmin(bin_rms), np.nanmax(bin_rms))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(16, 16), dpi=300)\n",
    "gs1 = gridspec.GridSpec(4, 1)\n",
    "gs1.update(wspace=0.025, hspace=0.0)\n",
    "\n",
    "ax1 = plt.subplot(gs1[0])\n",
    "ax1.scatter(bin_days[:-1], bin_rate, marker='.', s=80, color=xkcdcolor['xkcd:cherry red'])\n",
    "ax1.set_ylabel(\"Count rate (2-12 keV)\", fontproperties=font_prop)\n",
    "# ax1.set_ylim(20000, 67000)\n",
    "ax1.tick_params(axis='x', labelsize=20, bottom=True, top=True,\n",
    "                labelbottom=False, labeltop=False, direction=\"in\")\n",
    "ax1.tick_params(axis='y', labelsize=20, left=True, right=True,\n",
    "                labelleft=True, labelright=False, direction=\"in\")\n",
    "ax1.tick_params(which='major', width=1.5, length=9, direction=\"in\")\n",
    "ax1.tick_params(which='minor', width=1.5, length=6, direction=\"in\")\n",
    "ax1.xaxis.set_major_formatter(ScalarFormatter())\n",
    "ax1.xaxis.set_minor_formatter(NullFormatter())\n",
    "ax1.xaxis.set_ticks_position('both')\n",
    "ax1.yaxis.set_minor_locator(MultipleLocator(100))\n",
    "ax1.yaxis.set_ticks_position('both')\n",
    "for axis in ['top', 'bottom', 'left', 'right']:\n",
    "    ax1.spines[axis].set_linewidth(1.5)\n",
    "ax1.xaxis.grid(True, which='major')\n",
    "ax1.yaxis.grid(True, which='major')\n",
    "ax1.set_axisbelow(True)\n",
    "\n",
    "ax2 = plt.subplot(gs1[1], sharex=ax1)\n",
    "ax2.scatter(bin_days[:-1], bin_hr, marker='.', s=80, color=xkcdcolor['xkcd:medium green'])\n",
    "ax2.set_ylabel(\"Hardness (4-12/1-3)\", fontproperties=font_prop)\n",
    "# ax2.set_ylim(0.27, 0.85)\n",
    "ax2.tick_params(axis='x', labelsize=20, bottom=True, top=True,\n",
    "                labelbottom=False, labeltop=False, direction=\"in\")\n",
    "ax2.tick_params(axis='y', labelsize=20, left=True, right=True,\n",
    "                labelleft=True, labelright=False, direction=\"in\")\n",
    "ax2.tick_params(which='major', width=1.5, length=9, direction=\"in\")\n",
    "ax2.tick_params(which='minor', width=1.5, length=6, direction=\"in\")\n",
    "ax2.xaxis.set_major_formatter(ScalarFormatter())\n",
    "ax2.xaxis.set_minor_formatter(NullFormatter())\n",
    "ax2.xaxis.set_ticks_position('both')\n",
    "ax2.yaxis.set_minor_locator(MultipleLocator(0.01))\n",
    "ax2.yaxis.set_ticks_position('both')\n",
    "for axis in ['top', 'bottom', 'left', 'right']:\n",
    "    ax2.spines[axis].set_linewidth(1.5)\n",
    "ax2.xaxis.grid(True, which='major')\n",
    "ax2.yaxis.grid(True, which='major')\n",
    "ax2.set_axisbelow(True)\n",
    "\n",
    "ax4 = plt.subplot(gs1[2], sharex=ax1)\n",
    "ax4.scatter(bin_days[:-1], bin_fe, marker='.', s=80, color=xkcdcolor['xkcd:orange'])\n",
    "# ax4.scatter(bin_days[:-1], bin_fe63, marker='.', s=80, color=xkcdcolor['xkcd:orange'])\n",
    "# ax4.scatter(bin_days[:-1], bin_fe64, marker='.', s=80, color=xkcdcolor['xkcd:green'])\n",
    "# ax4.scatter(bin_days[:-1], bin_fe65, marker='.', s=80, color=xkcdcolor['xkcd:blue'])\n",
    "# ax4.scatter(bin_days[:-1], bin_fe66, marker='.', s=80, color=xkcdcolor['xkcd:purple'])\n",
    "ax4.set_ylim(0.8, 1.4)\n",
    "ax4.set_ylabel(\"Iron line flux\", fontproperties=font_prop)\n",
    "ax4.tick_params(axis='x', labelsize=20, bottom=True, top=True,\n",
    "                labelbottom=False, labeltop=False, direction=\"in\")\n",
    "ax4.tick_params(axis='y', labelsize=20, left=True, right=True,\n",
    "                labelleft=True, labelright=False, direction=\"in\")\n",
    "ax4.tick_params(which='major', width=1.5, length=9, direction=\"in\")\n",
    "ax4.tick_params(which='minor', width=1.5, length=6, direction=\"in\")\n",
    "ax4.xaxis.set_ticks_position('both')\n",
    "ax4.yaxis.set_minor_locator(MultipleLocator(0.025))\n",
    "ax4.yaxis.set_major_locator(MultipleLocator(0.1))\n",
    "ax4.yaxis.set_ticks_position('both')\n",
    "for axis in ['top', 'bottom', 'left', 'right']:\n",
    "    ax4.spines[axis].set_linewidth(1.5)\n",
    "ax4.xaxis.grid(True, which='major')\n",
    "ax4.yaxis.grid(True, which='major')\n",
    "ax4.set_axisbelow(True)\n",
    "\n",
    "ax3 = plt.subplot(gs1[3], sharex=ax1)\n",
    "ax3.scatter(bin_days[:-1], bin_rms, marker='.', s=80, color=xkcdcolor['xkcd:ultramarine'])\n",
    "ax3.set_xlabel(\"Days since %s (UTC)\" % isostart, fontproperties=font_prop)\n",
    "ax3.set_ylabel(\"rms (1-15 Hz)\", fontproperties=font_prop)\n",
    "# ax3.set_ylim(0.00, 0.16)\n",
    "ax3.set_xlim(0,np.ceil(times[-1]))\n",
    "# ax3.xaxis.set_minor_locator(MultipleLocator(1))\n",
    "# ax3.set_xlim(0, 8e5)\n",
    "ax3.xaxis.set_major_locator(MultipleLocator(5))\n",
    "ax3.xaxis.set_minor_locator(MultipleLocator(1))\n",
    "ax3.xaxis.set_major_formatter(ScalarFormatter())\n",
    "ax3.tick_params(axis='x', labelsize=20, bottom=True, top=True,\n",
    "                labelbottom=True, labeltop=False, direction=\"in\")\n",
    "ax3.tick_params(axis='y', labelsize=20, left=True, right=True,\n",
    "                labelleft=True, labelright=False, direction=\"in\")\n",
    "ax3.tick_params(which='major', width=1.5, length=9, direction=\"in\")\n",
    "ax3.tick_params(which='minor', width=1.5, length=6, direction=\"in\")\n",
    "ax3.xaxis.set_minor_formatter(NullFormatter())\n",
    "ax3.xaxis.set_ticks_position('both')\n",
    "ax3.yaxis.set_minor_locator(MultipleLocator(0.01))\n",
    "ax3.yaxis.set_major_locator(MultipleLocator(0.05))\n",
    "ax3.yaxis.set_ticks_position('both')\n",
    "for axis in ['top', 'bottom', 'left', 'right']:\n",
    "    ax3.spines[axis].set_linewidth(1.5)\n",
    "ax3.xaxis.grid(True, which='major')\n",
    "ax3.yaxis.grid(True, which='major')\n",
    "ax3.set_axisbelow(True)\n",
    "\n",
    "\n",
    "plt.savefig(\"%s/%s_outburst_lcs.pdf\" % (out_dir, obj_prefix))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "even_days = [i for i in range(0,int(np.ceil(times[-1])),2)]\n",
    "# print(even_days)\n",
    "even_idx = []\n",
    "for the_day in even_days:\n",
    "    temp, the_idx = find_nearest(bin_days, the_day)\n",
    "#     print(the_day, temp, the_idx)\n",
    "    even_idx.append(the_idx)\n",
    "# print(the_idx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(9, 6.75), dpi=300, tight_layout=True)\n",
    "the_cbar = ax.scatter(bin_rms, bin_hr, marker='o', s=80, c=np.arange(len(bin_hr)), cmap='magma', edgecolor='black')\n",
    "# ax.set_xlim(0.105, 0.20)\n",
    "# ax.set_ylim(0.27, 0.85)\n",
    "ax.set_xlabel(r\"rms (1$-$15$\\,$Hz)\", fontproperties=font_prop)\n",
    "ax.set_ylabel(r\"Hardness (4$-$12$\\,$keV$\\,/\\,$1$-$3$\\,$keV)\", fontproperties=font_prop)\n",
    "ax.tick_params(axis='x', labelsize=20, bottom=True, top=True,\n",
    "                labelbottom=True, labeltop=False, direction=\"in\")\n",
    "ax.tick_params(axis='y', labelsize=20, left=True, right=True,\n",
    "                labelleft=True, labelright=False, direction=\"in\")\n",
    "ax.tick_params(which='major', width=1.5, length=9, direction=\"in\")\n",
    "ax.tick_params(which='minor', width=1.5, length=6, direction=\"in\")\n",
    "ax.xaxis.set_minor_locator(MultipleLocator(0.01))\n",
    "ax.yaxis.set_minor_locator(MultipleLocator(0.01))\n",
    "ax.xaxis.set_major_formatter(ScalarFormatter())\n",
    "ax.xaxis.set_minor_formatter(NullFormatter())\n",
    "ax.xaxis.set_ticks_position('both')\n",
    "ax.yaxis.set_ticks_position('both')\n",
    "for axis in ['top', 'bottom', 'left', 'right']:\n",
    "    ax.spines[axis].set_linewidth(1.5)\n",
    "ax.xaxis.grid(True, which='major')\n",
    "ax.yaxis.grid(True, which='major')\n",
    "ax.set_axisbelow(True)\n",
    "cbar=plt.colorbar(the_cbar, label=\"Days since %s (UTC)\" % isostart, orientation=\"horizontal\")\n",
    "cbar.set_ticks(even_idx)\n",
    "cbar.set_ticklabels(even_days)\n",
    "plt.savefig(\"%s/%s_rms-hr.pdf\" % (out_dir, obj_prefix))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(9, 6.75), dpi=300, tight_layout=True)\n",
    "the_cbar = ax.scatter(bin_hr, bin_rate, marker='o', s=80, c=np.arange(len(bin_hr)), cmap='magma', edgecolor='black')\n",
    "ax.set_xlabel(r\"Hardness (4$-$12$\\,$keV$\\,/\\,$1$-$3$\\,$keV)\", fontproperties=font_prop)\n",
    "ax.set_ylabel(r\"Counts/s (2$-$12$\\,$keV)\", fontproperties=font_prop)\n",
    "ax.tick_params(axis='x', labelsize=20, bottom=True, top=True,\n",
    "                labelbottom=True, labeltop=False, direction=\"in\")\n",
    "ax.set_yscale('log')\n",
    "ax.set_xscale('log')\n",
    "ax.tick_params(axis='y', labelsize=20, left=True, right=True,\n",
    "                labelleft=True, labelright=False, direction=\"in\")\n",
    "ax.tick_params(which='major', width=1.5, length=9, direction=\"in\")\n",
    "ax.tick_params(which='minor', width=1.5, length=6, direction=\"in\")\n",
    "# ax.set_xlim(0.026, 0.14)\n",
    "# ax.xaxis.set_major_locator(MultipleLocator(0.02))\n",
    "# ax.xaxis.set_minor_locator(MultipleLocator(0.01))\n",
    "# ax.xaxis.set_major_formatter(ScalarFormatter())\n",
    "# ax.xaxis.set_minor_formatter(NullFormatter())\n",
    "# ax.yaxis.set_minor_formatter(NullFormatter())\n",
    "# ax.yaxis.set_major_locator(FixedLocator([300, 500, 1000, 1500]))\n",
    "# ax.set_yticklabels([\"300\", \"500\", \"1000\", \"1500\"])\n",
    "ax.xaxis.set_ticks_position('both')\n",
    "ax.yaxis.set_ticks_position('both')\n",
    "for axis in ['top', 'bottom', 'left', 'right']:\n",
    "    ax.spines[axis].set_linewidth(1.5)\n",
    "ax.xaxis.grid(True, which='major')\n",
    "ax.yaxis.grid(True, which='major')\n",
    "ax.set_axisbelow(True)\n",
    "cbar=plt.colorbar(the_cbar, label=\"Days since %s (UTC)\" % isostart, orientation=\"horizontal\")\n",
    "cbar.set_ticks(even_idx)\n",
    "cbar.set_ticklabels(even_days)\n",
    "plt.savefig(\"%s/%s_hid.pdf\" % (out_dir, obj_prefix))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ## Filtering seg_tab\n",
    "# time_mask = seg_tab['start_time'] <= 119138417.0\n",
    "# filt_tab = seg_tab[time_mask]\n",
    "# # print(filt_tab)\n",
    "# hr_mask = (filt_tab['hr'] < 0.031) & (filt_tab['hr'] > 0.021)\n",
    "# # hr_mask = (filt_tab['hr'] > 0.031) & (filt_tab['hr'] < 0.05)\n",
    "# filt_tab = filt_tab[hr_mask]\n",
    "# # print(filt_tab)\n",
    "# rms_mask = filt_tab['rms'] < 0.041\n",
    "# filt_tab = filt_tab[rms_mask]\n",
    "# # print(filt_tab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ## These ones have the criteria I want from the cell above\n",
    "# good_obsIDs = unique(filt_tab, keys='obsID')['obsID']\n",
    "# print(\"%d good obsIDs\" % len(good_obsIDs))\n",
    "# print(good_obsIDs)"
   ]
  },
  {
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}