{
"cells": [
{
"cell_type": "markdown",
"id": "german-smith",
"metadata": {},
"source": [
"**NICER Tutorial**\n",
"\n",
"**By David A. Espinoza-Galeas**\n",
"\n",
"- This is a very introductory NICER tutorial with the minimum tasks you will need to start to work with NICER observations in SciServer.\n",
"- The tutorial is intended for people that is not familiar with NICER, maybe students or researchers trying to give a try to NICER observations without getting too deep in details.\n",
"- To get started in SciServer please go to the nice tutorial at the HEASARC: https://heasarc.gsfc.nasa.gov/docs/sciserver/\n",
"- Any command line with a \"!\" before is a line that you can execute in the regular terminal at SciServer.\n",
"- Feel free to share if you want.\n",
"- If you have question or comments related with this tutorial write to 84espinozagaleas@cua.edu .\n",
"\n",
"**I hope you enjoy it!**"
]
},
{
"cell_type": "markdown",
"id": "continent-accountability",
"metadata": {},
"source": [
"# Importing Modules"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "tested-optimum",
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"import sys\n",
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"import xspec\n",
"\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"\n",
"from astropy.io import fits\n",
"import astropy.time as Time\n",
"import astropy.units as u"
]
},
{
"cell_type": "markdown",
"id": "single-penny",
"metadata": {},
"source": [
"# Installing background estimators\n",
"\n",
"- Most of this info is in the NICER background webpage at the HEASARC: https://heasarc.gsfc.nasa.gov/docs/nicer/tools/nicer_bkg_est_tools.html\n",
"\n",
"## Space Weather background estimator (SW)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "parallel-helicopter",
"metadata": {},
"outputs": [],
"source": [
"# Download the nicer_bkg_estimator module\n",
"!wget https://heasarc.gsfc.nasa.gov/docs/nicer/tools/nicer_bkg_estimator_v0p6.tar.gz"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "rubber-inclusion",
"metadata": {},
"outputs": [],
"source": [
"# Untar the file\n",
"!tar zxvf nicer_bkg_estimator_v0p6.tar.gz"
]
},
{
"cell_type": "markdown",
"id": "complimentary-leone",
"metadata": {},
"source": [
"- You can leave the folder in the same folder of your Jupyter
notebook or in your ${\\tt \\$PYTHONPATH}$ (as documentation suggests)\n",
"\n",
"## NICER background generator (3C50)\n",
"\n",
"**Step 1: Get the reference model file**"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "suspected-installation",
"metadata": {},
"outputs": [],
"source": [
"# Download Model reference\n",
"!wget https://heasarc.gsfc.nasa.gov/docs/nicer/tools/bg_models_3C50.tar.gz"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "expired-constitution",
"metadata": {},
"outputs": [],
"source": [
"# Untar files\n",
"!tar zxfv bg_models_3C50.tar.gz"
]
},
{
"cell_type": "markdown",
"id": "eligible-scanner",
"metadata": {},
"source": [
"- You can put the folder any place, just make sure you remember it.\n",
"\n",
"**Step 2: Get and install the software**"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "interior-pioneer",
"metadata": {},
"outputs": [],
"source": [
"# Where is heasoft nicer folder\n",
"!ls $HEADAS/nicer"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "smart-finder",
"metadata": {},
"outputs": [],
"source": [
"# If there is not \"tasks\" folder create one in the right place\n",
"!mkdir $HEADAS/nicer/x86_64-pc-linux-gnu-libc2.27/tasks"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "concerned-regulation",
"metadata": {},
"outputs": [],
"source": [
"# We are going to move, but first save your working directory\n",
"work_dir = !pwd\n",
"\n",
"# Move to the folder using the \"os\" module\n",
"os.chdir('/opt/heasoft/nicer/x86_64-pc-linux-gnu-libc2.27/tasks')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "shared-burns",
"metadata": {},
"outputs": [],
"source": [
"#Download the software\n",
"!wget https://heasarc.gsfc.nasa.gov/docs/nicer/tools/nibackgen3C50_v6f.tar.gz"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "occupational-accident",
"metadata": {},
"outputs": [],
"source": [
"#untar files\n",
"!tar zxvf $HEADAS/nicer/x86_64-pc-linux-gnu-libc2.27/tasks/nibackgen3C50_v6f.tar.gz"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "genuine-amazon",
"metadata": {},
"outputs": [],
"source": [
"# Move to the nibackgen3C50 folder using the \"os\" module\n",
"os.chdir('/opt/heasoft/nicer/x86_64-pc-linux-gnu-libc2.27/tasks/nibackgen3C50/')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "conceptual-carpet",
"metadata": {},
"outputs": [],
"source": [
"# Run the hmake\n",
"!hmake &> hmake.log"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "original-discipline",
"metadata": {},
"outputs": [],
"source": [
"# Install\n",
"!hmake install &> hmake_install.log"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "severe-lightweight",
"metadata": {},
"outputs": [],
"source": [
"# Comeback to your working directory \n",
"os.chdir(work_dir[0])"
]
},
{
"cell_type": "markdown",
"id": "direct-update",
"metadata": {},
"source": [
"# Download Observations\n",
"\n",
"- If you have the OBS ID and the month where the observations was allocated you can just copy the obervation folder in to your working directory."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "stuffed-costume",
"metadata": {},
"outputs": [],
"source": [
"!cp -a /home/idies/workspace/headata/FTP/nicer/data/obs/2019_05//2612011001/ ./"
]
},
{
"cell_type": "markdown",
"id": "hairy-moisture",
"metadata": {},
"source": [
"# Procesing NICER observation\n",
"\n",
"- First apply the recommended calibration process.\n",
"\n",
"- If you are working in a ${\\tt Jupyter}$ notebook I highly encourage you to dump the termial output into a log file.\n",
"\n",
"- Give it some time, can take a while."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "circular-layer",
"metadata": {},
"outputs": [],
"source": [
"!nicerl2 indir=\"./2612011001\" clobber=\"yes\" cor_range=\"1-15\" &> nicerl2.log"
]
},
{
"cell_type": "markdown",
"id": "organizational-replication",
"metadata": {},
"source": [
"- Create a ${\\tt mkf2}$ since you are going to need it to apply SW background estimator"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "organizational-forward",
"metadata": {},
"outputs": [],
"source": [
"!niprefilter2 indir=\"./2612011001\" infile=\"./2612011001/auxil/ni2612011001.mkf\" \\\n",
"outfile=\"./2612011001/auxil/ni2612011001.mkf2\" clobber=\"yes\" &> niprefilter2.log # &> /dev/null"
]
},
{
"cell_type": "markdown",
"id": "introductory-thesis",
"metadata": {},
"source": [
"NOTE:\n",
"- The best way to run any heasoft
tool in Python
is using the ${\\tt heasoftpy}$ module.\n",
"- You can run ${\\tt heasoftpy}$ in a IPython
terminal, but unfortunately at this moment, it is not working in Jupyter
notebooks at SciServer."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "accessible-guatemala",
"metadata": {},
"outputs": [],
"source": [
"# If you want to try the heasoftpy uncomment the next lines\n",
"#import heasoftpy as hsp\n",
"#res = hsp.nicerl2(indir=\"./2612011001\", clobber=\"yes\", cor_range=\"1-15\")\n",
"#print(res.stdout)"
]
},
{
"cell_type": "markdown",
"id": "optional-burden",
"metadata": {},
"source": [
"# Creating an Image, Lightcurve, and Spectrum\n",
"\n",
"- Unfortunately SciServer lack the capabilty to make you interact with a ${\\tt WX}$ server to watch LC and do cursor time filter with ${\\tt Xselect}$.\n",
"\n",
"- But you can write a ${\\tt Xselect}$ script to create products, take a look, and comeback to ${\\tt Xselect}$ to make a keyborad time filter. \n",
"\n",
"- For more about keyboard time filter see: https://heasarc.gsfc.nasa.gov/docs/asca/abc/node8.html#SECTION00852000000000000000"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "desirable-ranch",
"metadata": {},
"outputs": [],
"source": [
"# Write the script\n",
"script_xselect = open('script_xselect.xsel', 'w')\n",
"script_xselect.write('session1\\n')\n",
"script_xselect.write('read events ./2612011001/xti/event_cl/ni2612011001_0mpu7_cl.evt\\n')\n",
"script_xselect.write('./\\n')\n",
"script_xselect.write('yes\\n')\n",
"script_xselect.write('extract all\\n')\n",
"script_xselect.write('save all test\\n')\n",
"script_xselect.write('exit\\n')\n",
"script_xselect.write('no\\n')\n",
"script_xselect.close()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "spoken-heater",
"metadata": {},
"outputs": [],
"source": [
"# Run the script\n",
"!xselect @ script_xselect.xsel"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "formed-background",
"metadata": {},
"outputs": [],
"source": [
"# Open and getting LC using astropy.fits\n",
"fits_lc = fits.open('test.lc')\n",
"lc_header = fits_lc[0].header\n",
"lc = fits_lc[1].data\n",
"date = lc_header['DATE']\n",
"obj = lc_header['OBJECT']\n",
"exposure = lc_header['EXPOSURE']\n",
"tstart = float(lc_header['TSTART'])\n",
"fits_lc.close()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "seeing-butter",
"metadata": {},
"outputs": [],
"source": [
"# Take a look at using matplotlib\n",
"fig, ax = plt.subplots(figsize=(8,8))\n",
"ax.plot(lc['TIME'], lc['RATE'])\n",
"ax.set_ylabel('Rates (cnts/s)', fontsize=16)\n",
"ax.set_xlabel('Time (seconds)', fontsize=16)\n",
"ax.set_title(obj+' Date='+date+' Exposure='+str(exposure))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "communist-anthropology",
"metadata": {},
"outputs": [],
"source": [
"# If you want to make a keyboard time filter here is an example\n",
"to = tstart+0.0\n",
"tf = tstart+1000.0\n",
"script_xselect = open('script_xselect.xsel', 'w')\n",
"script_xselect.write('session1\\n')\n",
"script_xselect.write('read events ./2612011001/xti/event_cl/ni2612011001_0mpu7_cl.evt\\n')\n",
"script_xselect.write('./\\n')\n",
"script_xselect.write('yes\\n')\n",
"script_xselect.write('filter time scc\\n')\n",
"script_xselect.write(str(to)+','+str(tf)+'\\n')\n",
"script_xselect.write('x\\n')\n",
"script_xselect.write('extract all\\n')\n",
"script_xselect.write('save all test2\\n')\n",
"script_xselect.write('exit\\n')\n",
"script_xselect.write('no\\n')\n",
"script_xselect.close()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "alike-herald",
"metadata": {},
"outputs": [],
"source": [
"# Run the script to apply filter\n",
"!xselect @ script_xselect.xsel"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "derived-happiness",
"metadata": {},
"outputs": [],
"source": [
"# Open a getting LC (again)\n",
"fits_lc2 = fits.open('test2.lc')\n",
"lc_header2 = fits_lc2[0].header\n",
"lc2 = fits_lc2[1].data\n",
"date2 = lc_header2['DATE']\n",
"obj2 = lc_header2['OBJECT']\n",
"exposure2 = lc_header2['EXPOSURE']\n",
"fits_lc2.close()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "plain-video",
"metadata": {},
"outputs": [],
"source": [
"# Taking a new look at\n",
"fig, ax = plt.subplots(figsize=(8,8))\n",
"ax.plot(lc2['TIME'], lc2['RATE'])\n",
"ax.set_ylabel('Rates (cnts/s)', fontsize=16)\n",
"ax.set_xlabel('Time (seconds)', fontsize=16)\n",
"ax.set_title(obj2+' Date='+date2+' Exposure='+str(exposure2))"
]
},
{
"cell_type": "markdown",
"id": "cross-rotation",
"metadata": {},
"source": [
"# Estimating background"
]
},
{
"cell_type": "markdown",
"id": "juvenile-guidance",
"metadata": {},
"source": [
"## Space Weather (SW) estimator"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "tropical-tuesday",
"metadata": {},
"outputs": [],
"source": [
"#reading space weather python module\n",
"from nicergof.bkg import bkg_estimator as bk"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "royal-heating",
"metadata": {},
"outputs": [],
"source": [
"# Creating mkf3\n",
"mkf2 = '2612011001/auxil/ni2612011001.mkf2'\n",
"bk.add_kp(mkf2, clobber=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "composite-banks",
"metadata": {},
"outputs": [],
"source": [
"# Estimating background\n",
"pha = 'test.pha'\n",
"mkf3 = '2612011001/auxil/ni2612011001.mkf3'\n",
"bevt = '/home/idies/workspace/headata/FTP/caldb/data/nicer/xti/pcf/30nov18targskc_enhanced.evt'\n",
"bkg_chan, bkgspectot, btotexpo = bk.mk_bkg_spec_evt(pha, mkf3, bevt=bevt)"
]
},
{
"cell_type": "markdown",
"id": "informed-cocktail",
"metadata": {},
"source": [
"## NICER Background Generator (3C50)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fundamental-membership",
"metadata": {},
"outputs": [],
"source": [
"!nibackgen3C50 rootdir=\"./\" obsid=\"2612011001\" bkgidxdir=\"./bg_models_3C50\" bkglibdir=\"./bg_models_3C50\" gainepoch=\"2020\""
]
},
{
"cell_type": "markdown",
"id": "asian-haiti",
"metadata": {},
"source": [
"# Reading spectrum with Xspec"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "demographic-management",
"metadata": {},
"outputs": [],
"source": [
"xspec.AllData.clear()\n",
"xspec.AllModels.clear()\n",
"\n",
"s1 = xspec.Spectrum('./test.pha')\n",
"s1.response = '/home/idies/workspace/headata/FTP/caldb/data/nicer/xti/cpf/rmf/nixtiref20170601v002.rmf'\n",
"s1.response.arf = '/home/idies/workspace/headata/FTP/caldb/data/nicer/xti/cpf/arf/nixtiaveonaxis20170601v004.arf'\n",
"#s1.background = './test_bkg.pha'\n",
"s1.background = './nibackgen3C50_bkg.pi'\n",
"date = s1.fileinfo('DATE-OBS')\n",
"exposure = s1.fileinfo('EXPOSURE')\n",
"xspec.Plot.device = '/null'\n",
"#Plot.xAxis = 'angstrom'\n",
"xspec.Plot.xAxis = 'keV'\n",
"xspec.Plot.add = 'True'\n",
"xspec.Plot.setRebin(5.0, 10)\n",
"xspec.Plot('data')\n",
"x1 = np.asarray(xspec.Plot.x(1))\n",
"y1 = np.asarray(xspec.Plot.y(1))\n",
"x_err1 = np.asarray(xspec.Plot.xErr(1))\n",
"y_err1 = np.asarray(xspec.Plot.yErr(1))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "pressed-holder",
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(15,5))\n",
"ax.plot(x1, y1, label=\"2612011001\")\n",
"ax.set_xlim(0.2, 15)\n",
"ax.set_ylim(1e-2,15.0)\n",
"ax.set_yscale('log')\n",
"ax.set_xscale('log')\n",
"ax.set_xlabel('Energy (keV)',fontsize=16)\n",
"ax.set_ylabel('Normalized Counts s$^{-1}$ keV$^{-1}$', fontsize=16)\n",
"#title('NICER Cycle 2 Eta Car Spectra',fontsize = 20)\n",
"plt.minorticks_on()\n",
"plt.tick_params(axis='both', which='both',right = True, top=True,direction='in')\n",
"plt.legend(fontsize=12, ncol=2, numpoints=1, loc=3)\n",
"#savefig('./Eta_Car_spectrum.png', dpi=300)"
]
},
{
"cell_type": "markdown",
"id": "lesser-suggestion",
"metadata": {},
"source": [
"# Fitting a Model"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bridal-cathedral",
"metadata": {},
"outputs": [],
"source": [
"# Reading a model with Xspec\n",
"model = xspec.Model('tbabs * bvapec')\n",
"model.TBabs.nH.values = [1.0]\n",
"model.bvapec.kT.values = [1.0]\n",
"model.bvapec.He.values = [1.0]\n",
"model.bvapec.C.values = [1.0]\n",
"model.bvapec.N.values = [1.0]\n",
"model.bvapec.O.values = [1.0]\n",
"model.bvapec.Ne.values = [1.0]\n",
"model.bvapec.Mg.values = [1.0]\n",
"model.bvapec.Al.values = [1.0]\n",
"model.bvapec.Si.values = [1.0]\n",
"model.bvapec.S.values = [1.0]\n",
"model.bvapec.Ar.values = [1.0]\n",
"model.bvapec.Ca.values = [1.0]\n",
"model.bvapec.Fe.values = [1.0]\n",
"model.bvapec.Ni.values = [1.0]\n",
"model.bvapec.Redshift.values = [0.0]\n",
"model.bvapec.Velocity.values = [200]\n",
"model.bvapec.norm.values = [1.0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "senior-baseline",
"metadata": {},
"outputs": [],
"source": [
"# Creating a pandas DataFrame to save parameters\n",
"df_nicer = pd.DataFrame()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "severe-pride",
"metadata": {},
"outputs": [],
"source": [
"# Performing the fit with some other commented options that you may be interested \n",
"\n",
"#xspec.XspecSettings.abund = 'wilm'\n",
"#xspec.Fit.statMethod = 'chi'\n",
"xspec.Fit.query = \"yes\"\n",
"xspec.Fit.perform()\n",
"xspec.Plot.device = \"/null\"\n",
"xspec.Plot.xAxis = \"keV\"\n",
"#xspec.Plot.add = 'True'\n",
"#Plot.xAxis = \"angstrom\"\n",
"#s.ignore('**-1.7 10.0-**')\n",
"#xspec.Plot.setRebin(3.0, 10)\n",
"xspec.Plot('data')\n",
"x = np.asarray(xspec.Plot.x()) \n",
"y = np.asarray(xspec.Plot.y())\n",
"x_err = np.asarray(xspec.Plot.xErr())\n",
"y_err = np.asarray(xspec.Plot.yErr())\n",
"mod = xspec.Plot.model()\n",
"xspec.AllModels.calcFlux(\"2.0 10.0\")\n",
"df_nicer.at['2612011001', 'FLUX_2.0-10.0'] = s1.flux[0]\n",
"df_nicer.at['2612011001', 'Fit_Statistic'] = xspec.Fit.statistic\n",
"df_nicer.at['2612011001', 'dof'] = xspec.Fit.dof\n",
"red_chi = xspec.Fit.statistic/xspec.Fit.dof\n",
"for c in model.componentNames:\n",
" for p in model.__getattribute__(c).parameterNames:\n",
" df_nicer.at['2612011001', c+'_'+p] = model.__getattribute__(c).__getattribute__(p).values[0]\n",
" df_nicer.at['2612011001', c+'_'+p+'_err'] = model.__getattribute__(c).__getattribute__(p).sigma\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "essential-hypothetical",
"metadata": {},
"outputs": [],
"source": [
"# Take a look to the dataframe with parameters for OBS ID 2612011001\n",
"df_nicer"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "following-nurse",
"metadata": {},
"outputs": [],
"source": [
"# Taking a look at the spectrum and the fitting \n",
"\n",
"fig, ax = plt.subplots(2, sharex=True, figsize=(15, 8), gridspec_kw={'height_ratios': [3, 1]})\n",
"fig.subplots_adjust(hspace=0)\n",
"ax[0].errorbar(x, y, xerr=x_err, yerr=y_err, fmt='.', alpha=0.7)\n",
"ax[0].plot(x,mod, label='$\\chi^{2}_{red} = $'+str(red_chi))\n",
"ax[0].legend()\n",
"ax[0].minorticks_on()\n",
"ax[0].tick_params(axis='both', which='both',right = True, top=True,direction='in')\n",
"ax[0].set_ylabel('normalized $counts$ $s^{-1}$ $chan^{-1}$')\n",
"ax[1].set_xlabel('Energy (keV)')\n",
"ax[0].set_xlim(0.4, 12)\n",
"ax[0].set_ylim(1e-2,17.0)\n",
"ax[0].set_yscale('log')\n",
"ax[0].set_xscale('log')\n",
"ax[1].plot(x, y-mod)\n",
"ax[1].set_ylabel('Residual')\n",
"plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.8 (Heasarc)",
"language": "python",
"name": "python3.8"
},
"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"
},
"toc-autonumbering": true,
"toc-showcode": false,
"toc-showmarkdowntxt": false,
"toc-showtags": false
},
"nbformat": 4,
"nbformat_minor": 5
}